Mercurial > hg > ltpda
view m-toolbox/classes/@matrix/mcmc.m @ 28:01b86b780ba7 database-connection-manager
Remove LTPDARepositoryManager implementation. Java code
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 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)