+
−% Linear fit with singular value decomposition
+
−% 
+
−% INPUT:
+
−%   - exps: is a N dimensional struct whose fields are fitdata
+
−%   (experimental data), fitbasis (basis for the fit) and params (cell
+
−%   array with the names of the parameters used in the experiment).
+
−%   - groundexps: a M dim struct containing results from on ground experiments.
+
−%   Fuields are: pos (parameter position number), value (on ground measurement
+
−%   result) and err (experimental error on the on ground measurement).
+
−%   - sThreshold it's a threshold for singular values. It is a
+
−%     number, typically 1. It will remove singular values larger
+
−%     than sThreshold which corresponds to removing svd parameters estimated
+
−%     with an error larger than sThreshold.
+
−% 
+
−% OUTPUT:
+
−%   a:      params values
+
−%   Ca:     fit covariance matrix for A
+
−%   Corra:  fit correlation matrix for A
+
−%   Vu:     is the complete conversion matrix
+
−%   Cbu:    is the new variables covariance matrix
+
−%   Fbu:    is the information matrix for the new variable
+
−%   mse:    is the fit Mean Square Error
+
−%   dof:    degrees of freedom for the global estimation
+
−%   ppm:    number of svd parameters per measurements, provides also the
+
−%   number of independent combinations of parameters per each singular
+
−%   measurement. The coefficients of the combinations are then stored in Vu
+
−% 
+
−% 21-07-2009 L Ferraioli
+
−%       CREATION
+
−% 
+
−% 
+
−%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
−% $Id: linfitsvd.m,v 1.8 2010/10/29 12:40:47 ingo Exp $
+
−%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
−function [a,Ca,Corra,Vu,bu,Cbu,Fbu,mse,dof,ppm] = linfitsvd(varargin)
+
−  
+
−  % get inputs
+
−  if nargin == 1 % no ground experiment
+
−    exps = varargin{1};
+
−    groundexps = [];
+
−    remerr = false;
+
−  elseif nargin == 2
+
−    exps = varargin{1};
+
−    if isnumeric(varargin{2})
+
−      sThreshold = varargin{2};
+
−      remerr = true;
+
−      groundexps = [];
+
−    else
+
−      groundexps = varargin{2};
+
−      remerr = false;
+
−    end
+
−  elseif nargin == 3
+
−    exps = varargin{1};
+
−    groundexps = varargin{2};
+
−    sThreshold = varargin{3}; % admitted threshold for singular values
+
−    remerr = true;
+
−    if ~isnumeric(sThreshold)
+
−      sThreshold = inf;
+
−    end
+
−  end
+
−  
+
−  % init collect results struct
+
−  
+
−  % init union matrices
+
−  Vu = [];
+
−  Fbu = [];
+
−  wFbu = [];
+
−  bu = [];
+
−%   eCbu = [];
+
−  Cbu = [];
+
−  ppm = []; % number of svd parameters per measurements
+
−  
+
−  N = numel(exps);
+
−  % run over input experiments
+
−  for ii=1:N
+
−    % get data of the experiments out of the input struct
+
−    y = exps(ii).fitdata;
+
−    X = exps(ii).fitbasis;
+
−    
+
−    % willing to work with columns
+
−    if size(y,1)<size(y,2)
+
−      y = y.';
+
−    end
+
−    if size(X,1)<size(X,2)
+
−      X = X.';
+
−    end
+
−    
+
−    % get svd of X
+
−    [U,W,V] = svd(X,0);
+
−    
+
−    % removing zero singular values
+
−    svals = diag(W);    
+
−    idx = svals==0;    
+
−    svals(idx) = [];
+
−    
+
−    % removing columns of U and V corresponding to zero svalues
+
−    U(:,idx) = [];
+
−    V(:,idx) = [];
+
−    
+
−    if remerr
+
−      % Find params combinations with error larger than 1
+
−      idx = svals<sThreshold;
+
−      svals(idx) = [];
+
−
+
−      % removing columns of U and V corresponding to params combinations with
+
−      % error larger than 1
+
−      U(:,idx) = [];
+
−      V(:,idx) = [];
+
−    end
+
−    
+
−    % Sanity check
+
−    % if svals is empty you are going to gain no information from the
+
−    % current data series
+
−    if ~isempty(svals) % go with the fit
+
−      % rebuild W
+
−      W = diag(svals);
+
−      
+
−      % update ppm with the number of parameters for the given experiment
+
−      ppm = [ppm; numel(svals)];
+
−
+
−      % get new basis for the fit
+
−      K = U*W;
+
−
+
−      % get expected covariance matrix for b, assuming white noise
+
−      Fb = W*W; % information matrix for parameters b
+
−
+
−      % get b params with errors, mse is the mean square error, Cb is the
+
−      % covariance matrix of fitted b. It should be equal to eCb.*mse
+
−      [b,stdxb,mseb,Cb] = lscov(K,y);
+
−
+
−      % information matrix weighted for fit results Cb = inv(wFb)
+
−      wFb = Fb./mseb;
+
−
+
−      % add params from present experiment
+
−      % get union matrices
+
−      Vu = [Vu;V.'];
+
−      Fbu = blkdiag(Fbu,Fb);
+
−      wFbu = blkdiag(wFbu,wFb);
+
−      Cbu = blkdiag(Cbu,Cb);
+
−      bu = [bu;b];
+
−    
+
−    else % no information, no fit
+
−      % update ppm with the number of parameters for the given experiment
+
−      ppm = [ppm; 0];
+
−    end
+
−      
+
−    
+
−    
+
−  end
+
−  
+
−  % insert values from on-ground measured parameters if applicable
+
−  if nargin == 2
+
−    
+
−    for jj = 1:numel(groundexps)
+
−
+
−      % get values
+
−      pos = groundexps(jj).pos;
+
−      val = groundexps(jj).value;
+
−      err = groundexps(jj).err;
+
−
+
−      tV = zeros(1,size(Vu,2));
+
−      tV(1,pos) = 1;
+
−      Vu = [Vu;tV];
+
−      bu = [bu;val];
+
−      Cbu = blkdiag(Cbu,err^2);
+
−      Fbu = blkdiag(Fbu,1/(err^2));
+
−      wFbu = blkdiag(wFbu,1/(err^2));
+
−      ppm = [ppm; 1];
+
−
+
−    end
+
−    
+
−  end
+
−    
+
−  % get results on physical parameters, solve the system with proper weights
+
−  a = lscov(Vu,bu,1./diag(Cbu));
+
−    
+
−  % get params covariance
+
−  Ca = inv(Vu'*wFbu*Vu);
+
−  
+
−  % get params correlation matrix
+
−  Corra = Ca;
+
−  for tt=1:size(Corra,1)
+
−    for hh=1:size(Corra,2)
+
−      Corra(tt,hh) = Ca(tt,hh)/(sqrt(Ca(tt,tt))*sqrt(Ca(hh,hh)));
+
−    end
+
−  end
+
−  
+
−  % get chi square assuming unitary variance white noise
+
−  N = numel(exps);
+
−  mse = 0;
+
−  tL = 0;
+
−  % run over input experiments
+
−  for ii=1:N
+
−    % get data of the experiments out of the input struct
+
−    y = exps(ii).fitdata;
+
−    X = exps(ii).fitbasis;
+
−    
+
−    mse = mse + sum((y - X*a).^2);
+
−    tL = tL + numel(y);
+
−  end
+
−  dof = tL-numel(a);
+
−  mse = mse./dof;
+
−
+
−end
+
−
+
−
+
−
+
−
+
−
+
−
+
−
+
−
+
−