Mercurial > hg > ltpda

comparison m-toolbox/classes/+utils/@math/mhsample_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 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2 %
3 % Metropolis algorithm
4 %
5 % M Nofrarias 15-06-09
6 %
7 % $Id: mhsample_td.m,v 1.1 2010/11/26 13:48:48 luigi Exp $
8 %
9 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
10
11
12 function smpl = mhsample_td(model,in,out,cov,number,limit,parnames,Tc,xi,xo,search,jumps,parplot,dbg_info,inNames,outNames,Noise,cutbefore,cutafter)
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 nparnames = length(parnames);
28
29
30 % switch depending on model class (initial likelihood)
31 switch class(model)
32 case 'matrix'
33 % loglk1 = utils.math.loglikehood_matrix(xo,in,out,nse,model,parnames);
34 % loglk1 = utils.math.loglikehood_matrix(xo,in,out,nse,model,parnames,inModel,outModel);
35 case 'smodel'
36 % loglk1 = utils.math.loglikehood(xo,in,out,nse,model,parnames);
37 case 'ssm'
38 loglk1 = utils.math.loglikehood_ssm_td(xo,in,out,parnames,model,inNames,outNames,Noise,'cutbefore',cutbefore,'cutafter',cutafter);
39 otherwise
40 error('### Model must be either from the ''smodel'' or the ''ssm'' class. Please check the inputs')
41 end
42
43 % accept the first sample if no search active
44 if ~search
45 smpl(1,:) = xo;
46 nacc = 1;
47 else
48 smpl = [];
49 nacc = 0;
50 end
51 nrej = 0;
52
53 utils.helper.msg(msg.IMPORTANT, 'Starting Monte Carlo sampling', mfilename('class'), mfilename);
54
55
56 hjump = 0;
57 % main loop
58 while(nacc<number)
59 % jumping criterion during search phase
60 % - 2-sigma by default
61 % - 10-sigma jumps at mod(10) samples
62 % - 100-sigma jumps at mod(25) samples
63 % - 1000-sigma jumps at mod(100) samples
64 if search
65 if nacc <= Tc(1)
66 if(mod(nacc,10) == 0 && mod(nacc,25) ~= 0 && mod(nacc,100) ~= 0 && hjump ~= 1)
67 hjump = 1;
68 xn = mvnrnd(xo,jumps(2)^2*cov);
69 elseif(mod(nacc,25) == 0 && mod(nacc,100) ~= 0 && hjump ~= 1)
70 hjump = 1;
71 xn = mvnrnd(xo,jumps(3)^2*cov);
72 elseif(mod(nacc,100) == 0 && hjump ~= 1)
73 hjump = 1;
74 xn = mvnrnd(xo,jumps(4)^2*cov);
75 % xn = xo + range/2+rand(size(xo));
76 else
77 hjump = 0;
78 xn = mvnrnd(xo,jumps(1)^2*cov);
79 end
80 else
81 xn = mvnrnd(xo,cov);
82 end
83 else
84 xn = mvnrnd(xo,cov);
85 end
86
87 % check if out of limits
88 if( any(xn < limit(1,:)) || any(xn > limit(2,:)))
89 logprior = inf;
90 loglk2 = inf;
91 betta = inf;
92 else
93 % compute annealing
94 if ~isempty(Tc)
95 if nacc <= Tc(1)
96 betta = 1/2 * 10^(-xi*(1-Tc(1)/Tc(2)));
97 elseif Tc(1) < nacc && nacc <= Tc(2)
98 betta = 1/2 * 10^(-xi*(1-nacc/Tc(2)));
99 else
100 betta = 1/2;
101 end
102 else
103 betta = 1/2;
104 end
105 % switch depending on model class (new likelihood)
106 switch class(model)
107 case 'matrix'
108 % loglk2 = utils.math.loglikehood_matrix(xn,in,out,nse,model,parnames,inModel,outModel);
109 case 'smodel'
110 % loglk2 = utils.math.loglikehood(xn,in,out,nse,model,parnames);
111 case 'ssm'
112 loglk2 = utils.math.loglikehood_ssm_td(xn,in,out,parnames,model,inNames,outNames,Noise,'cutbefore',cutbefore,'cutafter',cutafter);
113 otherwise
114 error('### Model must be either from the ''smodel'' or the ''ssm'' class. Please check the inputs')
115 end
116 end
117
118 % likelihood ratio
119 logr = betta*(loglk2 - loglk1) ;
120 % decide if sample is accepted or not
121 if logr < 0
122 xo = xn;
123 nacc = nacc+1;
124 smpl(nacc,:) = xn;
125 loglk1 =loglk2;
126 if dbg_info
127 utils.helper.msg(msg.IMPORTANT, sprintf('acc.\t loglik: %d -> %d betta: %d ratio: %d',loglk1,loglk2,betta,logr), mfilename('class'), mfilename);
128 end
129 elseif rand(1) > (1 - exp(-logr))
130 xo = xn;
131 nacc = nacc+1;
132 smpl(nacc,:) = xn;
133 loglk1 =loglk2;
134 if dbg_info
135 utils.helper.msg(msg.IMPORTANT, sprintf('acc.\t loglik: %d -> %d betta: %d ratio: %d',loglk1,loglk2,betta,logr), mfilename('class'), mfilename);
136 end
137 else
138 nrej = nrej+1;
139 if dbg_info
140 utils.helper.msg(msg.IMPORTANT, sprintf('rej.\t loglik: %d -> %d betta: %d ratio: %d',loglk1,loglk2,betta,logr), mfilename('class'), mfilename);
141 end
142 end
143
144 % display and save
145 if(mod(nacc,100) == 0 && nacc ~= oldacc)
146 updacc = nacc-oldacc;
147 updrej = nrej-oldrej;
148 ratio(nacc/10,:) = updacc/(updacc+updrej);
149 utils.helper.msg(msg.IMPORTANT, sprintf('accepted: %d rejected : %d acc. rate: %4.2f',nacc,updrej,ratio(end)), mfilename('class'), mfilename);
150 for i = 1:numel(parnames)
151 fprintf('### Parameters: %s = %d \n',parnames{i},xn(i))
152 end
153 oldacc = nacc;
154 oldrej = nrej;
155 save('chain.txt','smpl','-ASCII')
156 save('acceptance.txt','ratio','-ASCII')
157 end
158
159 % plot
160 if ((mod(nacc,100) == 0) && (nacc ~= 0) && ~isempty(parplot))
161 for i = 1:numel(parplot)
162 figure(parplot(i))
163 plot(smpl(:,parplot(i)))
164 end
165 end
166
167 end
168 end