Mercurial > hg > ltpda

comparison m-toolbox/classes/+utils/@math/mhsample.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 % Metropolis algorithm
4 %
5 % M Nofrarias 15-06-09
6 %
7 % $Id: mhsample.m,v 1.19 2011/11/16 08:52:50 nikos Exp $
8 % prior
9 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
10
11
12 function [smpl smplr] = mhsample(model,in,out,nse,cov,number,limit,param,Tc,xi,xo,search,jumps,parplot,dbg_info,inNames,outNames,fpars,anneal,SNR0,DeltaL,inModel,outModel)
13
14 import utils.const.*
15
16 % compute mean and range
17 mn = mean(limit);
18 range = diff(limit);
19 % initialize
20 acc = [];
21 nacc = 1;
22 nrej = 1;
23 loop = 0;
24 oldacc = 0;
25 oldrej = 0;
26 chgsc = 0;
27 nparam = length(param);
28 delta = 1;
29 SNR = 0;
30 heat = xi;
31
32 % plist to pass in to the sampler for use in bode
33 spl = plist('outputs', outNames, ...
34 'inputs', inNames, ...
35 'reorganize', false,...
36 'numeric output',true);
37
38 % check if using priors
39 nparams = size(xo);
40 if isempty(fpars)
41 fpars = zeros(nparams(2),3);
42 end
43
44
45
46 % switch depending on model class (initial likelihood)
47 switch class(model)
48 case 'matrix'
49 [loglk1 SNR] = utils.math.loglikelihood_matrix(xo,in,out,nse,model,param,inModel,outModel);
50 case 'smodel'
51 loglk1 = utils.math.loglikelihood(xo,in,out,nse,model,param);
52 case 'ssm'
53 Amats = model.amats;
54 Bmats = model.bmats;
55 Cmats = model.cmats;
56 Dmats = model.dmats;
57 loglk1 = utils.math.loglikelihood_ssm(xo,in,out,nse,model,param,inNames,outNames, spl, Amats, Bmats, Cmats, Dmats);
58 otherwise
59 error('### Model must be either from the ''smodel'' or the ''ssm'' class. Please check the inputs')
60 end
61
62 % accept the first sample if no search active
63 if ~search
64 smpl(1,:) = xo;
65 nacc = 1;
66 else
67 smpl = [];
68 nacc = 0;
69 end
70 nrej = 0;
71
72 utils.helper.msg(msg.IMPORTANT, 'Starting Monte Carlo sampling', mfilename('class'), mfilename);
73
74 % compute prior for the initial sample of the chain
75 logprior1 = logPriors(xo,fpars);
76 % if initial guess is far away, this contidion sets our logprior
77 % to zero and takes account only the loglikelihood for the
78 % computation of logratio. From this initial point it starts
79 % searching for the small area of acceptance.
80 if logprior1 == -Inf;
81 logprior1 = 0;
82 end
83
84 T = Tc(1);
85 hjump = 0;
86 % main loop
87 while(nacc<number)
88 % jumping criterion during search phase
89 % - 2-sigma by default
90 % - 10-sigma jumps at mod(10) samples
91 % - 100-sigma jumps at mod(25) samples
92 % - 1000-sigma jumps at mod(100) samples
93
94
95 if search
96 if nacc <= Tc(1)
97 if(mod(nacc,10) == 0 && mod(nacc,25) ~= 0 && mod(nacc,100) ~= 0 && hjump ~= 1)
98 hjump = 1;
99 xn = mvnrnd(xo,jumps(2)^2*cov);
100 elseif(mod(nacc,20) == 0 && mod(nacc,100) ~= 0 && hjump ~= 1)
101 hjump = 1;
102 xn = mvnrnd(xo,jumps(3)^2*cov);
103 elseif(mod(nacc,50) == 0 && hjump ~= 1)
104 hjump = 1;
105 xn = mvnrnd(xo,jumps(4)^2*cov);
106 % xn = xo + range/2+rand(size(xo));
107 else
108 hjump = 0;
109 xn = mvnrnd(xo,jumps(1)^2*cov);
110 end
111 else
112 xn = mvnrnd(xo,cov); % changed that too
113 end
114 else
115 xn = mvnrnd(xo,cov);
116 end
117
118
119 % compute prior probability
120 logprior2 = logPriors(xn,fpars);
121
122 %check if out of limits
123 if( any(xn < limit(1,:)) || any(xn > limit(2,:)))
124 logprior2 = inf;
125 loglk2 = inf;
126 betta = inf;
127 %heat = inf;
128
129 % This condition save us computation time.
130 elseif logprior2 == -Inf
131 loglk2 = inf;
132 betta = inf;
133 else
134 % switch depending on model class (new likelihood)
135 switch class(model)
136 case 'matrix'
137 [loglk2 SNR] = utils.math.loglikelihood_matrix(xn,in,out,nse,model,param,inModel,outModel);
138 case 'smodel'
139 loglk2 = utils.math.loglikelihood(xn,in,out,nse,model,param);
140 case 'ssm'
141 [loglk2 SNR] = utils.math.loglikelihood_ssm(xn,in,out,nse,model,param,inNames,outNames, spl, Amats, Bmats, Cmats, Dmats);
142 otherwise
143 error('### Model must be either from the ''smodel'' or the ''ssm'' class. Please check the inputs')
144 end
145
146
147 % compute annealing
148 if ~isempty(Tc)
149 if nacc <= Tc(1)
150
151 % compute heat factor
152 switch anneal
153 case 'simul'
154 delta = 1;
155 case 'thermo'
156 if (0 <= SNR(1) && SNR(1) <= SNR0)
157 delta = 1;
158 elseif (SNR(1) > SNR0)
159 delta = (SNR(1)/SNR0)^2;
160 end
161 case 'simple'
162 if (nacc > 10 && mod(nacc,10) == 0)
163 deltalogp = std(smpl(nacc-10:nacc,1));
164 if deltalogp <= DeltaL(1)
165 T = T + DeltaL(3);
166 elseif deltalogp >= DeltaL(2)
167 T = T - DeltaL(4);
168 end
169 delta = Tc(1)/T;
170 end
171 heat = xi*delta;
172 end
173
174 betta = 1/2 * 10^(-heat*(1-Tc(1)/Tc(2)));
175 elseif Tc(1) < nacc && nacc <= Tc(2)
176 betta = 1/2 * 10^(-heat*(1-nacc/Tc(2)));
177 else
178 betta = 1/2;
179 end
180 else
181 betta = 1/2;
182 end
183
184 end % here we are
185
186 % likelihood ratio
187 logr = betta*(logprior2 + loglk2 - loglk1 - logprior1) ;
188 % decide if sample is accepted or not
189 if logr < 0
190 xo = xn;
191 nacc = nacc+1;
192 sumsamples = nacc + nrej;
193 smpl(nacc,1) = loglk2;
194 smpl(nacc,2:size(xn,2)+1) = xn;
195 smplr(sumsamples,1:size(xn,2)) = xn;
196 smplr(sumsamples,size(xn,2)+1) = 1;
197 smplr(sumsamples,size(xn,2)+2) = loglk2;
198 loglk1 =loglk2;
199 logprior1=logprior2;
200 if dbg_info
201 %utils.helper.msg(msg.IMPORTANT, sprintf('acc.\t loglik: %d -> %d priors: %d -> %d betta: %d ratio: %d',loglk1,loglk2,logprior1,logprior2,betta,logr), mfilename('class'), mfilename);
202 utils.helper.msg(msg.IMPORTANT, sprintf('acccount: %d SNR: %d %d Heat: %d',nacc,SNR(1),SNR(2),heat));
203 end
204 elseif rand(1) > (1 - exp(-logr))
205 xo = xn;
206 nacc = nacc+1;
207 sumsamples = nacc + nrej;
208 smpl(nacc,1) = loglk2;
209 smpl(nacc,2:size(xn,2)+1) = xn;
210 smplr(sumsamples,1:size(xn,2)) = xn;
211 smplr(sumsamples,size(xn,2)+1) = 1;
212 smplr(sumsamples,size(xn,2)+2) = loglk2;
213 loglk1 =loglk2;
214 logprior1=logprior2;
215 if dbg_info
216 %utils.helper.msg(msg.IMPORTANT, sprintf('acc.\t loglik: %d -> %d priors: %d -> %d betta: %d ratio: %d',loglk1,loglk2,logprior1,logprior2,betta,logr), mfilename('class'), mfilename);
217 utils.helper.msg(msg.IMPORTANT, sprintf('acccount: %d SNR: %d %d Heat: %d',nacc,SNR(1),SNR(2),heat));
218 end
219 else
220 nrej = nrej+1;
221 sumsamples = nacc + nrej;
222 smplr(sumsamples,1:size(xn,2)) = xn;
223 smplr(sumsamples,size(xn,2)+1) = 0;
224 smplr(sumsamples,size(xn,2)+2) = loglk2;
225 if dbg_info
226 %utils.helper.msg(msg.IMPORTANT, sprintf('rej.\t loglik: %d -> %d priors: %d -> %d betta: %d ratio: %d',loglk1,loglk2,logprior1,logprior2,betta,logr), mfilename('class'), mfilename);
227 utils.helper.msg(msg.IMPORTANT, sprintf('rejcount: %d SNR: %d %d Heat: %d',nrej,SNR(1),SNR(2),heat));
228 end
229
230 end
231
232 % display and save
233 if(mod(nacc,100) == 0 && nacc ~= oldacc)
234 updacc = nacc-oldacc;
235 updrej = nrej-oldrej;
236 ratio(nacc/10,:) = updacc/(updacc+updrej);
237 utils.helper.msg(msg.IMPORTANT, sprintf('accepted: %d rejected : %d acc. rate: %4.2f',nacc,updrej,ratio(end)), mfilename('class'), mfilename);
238 for i = 1:numel(param)
239 fprintf('### Parameters: %s = %d \n',param{i},xn(i))
240 end
241 oldacc = nacc;
242 oldrej = nrej;
243 save('chain.txt','smpl','-ASCII')
244 save('acceptance.txt','ratio','-ASCII')
245 end
246
247 % plot
248 if ((mod(sumsamples,100) == 0) && (nacc ~= 0) && ~isempty(parplot))
249 for i = 1:numel(parplot)
250 figure(parplot(i))
251 plot(smpl(:,parplot(i)))
252 end
253 figure(i+1)
254 plot(smplr(:,2))
255 end
256
257 end
258 end
259
260 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
261 % compute logpriors
262 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
263 function logprior = logPriors(parVect,fitparams)
264
265 D = size(parVect);
266
267 for ii=1:D(2)
268 prior(ii) = normpdf(parVect(ii),fitparams(ii,1),fitparams(ii,2))/fitparams(ii,3);
269 % checking if priors are used for this run
270 if isnan(prior(ii))
271 prior(ii) = 1;
272 end
273 end
274 prior = log(prior);
275 logprior = sum(prior);
276
277 % assuming that priors are independent, then
278 % p(x|8n)=p(x|81)p(x|82)p(x|83)...p(x|8n)=p1xp2xp3...xpn
279 % and log(p(x|8n)) = logp1+logp2+...+logpn
280
281 end
282
283