line source
+ − %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ − %
+ − % modelselect.m - method to compute the Bayes Factor using
+ − % reverse jump mcmc algorithm
+ − %
+ − % call - Bxy = modelselect(out,pl)
+ − %
+ − % inputs - out: matrix objects with measured outputs
+ − % pl: parameter list
+ − %
+ − % outputs - Bxy: an array of AOs containing the evolution
+ − % of the Bayes factors. (comparing each model
+ − % with each other)
+ − %
+ − %
+ − % <a href="matlab:utils.helper.displayMethodInfo('matrix','modelselect')">Parameters Description</a>
+ − %
+ − % N. Karnesis 27/09/2011
+ − %
+ − %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ −
+ − function varargout = modelselect(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');
+ − param = find(pl,'FitParams');
+ − cvar = find(pl, 'cov');
+ − rng = find(pl,'range');
+ − flim = find(pl,'frequencies');
+ − fsout = find(pl,'fsout');
+ − search = find(pl,'search');
+ − xo = find(pl,'x0');
+ − mtxin = find(pl,'input');
+ − mdlin = find(pl,'models');
+ − 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,'inModel1');
+ − outNames = find(pl,'outNames');
+ − inNames = find(pl,'inNames');
+ − anneal = find(pl,'anneal');
+ − DeltaL = find(pl,'DeltaL');
+ − SNR0 = find(pl,('SNR0'));
+ −
+ − nummodels = numel(mdlin);
+ −
+ − % Decide on a deep copy or a modify
+ − in = copy(mtxin, nargout);
+ − out = copy(mtxs, nargout);
+ − for k = 1:nummodels
+ − mdl(k) = copy(mdlin{1,k}, nargout);
+ − nparam(k) = numel(param{1,k});
+ − % lighten the model
+ − mdl(k).clearHistory;
+ − if isa(mdl(k), 'ssm')
+ − mdl(k).clearNumParams;
+ − mdl(k).clearAllUnits;
+ − mdl(k).params.simplify;
+ − end
+ − end
+ −
+ −
+ −
+ − % Check input parameters
+ − if isempty(rng)
+ − error('### Please define a search range ''range''');
+ − end
+ − if isempty(param)
+ − error('### Please define the parameters ''param''');
+ − end
+ − if size(mdl(1)) ~= 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 jj = 1:nummodels
+ − for i = 1:numel(param{1,jj})
+ − rang{1,jj}(:,i) = rng{1,jj}{i}';
+ − end
+ − end
+ − % Get number of experiments
+ − nexp = numel(in);
+ −
+ − switch class(mdl)
+ − case 'matrix'
+ − % Check model sizes
+ − for k = 1:nummodels
+ − if (isempty(outModel) && isempty(inModel))
+ − if (~(numel(in(1).objs) == mdl(k).ncols) || ~(numel(out(1).objs) == mdl(k).nrows))
+ − error('Check model or input/output sizes');
+ − end
+ − elseif isempty(inModel)
+ − if (~(numel(in(1).objs) == mdl(k).ncols) || ~(size(outModel,2) == mdl(k).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(k).ncols) || ~(numel(out(1).objs) == mdl(k).nrows))
+ − error('Check model or input/output sizes');
+ − end
+ − else
+ − if (~(numel(in(1).objs) == size(inModel,2)) || ~(size(inModel,1) == mdl(k).ncols) || ~(size(outModel,2) == mdl(k).nrows) || ~(numel(out(1).objs) == size(outModel,1)))
+ − error('Check model or input/output sizes');
+ − end
+ − 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]));
+ −
+ − elseif((numel(out(1).objs) == 1) && (numel(in(1).objs) == 1))
+ − % 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 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);
+ −
+ − else
+ − error('### Sorry, implemented only for 2 inputs - 2 outputs and 1 input - 1 output only...')
+ − end
+ −
+ − switch class(mdl)
+ −
+ − case 'matrix'
+ − for jj = 1:nummodels
+ − for i = 1:numel(mdl(jj).objs)
+ − if (mdlFreqDependent)
+ − % set Xvals
+ − mdl(jj).objs(i).setXvals(freqs);
+ − % set alias
+ − mdl(jj).objs(i).assignalias(mdl(jj).objs(i),plist('xvals',freqs));
+ − else
+ − mdl(jj).objs(i).setXvals(1);
+ − end
+ − end
+ − end
+ − case 'ssm'
+ − if k < 2
+ − for jj = 1:nummodels
+ − mdl(jj).clearNumParams;
+ − spl = plist('set', 'for bode', ...
+ − 'outputs', outNames, ...
+ − 'inputs', inNames);
+ − % first optimise our model for the case in hand
+ − mdl(jj).reorganize(spl);
+ − % make it lighter
+ − mdl(jj).optimiseForFitting();
+ − end
+ − end
+ − end
+ − % add lines to the matrix of models. Each line corresponds to one
+ − % experiment
+ − mmdl(k,:) = mdl;
+ − end
+ −
+ − [Bxy LogLambda] = utils.math.rjsample(mmdl,fin,fout,mnse,cvar,N,rang,param,Tc,xi,xo,search,jumps,parplot,debug,inNames,outNames,anneal,SNR0,DeltaL,inModel,outModel);
+ −
+ − % Set output
+ − numfactors = 0;
+ − for gg = 1:(nummodels)
+ − for dd = gg+1:(nummodels)
+ − numfactors = numfactors + 1;
+ − BayesFactor(numfactors) = ao(Bxy{gg,dd});
+ − BayesFactor(numfactors).setName(sprintf('B%d%d',gg,dd));
+ − BayesFactor(numfactors).procinfo = plist('Likelihoods',LogLambda);
+ − end
+ − end
+ −
+ − output = {BayesFactor};
+ − varargout = output(1:nargout);
+ −
+ −
+ − 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.33 2011/07/31 17:38:34 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','An array input of models.'}, paramValue.EMPTY_STRING);
+ − pl.append(p);
+ −
+ − % Param
+ − p = param({'FitParams','A cell array of evaluated parameters for each model.'}, 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','Cell array containing the covariances of the gaussian jumping distribution for each model.'}, 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);
+ −
+ − % 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);
+ −
+ − % initial values
+ − p = param({'x0','The proposed initial values (cell array again).'}, paramValue.EMPTY_DOUBLE);
+ − pl.append(p);
+ −
+ − % ranges
+ − p = param({'range','ranges'}, 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','If set equal to one, the evolution of the Bayes factor is plotted every 500 steps.'}, 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);
+ −
+ − % outModel
+ − p = param({'outModel','Output model. 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
+ −
+ −