Mercurial > hg > ltpda

comparison m-toolbox/classes/@matrix/mcmc.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 % MCMC estimates paramters using a Monte Carlo Markov Chain.
2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3 %
4 % DESCRIPTION: MCMC estimate the parameters of a given model given
5 % inputs, outputs and noise using a Metropolis algorithm.
6 %
7 % CALL: [b smplr] = mcmc(out,pl)
8 %
9 % INPUTS: out - matrix objects with measured outputs
10 % pl - parameter list
11 %
12 % OUTPUTS: b - pest object contatining estimate information
13 % smplr - matrix containing info about the rejected points
14 %
15 % <a href="matlab:utils.helper.displayMethodInfo('matrix', 'mcmc')">Parameters Description</a>
16 %
17 % VERSION: $Id: mcmc.m,v 1.35 2011/11/16 15:21:13 nikos Exp $
18 %
19 % References: M Nofrarias et al. Phys. Rev. D 82, 122002 (2010)
20 %
21 %
22 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
23
24 % TODO: multiple chain option not implemented yet
25 % metropolis/hastings not implemented
26 % empty initial values not implemented
27
28
29 function varargout = mcmc(varargin)
30
31 % Check if this is a call for parameters
32 if utils.helper.isinfocall(varargin{:})
33 varargout{1} = getInfo(varargin{3});
34 return
35 end
36
37 import utils.const.*
38 utils.helper.msg(msg.PROC3, 'running %s/%s', mfilename('class'), mfilename);
39
40 % Method can not be used as a modifier
41 if nargout == 0
42 error('### mcmc cannot be used as a modifier. Please give an output variable.');
43 end
44
45 % Collect input variable names
46 in_names = cell(size(varargin));
47 for ii = 1:nargin,in_names{ii} = inputname(ii);end
48
49 % Collect all AOs smodels and plists
50 [mtxs, mtxs_invars] = utils.helper.collect_objects(varargin(:), 'matrix', in_names);
51 pl = utils.helper.collect_objects(varargin(:), 'plist', in_names);
52
53 % Combine plists
54 pl = parse(pl, getDefaultPlist);
55
56 % Get parameters
57 N = find(pl, 'N');
58 Tc = find(pl,'Tc');
59 xi = find(pl,'heat');
60 cvar = find(pl, 'cov');
61 rng = find(pl,'range');
62 flim = find(pl,'frequencies');
63 fsout = find(pl,'fsout');
64 search = find(pl,'search');
65 simplex = find(pl,'simplex');
66 mhsampleTrue = find(pl,'mhsample');
67 xo = find(pl,'x0');
68 mtxin = find(pl,'input');
69 mdlin = find(pl,'model');
70 mdlFreqDependent = find(pl,'modelFreqDependent');
71 nse = find(pl,'noise');
72 jumps = find(pl,'jumps');
73 parplot = find(pl,'plot');
74 debug = find(pl,'debug');
75 outModel = find(pl,'outModel');
76 inModel = find(pl,'inModel');
77 outNames = find(pl,'outNames');
78 inNames = find(pl,'inNames');
79 fpars = find(pl,'prior');
80 anneal = find(pl,'anneal');
81 DeltaL = find(pl,'DeltaL');
82 SNR0 = find(pl,('SNR0'));
83
84 % Decide on a deep copy or a modify
85 in = copy(mtxin, nargout);
86 out = copy(mtxs, nargout);
87 mdl = copy(mdlin, nargout);
88
89 % checking if input or noise are aos
90 if isa(in, 'ao')
91 in = matrix(in,plist('shape',[size(in,1) size(in,2)]));
92 end
93 if isa(nse, 'ao')
94 nse = matrix(nse,plist('shape',[size(nse,1) size(nse,2)]));
95 end
96
97 % lighten the model
98 mdl.clearHistory;
99 if isa(mdl, 'ssm')
100 mdl.clearNumParams;
101 mdl.clearAllUnits;
102 mdl.params.simplify;
103 end
104
105 % Check input parameters
106 if isempty(rng)
107 error('### Please define a search range ''range''');
108 end
109 param = find(pl,'FitParams');
110 if isempty(param)
111 error('### Please define the parameters ''param''');
112 end
113 nparam = numel(param);
114 if size(mdl) ~= size(nse)
115 error('### Parameters ''model'' and ''noise'' must be the same dimension');
116 end
117 if size(in) ~= size(out)
118 error('### Number of input and output experiments must be the same');
119 end
120
121
122 % Get range for parameters
123 for i = 1:nparam
124 rang(:,i) = rng{i};
125 end
126
127 % Get number of experiments
128 nexp = numel(in);
129
130 switch class(mdl)
131 case 'matrix'
132 % Check model sizes
133 if (isempty(outModel) && isempty(inModel))
134 if (~(numel(in(1).objs) == mdl.ncols) || ~(numel(out(1).objs) == mdl.nrows))
135 error('Check model or input/output sizes');
136 end
137 elseif isempty(inModel)
138 if (~(numel(in(1).objs) == mdl.ncols) || ~(size(outModel,2) == mdl.nrows) || ~(size(outModel,1) == numel(out(1).objs)))
139 error('Check model or input/output sizes');
140 end
141 elseif isempty(outModel)
142 if (~(numel(in(1).objs) == size(inModel,2)) || ~(size(inModel,1) == mdl.ncols) || ~(numel(out(1).objs) == mdl.nrows))
143 error('Check model or input/output sizes');
144 end
145 else
146 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)))
147 error('Check model or input/output sizes');
148 end
149 end
150 case 'ssm'
151 inNames = find(pl,'inNames');
152 outNames = find(pl,'outNames');
153 if( (numel(inNames) ~= numel(in(1).objs)) || numel(outNames) ~= numel(out(1).objs))
154 error('Check model inNames and outNames, they do not match with the input objects')
155 end
156 otherwise
157 error('### Model must be either from the ''matrix'' or the ''ssm'' class. Please check the inputs')
158 end
159
160
161
162 %%%%%%%%%%%%%%%%%%%%%%%%%%%% Frequency domain pre-process %%%%%%%%%%%%%%%%%%%%%%%%
163
164 utils.helper.msg(msg.IMPORTANT, 'Computing frequency domain objects', mfilename('class'), mfilename);
165 for k = 1:nexp
166
167 if((numel(out(1).objs) == 2) && (numel(in(1).objs) == 2)) % previous code done by miquel
168 % optional resampling (not matrix/resample method yet)
169 if ~isempty(fsout)
170 for i = 1:numel(in(k).objs(:))
171 in(k).objs(i).resample(plist('fsout',fsout));
172 out(k).objs(i).resample(plist('fsout',fsout));
173 end
174 end
175
176 % fft
177 fin(k) = fft(in(k));
178 fout(k) = fft(out(k));
179
180 Nsamples = length(fin(k).getObjectAtIndex(1).x);
181 fs = fin(k).getObjectAtIndex(1).fs;
182
183
184 if (isempty(outModel))
185 % Get rid of fft f =0, reduce frequency range if needed
186 if ~isempty(flim)
187 fin(k) = split(fin(k),plist('frequencies',[flim(1) flim(end)]));
188 fout(k) = split(fout(k),plist('frequencies',[flim(1) flim(end)]));
189 end
190 elseif (~isempty(outModel))
191 % Get rid of fft f =0, reduce frequency range if needed
192 if ~isempty(flim)
193 fin(k) = split(fin(k),plist('frequencies',[flim(1) flim(end)]));
194 fout(k) = split(fout(k),plist('frequencies',[flim(1) flim(end)]));
195 if(~isempty(outModel))
196 for lll=1:size(outModel,1)
197 for kkk=1:size(outModel,2)
198 outModel(lll,kkk) = split(outModel(lll,kkk),plist('frequencies',[flim(1) flim(end)]));
199 end
200 end
201 end
202 if(~isempty(inModel))
203 inModel = split(inModel,plist('frequencies',[flim(1) flim(end)]));
204 end
205 end
206 end
207
208 % use signal fft to get frequency vector. Take into account signal
209 % could be empty or set to zero
210
211 % 1st channel
212 if all(fin(k).getObjectAtIndex(1,1).y == 0) || isempty(fin(k).getObjectAtIndex(1,1).y)
213 i1 = ao(plist('type','fsdata','xvals',0,'yvals',0));
214 % rebuild input with new zeroed signal
215 fin(k) = matrix(i1,fin(k).getObjectAtIndex(2,1),plist('shape',[2 1]));
216 else
217 freqs = fin(k).getObjectAtIndex(1,1).x;
218 end
219 % 2nd channel
220 if all(fin(k).getObjectAtIndex(2,1).y == 0) || isempty(fin(k).getObjectAtIndex(2,1).y)
221 i2 = ao(plist('type','fsdata','xvals',0,'yvals',0));
222 % rebuild input with new zeroed signal
223 fin(k) = matrix(fin(k).getObjectAtIndex(1,1),i2,plist('shape',[2 1]));
224 else
225 freqs = fin(k).getObjectAtIndex(2,1).x;
226 end
227
228 % Compute psd
229 n1 = Nsamples*fs/2*psd(nse(k).getObjectAtIndex(1,1), pl);
230 n2 = Nsamples*fs/2*psd(nse(k).getObjectAtIndex(2,1), pl);
231 n12 = Nsamples*fs/2*cpsd(nse(k).getObjectAtIndex(1,1),nse(k).getObjectAtIndex(2,1), pl);
232
233 % interpolate noise to fft frequencies
234 S11 = interp(n1,plist('vertices',freqs));
235 S12 = interp(n12,plist('vertices',freqs));
236 S22 = interp(n2,plist('vertices',freqs));
237 S21 = conj(S12);
238 % build cross-spectrum matrix
239 mnse(k) = matrix(S11,S21,S12,S22,plist('shape',[2 2]));
240
241 switch class(mdl)
242 case 'matrix'
243 for i = 1:numel(mdl.objs)
244 if (mdlFreqDependent)
245 % set Xvals
246 mdl.objs(i).setXvals(freqs);
247 % set alias
248 mdl.objs(i).assignalias(mdl.objs(i),plist('xvals',freqs));
249 else
250 mdl.objs(i).setXvals(1);
251 end
252 end
253 case 'ssm'
254 if k<2
255 mdl.clearNumParams;
256 spl = plist('set', 'for bode', ...
257 'outputs', outNames, ...
258 'inputs', inNames);
259 % first optimise our model for the case in hand
260 mdl.reorganize(spl);
261 % make it lighter
262
263 %for k = 1:nexp
264
265 mdl(k).optimiseForFitting();
266 %end
267 end
268 end
269 % set Xvals for model (frequency vector could change for each experiment)
270
271 mmdl(k) = mdl;
272
273 elseif ((numel(out(1).objs) == 3) && (numel(in(1).objs) == 4))
274
275 % here we are implementing only the magnetic case
276 % We have 4 inputs (the 4 conformator waveforms of the magnetic
277 % analysis and
278 % 3 outputs (that correspond to the IFO.x12 and IFO.ETA1 and
279 % IFO.PHI1
280
281
282 % fft
283 fin(k) = fft(in(k));
284 fout(k) = fft(out(k));
285 Nsamples = length(fin(k).getObjectAtIndex(1,1).x);
286 fs = fin(k).getObjectAtIndex(1).fs;
287
288 % Get rid of fft f =0, reduce frequency range if needed
289 if ~isempty(flim)
290 fin(k) = split(fin(k),plist('frequencies',[flim(1) flim(end)]));
291 fout(k) = split(fout(k),plist('frequencies',[flim(1) flim(end)]));
292
293 if(~isempty(outModel))
294 for lll=1:size(outModel,1)
295 for kkk=1:size(outModel,2)
296 outModel(lll,kkk) = split(outModel(lll,kkk),plist('frequencies',[flim(1) flim(end)]));
297 end
298 end
299 end
300 if(~isempty(inModel))
301 inModel = split(inModel,plist('frequencies',[flim(1) flim(end)]));
302 end
303 end
304
305 freqs = fin(k).getObjectAtIndex(1,1).x;
306
307 % Build noise model
308 for ii = 1:numel(out.objs)
309 for jj = ii:numel(out.objs)
310 % Compute psd
311 if (ii==jj)
312 n(ii,jj) = Nsamples*fs/2*psd(nse(k).getObjectAtIndex(ii), pl);
313 S(ii,jj) = interp(n(ii,jj),plist('vertices',freqs,'method','linear'));
314 else
315 n(ii,jj) = Nsamples*fs/2*cpsd(nse(k).getObjectAtIndex(ii),nse(k).getObjectAtIndex(jj),pl);
316 S(ii,jj) = interp(n(ii,jj),plist('vertices',freqs,'method','linear'));
317 S(jj,ii) = conj(S(ii,jj));
318 end
319 end
320 end
321 % build cross-spectrum matrix
322 mnse(k) = matrix(S,plist('shape',[numel(out.objs) numel(out.objs)]));
323
324 % Check model sizes
325
326 switch class(mdl)
327 case 'matrix'
328 % set Xvals for model (frequency vector could change for each experiment)
329 for i = 1:numel(mdl.objs)
330 % set Xvals
331 if (mdlFreqDependent)
332 mdl.objs(i).setXvals(freqs);
333 else
334 mdl.objs(i).setXvals(1);
335 end
336 end
337 case 'ssm'
338 spl = plist('set', 'for bode', ...
339 'outputs', outNames, ...
340 'inputs', inNames);
341 % first optimise our model for the case in hand
342 mdl.reorganize(spl);
343 % make it lighter
344 for k = 1:nexp
345 mmdl(k).optimiseForFitting();
346 end
347 end
348 mmdl(k) = mdl;
349
350 elseif ((numel(out(1).objs) == 1) && (numel(in(1).objs) == 1))
351
352 % case 1 input, 1 output used mostly for passing
353 % ao objects from ao.mcmc.
354
355 % optional resampling (not matrix/resample method yet)
356 if ~isempty(fsout)
357 in(k).objs(1).resample(plist('fsout',fsout));
358 out(k).objs(1).resample(plist('fsout',fsout));
359 end
360
361 % fft
362 fin(k) = fft(in(k));
363 fout(k) = fft(out(k));
364
365 Nsamples = length(fin(k).getObjectAtIndex(1).x);
366 fs = fin(k).getObjectAtIndex(1).fs;
367
368
369 if (isempty(outModel))
370 % Get rid of fft f =0, reduce frequency range if needed
371 if ~isempty(flim)
372 fin(k) = split(fin(k),plist('frequencies',[flim(1) flim(end)]));
373 fout(k) = split(fout(k),plist('frequencies',[flim(1) flim(end)]));
374 end
375 elseif (~isempty(outModel))
376 % Get rid of fft f =0, reduce frequency range if needed
377 if ~isempty(flim)
378 fin(k) = split(fin(k),plist('frequencies',[flim(1) flim(end)]));
379 fout(k) = split(fout(k),plist('frequencies',[flim(1) flim(end)]));
380 if(~isempty(outModel))
381 for lll=1:size(outModel,1)
382 for kkk=1:size(outModel,2)
383 outModel(lll,kkk) = split(outModel(lll,kkk),plist('frequencies',[flim(1) flim(end)]));
384 end
385 end
386 end
387 if(~isempty(inModel))
388 inModel = split(inModel,plist('frequencies',[flim(1) flim(end)]));
389 end
390 end
391 end
392
393 % use signal fft to get frequency vector. Take into account signal
394 % could be empty or set to zero
395
396 % 1st and only channel
397 if all(fin(k).getObjectAtIndex(1,1).y == 0) || isempty(fin(k).getObjectAtIndex(1,1).y)
398 i1 = ao(plist('type','fsdata','xvals',0,'yvals',0));
399 % rebuild input with new zeroed signal
400 fin(k) = matrix(i1,fin(k).getObjectAtIndex(2,1),plist('shape',[2 1]));
401 else
402 freqs = fin(k).getObjectAtIndex(1,1).x;
403 end
404
405 % Compute psd
406 n1 = Nsamples*fs/2*psd(nse(k).getObjectAtIndex(1,1), pl);
407
408 % interpolate noise to fft frequencies
409 S11 = interp(n1,plist('vertices',freqs));
410
411 % build cross-spectrum matrix
412 mnse(k) = matrix(S11);
413
414 switch class(mdl)
415 case 'matrix'
416 for i = 1:numel(mdl.objs)
417 if (mdlFreqDependent)
418 % set Xvals
419 mdl.objs(i).setXvals(freqs);
420 % set alias
421 mdl.objs(i).assignalias(mdl.objs(i),plist('xvals',freqs));
422 else
423 mdl.objs(i).setXvals(1);
424 end
425 end
426 case 'ssm'
427 if k<2
428 mdl.clearNumParams;
429 spl = plist('set', 'for bode', ...
430 'outputs', outNames, ...
431 'inputs', inNames);
432 % first optimise our model for the case in hand
433 mdl.reorganize(spl);
434 % make it lighter
435 mdl(k).optimiseForFitting();
436 end
437 end
438 % set Xvals for model (frequency vector could change for each experiment)
439
440 mmdl(k) = mdl;
441
442
443 else
444 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)');
445 end
446 end
447
448 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
449
450 % do simplex
451 if simplex
452 if isempty(xo)
453 error('### Simplex needs a starting guess. Please input a ''x0''.');
454 else
455 switch class(mmdl)
456 case 'matrix'
457 xo = fminsearch(@(x) utils.math.loglikelihood_matrix(x,fin,fout,mnse,mmdl,param,inModel,outModel),xo,optimset('Display','iter'));
458 % case 'smodel'
459 % xo = fminsearch(@(x) utils.math.loglikelihood(x,in,out,nse,mmdl,param),xo,optimset('Display','iter'));
460 case 'ssm'
461 Amats = mdl.amats;
462 Bmats = mdl.bmats;
463 Cmats = mdl.cmats;
464 Dmats = mdl.dmats;
465 xo = fminsearch(@(x) utils.math.loglikelihood_ssm(x,fin,fout,mnse,mmdl,param,inNames,outNames,spl,Amats,Bmats,Cmats,Dmats),xo,optimset('Display','iter'));
466 otherwise
467 error('### Model must be either from the ''smodel'' or the ''ssm'' class. Please check the inputs')
468 end
469
470 for i = 1:numel(param)
471 fprintf('### Simplex estimate: %s = %d \n',param{i},xo(i))
472 save('parameters_simplex.txt','xo','-ASCII')
473 end
474 end
475
476 p = pest(xo);
477 p.setName('SimplexEstimate');
478 p.setNames(param{:});
479 p.setModels(mmdl);
480
481
482 end
483
484 if mhsampleTrue
485
486 % Sampling and saving rejected points. smplr is a matrix of
487 % the size of (# of samples)x(# of params + 2) wich contains the whole
488 % chain for each parameter, a column of zeros and ones (rejected or
489 % accepted point) and the value of the loglikelihood.
490 [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);
491
492 % statistics of the chain
493 if isempty(Tc)
494 initial =1;
495 else
496 initial = Tc(2)+1;
497 end
498 mn = mean(smpl(initial:end,:));
499 cv = cov(smpl(initial:end,:));
500
501 for i = 1:nparam
502 for j = 1:nparam
503 cr(i,j) = cv(i,j)/sqrt(cv(i,i)*cv(j,j));
504 end
505 end
506
507 % compute dof
508 ndata = 0;
509 for i = 1:length(mtxs)
510 for j = 1:length(mtxs(1).objs)
511 ndata = ndata + mtxs(i).objs(j).len;
512 end
513 end
514 dof = ndata-nparam;
515
516 % create pest output
517 p = pest(mn);
518 p.setName('mcmc');
519 p.setNames(param{:});
520 % add statistical info
521 p.setCov(cv);
522 p.setCorr(cr);
523 p.setDy(sqrt(diag(cv)));
524 p.setChain(smpl);
525 %p.computePdf;
526 p.setDof(dof);
527 p.setModels(mdl);
528 % set history
529 %p.addHistory(getInfo('None'), pl, mtxs_invars(:), [out(:).hist mdl(:).hist]);
530 end
531 % Set output
532 % Implemented it this way in order to leave the pest class as it is.
533 output = {p,smplr};
534 %varargout{1} = p;
535 varargout = output(1:nargout);
536
537 end
538
539 % %--------------------------------------------------------------------------
540 % % Compute priors means and sigmas
541 % %--------------------------------------------------------------------------
542 % function fit = fitPrior(prior,nparam)
543 % fit = [];
544 % for ii=1:nparam
545 % [muhat,sigmahat] = normfit(prior(:,2*ii-1),prior(:,2*ii));
546 % fit = [fit ,[muhat,sigmahat]'];
547 % end
548 % end
549
550 %--------------------------------------------------------------------------
551 % Get Info Object
552 %--------------------------------------------------------------------------
553 function ii = getInfo(varargin)
554 if nargin == 1 && strcmpi(varargin{1}, 'None')
555 sets = {};
556 pl = [];
557 else
558 sets = {'Default'};
559 pl = getDefaultPlist;
560 end
561 % Build info object
562 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);
563 end
564
565 %--------------------------------------------------------------------------
566 % Get Default Plist
567 %--------------------------------------------------------------------------
568 function plout = getDefaultPlist()
569 persistent pl;
570 if exist('pl', 'var')==0 || isempty(pl)
571 pl = buildplist();
572 end
573 plout = pl;
574 end
575
576 function pl = buildplist()
577 pl = plist();
578
579 % % inNames
580 % p = param({'inNames','Input names. Used for ssm models'}, paramValue.EMPTY_STRING);
581 % pl.append(p);
582 %
583 % % outNames
584 % p = param({'outNames','Output names. Usde for ssm models'}, paramValue.EMPTY_STRING);
585 % pl.append(p);
586 %
587 % Model
588 % p = param({'model','A matrix array of models.'}, paramValue.EMPTY_STRING);
589 % pl.append(p);
590
591 % Param
592 p = param({'FitParams','A cell array of evaluated parameters.'}, paramValue.EMPTY_DOUBLE);
593 pl.append(p);
594
595 % Input
596 p = param({'input','A matrix array of input signals.'}, paramValue.EMPTY_STRING);
597 pl.append(p);
598
599 pl = plist.MCH_FIT_PLIST;
600
601 % N
602 p = param({'N','number of samples of the chain.'}, paramValue.DOUBLE_VALUE(1000));
603 pl.append(p);
604
605 % Sigma
606 p = param({'cov','covariance of the gaussian jumping distribution.'}, paramValue.DOUBLE_VALUE(1e-4));
607 pl.append(p);
608
609 % Noise
610 p = param({'noise','A matrix array of noise spectrum (PSD) used to compute the likelihood.'}, paramValue.EMPTY_STRING);
611 pl.append(p);
612
613 % Search
614 p = param({'modelFreqDependent','Set to true to use frequency dependent s models, set to false when using constant models'}, paramValue.TRUE_FALSE);
615 pl.append(p);
616
617 % Search
618 p = param({'search','Set to true to use bigger jumps in parameter space during annealing and cool down.'}, paramValue.TRUE_FALSE);
619 pl.append(p);
620
621 % Frequencies
622 p = param({'frequencies','Range of frequencies where the analysis is performed. If an array, only first and last are used'}, paramValue.EMPTY_DOUBLE);
623 pl.append(p);
624
625 % Resample
626 p = param({'fsout','Desired sampling frequency to resample the input time series'}, paramValue.EMPTY_DOUBLE);
627 pl.append(p);
628
629 % Simplex
630 p = param({'simplex','Set to true to perform a simplex search to find the starting parameters of the MCMC chain.'}, paramValue.TRUE_FALSE);
631 pl.append(p);
632
633 % mhsampleTrue
634 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);
635 pl.append(p);
636
637 % heat
638 p = param({'heat','The heat index flattening likelihood surface during annealing.'}, paramValue.DOUBLE_VALUE(1));
639 pl.append(p);
640
641 % Tc
642 p = param({'Tc','An array of two values setting the initial and final value for the cooling down.'}, paramValue.EMPTY_STRING);
643 pl.append(p);
644
645 % heat
646 p = param({'x0','The proposed initial values.'}, paramValue.EMPTY_DOUBLE);
647 pl.append(p);
648
649 % jumps
650 p = param({'jumps','An array of four numbers setting the rescaling of the covariance matrix during the search phase.',...
651 'The first value is the one applied by default, the following thhree apply just when the chain sample is',...
652 'mod(10), mod(25) and mod(100) respectively.'}, paramValue.EMPTY_DOUBLE);
653 pl.append(p);
654
655 % heat
656 p = param({'plot','Select indexes of the parameters to be plotted.'}, paramValue.EMPTY_DOUBLE);
657 pl.append(p);
658
659 % debug
660 p = param({'debug','Set to true to get debug information of the MCMC process.'}, paramValue.FALSE_TRUE);
661 pl.append(p);
662
663 % inModel
664 p = param({'inModel','Input model. Still under test'}, paramValue.EMPTY_STRING);
665 pl.append(p);
666
667 % outModel
668 p = param({'outModel','Output model. Still under test'}, paramValue.EMPTY_STRING);
669 pl.append(p);
670
671 p = param({'prior','Mean, sigma and normalization factor for priors. Still under test'}, paramValue.EMPTY_STRING);
672 pl.append(p);
673
674 p = param({'anneal',['Choose type of annealing during sampling. Default value is ',...
675 'simulated annealing. Choose "thermo" for annealing with a thermostat.',...
676 ' SNR is computed and if it is larger than a fixed value SNR0 (provided also in the plist), ',...
677 'then the chains are heated by a factor of (SNR(1)/SNR0)^2. Choosing "simple" ',...
678 'the deviation of the loglikelihood of every 10 points in the chains is stored. If this deviation ',...
679 'is larger or smaller than two fixed values the chains are cooled or heated respectively.']}, {1, {'simul','thermo', 'simple'}, paramValue.SINGLE});
680 pl.append(p);
681
682 p = param({'SNR0','Fixed value for thermostated annealing.'}, {1, {200}, paramValue.OPTIONAL});
683 pl.append(p);
684
685 p = param({'DeltaL',['Deviation of Loglikelihood for 10 points of the chains. Used for',...
686 'the "simple" choice of annealing with a thermostat.']}, {1, {[100 600 2 3]}, paramValue.OPTIONAL});
687 pl.append(p);
688
689 end
690
691
692 % PARAMETERS: J - number of chains to compute.
693 % sigma - dispersion of the jumping distribution.
694 % range - range where the parameteters are sampled.
695 % param - a cell array of evaluated parameters (from the
696 % smodel).
697 % Tc - an array with two values (default: 0, i.e. no annealing).
698 % heat - heat index flattening likelihood surface (default: 1)
699 % x0 - proposed initial values, if empty selected randomly.
700 % (default: empty)