line source
+ − % combineExps combine the results of different parameter estimation
+ − % experimets and give the final joint estimate
+ − %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ − %
+ − % DESCRIPTION: combineExps combine different parameter estimates into a
+ − % single joint estimate. Joint covariance is also computed.
+ − %
+ − % CALL: obj = combineExps(objs);
+ − % obj = combineExps(objs);
+ − %
+ − % INPUTS: obj - can be a vector
+ − %
+ − % <a href="matlab:utils.helper.displayMethodInfo('pest', 'combineExps')">Parameters Description</a>
+ − %
+ − % VERSION: $Id: combineExps.m,v 1.9 2011/04/08 08:56:25 hewitson Exp $
+ − %
+ − % HISTORY: 10-12-2010 G. Congedo
+ − %
+ − %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ −
+ − function varargout = combineExps(varargin)
+ −
+ − %%% Check if this is a call for parameters
+ − if utils.helper.isinfocall(varargin{:})
+ − varargout{1} = getInfo(varargin{3});
+ − return
+ − end
+ −
+ − import utils.const.*
+ − utils.helper.msg(msg.PROC3, 'running %s/%s', mfilename('class'), mfilename);
+ −
+ − % Collect input variable names
+ − in_names = cell(size(varargin));
+ − for ii = 1:nargin,in_names{ii} = inputname(ii);end
+ −
+ − % Collect all AOs and plists
+ − [pests, pest_invars] = utils.helper.collect_objects(varargin(:), 'pest', in_names);
+ − pl = utils.helper.collect_objects(varargin(:), 'plist', in_names);
+ −
+ −
+ − % combine plists
+ − % pl = parse(pl, getDefaultPlist());
+ −
+ − if nargout == 0
+ − error('### combineExps cannot be used as a modifier. Please give an output variable.');
+ − end
+ −
+ − if ~all(isa(pests, 'pest'))
+ − error('### combineExps must be only applied to pest objects.');
+ − end
+ −
+ − N = numel(pests);
+ −
+ − I = cell(N,1);
+ − C = cell(N,1);
+ − p = cell(N,1);
+ − pnames = cell(N,1);
+ − chi2 = zeros(N,1);
+ − dof = zeros(N,1);
+ −
+ − % Iterate over all input pest
+ − for ii=1:N
+ −
+ − counter = sprintf('%u',ii);
+ −
+ − % Extract values
+ − pnames{ii} = pests(ii).names;
+ − p{ii} = pests(ii).y;
+ − chi2(ii) = pests(ii).chi2;
+ − if isempty(chi2(ii))
+ − error(['### Could not find chi2 value for input pest #' counter]);
+ − end
+ − dof(ii) = pests(ii).dof;
+ − if isempty(dof(ii))
+ − error(['### Could not find dof value for input pest #' counter]);
+ − end
+ − I{ii} = pests(ii).procinfo.find('infomatrix');
+ − C{ii} = pests(ii).cov;
+ − if isempty(C{ii})
+ − C{ii} = inv(I{ii});
+ − elseif isempty(I{ii})
+ − I{ii} = inv(C{ii});
+ − elseif isempty(I{ii}) && isempty(C{ii});
+ − error(['### Could not find covariance matrix for input pest #' counter]);
+ − end
+ −
+ − end
+ −
+ − % Now redefine all quantities to handle the general case when we have
+ − % different sizes
+ − Inew = cell(N,1);
+ − pnew = cell(N,1);
+ − pnamesAll = pnames{1};
+ − dofnew = zeros(N,1);
+ − chi2new = zeros(N,1);
+ − for ii=1:N
+ − pnamesAll = union(pnamesAll,pnames{ii});
+ − end
+ − Np = numel(pnamesAll);
+ − for ii=1:N
+ − Inew{ii} = zeros(Np);
+ − pnew{ii} = zeros(Np,1);
+ − % dofnew{ii} = zeros(Np,1);
+ − % information matrix
+ − for kk=1:Np
+ − for ll=1:Np
+ − ixk = strcmp(pnamesAll{kk},pnames{ii});
+ − ixl = strcmp(pnamesAll{ll},pnames{ii});
+ − ixk = find(ixk);
+ − ixl = find(ixl);
+ − if ixk~=0 & ixl~=0
+ − Inew{ii}(kk,ll) = I{ii}(ixk,ixl);
+ − end
+ − end
+ − end
+ − % param vector
+ − for kk=1:Np
+ − ix = strcmp(pnamesAll{kk},pnames{ii});
+ − ix = find(ix);
+ − if ix~=0
+ − pnew{ii}(kk) = p{ii}(ix);
+ − end
+ − end
+ − % dof & chi2
+ − dofnew(ii) = dof(ii) + numel(p{ii}) - Np;
+ − chi2new(ii) = chi2(ii) * dof(ii) / dofnew(ii);
+ −
+ − % for kk=1:Np
+ − % ix = strcmp(pnamesAll,pnames{ii}{kk});
+ − % ix = find(ix);
+ − % arr = I{ii}(kk,:);
+ − % % (1:kk) = arr(1:ix)
+ − % % Inew{ii}(ix,:) = [0 arr(kk:end)];
+ − % % pnew{ii}(kk) = p{ii}(ix);
+ − % end
+ − end
+ − I = Inew;
+ − p = pnew;
+ − dof = dofnew;
+ − chi2 = chi2new;
+ −
+ − % Compute joint information matrix
+ − II = zeros(size(I{1}));
+ − for ii=1:N
+ − II = II+I{ii};
+ − end
+ − % CC = inv(II);
+ −
+ − % Compute joint parameter estimate
+ − pp = zeros(size(p{1}));
+ − for ii=1:N
+ − pp = pp+I{ii}*p{ii};
+ − end
+ − pp = II\pp;
+ −
+ − % Compute joint errors and correlation matrix
+ − CC = inv(II);
+ − dp = sqrt(diag(CC));
+ − corr = zeros(size(CC));
+ − for ll=1:size(CC,1)
+ − for mm=1:size(CC,2);
+ − corr(ll,mm) = CC(ll,mm)/dp(ll)/dp(mm);
+ − end
+ − end
+ −
+ − % Compute joint chi2 and dof
+ − DOF = sum(dof);
+ − CHI2 = chi2'*dof/DOF;
+ −
+ − % Output pest
+ − out = pest();
+ − out = out.setNames(pnamesAll);
+ − % out = out.setYunits(pests(1).yunits);
+ − out = out.setY(pp);
+ − out = out.setCov(CC);
+ − out = out.setDy(dp);
+ − out = out.setCorr(corr);
+ − out = out.setChi2(CHI2);
+ − out = out.setDof(DOF);
+ − name = pests(1).name;
+ − if N>1
+ − for ii=2:N
+ − name = [name ',' pests(ii).name];
+ − end
+ − end
+ − out = out.setName(['combineExps(' name ')']);
+ −
+ − out.procinfo = plist('infomatrix',II);
+ − out.addHistory(getInfo('None'), pl, pest_invars(:), [pests(:).hist]);
+ −
+ − % Set outputs
+ − if nargout > 0
+ − varargout{1} = out;
+ − end
+ −
+ − end
+ −
+ −
+ − %--------------------------------------------------------------------------
+ − % Get Info Object
+ − %--------------------------------------------------------------------------
+ − function ii = getInfo(varargin)
+ − if nargin == 1 && strcmpi(varargin{1}, 'None')
+ − sets = {};
+ − pl = [];
+ − else
+ − sets = {'Default'};
+ − pl = getDefaultPlist;
+ − end
+ − % Build info object
+ − ii = minfo(mfilename, 'pest', 'ltpda', utils.const.categories.helper, '$Id: combineExps.m,v 1.9 2011/04/08 08:56:25 hewitson Exp $', sets, pl);
+ − end
+ −
+ − %--------------------------------------------------------------------------
+ − % Get Default Plist
+ − %--------------------------------------------------------------------------
+ − function plout = getDefaultPlist()
+ − persistent pl;
+ − if exist('pl', 'var')==0 || isempty(pl)
+ − pl = buildplist();
+ − end
+ − plout = pl;
+ − end
+ −
+ − function plo = buildplist()
+ − plo = plist.EMPTY_PLIST;
+ − end
+ −