0
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 1 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 2 %
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 3 % modelselect.m - method to compute the Bayes Factor using
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 4 % reverse jump mcmc algorithm
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 5 %
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 6 % call - Bxy = modelselect(out,pl)
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 7 %
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 8 % inputs - out: matrix objects with measured outputs
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 9 % pl: parameter list
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 10 %
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 11 % outputs - Bxy: an array of AOs containing the evolution
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 12 % of the Bayes factors. (comparing each model
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 13 % with each other)
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 14 %
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 15 %
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 16 % <a href="matlab:utils.helper.displayMethodInfo('matrix','modelselect')">Parameters Description</a>
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 17 %
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 18 % N. Karnesis 27/09/2011
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 19 %
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 20 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 21
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 22 function varargout = modelselect(varargin)
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 23
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 24 % Check if this is a call for parameters
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 25 if utils.helper.isinfocall(varargin{:})
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 26 varargout{1} = getInfo(varargin{3});
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 27 return
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 28 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 29
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 30 import utils.const.*
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 31 utils.helper.msg(msg.PROC3, 'running %s/%s', mfilename('class'), mfilename);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 32
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 33 % Method can not be used as a modifier
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 34 if nargout == 0
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 35 error('### mcmc cannot be used as a modifier. Please give an output variable.');
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 36 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 37
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 38 % Collect input variable names
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 39 in_names = cell(size(varargin));
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 40 for ii = 1:nargin,in_names{ii} = inputname(ii);end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 41
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 42 % Collect all AOs smodels and plists
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 43 [mtxs, mtxs_invars] = utils.helper.collect_objects(varargin(:), 'matrix', in_names);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 44 pl = utils.helper.collect_objects(varargin(:), 'plist', in_names);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 45
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 46 % Combine plists
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 47 pl = parse(pl, getDefaultPlist);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 48
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 49 % Get parameters
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 50 N = find(pl, 'N');
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 51 Tc = find(pl,'Tc');
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 52 xi = find(pl,'heat');
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 53 param = find(pl,'FitParams');
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 54 cvar = find(pl, 'cov');
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 55 rng = find(pl,'range');
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 56 flim = find(pl,'frequencies');
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 57 fsout = find(pl,'fsout');
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 58 search = find(pl,'search');
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 59 xo = find(pl,'x0');
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 60 mtxin = find(pl,'input');
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 61 mdlin = find(pl,'models');
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 62 mdlFreqDependent = find(pl,'modelFreqDependent');
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 63 nse = find(pl,'noise');
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 64 jumps = find(pl,'jumps');
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 65 parplot = find(pl,'plot');
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 66 debug = find(pl,'debug');
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 67 outModel = find(pl,'outModel');
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 68 inModel = find(pl,'inModel1');
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 69 outNames = find(pl,'outNames');
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 70 inNames = find(pl,'inNames');
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 71 anneal = find(pl,'anneal');
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 72 DeltaL = find(pl,'DeltaL');
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 73 SNR0 = find(pl,('SNR0'));
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 74
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 75 nummodels = numel(mdlin);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 76
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 77 % Decide on a deep copy or a modify
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 78 in = copy(mtxin, nargout);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 79 out = copy(mtxs, nargout);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 80 for k = 1:nummodels
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 81 mdl(k) = copy(mdlin{1,k}, nargout);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 82 nparam(k) = numel(param{1,k});
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 83 % lighten the model
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 84 mdl(k).clearHistory;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 85 if isa(mdl(k), 'ssm')
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 86 mdl(k).clearNumParams;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 87 mdl(k).clearAllUnits;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 88 mdl(k).params.simplify;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 89 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 90 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 91
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 92
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 93
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 94 % Check input parameters
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 95 if isempty(rng)
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 96 error('### Please define a search range ''range''');
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 97 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 98 if isempty(param)
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 99 error('### Please define the parameters ''param''');
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 100 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 101 if size(mdl(1)) ~= size(nse)
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 102 error('### Parameters ''model'' and ''noise'' must be the same dimension');
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 103 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 104 if size(in) ~= size(out)
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 105 error('### Number of input and output experiments must be the same');
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 106 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 107
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 108
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 109 % Get range for parameters
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 110 for jj = 1:nummodels
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 111 for i = 1:numel(param{1,jj})
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 112 rang{1,jj}(:,i) = rng{1,jj}{i}';
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 113 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 114 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 115 % Get number of experiments
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 116 nexp = numel(in);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 117
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 118 switch class(mdl)
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 119 case 'matrix'
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 120 % Check model sizes
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 121 for k = 1:nummodels
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 122 if (isempty(outModel) && isempty(inModel))
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 123 if (~(numel(in(1).objs) == mdl(k).ncols) || ~(numel(out(1).objs) == mdl(k).nrows))
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 124 error('Check model or input/output sizes');
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 125 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 126 elseif isempty(inModel)
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 127 if (~(numel(in(1).objs) == mdl(k).ncols) || ~(size(outModel,2) == mdl(k).nrows) || ~(size(outModel,1) == numel(out(1).objs)))
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 128 error('Check model or input/output sizes');
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 129 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 130 elseif isempty(outModel)
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 131 if (~(numel(in(1).objs) == size(inModel,2)) || ~(size(inModel,1) == mdl(k).ncols) || ~(numel(out(1).objs) == mdl(k).nrows))
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 132 error('Check model or input/output sizes');
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 133 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 134 else
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 135 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)))
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 136 error('Check model or input/output sizes');
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 137 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 138 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 139 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 140 case 'ssm'
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 141 inNames = find(pl,'inNames');
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 142 outNames = find(pl,'outNames');
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 143 if( (numel(inNames) ~= numel(in(1).objs)) || numel(outNames) ~= numel(out(1).objs))
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 144 error('Check model inNames and outNames, they do not match with the input objects')
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 145 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 146 otherwise
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 147 error('### Model must be either from the ''matrix'' or the ''ssm'' class. Please check the inputs')
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 148 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 149
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 150
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 151
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 152 %%%%%%%%%%%%%%%%%%%%%%%%%%%% Frequency domain pre-process %%%%%%%%%%%%%%%%%%%%%%%%
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 153
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 154 utils.helper.msg(msg.IMPORTANT, 'Computing frequency domain objects', mfilename('class'), mfilename);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 155 for k = 1:nexp
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 156
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 157 if((numel(out(1).objs) == 2) && (numel(in(1).objs) == 2)) % previous code done by miquel
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 158 % optional resampling (not matrix/resample method yet)
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 159 if ~isempty(fsout)
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 160 for i = 1:numel(in(k).objs(:))
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 161 in(k).objs(i).resample(plist('fsout',fsout));
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 162 out(k).objs(i).resample(plist('fsout',fsout));
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 163 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 164 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 165
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 166 % fft
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 167 fin(k) = fft(in(k));
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 168 fout(k) = fft(out(k));
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 169
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 170 Nsamples = length(fin(k).getObjectAtIndex(1).x);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 171 fs = fin(k).getObjectAtIndex(1).fs;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 172
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 173
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 174 if (isempty(outModel))
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 175 % Get rid of fft f =0, reduce frequency range if needed
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 176 if ~isempty(flim)
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 177 fin(k) = split(fin(k),plist('frequencies',[flim(1) flim(end)]));
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 178 fout(k) = split(fout(k),plist('frequencies',[flim(1) flim(end)]));
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 179 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 180 elseif (~isempty(outModel))
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 181 % Get rid of fft f =0, reduce frequency range if needed
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 182 if ~isempty(flim)
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 183 fin(k) = split(fin(k),plist('frequencies',[flim(1) flim(end)]));
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 184 fout(k) = split(fout(k),plist('frequencies',[flim(1) flim(end)]));
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 185 if(~isempty(outModel))
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 186 for lll=1:size(outModel,1)
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 187 for kkk=1:size(outModel,2)
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 188 outModel(lll,kkk) = split(outModel(lll,kkk),plist('frequencies',[flim(1) flim(end)]));
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 189 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 190 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 191 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 192 if(~isempty(inModel))
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 193 inModel = split(inModel,plist('frequencies',[flim(1) flim(end)]));
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 194 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 195 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 196 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 197
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 198 % use signal fft to get frequency vector. Take into account signal
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 199 % could be empty or set to zero
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 200
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 201 % 1st channel
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 202 if all(fin(k).getObjectAtIndex(1,1).y == 0) || isempty(fin(k).getObjectAtIndex(1,1).y)
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 203 i1 = ao(plist('type','fsdata','xvals',0,'yvals',0));
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 204 % rebuild input with new zeroed signal
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 205 fin(k) = matrix(i1,fin(k).getObjectAtIndex(2,1),plist('shape',[2 1]));
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 206 else
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 207 freqs = fin(k).getObjectAtIndex(1,1).x;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 208 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 209 % 2nd channel
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 210 if all(fin(k).getObjectAtIndex(2,1).y == 0) || isempty(fin(k).getObjectAtIndex(2,1).y)
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 211 i2 = ao(plist('type','fsdata','xvals',0,'yvals',0));
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 212 % rebuild input with new zeroed signal
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 213 fin(k) = matrix(fin(k).getObjectAtIndex(1,1),i2,plist('shape',[2 1]));
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 214 else
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 215 freqs = fin(k).getObjectAtIndex(2,1).x;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 216 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 217
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 218 % Compute psd
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 219 n1 = Nsamples*fs/2*psd(nse(k).getObjectAtIndex(1,1), pl);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 220 n2 = Nsamples*fs/2*psd(nse(k).getObjectAtIndex(2,1), pl);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 221 n12 = Nsamples*fs/2*cpsd(nse(k).getObjectAtIndex(1,1),nse(k).getObjectAtIndex(2,1), pl);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 222
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 223 % interpolate noise to fft frequencies
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 224 S11 = interp(n1,plist('vertices',freqs));
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 225 S12 = interp(n12,plist('vertices',freqs));
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 226 S22 = interp(n2,plist('vertices',freqs));
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 227 S21 = conj(S12);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 228 % build cross-spectrum matrix
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 229 mnse(k) = matrix(S11,S21,S12,S22,plist('shape',[2 2]));
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 230
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 231 elseif((numel(out(1).objs) == 1) && (numel(in(1).objs) == 1))
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 232 % optional resampling (not matrix/resample method yet)
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 233 if ~isempty(fsout)
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 234 for i = 1:numel(in(k).objs(:))
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 235 in(k).objs(i).resample(plist('fsout',fsout));
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 236 out(k).objs(i).resample(plist('fsout',fsout));
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 237 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 238 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 239
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 240 % fft
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 241 fin(k) = fft(in(k));
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 242 fout(k) = fft(out(k));
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 243
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 244 Nsamples = length(fin(k).getObjectAtIndex(1).x);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 245 fs = fin(k).getObjectAtIndex(1).fs;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 246
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 247
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 248 if (isempty(outModel))
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 249 % Get rid of fft f =0, reduce frequency range if needed
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 250 if ~isempty(flim)
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 251 fin(k) = split(fin(k),plist('frequencies',[flim(1) flim(end)]));
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 252 fout(k) = split(fout(k),plist('frequencies',[flim(1) flim(end)]));
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 253 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 254 elseif (~isempty(outModel))
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 255 % Get rid of fft f =0, reduce frequency range if needed
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 256 if ~isempty(flim)
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 257 fin(k) = split(fin(k),plist('frequencies',[flim(1) flim(end)]));
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 258 fout(k) = split(fout(k),plist('frequencies',[flim(1) flim(end)]));
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 259 if(~isempty(outModel))
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 260 for lll=1:size(outModel,1)
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 261 for kkk=1:size(outModel,2)
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 262 outModel(lll,kkk) = split(outModel(lll,kkk),plist('frequencies',[flim(1) flim(end)]));
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 263 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 264 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 265 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 266 if(~isempty(inModel))
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 267 inModel = split(inModel,plist('frequencies',[flim(1) flim(end)]));
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 268 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 269 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 270 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 271
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 272 % use signal fft to get frequency vector. Take into account signal
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 273 % could be empty or set to zero
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 274
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 275 % 1st and only channel
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 276 if all(fin(k).getObjectAtIndex(1,1).y == 0) || isempty(fin(k).getObjectAtIndex(1,1).y)
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 277 i1 = ao(plist('type','fsdata','xvals',0,'yvals',0));
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 278 % rebuild input with new zeroed signal
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 279 fin(k) = matrix(i1,fin(k).getObjectAtIndex(2,1),plist('shape',[2 1]));
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 280 else
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 281 freqs = fin(k).getObjectAtIndex(1,1).x;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 282 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 283
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 284 % Compute psd
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 285 n1 = Nsamples*fs/2*psd(nse(k).getObjectAtIndex(1,1), pl);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 286
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 287 % interpolate noise to fft frequencies
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 288 S11 = interp(n1,plist('vertices',freqs));
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 289
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 290 % build cross-spectrum matrix
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 291 mnse(k) = matrix(S11);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 292
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 293 else
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 294 error('### Sorry, implemented only for 2 inputs - 2 outputs and 1 input - 1 output only...')
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 295 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 296
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 297 switch class(mdl)
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 298
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 299 case 'matrix'
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 300 for jj = 1:nummodels
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 301 for i = 1:numel(mdl(jj).objs)
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 302 if (mdlFreqDependent)
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 303 % set Xvals
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 304 mdl(jj).objs(i).setXvals(freqs);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 305 % set alias
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 306 mdl(jj).objs(i).assignalias(mdl(jj).objs(i),plist('xvals',freqs));
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 307 else
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 308 mdl(jj).objs(i).setXvals(1);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 309 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 310 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 311 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 312 case 'ssm'
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 313 if k < 2
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 314 for jj = 1:nummodels
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 315 mdl(jj).clearNumParams;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 316 spl = plist('set', 'for bode', ...
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 317 'outputs', outNames, ...
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 318 'inputs', inNames);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 319 % first optimise our model for the case in hand
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 320 mdl(jj).reorganize(spl);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 321 % make it lighter
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 322 mdl(jj).optimiseForFitting();
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 323 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 324 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 325 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 326 % add lines to the matrix of models. Each line corresponds to one
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 327 % experiment
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 328 mmdl(k,:) = mdl;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 329 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 330
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 331 [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);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 332
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 333 % Set output
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 334 numfactors = 0;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 335 for gg = 1:(nummodels)
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 336 for dd = gg+1:(nummodels)
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 337 numfactors = numfactors + 1;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 338 BayesFactor(numfactors) = ao(Bxy{gg,dd});
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 339 BayesFactor(numfactors).setName(sprintf('B%d%d',gg,dd));
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 340 BayesFactor(numfactors).procinfo = plist('Likelihoods',LogLambda);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 341 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 342 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 343
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 344 output = {BayesFactor};
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 345 varargout = output(1:nargout);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 346
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 347
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 348 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 349
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 350
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 351 %--------------------------------------------------------------------------
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 352 % Get Info Object
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 353 %--------------------------------------------------------------------------
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 354 function ii = getInfo(varargin)
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 355 if nargin == 1 && strcmpi(varargin{1}, 'None')
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 356 sets = {};
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 357 pl = [];
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 358 else
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 359 sets = {'Default'};
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 360 pl = getDefaultPlist;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 361 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 362 % Build info object
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 363 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);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 364 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 365
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 366 %--------------------------------------------------------------------------
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 367 % Get Default Plist
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 368 %--------------------------------------------------------------------------
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 369 function plout = getDefaultPlist()
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 370 persistent pl;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 371 if exist('pl', 'var')==0 || isempty(pl)
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 372 pl = buildplist();
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 373 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 374 plout = pl;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 375 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 376
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 377 function pl = buildplist()
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 378 pl = plist();
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 379
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 380 % inNames
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 381 p = param({'inNames','Input names. Used for ssm models'}, paramValue.EMPTY_STRING);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 382 pl.append(p);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 383
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 384 % outNames
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 385 p = param({'outNames','Output names. Usde for ssm models'}, paramValue.EMPTY_STRING);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 386 pl.append(p);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 387
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 388 %Model
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 389 p = param({'model','An array input of models.'}, paramValue.EMPTY_STRING);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 390 pl.append(p);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 391
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 392 % Param
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 393 p = param({'FitParams','A cell array of evaluated parameters for each model.'}, paramValue.EMPTY_DOUBLE);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 394 pl.append(p);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 395
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 396 % Input
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 397 p = param({'input','A matrix array of input signals.'}, paramValue.EMPTY_STRING);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 398 pl.append(p);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 399
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 400 pl = plist.MCH_FIT_PLIST;
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 401
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 402 % N
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 403 p = param({'N','number of samples of the chain.'}, paramValue.DOUBLE_VALUE(1000));
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 404 pl.append(p);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 405
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 406 % Sigma
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 407 p = param({'cov','Cell array containing the covariances of the gaussian jumping distribution for each model.'}, paramValue.DOUBLE_VALUE(1e-4));
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 408 pl.append(p);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 409
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 410 % Noise
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 411 p = param({'noise','A matrix array of noise spectrum (PSD) used to compute the likelihood.'}, paramValue.EMPTY_STRING);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 412 pl.append(p);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 413
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 414 % Search
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 415 p = param({'modelFreqDependent','Set to true to use frequency dependent s models, set to false when using constant models'}, paramValue.TRUE_FALSE);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 416 pl.append(p);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 417
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 418 % Search
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 419 p = param({'search','Set to true to use bigger jumps in parameter space during annealing and cool down.'}, paramValue.TRUE_FALSE);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 420 pl.append(p);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 421
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 422 % Frequencies
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 423 p = param({'frequencies','Range of frequencies where the analysis is performed. If an array, only first and last are used'}, paramValue.EMPTY_DOUBLE);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 424 pl.append(p);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 425
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 426 % Resample
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 427 p = param({'fsout','Desired sampling frequency to resample the input time series'}, paramValue.EMPTY_DOUBLE);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 428 pl.append(p);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 429
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 430 % heat
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 431 p = param({'heat','The heat index flattening likelihood surface during annealing.'}, paramValue.DOUBLE_VALUE(1));
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 432 pl.append(p);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 433
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 434 % Tc
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 435 p = param({'Tc','An array of two values setting the initial and final value for the cooling down.'}, paramValue.EMPTY_STRING);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 436 pl.append(p);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 437
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 438 % initial values
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 439 p = param({'x0','The proposed initial values (cell array again).'}, paramValue.EMPTY_DOUBLE);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 440 pl.append(p);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 441
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 442 % ranges
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 443 p = param({'range','ranges'}, paramValue.EMPTY_DOUBLE);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 444 pl.append(p);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 445
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 446 % jumps
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 447 p = param({'jumps','An array of four numbers setting the rescaling of the covariance matrix during the search phase.',...
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 448 'The first value is the one applied by default, the following thhree apply just when the chain sample is',...
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 449 'mod(10), mod(25) and mod(100) respectively.'}, paramValue.EMPTY_DOUBLE);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 450 pl.append(p);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 451
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 452 % heat
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 453 p = param({'plot','If set equal to one, the evolution of the Bayes factor is plotted every 500 steps.'}, paramValue.EMPTY_DOUBLE);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 454 pl.append(p);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 455
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 456 % debug
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 457 p = param({'debug','Set to true to get debug information of the MCMC process.'}, paramValue.FALSE_TRUE);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 458 pl.append(p);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 459
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 460 % outModel
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 461 p = param({'outModel','Output model. Still under test'}, paramValue.EMPTY_STRING);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 462 pl.append(p);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 463
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 464 p = param({'anneal',['Choose type of annealing during sampling. Default value is ',...
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 465 'simulated annealing. Choose "thermo" for annealing with a thermostat.',...
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 466 ' SNR is computed and if it is larger than a fixed value SNR0 (provided also in the plist), ',...
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 467 'then the chains are heated by a factor of (SNR(1)/SNR0)^2. Choosing "simple" ',...
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 468 'the deviation of the loglikelihood of every 10 points in the chains is stored. If this deviation ',...
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 469 'is larger or smaller than two fixed values the chains are cooled or heated respectively.']}, {1, {'simul','thermo', 'simple'}, paramValue.SINGLE});
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 470 pl.append(p);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 471
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 472 p = param({'SNR0','Fixed value for thermostated annealing.'}, {1, {200}, paramValue.OPTIONAL});
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 473 pl.append(p);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 474
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 475 p = param({'DeltaL',['Deviation of Loglikelihood for 10 points of the chains. Used for',...
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 476 'the "simple" choice of annealing with a thermostat.']}, {1, {[100 600 2 3]}, paramValue.OPTIONAL});
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 477 pl.append(p);
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 478
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 479
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 480 end
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 481
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
+ − 482