Mercurial > hg > ltpda

comparison m-toolbox/classes/@ao/mcmc_td.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_TD 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 %
8 % CALL: b = mcmc(out,pl)
9 %
10 % INPUTS:
11 % out - analysis objects with measured outputs
12 % pl - parameter list
13 %
14 % OUTPUTS: b - pest object contatining estimate information
15 %
16 % <a href="matlab:utils.helper.displayMethodInfo('ao', 'mcmc_td')">Parameters Description</a>
17 %
18 % VERSION: $Id: mcmc_td.m,v 1.6 2011/04/08 08:56:12 hewitson Exp $
19 %
20 % References: "Catching supermassive black holes binaries without a net"
21 % N.J. Cornish, E.K. Porter, Phys.Rev.D 75, 021301, 2007
22 %
23 % TODO: multiple chain option not implemented yet
24 % metropolis/hastings not implemented
25 % empty initial values not implemented
26 %
27 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
28
29 function varargout = mcmc_td(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('### metropolis2D 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 bs = copy(as, nargout);
55
56 out = bs;
57
58
59 % Combine plists
60 pl = parse(pl, getDefaultPlist);
61
62
63 % Get parameters
64 N = find(pl, 'N');
65 Tc = find(pl,'Tc');
66 xi = find(pl,'heat');
67 cvar = find(pl, 'cov');
68 rng = find(pl,'range');
69 search = find(pl,'search');
70 simplex = find(pl,'simplex');
71 x0 = find(pl,'x0');
72 mdl = find(pl,'model');
73 %nse = find(pl,'noise');
74 jumps = find(pl,'jumps');
75 parplot = find(pl,'plot');
76 debug = find(pl,'debug');
77
78 in = find(pl,'Input');
79 inNames = find(pl,'inNames');
80 outNames = find(pl,'outNames');
81 Noise = find(pl,'Noise');
82 parnames = find(pl,'parnames');
83 cutbefore = find(pl,'cutbefore');
84 cutafter = find(pl,'cutafter');
85
86 % Check input parameters
87 if isempty(rng)
88 error('### Please define a search range ''range''');
89 end
90
91 if isempty(parnames)
92 error('### Please define the parameters ''parnames''');
93 end
94 nparam = numel(parnames);
95
96
97 % Get range for parameters
98 for i = 1:nparam
99 range(:,i) = rng{i};
100 end
101
102
103
104 % do simplex
105 if simplex
106 if isempty(x0)
107 error('### Simplex needs a starting guess. Please input a ''x0''.');
108 else
109 switch class(mdl)
110 case 'smodel'
111 %xo = fminsearch(@(x) utils.math.loglikehood(x,in,out,nse,mdl,param),xo,optimset('Display','iter'));
112 case 'ssm'
113 x0 = fminsearch(@(x)utils.math.loglikehood_ssm_td(x,in,out,parnames,mdl,inNames,outNames,Noise,'cutbefore',cutbefore,'cutafter',cutafter),x0,optimset('Display','iter'));
114 otherwise
115 error('### Model must be either from the ''smodel'' or the ''ssm'' class. Please check the inputs')
116 end
117
118 for i = 1:numel(parnames)
119 fprintf('### Simplex estimate: %s = %d \n',parnames{i},x0(i))
120 %save('parameters_simplex.txt','x0','-ASCII')
121 end
122 end
123 end
124
125 % sample distribution
126 smpl = utils.math.mhsample_td(mdl,in,out,cvar,N,range,parnames,Tc,xi,x0,search,jumps,parplot,debug,inNames,outNames,Noise,cutbefore,cutafter);
127
128 % statistics of the chain
129 if isempty(Tc)
130 initial =1;
131 else
132 initial = Tc(2)+1;
133 end
134
135 mn = mean(smpl(initial:end,:));
136 cv = cov(smpl(initial:end,:));
137
138 cr = utils.math.cov2corr(cv);
139
140 % compute pdf
141 for i = 1:nparam
142 [count,bin] = hist(smpl(initial:end,i));
143 end
144
145 % create pest output
146 p = pest(mn);
147 p.setName('metropolis2D');
148 p.setNames(parnames{:});
149 % add statistical info
150 p.setCov(cv);
151 p.setCorr(cr);
152 p.setDy(sqrt(diag(cv))); % need to change if non-gaussian
153 p.setPdf([reshape(bin, 1, []), reshape(count, 1, [])]); % not yet able to store bins
154 p.setChain(smpl);
155 p.setModels(mdl);
156 % set history
157 p.addHistory(getInfo('None'), pl, [ao_invars(:)], [bs(:).hist mdl(:).hist]);
158
159 % Set output
160 varargout{1} = p;
161
162 end
163
164
165
166 %--------------------------------------------------------------------------
167 % Get Info Object
168 %--------------------------------------------------------------------------
169 function ii = getInfo(varargin)
170 if nargin == 1 && strcmpi(varargin{1}, 'None')
171 sets = {};
172 pl = [];
173 else
174 sets = {'Default'};
175 pl = getDefaultPlist;
176 end
177 % Build info object
178 ii = minfo(mfilename, 'ao', 'ltpda', utils.const.categories.sigproc, '$Id: mcmc_td.m,v 1.6 2011/04/08 08:56:12 hewitson Exp $', sets, pl);
179 end
180
181 %--------------------------------------------------------------------------
182 % Get Default Plist
183 %--------------------------------------------------------------------------
184 function plout = getDefaultPlist()
185 persistent pl;
186 if exist('pl', 'var')==0 || isempty(pl)
187 pl = buildplist();
188 end
189 plout = pl;
190 end
191
192 function pl = buildplist()
193 pl = plist();
194
195 % N
196 p = param({'N','number of samples of the chain.'}, paramValue.DOUBLE_VALUE(1000));
197 pl.append(p);
198
199 % Sigma
200 p = param({'cov','covariance of the gaussian jumping distribution.'}, paramValue.DOUBLE_VALUE(1e-4));
201 pl.append(p);
202
203 % Param
204 p = param({'parnames','A cell array of evaluated parameters.'}, paramValue.EMPTY_CELL);
205 pl.append(p);
206
207
208 % Model
209 p = param({'model','The model for the system, at the moment only ssm is supported'}, paramValue.EMPTY_STRING);
210 pl.append(p);
211
212 % Search
213 p = param({'search','Set to true to use bigger jumps in parameter space during annealing and cool down.'}, paramValue.TRUE_FALSE);
214 pl.append(p);
215
216 % Simplex
217 p = param({'simplex','Set to true to perform a simplex search to find the starting parameters of the MCMC chain.'}, paramValue.TRUE_FALSE);
218 pl.append(p);
219
220 % heat
221 p = param({'heat','The heat index flattening likelihood surface during annealing.'}, paramValue.DOUBLE_VALUE(1));
222 pl.append(p);
223
224 % Tc
225 p = param({'Tc','An array of two values setting the initial and final value for the cooling down.'}, paramValue.EMPTY_STRING);
226 pl.append(p);
227
228 % x0
229 p = param({'x0','The proposed initial values.'}, paramValue.EMPTY_DOUBLE);
230 pl.append(p);
231
232 % jumps
233 p = param({'jumps','An array of four numbers setting the rescaling of the covariance matrix during the search phase.',...
234 'The first value is the one applied by default, the following thhree apply just when the chain sample is',...
235 'mod(10), mod(25) and mod(100) respectively.'}, paramValue.EMPTY_DOUBLE);
236 pl.append(p);
237
238 % plot
239 p = param({'plot','Select indexes of the parameters to be plotted.'}, paramValue.EMPTY_DOUBLE);
240 pl.append(p);
241
242 % debug
243 p = param({'debug','Set to true to get debug information of the MCMC process.'}, paramValue.FALSE_TRUE);
244 pl.append(p);
245
246 % Input
247 p = param({'Input','An array of analysis object containing the input signals'},paramValue.EMPTY_DOUBLE);
248 pl.append(p);
249
250 % inNames
251 p = param({'inNames','Cell array of input port names'}, paramValue.EMPTY_CELL);
252 pl.append(p);
253
254 % outNames
255 p = param({'outNames','Cell array of output port names'}, paramValue.EMPTY_CELL);
256 pl.append(p);
257
258
259 % Noise
260 p = param({'Noise','An array of analysis objects containg noise series'}, paramValue.EMPTY_DOUBLE);
261 pl.append(p);
262
263 % cutbefore
264 p = param({'cutbefore','the data samples to cut at the starting of the series'}, paramValue.EMPTY_DOUBLE);
265 pl.append(p);
266
267 % cutafter
268 p = param({'cutafter','the data samples to cut at the ending of the series'}, paramValue.EMPTY_DOUBLE);
269 pl.append(p);
270
271
272 end
273
274
275