%% 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)