Mercurial > hg > ltpda

comparison m-toolbox/classes/@ao/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 % It handles (1 input / 1 output) systems, (2 input / 1 output) systems,
7 % and (2 input / 2 output) systems.
8 %
9 % CALL: b = mcmc(in,out,pl)
10 %
11 % INPUTS: out - analysis objects with measured outputs
12 %
13 % pl - parameter list
14 %
15 % OUTPUTS: b - pest object contatining estimate information
16 %
17 % <a href="matlab:utils.helper.displayMethodInfo('ao', 'mcmc')">Parameters Description</a>
18 %
19 % VERSION: $Id: mcmc.m,v 1.24 2011/11/16 15:21:13 nikos Exp $
20 %
21 % References: "Catching supermassive black holes binaries without a net"
22 % N.J. Cornish, E.K. Porter, Phys.Rev.D 75, 021301, 2007
23 %
24 % TODO: multiple chain option not implemented yet
25 % metropolis/hastings not implemented
26 %
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 [as, ao_invars] = utils.helper.collect_objects(varargin(:), 'ao', in_names);
51 pl = utils.helper.collect_objects(varargin(:), 'plist', in_names);
52
53 % Decide on a deep copy or a modify
54 out = copy(as, nargout);
55
56 % Combine plists
57 pl = parse(pl, getDefaultPlist);
58
59
60 % Get parameters from plist
61 mtxin = find(pl,'input');
62 noi = find(pl,'noise');
63 mdlin = find(pl,'model');
64
65 % Decide on a deep copy or a modify
66 in = copy(mtxin, nargout);
67 mdl = copy(mdlin, nargout);
68
69 numexp = numel(in(1,:));
70 numchannels_in = numel(in(:,1));
71 numchannels_out = numel(out(:,1));
72
73 if (numel(in) ~= numel(out))
74 error('### Dimensions of input and output must be equal.');
75 end
76
77 % Redefining the ao inputs as matrix objects
78 for ii = 1:numexp
79 matin(ii) = matrix(in(:,ii),plist('shape',[numchannels_in 1]));
80 matout(ii) = matrix(out(:,ii),plist('shape',[numchannels_in 1]));
81 matnoise(ii) = matrix(noi(:,ii),plist('shape',[numchannels_in 1]));
82 end
83
84 switch class(mdl)
85 case 'smodel'
86 matmodel = matrix(mdl);
87 case 'ssm'
88 matmodel = mdl;
89 case 'matrix'
90 matmodel = mdl;
91 end
92
93 pl = pset(pl,'model',matmodel,'input',matin,'noise',matnoise);
94 mcmc_obj = mcmc(matout,pl);
95
96
97 % Set output
98 varargout{1} = mcmc_obj;
99
100 end
101
102
103
104 %--------------------------------------------------------------------------
105 % Get Info Object
106 %--------------------------------------------------------------------------
107 function ii = getInfo(varargin)
108 if nargin == 1 && strcmpi(varargin{1}, 'None')
109 sets = {};
110 pl = [];
111 else
112 sets = {'Default'};
113 pl = getDefaultPlist;
114 end
115 % Build info object
116 ii = minfo(mfilename, 'ao', 'ltpda', utils.const.categories.sigproc, '$Id: mcmc.m,v 1.24 2011/11/16 15:21:13 nikos Exp $', sets, pl);
117 end
118
119 %--------------------------------------------------------------------------
120 % Get Default Plist
121 %--------------------------------------------------------------------------
122 function plout = getDefaultPlist()
123 persistent pl;
124 if exist('pl', 'var')==0 || isempty(pl)
125 pl = buildplist();
126 end
127 plout = pl;
128 end
129
130 function pl = buildplist()
131 pl = plist();
132
133 % N
134 p = param({'N','number of samples of the chain.'}, paramValue.DOUBLE_VALUE(1000));
135 pl.append(p);
136
137 % Sigma
138 p = param({'cov','covariance of the gaussian jumping distribution.'}, paramValue.DOUBLE_VALUE(1e-4));
139 pl.append(p);
140
141 % Param
142 p = param({'Fitparams','A cell array of evaluated parameters.'}, paramValue.EMPTY_DOUBLE);
143 pl.append(p);
144
145 % Input
146 p = param({'input','A matrix array of input signals.'}, paramValue.EMPTY_STRING);
147 pl.append(p);
148
149 % Frequencies
150 p = param({'frequencies','Range of frequencies where the analysis is performed. If an array, only first and last are used'}, paramValue.EMPTY_DOUBLE);
151 pl.append(p);
152
153 % Noise
154 p = param({'noise','An array of noise spectrum (PSD) used to compute the likelihood.<ul>', ...
155 '<li>1 channel - Input one object</li>', ...
156 '<li>2 channel - Input four objects defining the power spectrum matrix </li>',...
157 '</ul>'}, paramValue.EMPTY_STRING);
158 pl.append(p);
159
160 % Model
161 p = param({'model','An array of models.',...
162 '<li>1 input / 1 output - Input one object</li>', ...
163 '<li>2 input / 1 output - Input two objects</li>', ...
164 '<li>2 channel (2 input/ 2 output) - Input four objects defining the transfer function matrix </li>',...
165 '</ul>'}, paramValue.EMPTY_STRING);
166 pl.append(p);
167
168 % Search
169 p = param({'search','Set to true to use bigger jumps in parameter space during annealing and cool down.'}, paramValue.TRUE_FALSE);
170 pl.append(p);
171
172 % Simplex
173 p = param({'simplex','Set to true to perform a simplex search to find the starting parameters of the MCMC chain.'}, paramValue.TRUE_FALSE);
174 pl.append(p);
175
176 % heat
177 p = param({'heat','The heat index flattening likelihood surface during annealing.'}, paramValue.DOUBLE_VALUE(1));
178 pl.append(p);
179
180 % Tc
181 p = param({'Tc','An array of two values setting the initial and final value for the cooling down.'}, paramValue.EMPTY_DOUBLE);
182 pl.append(p);
183
184 % x0
185 p = param({'x0','The proposed initial values.'}, paramValue.EMPTY_DOUBLE);
186 pl.append(p);
187
188 % jumps
189 p = param({'jumps','An array of four numbers setting the rescaling of the covariance matrix during the search phase.',...
190 'The first value is the one applied by default, the following thhree apply just when the chain sample is',...
191 'mod(10), mod(25) and mod(100) respectively.'}, paramValue.EMPTY_DOUBLE);
192 pl.append(p);
193
194 % plot
195 p = param({'plot','Select indexes of the parameters to be plotted.'}, paramValue.EMPTY_DOUBLE);
196 pl.append(p);
197
198 % debug
199 p = param({'debug','Set to true to get debug information of the MCMC process.'}, paramValue.FALSE_TRUE);
200 pl.append(p);
201
202 p = param({'prior','Mean, sigma and normalization factor for priors. Still under test'}, paramValue.EMPTY_STRING);
203 pl.append(p);
204
205 p = param({'anneal',['Choose type of annealing during sampling. Default value is ',...
206 'simulated annealing. Choose "thermo" for annealing with a thermostat.',...
207 ' SNR is computed and if it is larger than a fixed value SNR0 (provided also in the plist), ',...
208 'then the chains are heated by a factor of (SNR(1)/SNR0)^2. Choosing "simple" ',...
209 'the deviation of the loglikelihood of every 10 points in the chains is stored. If this deviation ',...
210 'is larger or smaller than two fixed values the chains are cooled or heated respectively.']}, {1, {'simul','thermo', 'simple'}, paramValue.SINGLE});
211 pl.append(p);
212
213 p = param({'SNR0','Fixed value for thermostated annealing.'}, {1, {200}, paramValue.OPTIONAL});
214 pl.append(p);
215
216 p = param({'DeltaL',['Deviation of Loglikelihood for 10 points of the chains. Used for',...
217 'the "simple" choice of annealing with a thermostat.']}, {1, {[100 600 2 3]}, paramValue.OPTIONAL});
218 pl.append(p);
219
220 end
221
222
223 % PARAMETERS: J - number of chains to compute.
224 % sigma - dispersion of the jumping distribution.
225 % range - range where the parameteters are sampled.
226 % param - a cell array of evaluated parameters (from the
227 % smodel).
228 % Tc - an array with two values (default: 0, i.e. no annealing).
229 % heat - heat index flattening likelihood surface (default: 1)
230 % x0 - proposed initial values, if empty selected randomly.
231 % (default: empty)
232
233