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.hxx
9 : Renamed to markovchain.h 2023
10 :
11 : */
12 : #ifndef _BNMIN1_MARKOVCHAIN_HXX__
13 : #define _BNMIN1_MARKOVCHAIN_HXX__
14 :
15 : #include <vector>
16 :
17 : #include <random>
18 :
19 : #include <functional> // was: <boost/function.hpp>
20 :
21 : namespace Minim {
22 :
23 : struct MCPoint2
24 : {
25 : /// The position in space of the point
26 : std::vector<double> x;
27 :
28 : /// The value of likelihood
29 : double l;
30 :
31 : /// The value of the prior
32 : double p;
33 : };
34 :
35 : /** \brief Base class for chains-type exploration algorithms
36 :
37 : */
38 : class ChainBase
39 : {
40 :
41 : public:
42 :
43 : /// Type of the integer random number generator
44 : typedef std::mt19937 irg_t;
45 :
46 : /// Vector type
47 : typedef std::vector<double> v_t;
48 :
49 : /// A function of a point in the chain
50 : typedef std::function< double (const v_t &x) > fx_t;
51 :
52 : private:
53 :
54 : /// The random number generator
55 : irg_t igen;
56 :
57 : /// The distributions
58 : std::normal_distribution<double> norm_dist;
59 :
60 : std::uniform_real_distribution<double> u_dist;
61 :
62 : protected:
63 :
64 : /// The likelihood function
65 : fx_t fLkl;
66 :
67 : /// The prior function
68 : fx_t fPr;
69 :
70 : /// The current position of the chain
71 : MCPoint2 c;
72 :
73 : /// The next proposed point
74 : MCPoint2 p;
75 :
76 : public:
77 :
78 : /// The normal distribution generator (between 0 and 1)
79 : double ngen();
80 :
81 : // Uniform random numbers between 0 and 1
82 :
83 : double u01gen();
84 :
85 : // ---------- Public data -----------------------------
86 :
87 : /// Number of parameters
88 : const size_t n;
89 :
90 : // ---------- Construction / Destruction --------------
91 :
92 : ChainBase(const v_t &ic,
93 : fx_t fLkl,
94 : fx_t fPr);
95 :
96 : virtual ~ChainBase();
97 :
98 : // ---------- Public interface --------------------------
99 :
100 : /** Reset the chain, discard existing data and set x as the
101 : starting point, optionally reset the seed (default 0 = do not reset) */
102 : void reset(const v_t &x, unsigned seed=0);
103 :
104 : /** \brief Propose x as the next point in chain
105 :
106 : \returns true if the point has been accepted
107 : */
108 : virtual bool propose(const v_t &x) =0;
109 :
110 :
111 : // ------------ Access type functions ---------------
112 :
113 : /// Return the current point of the chain
114 77873424 : const std::vector<double> &gcx(void) const
115 : {
116 77873424 : return c.x;
117 : }
118 :
119 : /// Return the current likelihood
120 771024 : double gcl(void) const
121 : {
122 771024 : return c.l;
123 : }
124 :
125 : };
126 :
127 : /** \brief Markov chain class
128 :
129 : Proposals depend only on the last point in the chain
130 : */
131 : class MarkovChain:
132 : public ChainBase
133 : {
134 : public:
135 :
136 : /// A function to return the acceptance probability
137 : typedef std::function< double (const MCPoint2 &c, const MCPoint2 &p) > fa_t;
138 :
139 : private:
140 :
141 : /// The acceptance function
142 : fa_t fAccept;
143 :
144 : public:
145 :
146 : // ---------- Construction / Destruction --------------
147 :
148 : MarkovChain(const v_t &ic,
149 : fx_t fLkl,
150 : fx_t fPr,
151 : fa_t fAccept);
152 :
153 : // ---------- Public interface --------------------------
154 :
155 :
156 : bool propose(const v_t &x);
157 :
158 : };
159 :
160 : /** Minimum information for constrained prior exploration
161 :
162 : Proposal depend on the first and the last points in the chain
163 : */
164 : class InitPntChain:
165 : public ChainBase
166 : {
167 : public:
168 :
169 : /** \brief Function type to return the acceptance probability
170 :
171 : \param f is the first point in the chain
172 : */
173 : typedef std::function< double (const MCPoint2 &f, const MCPoint2 &c, const MCPoint2 &p) > fa_t;
174 :
175 : private:
176 :
177 : /// First point in the chain
178 : MCPoint2 f;
179 :
180 : /// The acceptance function
181 : fa_t fAccept;
182 :
183 : public:
184 :
185 : // ---------- Construction / Destruction --------------
186 :
187 : InitPntChain(const v_t &ic,
188 : fx_t fLkl,
189 : fx_t fPr,
190 : fa_t fAccept);
191 :
192 : // ---------- Public interface --------------------------
193 :
194 : /** Propose x as the next point in chain
195 : */
196 : bool propose(const v_t &x);
197 :
198 : void reset(const v_t &x, unsigned seed=0);
199 :
200 : };
201 :
202 : /** A chain which depends on current point and a likelihood value
203 : defined at the reset stage
204 : */
205 : class ILklChain:
206 : public ChainBase
207 : {
208 : public:
209 :
210 : /** \brief Function type to return the acceptance probability
211 :
212 : \param L is the likelihood supplied to reset
213 : */
214 : typedef std::function< double (double L, const MCPoint2 &c, const MCPoint2 &p) > fa_t;
215 :
216 : private:
217 :
218 : /// The initial likelihood
219 : double L;
220 :
221 : /// The acceptance function
222 : fa_t fAccept;
223 :
224 : public:
225 :
226 : // ---------- Construction / Destruction --------------
227 :
228 : ILklChain(const v_t &ic,
229 : fx_t fLkl,
230 : fx_t fPr,
231 : fa_t fAccept);
232 :
233 : // ---------- Public interface --------------------------
234 :
235 : /** Propose x as the next point in chain
236 : */
237 : bool propose(const v_t &x);
238 :
239 : void reset(const v_t &x,
240 : double L, unsigned seed=0);
241 :
242 : };
243 :
244 : bool normProp(ChainBase &c,
245 : const std::vector<double> &sigma);
246 :
247 : /** Propose in dimension i only
248 : */
249 : bool normProp(ChainBase &c,
250 : size_t i,
251 : double s);
252 :
253 : /** Propose based on a combination of eigenvectors
254 : */
255 : bool eigenProp(ChainBase &c,
256 : const std::vector<double> &eigvals,
257 : const std::vector<double> &eigvects
258 : );
259 :
260 :
261 : /** \brief The standard metropolis acceptance function
262 : */
263 : double metropolis(const MCPoint2 &c,
264 : const MCPoint2 &p);
265 :
266 : /** \brief The nested sampling acceptance function
267 : */
268 : double constrPrior(const MCPoint2 &c,
269 : const MCPoint2 &p);
270 :
271 :
272 : double constrPriorP(const MCPoint2 &f,
273 : const MCPoint2 &c,
274 : const MCPoint2 &p);
275 :
276 : /** Accept proportionally to prior mass but only if likelihood is
277 : greater than L
278 : */
279 : double constrPriorL(double L,
280 : const MCPoint2 &c,
281 : const MCPoint2 &p);
282 :
283 :
284 :
285 :
286 : }
287 :
288 : #endif
|
