ltpda: m-toolbox/classes/@pest/combineExps.m comparison

comparison m-toolbox/classes/@pest/combineExps.m @ 0:f0afece42f48

Import.

author	Daniele Nicolodi <nicolodi@science.unitn.it>
date	Wed, 23 Nov 2011 19:22:13 +0100
parents
children

comparison

equal deleted inserted replaced

--1:000000000000
+:f0afece42f48
+% 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

Mercurial > hg > ltpda

comparison m-toolbox/classes/@pest/combineExps.m @ 0:f0afece42f48