Mercurial > hg > ltpda
view m-toolbox/classes/@ao/mcmc_td.m @ 30:317b5f447f3e database-connection-manager
Update workspaceBrowser
author | Daniele Nicolodi <nicolodi@science.unitn.it> |
---|---|
date | Mon, 05 Dec 2011 16:20:06 +0100 |
parents | f0afece42f48 |
children |
line wrap: on
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