%% Basis functions and matrix inversion
1*[1 0 0] + 2*[0 1 0] + 3*[0 0 1]
1*[1 1 1] + 1*[0 1 0] + 2*[0 0 1]
-1*[1 -1 0] + -1*[0 1 -1] + 2*[1 1 1]
% creating a linear system
bf1 = [1 -1 0];
bf2 = [0 1 -1];
bf3 = [1 1 1];
% You may call this the "design matrix"
BF = [bf1', bf2', bf3'];
data = [1 2 3]';
% solve equation
sol = inv(BF)*data
% check if it's really a solution
BF*sol
% check if inverse does what it's supposed to do
inv(BF)*BF
BF*inv(BF)
% note: matrix inversion is not element-wise division
inv(BF)
BF.^-1
% you can also use
y = BF \ data
y = pinv(BF)*data
% let's introduce collinearity and see what happens
BF(:,1) = [1 1 1]'
% let's "fix" the collinearity
BF(1,3) = 1 + 1e-12
inv(BF)
% the inverse matrix has huge values - this may amplify errors in the data
%% parameter estimation - forward model
x = 0:0.1:8*pi; % abscissa values
figure;
y1 = x'; % straight line (note transpose)
plot(x, y1);
figure;
y2 = sin(x)'; % sine curve
plot(x, y2);
figure;
y3 = exp(-x/2)'; % exponential function
plot(x, y3);
% Generate data with some noise
figure;
y = y1/5 + y2 + 5*y3 + 0.5*randn(length(y1),1);
plot(x, y);
whos
%% Estimating the parameters using minimum least-squares (MNLS)
% creating the design matrix
M = [y1 y2 y3];
whos
% note that in this case there are many more data points than parameters
% i.e. the problem is over-determined
% check correlation among basis functions
cM = corr(M)
imagesc(cM)
colorbar
% (pseudo)inverting the design matrix
Minv = pinv(M);
whos
% getting the solution
b = Minv*y
% checking how solution predicts the data
figure;
ypred = b(1)*y1 + b(2)*y2 + b(3)*y3;
plot(x, ypred)
% compare with "ground truth"
hold on; % don't overwrite figure
yreal = y1/5 + y2 + 5*y3;
plot(x, yreal, 'r');
% plot the difference between measured and predicted data
figure;
ydiff = y - ypred;
plot(x, ydiff);
% compute your very own "pseudoinverse"
mypinv = inv(M'*M)*M';
% check if the results is the same
b2 = mypinv*y;
% it could have been so easy:
b3 = regress(y, M)