Attachment 'GLM1_script.m'
Download 1 %% Basis functions and matrix inversion
2
3 1*[1 0 0] + 2*[0 1 0] + 3*[0 0 1]
4
5 1*[1 1 1] + 1*[0 1 0] + 2*[0 0 1]
6
7 -1*[1 -1 0] + -1*[0 1 -1] + 2*[1 1 1]
8
9 % creating a linear system
10 bf1 = [1 -1 0];
11 bf2 = [0 1 -1];
12 bf3 = [1 1 1];
13
14 % You may call this the "design matrix"
15 BF = [bf1', bf2', bf3'];
16
17 data = [1 2 3]';
18
19 % solve equation
20 sol = inv(BF)*data
21
22 % check if it's really a solution
23 BF*sol
24
25 % check if inverse does what it's supposed to do
26 inv(BF)*BF
27 BF*inv(BF)
28
29 % note: matrix inversion is not element-wise division
30 inv(BF)
31 BF.^-1
32
33 % you can also use
34 y = BF \ data
35 y = pinv(BF)*data
36
37 % let's introduce collinearity and see what happens
38 BF(:,1) = [1 1 1]'
39
40 % let's "fix" the collinearity
41 BF(1,3) = 1 + 1e-12
42
43 inv(BF)
44 % the inverse matrix has huge values - this may amplify errors in the data
45
46 %% parameter estimation - forward model
47 x = 0:0.1:8*pi; % abscissa values
48
49 figure;
50 y1 = x'; % straight line (note transpose)
51 plot(x, y1);
52
53 figure;
54 y2 = sin(x)'; % sine curve
55 plot(x, y2);
56
57 figure;
58 y3 = exp(-x/2)'; % exponential function
59 plot(x, y3);
60
61 % Generate data with some noise
62 figure;
63 y = y1/5 + y2 + 5*y3 + 0.5*randn(length(y1),1);
64 plot(x, y);
65
66 whos
67 %% Estimating the parameters using minimum least-squares (MNLS)
68 % creating the design matrix
69 M = [y1 y2 y3];
70 whos
71 % note that in this case there are many more data points than parameters
72 % i.e. the problem is over-determined
73
74 % check correlation among basis functions
75 cM = corr(M)
76 imagesc(cM)
77 colorbar
78
79 % (pseudo)inverting the design matrix
80 Minv = pinv(M);
81 whos
82
83 % getting the solution
84 b = Minv*y
85
86 % checking how solution predicts the data
87 figure;
88 ypred = b(1)*y1 + b(2)*y2 + b(3)*y3;
89 plot(x, ypred)
90
91 % compare with "ground truth"
92 hold on; % don't overwrite figure
93 yreal = y1/5 + y2 + 5*y3;
94 plot(x, yreal, 'r');
95
96 % plot the difference between measured and predicted data
97 figure;
98 ydiff = y - ypred;
99 plot(x, ydiff);
100
101 % compute your very own "pseudoinverse"
102 mypinv = inv(M'*M)*M';
103
104 % check if the results is the same
105 b2 = mypinv*y;
106
107 % it could have been so easy:
108 b3 = regress(y, M)
Attached Files
To refer to attachments on a page, use attachment:filename, as shown below in the list of files. Do NOT use the URL of the [get] link, since this is subject to change and can break easily.You are not allowed to attach a file to this page.