--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/m-toolbox/classes/@ao/mcmc_td.m	Wed Nov 23 19:22:13 2011 +0100
@@ -0,0 +1,275 @@
+% MCMC_TD estimates paramters using a Monte Carlo Markov Chain.
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%
+% DESCRIPTION: MCMC estimate the parameters of a given model given
+%              inputs, outputs and noise using a Metropolis algorithm.
+%              
+%
+% CALL:        b = mcmc(out,pl)
+%
+% INPUTS:      
+%              out     - analysis objects with measured outputs
+%              pl      - parameter list
+%
+% OUTPUTS:     b   - pest object contatining estimate information
+%
+% <a href="matlab:utils.helper.displayMethodInfo('ao', 'mcmc_td')">Parameters Description</a>
+%
+% VERSION:    $Id: mcmc_td.m,v 1.6 2011/04/08 08:56:12 hewitson Exp $
+%
+% References:  "Catching supermassive black holes binaries without a net"
+%              N.J. Cornish, E.K. Porter, Phys.Rev.D 75, 021301, 2007
+%
+% TODO: multiple chain option not implemented yet
+%       metropolis/hastings not implemented
+%       empty initial values not implemented
+%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+function varargout = mcmc_td(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);
+  
+  % Method can not be used as a modifier
+  if nargout == 0
+    error('### metropolis2D cannot be used as a modifier. Please give an output variable.');
+  end
+  
+  % Collect input variable names
+  in_names = cell(size(varargin));
+  for ii = 1:nargin,in_names{ii} = inputname(ii);end
+  
+  % Collect all AOs smodels and plists
+  [as, ao_invars] = utils.helper.collect_objects(varargin(:), 'ao', in_names);
+  pl              = utils.helper.collect_objects(varargin(:), 'plist', in_names);
+  
+  % Decide on a deep copy or a modify
+  bs = copy(as, nargout);
+  
+  out = bs;
+  
+  
+  % Combine plists
+  pl = parse(pl, getDefaultPlist);
+  
+  
+  % Get parameters
+  N = find(pl, 'N');
+  Tc = find(pl,'Tc');
+  xi = find(pl,'heat');
+  cvar  = find(pl, 'cov');
+  rng = find(pl,'range');
+  search = find(pl,'search');
+  simplex = find(pl,'simplex');
+  x0 = find(pl,'x0');
+  mdl = find(pl,'model');
+  %nse = find(pl,'noise');
+  jumps = find(pl,'jumps');
+  parplot = find(pl,'plot');
+  debug = find(pl,'debug');
+  
+  in = find(pl,'Input');
+  inNames = find(pl,'inNames');
+  outNames = find(pl,'outNames');
+  Noise = find(pl,'Noise');
+  parnames = find(pl,'parnames');
+  cutbefore = find(pl,'cutbefore');
+  cutafter = find(pl,'cutafter');
+  
+  % Check input parameters
+  if isempty(rng)
+    error('### Please define a search range ''range''');
+  end
+
+  if isempty(parnames)
+    error('### Please define the parameters ''parnames''');
+  end
+  nparam = numel(parnames);
+
+  
+  % Get range for parameters
+  for i = 1:nparam
+    range(:,i) = rng{i};
+  end
+  
+  
+  
+  % do simplex
+  if simplex
+    if isempty(x0)
+      error('### Simplex needs a starting guess. Please input a ''x0''.');
+    else
+      switch class(mdl)
+        case 'smodel'
+          %xo = fminsearch(@(x) utils.math.loglikehood(x,in,out,nse,mdl,param),xo,optimset('Display','iter'));
+        case 'ssm'
+          x0 = fminsearch(@(x)utils.math.loglikehood_ssm_td(x,in,out,parnames,mdl,inNames,outNames,Noise,'cutbefore',cutbefore,'cutafter',cutafter),x0,optimset('Display','iter'));
+        otherwise
+          error('### Model must be either from the ''smodel'' or the ''ssm'' class. Please check the inputs')
+      end
+      
+      for i = 1:numel(parnames)
+        fprintf('###  Simplex estimate: %s = %d \n',parnames{i},x0(i))
+        %save('parameters_simplex.txt','x0','-ASCII')
+      end
+    end
+  end
+  
+  % sample distribution
+  smpl = utils.math.mhsample_td(mdl,in,out,cvar,N,range,parnames,Tc,xi,x0,search,jumps,parplot,debug,inNames,outNames,Noise,cutbefore,cutafter);
+  
+  % statistics of the chain
+  if isempty(Tc)
+    initial =1;
+  else
+    initial = Tc(2)+1;
+  end
+  
+  mn = mean(smpl(initial:end,:));
+  cv = cov(smpl(initial:end,:));
+  
+  cr = utils.math.cov2corr(cv);
+  
+  % compute pdf
+  for i = 1:nparam
+    [count,bin] = hist(smpl(initial:end,i));
+  end
+  
+  % create pest output
+  p = pest(mn);
+  p.setName('metropolis2D');
+  p.setNames(parnames{:});
+  % add statistical info
+  p.setCov(cv);
+  p.setCorr(cr);
+  p.setDy(sqrt(diag(cv))); % need to change if non-gaussian
+  p.setPdf([reshape(bin, 1, []), reshape(count, 1, [])]); % not yet able to store bins
+  p.setChain(smpl);
+  p.setModels(mdl);
+  % set history
+  p.addHistory(getInfo('None'), pl, [ao_invars(:)], [bs(:).hist mdl(:).hist]);
+  
+  % Set output
+  varargout{1} = p;
+  
+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, 'ao', 'ltpda', utils.const.categories.sigproc, '$Id: mcmc_td.m,v 1.6 2011/04/08 08:56:12 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 pl = buildplist()
+  pl = plist();
+  
+  % N
+  p = param({'N','number of samples of the chain.'}, paramValue.DOUBLE_VALUE(1000));
+  pl.append(p);
+  
+  % Sigma
+  p = param({'cov','covariance of the gaussian jumping distribution.'}, paramValue.DOUBLE_VALUE(1e-4));
+  pl.append(p);
+  
+  % Param
+  p = param({'parnames','A cell array of evaluated parameters.'}, paramValue.EMPTY_CELL);
+  pl.append(p);
+  
+  
+  % Model
+  p = param({'model','The model for the system, at the moment only ssm is supported'}, paramValue.EMPTY_STRING);
+  pl.append(p);
+  
+  % Search
+  p = param({'search','Set to true to use bigger jumps in parameter space during annealing and cool down.'}, paramValue.TRUE_FALSE);
+  pl.append(p);
+  
+  % Simplex
+  p = param({'simplex','Set to true to perform a simplex search to find the starting parameters of the MCMC chain.'}, paramValue.TRUE_FALSE);
+  pl.append(p);
+  
+  % heat
+  p = param({'heat','The heat index flattening likelihood surface during annealing.'}, paramValue.DOUBLE_VALUE(1));
+  pl.append(p);
+  
+  % Tc
+  p = param({'Tc','An array of two values setting the initial and final value for the cooling down.'}, paramValue.EMPTY_STRING);
+  pl.append(p);
+  
+  % x0
+  p = param({'x0','The proposed initial values.'}, paramValue.EMPTY_DOUBLE);
+  pl.append(p);
+  
+  % jumps
+  p = param({'jumps','An array of four numbers setting the rescaling of the covariance matrix during the search phase.',...
+    'The first value is the one applied by default, the following thhree apply just when the chain sample is',...
+    'mod(10), mod(25) and mod(100) respectively.'}, paramValue.EMPTY_DOUBLE);
+  pl.append(p);
+  
+  % plot
+  p = param({'plot','Select indexes of the parameters to be plotted.'}, paramValue.EMPTY_DOUBLE);
+  pl.append(p);
+  
+  % debug
+  p = param({'debug','Set to true to get debug information of the MCMC process.'}, paramValue.FALSE_TRUE);
+  pl.append(p);
+  
+  % Input
+  p = param({'Input','An array of analysis object containing the input signals'},paramValue.EMPTY_DOUBLE);
+  pl.append(p);
+  
+  % inNames
+  p = param({'inNames','Cell array of input port names'}, paramValue.EMPTY_CELL);
+  pl.append(p);
+  
+  % outNames
+  p = param({'outNames','Cell array of output port names'}, paramValue.EMPTY_CELL);
+  pl.append(p);
+  
+  
+  % Noise
+  p = param({'Noise','An array of analysis objects containg noise series'}, paramValue.EMPTY_DOUBLE);
+  pl.append(p);
+  
+  % cutbefore
+  p = param({'cutbefore','the data samples to cut at the starting of the series'}, paramValue.EMPTY_DOUBLE);
+  pl.append(p);
+  
+  % cutafter
+  p = param({'cutafter','the data samples to cut at the ending of the series'}, paramValue.EMPTY_DOUBLE);
+  pl.append(p);
+  
+  
+end
+
+
+