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