line source
+ − % MCMC 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 smplr] = mcmc(out,pl)
+ − %
+ − % INPUTS: out - matrix objects with measured outputs
+ − % pl - parameter list
+ − %
+ − % OUTPUTS: b - pest object contatining estimate information
+ − % smplr - matrix containing info about the rejected points
+ − %
+ − % <a href="matlab:utils.helper.displayMethodInfo('matrix', 'mcmc')">Parameters Description</a>
+ − %
+ − % VERSION: $Id: mcmc.m,v 1.35 2011/11/16 15:21:13 nikos Exp $
+ − %
+ − % References: M Nofrarias et al. Phys. Rev. D 82, 122002 (2010)
+ − %
+ − %
+ − %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ −
+ − % TODO: multiple chain option not implemented yet
+ − % metropolis/hastings not implemented
+ − % empty initial values not implemented
+ −
+ −
+ − function varargout = mcmc(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('### mcmc 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
+ − [mtxs, mtxs_invars] = utils.helper.collect_objects(varargin(:), 'matrix', in_names);
+ − pl = utils.helper.collect_objects(varargin(:), 'plist', in_names);
+ −
+ − % 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');
+ − flim = find(pl,'frequencies');
+ − fsout = find(pl,'fsout');
+ − search = find(pl,'search');
+ − simplex = find(pl,'simplex');
+ − mhsampleTrue = find(pl,'mhsample');
+ − xo = find(pl,'x0');
+ − mtxin = find(pl,'input');
+ − mdlin = find(pl,'model');
+ − mdlFreqDependent = find(pl,'modelFreqDependent');
+ − nse = find(pl,'noise');
+ − jumps = find(pl,'jumps');
+ − parplot = find(pl,'plot');
+ − debug = find(pl,'debug');
+ − outModel = find(pl,'outModel');
+ − inModel = find(pl,'inModel');
+ − outNames = find(pl,'outNames');
+ − inNames = find(pl,'inNames');
+ − fpars = find(pl,'prior');
+ − anneal = find(pl,'anneal');
+ − DeltaL = find(pl,'DeltaL');
+ − SNR0 = find(pl,('SNR0'));
+ −
+ − % Decide on a deep copy or a modify
+ − in = copy(mtxin, nargout);
+ − out = copy(mtxs, nargout);
+ − mdl = copy(mdlin, nargout);
+ −
+ − % checking if input or noise are aos
+ − if isa(in, 'ao')
+ − in = matrix(in,plist('shape',[size(in,1) size(in,2)]));
+ − end
+ − if isa(nse, 'ao')
+ − nse = matrix(nse,plist('shape',[size(nse,1) size(nse,2)]));
+ − end
+ −
+ − % lighten the model
+ − mdl.clearHistory;
+ − if isa(mdl, 'ssm')
+ − mdl.clearNumParams;
+ − mdl.clearAllUnits;
+ − mdl.params.simplify;
+ − end
+ −
+ − % Check input parameters
+ − if isempty(rng)
+ − error('### Please define a search range ''range''');
+ − end
+ − param = find(pl,'FitParams');
+ − if isempty(param)
+ − error('### Please define the parameters ''param''');
+ − end
+ − nparam = numel(param);
+ − if size(mdl) ~= size(nse)
+ − error('### Parameters ''model'' and ''noise'' must be the same dimension');
+ − end
+ − if size(in) ~= size(out)
+ − error('### Number of input and output experiments must be the same');
+ − end
+ −
+ −
+ − % Get range for parameters
+ − for i = 1:nparam
+ − rang(:,i) = rng{i};
+ − end
+ −
+ − % Get number of experiments
+ − nexp = numel(in);
+ −
+ − switch class(mdl)
+ − case 'matrix'
+ − % Check model sizes
+ − if (isempty(outModel) && isempty(inModel))
+ − if (~(numel(in(1).objs) == mdl.ncols) || ~(numel(out(1).objs) == mdl.nrows))
+ − error('Check model or input/output sizes');
+ − end
+ − elseif isempty(inModel)
+ − if (~(numel(in(1).objs) == mdl.ncols) || ~(size(outModel,2) == mdl.nrows) || ~(size(outModel,1) == numel(out(1).objs)))
+ − error('Check model or input/output sizes');
+ − end
+ − elseif isempty(outModel)
+ − if (~(numel(in(1).objs) == size(inModel,2)) || ~(size(inModel,1) == mdl.ncols) || ~(numel(out(1).objs) == mdl.nrows))
+ − error('Check model or input/output sizes');
+ − end
+ − else
+ − if (~(numel(in(1).objs) == size(inModel,2)) || ~(size(inModel,1) == mdl.ncols) || ~(size(outModel,2) == mdl.nrows) || ~(numel(out(1).objs) == size(outModel,1)))
+ − error('Check model or input/output sizes');
+ − end
+ − end
+ − case 'ssm'
+ − inNames = find(pl,'inNames');
+ − outNames = find(pl,'outNames');
+ − if( (numel(inNames) ~= numel(in(1).objs)) || numel(outNames) ~= numel(out(1).objs))
+ − error('Check model inNames and outNames, they do not match with the input objects')
+ − end
+ − otherwise
+ − error('### Model must be either from the ''matrix'' or the ''ssm'' class. Please check the inputs')
+ − end
+ −
+ −
+ −
+ − %%%%%%%%%%%%%%%%%%%%%%%%%%%% Frequency domain pre-process %%%%%%%%%%%%%%%%%%%%%%%%
+ −
+ − utils.helper.msg(msg.IMPORTANT, 'Computing frequency domain objects', mfilename('class'), mfilename);
+ − for k = 1:nexp
+ −
+ − if((numel(out(1).objs) == 2) && (numel(in(1).objs) == 2)) % previous code done by miquel
+ − % optional resampling (not matrix/resample method yet)
+ − if ~isempty(fsout)
+ − for i = 1:numel(in(k).objs(:))
+ − in(k).objs(i).resample(plist('fsout',fsout));
+ − out(k).objs(i).resample(plist('fsout',fsout));
+ − end
+ − end
+ −
+ − % fft
+ − fin(k) = fft(in(k));
+ − fout(k) = fft(out(k));
+ −
+ − Nsamples = length(fin(k).getObjectAtIndex(1).x);
+ − fs = fin(k).getObjectAtIndex(1).fs;
+ −
+ −
+ − if (isempty(outModel))
+ − % Get rid of fft f =0, reduce frequency range if needed
+ − if ~isempty(flim)
+ − fin(k) = split(fin(k),plist('frequencies',[flim(1) flim(end)]));
+ − fout(k) = split(fout(k),plist('frequencies',[flim(1) flim(end)]));
+ − end
+ − elseif (~isempty(outModel))
+ − % Get rid of fft f =0, reduce frequency range if needed
+ − if ~isempty(flim)
+ − fin(k) = split(fin(k),plist('frequencies',[flim(1) flim(end)]));
+ − fout(k) = split(fout(k),plist('frequencies',[flim(1) flim(end)]));
+ − if(~isempty(outModel))
+ − for lll=1:size(outModel,1)
+ − for kkk=1:size(outModel,2)
+ − outModel(lll,kkk) = split(outModel(lll,kkk),plist('frequencies',[flim(1) flim(end)]));
+ − end
+ − end
+ − end
+ − if(~isempty(inModel))
+ − inModel = split(inModel,plist('frequencies',[flim(1) flim(end)]));
+ − end
+ − end
+ − end
+ −
+ − % use signal fft to get frequency vector. Take into account signal
+ − % could be empty or set to zero
+ −
+ − % 1st channel
+ − if all(fin(k).getObjectAtIndex(1,1).y == 0) || isempty(fin(k).getObjectAtIndex(1,1).y)
+ − i1 = ao(plist('type','fsdata','xvals',0,'yvals',0));
+ − % rebuild input with new zeroed signal
+ − fin(k) = matrix(i1,fin(k).getObjectAtIndex(2,1),plist('shape',[2 1]));
+ − else
+ − freqs = fin(k).getObjectAtIndex(1,1).x;
+ − end
+ − % 2nd channel
+ − if all(fin(k).getObjectAtIndex(2,1).y == 0) || isempty(fin(k).getObjectAtIndex(2,1).y)
+ − i2 = ao(plist('type','fsdata','xvals',0,'yvals',0));
+ − % rebuild input with new zeroed signal
+ − fin(k) = matrix(fin(k).getObjectAtIndex(1,1),i2,plist('shape',[2 1]));
+ − else
+ − freqs = fin(k).getObjectAtIndex(2,1).x;
+ − end
+ −
+ − % Compute psd
+ − n1 = Nsamples*fs/2*psd(nse(k).getObjectAtIndex(1,1), pl);
+ − n2 = Nsamples*fs/2*psd(nse(k).getObjectAtIndex(2,1), pl);
+ − n12 = Nsamples*fs/2*cpsd(nse(k).getObjectAtIndex(1,1),nse(k).getObjectAtIndex(2,1), pl);
+ −
+ − % interpolate noise to fft frequencies
+ − S11 = interp(n1,plist('vertices',freqs));
+ − S12 = interp(n12,plist('vertices',freqs));
+ − S22 = interp(n2,plist('vertices',freqs));
+ − S21 = conj(S12);
+ − % build cross-spectrum matrix
+ − mnse(k) = matrix(S11,S21,S12,S22,plist('shape',[2 2]));
+ −
+ − switch class(mdl)
+ − case 'matrix'
+ − for i = 1:numel(mdl.objs)
+ − if (mdlFreqDependent)
+ − % set Xvals
+ − mdl.objs(i).setXvals(freqs);
+ − % set alias
+ − mdl.objs(i).assignalias(mdl.objs(i),plist('xvals',freqs));
+ − else
+ − mdl.objs(i).setXvals(1);
+ − end
+ − end
+ − case 'ssm'
+ − if k<2
+ − mdl.clearNumParams;
+ − spl = plist('set', 'for bode', ...
+ − 'outputs', outNames, ...
+ − 'inputs', inNames);
+ − % first optimise our model for the case in hand
+ − mdl.reorganize(spl);
+ − % make it lighter
+ −
+ − %for k = 1:nexp
+ −
+ − mdl(k).optimiseForFitting();
+ − %end
+ − end
+ − end
+ − % set Xvals for model (frequency vector could change for each experiment)
+ −
+ − mmdl(k) = mdl;
+ −
+ − elseif ((numel(out(1).objs) == 3) && (numel(in(1).objs) == 4))
+ −
+ − % here we are implementing only the magnetic case
+ − % We have 4 inputs (the 4 conformator waveforms of the magnetic
+ − % analysis and
+ − % 3 outputs (that correspond to the IFO.x12 and IFO.ETA1 and
+ − % IFO.PHI1
+ −
+ −
+ − % fft
+ − fin(k) = fft(in(k));
+ − fout(k) = fft(out(k));
+ − Nsamples = length(fin(k).getObjectAtIndex(1,1).x);
+ − fs = fin(k).getObjectAtIndex(1).fs;
+ −
+ − % Get rid of fft f =0, reduce frequency range if needed
+ − if ~isempty(flim)
+ − fin(k) = split(fin(k),plist('frequencies',[flim(1) flim(end)]));
+ − fout(k) = split(fout(k),plist('frequencies',[flim(1) flim(end)]));
+ −
+ − if(~isempty(outModel))
+ − for lll=1:size(outModel,1)
+ − for kkk=1:size(outModel,2)
+ − outModel(lll,kkk) = split(outModel(lll,kkk),plist('frequencies',[flim(1) flim(end)]));
+ − end
+ − end
+ − end
+ − if(~isempty(inModel))
+ − inModel = split(inModel,plist('frequencies',[flim(1) flim(end)]));
+ − end
+ − end
+ −
+ − freqs = fin(k).getObjectAtIndex(1,1).x;
+ −
+ − % Build noise model
+ − for ii = 1:numel(out.objs)
+ − for jj = ii:numel(out.objs)
+ − % Compute psd
+ − if (ii==jj)
+ − n(ii,jj) = Nsamples*fs/2*psd(nse(k).getObjectAtIndex(ii), pl);
+ − S(ii,jj) = interp(n(ii,jj),plist('vertices',freqs,'method','linear'));
+ − else
+ − n(ii,jj) = Nsamples*fs/2*cpsd(nse(k).getObjectAtIndex(ii),nse(k).getObjectAtIndex(jj),pl);
+ − S(ii,jj) = interp(n(ii,jj),plist('vertices',freqs,'method','linear'));
+ − S(jj,ii) = conj(S(ii,jj));
+ − end
+ − end
+ − end
+ − % build cross-spectrum matrix
+ − mnse(k) = matrix(S,plist('shape',[numel(out.objs) numel(out.objs)]));
+ −
+ − % Check model sizes
+ −
+ − switch class(mdl)
+ − case 'matrix'
+ − % set Xvals for model (frequency vector could change for each experiment)
+ − for i = 1:numel(mdl.objs)
+ − % set Xvals
+ − if (mdlFreqDependent)
+ − mdl.objs(i).setXvals(freqs);
+ − else
+ − mdl.objs(i).setXvals(1);
+ − end
+ − end
+ − case 'ssm'
+ − spl = plist('set', 'for bode', ...
+ − 'outputs', outNames, ...
+ − 'inputs', inNames);
+ − % first optimise our model for the case in hand
+ − mdl.reorganize(spl);
+ − % make it lighter
+ − for k = 1:nexp
+ − mmdl(k).optimiseForFitting();
+ − end
+ − end
+ − mmdl(k) = mdl;
+ −
+ − elseif ((numel(out(1).objs) == 1) && (numel(in(1).objs) == 1))
+ −
+ − % case 1 input, 1 output used mostly for passing
+ − % ao objects from ao.mcmc.
+ −
+ − % optional resampling (not matrix/resample method yet)
+ − if ~isempty(fsout)
+ − in(k).objs(1).resample(plist('fsout',fsout));
+ − out(k).objs(1).resample(plist('fsout',fsout));
+ − end
+ −
+ − % fft
+ − fin(k) = fft(in(k));
+ − fout(k) = fft(out(k));
+ −
+ − Nsamples = length(fin(k).getObjectAtIndex(1).x);
+ − fs = fin(k).getObjectAtIndex(1).fs;
+ −
+ −
+ − if (isempty(outModel))
+ − % Get rid of fft f =0, reduce frequency range if needed
+ − if ~isempty(flim)
+ − fin(k) = split(fin(k),plist('frequencies',[flim(1) flim(end)]));
+ − fout(k) = split(fout(k),plist('frequencies',[flim(1) flim(end)]));
+ − end
+ − elseif (~isempty(outModel))
+ − % Get rid of fft f =0, reduce frequency range if needed
+ − if ~isempty(flim)
+ − fin(k) = split(fin(k),plist('frequencies',[flim(1) flim(end)]));
+ − fout(k) = split(fout(k),plist('frequencies',[flim(1) flim(end)]));
+ − if(~isempty(outModel))
+ − for lll=1:size(outModel,1)
+ − for kkk=1:size(outModel,2)
+ − outModel(lll,kkk) = split(outModel(lll,kkk),plist('frequencies',[flim(1) flim(end)]));
+ − end
+ − end
+ − end
+ − if(~isempty(inModel))
+ − inModel = split(inModel,plist('frequencies',[flim(1) flim(end)]));
+ − end
+ − end
+ − end
+ −
+ − % use signal fft to get frequency vector. Take into account signal
+ − % could be empty or set to zero
+ −
+ − % 1st and only channel
+ − if all(fin(k).getObjectAtIndex(1,1).y == 0) || isempty(fin(k).getObjectAtIndex(1,1).y)
+ − i1 = ao(plist('type','fsdata','xvals',0,'yvals',0));
+ − % rebuild input with new zeroed signal
+ − fin(k) = matrix(i1,fin(k).getObjectAtIndex(2,1),plist('shape',[2 1]));
+ − else
+ − freqs = fin(k).getObjectAtIndex(1,1).x;
+ − end
+ −
+ − % Compute psd
+ − n1 = Nsamples*fs/2*psd(nse(k).getObjectAtIndex(1,1), pl);
+ −
+ − % interpolate noise to fft frequencies
+ − S11 = interp(n1,plist('vertices',freqs));
+ −
+ − % build cross-spectrum matrix
+ − mnse(k) = matrix(S11);
+ −
+ − switch class(mdl)
+ − case 'matrix'
+ − for i = 1:numel(mdl.objs)
+ − if (mdlFreqDependent)
+ − % set Xvals
+ − mdl.objs(i).setXvals(freqs);
+ − % set alias
+ − mdl.objs(i).assignalias(mdl.objs(i),plist('xvals',freqs));
+ − else
+ − mdl.objs(i).setXvals(1);
+ − end
+ − end
+ − case 'ssm'
+ − if k<2
+ − mdl.clearNumParams;
+ − spl = plist('set', 'for bode', ...
+ − 'outputs', outNames, ...
+ − 'inputs', inNames);
+ − % first optimise our model for the case in hand
+ − mdl.reorganize(spl);
+ − % make it lighter
+ − mdl(k).optimiseForFitting();
+ − end
+ − end
+ − % set Xvals for model (frequency vector could change for each experiment)
+ −
+ − mmdl(k) = mdl;
+ −
+ −
+ − else
+ − error('Implemented cases: 1 inputs / 1output, 2 inputs / 2outputs (TN3045 analysis), and 4 inputs / 3 outpus (magnetic complete analysis model. Other cases have not been implemented yet. Sorry for the inconvenience)');
+ − end
+ − end
+ −
+ − %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ −
+ − % do simplex
+ − if simplex
+ − if isempty(xo)
+ − error('### Simplex needs a starting guess. Please input a ''x0''.');
+ − else
+ − switch class(mmdl)
+ − case 'matrix'
+ − xo = fminsearch(@(x) utils.math.loglikelihood_matrix(x,fin,fout,mnse,mmdl,param,inModel,outModel),xo,optimset('Display','iter'));
+ − % case 'smodel'
+ − % xo = fminsearch(@(x) utils.math.loglikelihood(x,in,out,nse,mmdl,param),xo,optimset('Display','iter'));
+ − case 'ssm'
+ − Amats = mdl.amats;
+ − Bmats = mdl.bmats;
+ − Cmats = mdl.cmats;
+ − Dmats = mdl.dmats;
+ − xo = fminsearch(@(x) utils.math.loglikelihood_ssm(x,fin,fout,mnse,mmdl,param,inNames,outNames,spl,Amats,Bmats,Cmats,Dmats),xo,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(param)
+ − fprintf('### Simplex estimate: %s = %d \n',param{i},xo(i))
+ − save('parameters_simplex.txt','xo','-ASCII')
+ − end
+ − end
+ −
+ − p = pest(xo);
+ − p.setName('SimplexEstimate');
+ − p.setNames(param{:});
+ − p.setModels(mmdl);
+ −
+ −
+ − end
+ −
+ − if mhsampleTrue
+ −
+ − % Sampling and saving rejected points. smplr is a matrix of
+ − % the size of (# of samples)x(# of params + 2) wich contains the whole
+ − % chain for each parameter, a column of zeros and ones (rejected or
+ − % accepted point) and the value of the loglikelihood.
+ − [smpl smplr] = utils.math.mhsample(mmdl,fin,fout,mnse,cvar,N,rang,param,Tc,xi,xo,search,jumps,parplot,debug,inNames,outNames,fpars,anneal,SNR0,DeltaL,inModel,outModel);
+ −
+ − % 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,:));
+ −
+ − for i = 1:nparam
+ − for j = 1:nparam
+ − cr(i,j) = cv(i,j)/sqrt(cv(i,i)*cv(j,j));
+ − end
+ − end
+ −
+ − % compute dof
+ − ndata = 0;
+ − for i = 1:length(mtxs)
+ − for j = 1:length(mtxs(1).objs)
+ − ndata = ndata + mtxs(i).objs(j).len;
+ − end
+ − end
+ − dof = ndata-nparam;
+ −
+ − % create pest output
+ − p = pest(mn);
+ − p.setName('mcmc');
+ − p.setNames(param{:});
+ − % add statistical info
+ − p.setCov(cv);
+ − p.setCorr(cr);
+ − p.setDy(sqrt(diag(cv)));
+ − p.setChain(smpl);
+ − %p.computePdf;
+ − p.setDof(dof);
+ − p.setModels(mdl);
+ − % set history
+ − %p.addHistory(getInfo('None'), pl, mtxs_invars(:), [out(:).hist mdl(:).hist]);
+ − end
+ − % Set output
+ − % Implemented it this way in order to leave the pest class as it is.
+ − output = {p,smplr};
+ − %varargout{1} = p;
+ − varargout = output(1:nargout);
+ −
+ − end
+ −
+ − % %--------------------------------------------------------------------------
+ − % % Compute priors means and sigmas
+ − % %--------------------------------------------------------------------------
+ − % function fit = fitPrior(prior,nparam)
+ − % fit = [];
+ − % for ii=1:nparam
+ − % [muhat,sigmahat] = normfit(prior(:,2*ii-1),prior(:,2*ii));
+ − % fit = [fit ,[muhat,sigmahat]'];
+ − % 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, 'ao', 'ltpda', utils.const.categories.sigproc, '$Id: mcmc.m,v 1.35 2011/11/16 15:21:13 nikos 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();
+ −
+ − % % inNames
+ − % p = param({'inNames','Input names. Used for ssm models'}, paramValue.EMPTY_STRING);
+ − % pl.append(p);
+ − %
+ − % % outNames
+ − % p = param({'outNames','Output names. Usde for ssm models'}, paramValue.EMPTY_STRING);
+ − % pl.append(p);
+ − %
+ − % Model
+ − % p = param({'model','A matrix array of models.'}, paramValue.EMPTY_STRING);
+ − % pl.append(p);
+ −
+ − % Param
+ − p = param({'FitParams','A cell array of evaluated parameters.'}, paramValue.EMPTY_DOUBLE);
+ − pl.append(p);
+ −
+ − % Input
+ − p = param({'input','A matrix array of input signals.'}, paramValue.EMPTY_STRING);
+ − pl.append(p);
+ −
+ − pl = plist.MCH_FIT_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);
+ −
+ − % Noise
+ − p = param({'noise','A matrix array of noise spectrum (PSD) used to compute the likelihood.'}, paramValue.EMPTY_STRING);
+ − pl.append(p);
+ −
+ − % Search
+ − p = param({'modelFreqDependent','Set to true to use frequency dependent s models, set to false when using constant models'}, paramValue.TRUE_FALSE);
+ − 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);
+ −
+ − % Frequencies
+ − p = param({'frequencies','Range of frequencies where the analysis is performed. If an array, only first and last are used'}, paramValue.EMPTY_DOUBLE);
+ − pl.append(p);
+ −
+ − % Resample
+ − p = param({'fsout','Desired sampling frequency to resample the input time series'}, paramValue.EMPTY_DOUBLE);
+ − 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);
+ −
+ − % mhsampleTrue
+ − p = param({'mhsample','Set to true to perform a mhsample search. This is set to true by default. Only to be set to false by the user if we does not want to perform the mcmc search'}, 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);
+ −
+ − % heat
+ − 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);
+ −
+ − % heat
+ − 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);
+ −
+ − % inModel
+ − p = param({'inModel','Input model. Still under test'}, paramValue.EMPTY_STRING);
+ − pl.append(p);
+ −
+ − % outModel
+ − p = param({'outModel','Output model. Still under test'}, paramValue.EMPTY_STRING);
+ − pl.append(p);
+ −
+ − p = param({'prior','Mean, sigma and normalization factor for priors. Still under test'}, paramValue.EMPTY_STRING);
+ − pl.append(p);
+ −
+ − p = param({'anneal',['Choose type of annealing during sampling. Default value is ',...
+ − 'simulated annealing. Choose "thermo" for annealing with a thermostat.',...
+ − ' SNR is computed and if it is larger than a fixed value SNR0 (provided also in the plist), ',...
+ − 'then the chains are heated by a factor of (SNR(1)/SNR0)^2. Choosing "simple" ',...
+ − 'the deviation of the loglikelihood of every 10 points in the chains is stored. If this deviation ',...
+ − 'is larger or smaller than two fixed values the chains are cooled or heated respectively.']}, {1, {'simul','thermo', 'simple'}, paramValue.SINGLE});
+ − pl.append(p);
+ −
+ − p = param({'SNR0','Fixed value for thermostated annealing.'}, {1, {200}, paramValue.OPTIONAL});
+ − pl.append(p);
+ −
+ − p = param({'DeltaL',['Deviation of Loglikelihood for 10 points of the chains. Used for',...
+ − 'the "simple" choice of annealing with a thermostat.']}, {1, {[100 600 2 3]}, paramValue.OPTIONAL});
+ − pl.append(p);
+ −
+ − end
+ −
+ −
+ − % PARAMETERS: J - number of chains to compute.
+ − % sigma - dispersion of the jumping distribution.
+ − % range - range where the parameteters are sampled.
+ − % param - a cell array of evaluated parameters (from the
+ − % smodel).
+ − % Tc - an array with two values (default: 0, i.e. no annealing).
+ − % heat - heat index flattening likelihood surface (default: 1)
+ − % x0 - proposed initial values, if empty selected randomly.
+ − % (default: empty)