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