line source
+ − % 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
+ −
+ −
+ −
+ −
+ −
+ −
+ −
+ −
+ −