Mercurial > hg > ltpda

comparison m-toolbox/classes/@matrix/modelselect.m @ 0:f0afece42f48

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