Mercurial > hg > ltpda
view 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 |
line wrap: on
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