Introduction
Introduction Statistics Contact Development Disclaimer Help
tMetHasLongstep4.m - cosmo - front and backend for Markov-Chain Monte Carlo inv…
git clone git://src.adamsgaard.dk/cosmo
Log
Files
Refs
README
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 defi…
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-variab…
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 w…
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…
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 acce…
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),f…
131 % disp(['sqrtCpost updated. Cond(sqrtCpost)=',num2str(cond(sqrtCpost…
132 iCpostUpdate = iCpostUpdate+1;
133 sqrtCposts{iCpostUpdate}=sqrtCpostpr;
134 end
135
136 [mprop,propaccepts(it)] = propLongstep(mcur,StepFact,fsSampling,fixed_…
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
You are viewing proxied material from mx1.adamsgaard.dk. The copyright of proxied material belongs to its original authors. Any comments or complaints in relation to proxied material should be directed to the original authors of the content concerned. Please see the disclaimer for more details.