tMetHasLongstep4.m - cosmo - front and backend for Markov-Chain Monte Carlo inversion of cosmogenic nuclide concentrations
(HTM) git clone git://src.adamsgaard.dk/cosmo
(DIR) Log
(DIR) Files
(DIR) Refs
(DIR) README
(DIR) LICENSE
---
tMetHasLongstep4.m (7389B)
---
1 %file: MetHasLongstep4.m
2 %'cosmolongstep' changed to 'CosmoLongsteps'.
3
4 %Developed from MetHasLongstep2
5 %Developed from MetHasLongstep
6 %But now including
7 % - linearized inspiration proposer using partial derivatives
8 % - computed either separately or
9 % - computed from from previous samples
10 %
11 %Metropolis-Hastings iterations
12 %Inside a function call
13 % d = g(m_prop,fixed_stuff)
14 %computes the data to be compared to the measured data, d_obs.
15 %Make sure to let function g be availble at run-time.
16 %Make sure that the paramater fixed_stuff (possibly a structure) is defined
17 %properly, including measurement geometry etc.
18 %INPUTS:
19 % m_start: First current model
20 % seed: if ~isempty(seed), setseed(seed), end Makes walks repeatable
21 % fixed-stuff: A variable, possible a structure which defines non-variable
22 % parameters for function g
23 %OUTPUTS:
24 % ms: Each column is a stored state. Note that only one out of N_skip
25 % states were stored.
26 % accepts: accepts(i)=1 if this proposal was accepted, else accepts(i)=0
27 % Qs: Qds(i) + Qps(i)
28 % Qds: Qds(i) = (d_obs-g(ms(i))'*inv(C_obs)*(d_obs-g(ms(i))
29 % Qps: Qps(i) = (m_prior - ms(i))'*inv(C_prior)*(m_prior-ms(i))
30 % lump: A number of additional diagnostics
31
32 % function [ms,accepts,Qs,Qds,lump]=MetHasLongstep3(...
33 % m_start,d_obs,seed,fixed_stuff);
34 function [ms,accepts,Qs,Qds,lump]=MetHasLongstep(...
35 m_start,seed,isBurnIn,fixed_stuff, ...
36 statusfile) %%%% ADDED BY ANDERS
37 d_obs = fixed_stuff.d_obs;
38 if isBurnIn
39 fsSampling = fixed_stuff.BurnIn;
40 else
41 fsSampling = fixed_stuff.Sampling;
42 end
43
44 %Initializations
45 diagCobs = fixed_stuff.ErrorStdObs.^2; %the variances
46 C_obs = diag(diagCobs);
47 iCobs = inv(C_obs);
48 StepFact = fsSampling.StepFactor;
49 sqrtCpostpr = []; %We have to bootstrap
50
51 if ~isempty(seed)
52 setseed(seed);
53 end
54
55 %Preallocations
56 M = length(m_start);
57 N_run = fsSampling.Nsamp;
58 N_skip = fsSampling.Nskip;
59 ms = NaN*ones(M,floor(N_run/N_skip));
60 Qs = NaN*ones(1,floor(N_run/N_skip));
61 Qds = Qs;
62 % Qps = Qs;
63
64 accepts = zeros(N_run,1);
65 propaccepts = accepts;
66
67 mprops = zeros(M,N_run);
68 Qprops = zeros(1,N_run);
69 Qdprops = Qprops;
70 % Qpprops = Qprops;
71
72 mcurs = zeros(M,N_run);
73 Qcurs = zeros(1,N_run);
74 Qdcurs = Qcurs;
75 % Qpcurs = Qcurs;
76
77 dprops = NaN(length(d_obs),N_run);
78 dcurs = dprops;
79
80 %Initializing mcur as m_start
81 mcur = m_start;
82 mprops(:,1)=m_start; %Just to make updateCpost have something to start with
83 [dcur,lumpprop] = g(mcur,fixed_stuff);
84 dd_obs = d_obs - dcur;
85 Qdcur = dd_obs'*iCobs*dd_obs;
86
87 % dm_pri = m_prior - mcur;
88 % Qpcur = dm_pri'*iCpri*dm_pri;
89
90 Qcur = Qdcur; %+Qpcur;
91 if strcmp(fixed_stuff.g_case,'CosmoLongsteps')
92 zsss{1}=lumpprop.zss;
93 tss{1}=lumpprop.ts;
94 ExposureTimeSinceNows{1}=lumpprop.ExposureTimeSinceNow;
95 c14Csss{1}=lumpprop.c14Css;
96 c10Besss{1}=lumpprop.c10Bess;
97 c26Alsss{1}=lumpprop.c26Alss;
98 c21Nesss{1}=lumpprop.c21Ness;
99
100 end
101 i_keep = 0; %index of the so far last stored model
102 iCpostUpdate = 0;
103
104 %................... The main iteration loop
105 for it=1:N_run
106 if rem(it,1000)==0,
107 [ElapsedTime,RemainingTime,TotalTime,EndTime,IterationsPerSecond,PrintString] = ...
108 RuntimeStatus(it,isBurnIn,fixed_stuff);
109 disp([num2str(it),': ',PrintString])
110
111 %%% INSERTED BY ANDERS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
112 fid = fopen(statusfile, 'w');
113 statusinfo = strcat(num2str(datestr(RemainingTime,13)), ...
114 ' remaining', ...
115 ' <!-- ', num2str(ElapsedTime/TotalTime*100., '%3.0f'), '-->');
116 fprintf(fid, statusinfo);
117 fclose(fid);
118 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
119 end
120
121 if rem(it,fsSampling.StepsPerFactorUpdate)==0
122 %update StepFact
123 acc=accepts((it-fsSampling.StepsPerFactorUpdate+1):it); %recent accepts
124 StepFact = UpdateStepFact(StepFact,acc,fsSampling);
125 disp(['StepFact updated to ',num2str(StepFact),'. It=',num2str(it)])
126 end
127 if rem(it,fsSampling.StepsPerCpostUpdate)==1
128 %update sqrtCpost already at the start of the first iterations
129 is = max(1,it-fsSampling.StepsPerCpostUpdate-1):it;
130 [sqrtCpostpr,Cpostpr,Gestpr]=updateCpost(mprops(:,is),dprops(:,is),fsSampling,fixed_stuff);
131 % disp(['sqrtCpost updated. Cond(sqrtCpost)=',num2str(cond(sqrtCpostpr)),'. It=',num2str(it)])
132 iCpostUpdate = iCpostUpdate+1;
133 sqrtCposts{iCpostUpdate}=sqrtCpostpr;
134 end
135
136 [mprop,propaccepts(it)] = propLongstep(mcur,StepFact,fsSampling,fixed_stuff,sqrtCpostpr); %Mads' proposer, takes care of hard limits
137 dmpropspr(:,it) = m2mpr(mprop,fixed_stuff) - m2mpr(mcur,fixed_stuff);
138
139 [dprop,lumpprop] = g(mprop,fixed_stuff);
140 dd_obs = d_obs - dprop;
141 Qdprop = dd_obs'*iCobs*dd_obs;
142 Qprop = Qdprop; %+Qpprop;
143
144 alpha = rand;
145
146 if alpha < exp(-0.5*(Qprop-Qcur))
147 if propaccepts(it)
148 accepts(it)=1;
149 end
150 mcur=mprop;
151 dcur=dprop;
152 Qdcur = Qdprop;
153 % Qpcur = Qpprop;
154 Qcur = Qprop;
155 end
156 if strcmp(fixed_stuff.g_case,'CosmoLongsteps')
157 if accepts(it)
158 zsss{it+1}=lumpprop.zss; %note: varying length
159 tss{it+1}=lumpprop.ts;
160 ExposureTimeSinceNows{it+1}=lumpprop.ExposureTimeSinceNow;
161 c14Csss{it+1}=lumpprop.c14Css;
162 c10Besss{it+1}=lumpprop.c10Bess;
163 c26Alsss{it+1}=lumpprop.c26Alss;
164 c21Nesss{it+1}=lumpprop.c21Ness;
165
166
167 else %mcur remains; take the previous model
168 zsss{it+1}=zsss{it};
169 tss{it+1}=tss{it};
170 ExposureTimeSinceNows{it+1}=ExposureTimeSinceNows{it};
171 c14Csss{it+1}=c14Csss{it};
172 c10Besss{it+1}=c10Besss{it};
173 c26Alsss{it+1}=c26Alsss{it};
174 c21Nesss{it+1}=c21Nesss{it};
175 end
176 end
177 dmcurspr(:,it) = dmpropspr(:,it)*accepts(it);
178
179 if mod(it,N_skip)==0
180 %pick out sample to keep
181 i_keep = i_keep+1;
182 ms(:,i_keep)=mcur;
183 Qds(i_keep)=Qdcur;
184 % Qps(i_keep)=Qpcur;
185 Qs(i_keep)=Qcur;
186 end
187
188 mcurs(:,it) = mcur;
189 dcurs(:,it) = dcur;
190 Qcurs(it) = Qcur;
191 Qdcurs(it) = Qdcur;
192 % Qpcurs(it) = Qpcur;
193
194 mprops(:,it) = mprop;
195 dprops(:,it) = dprop;
196 Qprops(:,it) = Qprop;
197 Qdprops(:,it) = Qdprop;
198 % Qpprops(:,it) = Qpprop;
199 StepFacts(it) = StepFact;
200 end
201
202 if nargout==5
203
204 lump.StepFacts = StepFacts;
205
206 lump.mcurs = mcurs;
207 lump.Qcurs = Qcurs;
208 lump.Qdcurs = Qdcurs;
209 lump.dmcurspr = dmcurspr;
210 % lump.Qpcurs = Qpcurs;
211
212 lump.mprops = mprops;
213 lump.Qprops = Qprops;
214 lump.Qdprops = Qdprops;
215 lump.dmpropspr = dmpropspr;
216 % lump.Qpprops = Qpprops;
217
218 if strcmp(fixed_stuff.g_case,'CosmoLongsteps')
219 lump.zsss = zsss;
220 lump.tss = tss;
221 lump.ExposureTimeSinceNows=ExposureTimeSinceNows;
222 lump.c14Csss=c14Csss;
223 lump.c10Besss=c10Besss;
224 lump.c26Alsss=c26Alsss;
225 lump.c21Nesss=c21Nesss;
226 end
227 end
228
229
230 %{
231 function u = prop(proposer_type,sqrtC)
232 switch proposer_type
233 case 1 % isotropic unit coordinate variances:
234 u = randn(size(sqrtC,1),1);
235 case 2 % box-shaped coordinate variances:
236 u = sqrt(12)*(rand(size(sqrtC,1))-0.5);
237 case {3,4} %prior proposer:
238 u = sqrtC*randn(size(m));
239 case 5 %proposer defined by external function
240 u = prop_external;
241 end
242