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)