Line data Source code
1 : /**
2 : Bojan Nikolic <bojan@bnikolic.co.uk>
3 : Initial version 2009
4 :
5 : This file is part of BNMin1 and is licensed under GNU General
6 : Public License version 2
7 :
8 : \file markovchain.cxx
9 : Renamed to markovchain.cc
10 :
11 : */
12 :
13 : #include "markovchain.h"
14 :
15 : namespace Minim {
16 :
17 :
18 0 : ChainBase::ChainBase(const v_t &ic,
19 : fx_t fLkl,
20 0 : fx_t fPr):
21 0 : fLkl(fLkl),
22 0 : fPr(fPr),
23 0 : n(ic.size()),
24 0 : norm_dist(0.,1.),
25 0 : u_dist(0.,1.)
26 : {
27 0 : reset(ic);
28 0 : }
29 :
30 0 : ChainBase::~ChainBase()
31 : {
32 0 : }
33 :
34 0 : double ChainBase::ngen(){
35 0 : return norm_dist(igen);
36 : }
37 :
38 0 : double ChainBase::u01gen(){
39 0 : return u_dist(igen);
40 : }
41 :
42 0 : void ChainBase::reset(const v_t &x, unsigned seed)
43 : {
44 0 : c.x=x;
45 0 : c.l=fLkl(x);
46 0 : c.p=fPr(x);
47 0 : if(seed != 0){
48 0 : igen.seed(seed); // re-seed the generator
49 : }
50 0 : }
51 :
52 0 : MarkovChain::MarkovChain(const v_t &ic,
53 : fx_t fLkl,
54 : fx_t fPr,
55 0 : fa_t fAccept):
56 : ChainBase(ic,
57 : fLkl,
58 : fPr),
59 0 : fAccept(fAccept)
60 : {
61 0 : }
62 :
63 0 : bool MarkovChain::propose(const std::vector<double> &x)
64 : {
65 0 : p.x=x;
66 0 : p.l=fLkl(x);
67 0 : p.p=fPr(x);
68 :
69 0 : const double aprob=fAccept(c,p);
70 :
71 0 : if (aprob>=1.0)
72 : {
73 0 : c=p;
74 0 : return true;
75 : }
76 : else
77 : {
78 0 : if(u01gen() < aprob)
79 : {
80 0 : c=p;
81 0 : return true;
82 : }
83 : }
84 0 : return false;
85 : }
86 :
87 0 : InitPntChain::InitPntChain(const v_t &ic,
88 : fx_t fLkl,
89 : fx_t fPr,
90 0 : fa_t fAccept):
91 : ChainBase(ic,
92 : fLkl,
93 : fPr),
94 0 : fAccept(fAccept)
95 : {
96 0 : f=c;
97 0 : }
98 :
99 0 : bool InitPntChain::propose(const v_t &x)
100 : {
101 0 : p.x=x;
102 0 : p.l=fLkl(x);
103 0 : p.p=fPr(x);
104 :
105 0 : const double aprob=fAccept(f,c,p);
106 :
107 0 : if (aprob>=1.0)
108 : {
109 0 : c=p;
110 0 : return true;
111 : }
112 : else
113 : {
114 0 : if(u01gen() < aprob)
115 : {
116 0 : c=p;
117 0 : return true;
118 : }
119 : }
120 0 : return false;
121 : }
122 :
123 0 : void InitPntChain::reset(const v_t &x, unsigned seed)
124 : {
125 0 : ChainBase::reset(x, seed);
126 0 : f=c;
127 0 : }
128 :
129 0 : ILklChain::ILklChain(const v_t &ic,
130 : fx_t fLkl,
131 : fx_t fPr,
132 0 : fa_t fAccept):
133 : ChainBase(ic,
134 : fLkl,
135 : fPr),
136 0 : fAccept(fAccept)
137 : {
138 0 : L=c.l;
139 0 : }
140 :
141 0 : bool ILklChain::propose(const v_t &x)
142 : {
143 0 : p.x=x;
144 0 : p.l=fLkl(x);
145 0 : p.p=fPr(x);
146 :
147 0 : const double aprob=fAccept(L,c,p);
148 :
149 0 : if (aprob>=1.0)
150 : {
151 0 : c=p;
152 0 : return true;
153 : }
154 : else
155 : {
156 0 : if(u01gen() < aprob)
157 : {
158 0 : c=p;
159 0 : return true;
160 : }
161 : }
162 0 : return false;
163 : }
164 :
165 0 : void ILklChain::reset(const v_t &x,
166 : double Ls, unsigned seed)
167 : {
168 0 : ChainBase::reset(x, seed);
169 0 : L=Ls;
170 0 : }
171 :
172 0 : bool normProp(ChainBase &c,
173 : const std::vector<double> &sigma)
174 : {
175 0 : const std::vector<double> & cx=c.gcx();
176 0 : std::vector<double> px(cx.size());
177 0 : for(size_t i=0; i<cx.size(); ++i)
178 : {
179 0 : px[i]= cx[i]+ sigma[i]*c.ngen();
180 : }
181 0 : return c.propose(px);
182 0 : }
183 :
184 0 : bool normProp(ChainBase &c,
185 : size_t i,
186 : double s)
187 : {
188 0 : const std::vector<double> &cx=c.gcx();
189 0 : std::vector<double> px(cx);
190 0 : px[i]+= s*c.ngen();
191 0 : return c.propose(px);
192 0 : }
193 :
194 0 : bool eigenProp(ChainBase &c,
195 : const std::vector<double> &eigvals,
196 : const std::vector<double> &eigvects
197 : )
198 : {
199 0 : const size_t n=eigvals.size();
200 0 : const std::vector<double> &cx=c.gcx();
201 0 : std::vector<double> px(cx);
202 0 : for(size_t i=0; i<n; ++i)
203 : {
204 0 : const double s=eigvals[i]*c.ngen();
205 0 : for(size_t j=0; j<n; ++j)
206 : {
207 0 : px[j]+= s*eigvects[i*n+j];
208 : }
209 : }
210 0 : return c.propose(px);
211 :
212 0 : }
213 :
214 0 : double metropolis(const MCPoint2 &c,
215 : const MCPoint2 &p)
216 : {
217 0 : if (p.l<c.l)
218 0 : return 1;
219 : else
220 0 : return exp(c.l-p.l);
221 : }
222 :
223 0 : double constrPrior(const MCPoint2 &c,
224 : const MCPoint2 &p)
225 : {
226 0 : if(p.l >= c.l)
227 0 : return 0;
228 : else
229 : {
230 0 : if (p.p<c.p)
231 : {
232 0 : return 1;
233 : }
234 : else
235 : {
236 0 : return exp(c.p-p.p);
237 : }
238 : }
239 : }
240 :
241 0 : double constrPriorP(const MCPoint2 &f,
242 : const MCPoint2 &c,
243 : const MCPoint2 &p)
244 : {
245 0 : if(p.l >= f.l)
246 0 : return 0;
247 : else
248 : {
249 0 : if (p.p<c.p)
250 : {
251 0 : return 1;
252 : }
253 : else
254 : {
255 0 : return exp(c.p-p.p);
256 : }
257 : }
258 : }
259 :
260 0 : double constrPriorL(double L,
261 : const MCPoint2 &c,
262 : const MCPoint2 &p)
263 : {
264 0 : if(p.l >= L)
265 0 : return 0;
266 : else
267 : {
268 0 : if (p.p<c.p)
269 : {
270 0 : return 1;
271 : }
272 : else
273 : {
274 0 : return exp(c.p-p.p);
275 : }
276 : }
277 : }
278 :
279 :
280 : }
281 :
282 :
|
