Theoretica
A C++ numerical and automatic mathematical library
Loading...
Searching...
No Matches
montecarlo.h
Go to the documentation of this file.
1
5
6#ifndef THEORETICA_MONTECARLO_H
7#define THEORETICA_MONTECARLO_H
8
9#include "./prng.h"
10#include "./quasirandom.h"
11#include "./sampling.h"
12#include "../algebra/algebra_types.h"
13
14
15namespace theoretica {
16
17
24 inline real integral_crude(
25 real_function f,
26 real a, real b,
27 PRNG& g, unsigned int N = 1000) {
28
29 real sum_y = 0;
30
31 for (unsigned int i = 0; i < N; ++i)
32 sum_y += f(rand_uniform(a, b, g));
33
34 return (b - a) * sum_y / static_cast<real>(N);
35 }
36
37
43 template<unsigned int S>
44 inline real integral_crude(
45 real(*f)(vec<real, S>), vec<vec2, S> extremes,
46 PRNG& g, unsigned int N = 1000) {
47
48 real sum_y = 0;
49
50 // Sample the function at random points in the integration region
51 for (unsigned int i = 0; i < N; ++i) {
52
53 vec<real, S> v;
54 for (unsigned int k = 0; k < S; ++k)
55 v[k] = rand_uniform(extremes[k][0], extremes[k][1], g);
56
57 sum_y += f(v);
58 }
59
60 // Volume of the integration region
61 real volume = 1;
62 for (unsigned int i = 0; i < S; ++i)
63 volume *= (extremes[i][1] - extremes[i][0]);
64
65 return volume * sum_y / static_cast<real>(N);
66 }
67
68
76 inline real integral_quasi_crude(
77 real_function f,
78 real a, real b,
79 unsigned int N = 1000) {
80
81 real sum_y = 0;
82
83 for (unsigned int i = 0; i < N; ++i)
84 sum_y += f(a + qrand_weyl(i) * (b - a));
85
86 return (b - a) * sum_y / static_cast<real>(N);
87 }
88
89
97 template<unsigned int S>
98 inline real integral_quasi_crude(
99 real(*f)(vec<real, S>), vec<vec2, S> extremes,
100 unsigned int N, vec<real, S> alpha) {
101
102 real sum_y = 0;
103
104 for (unsigned int i = 0; i < N; ++i) {
105
106 vec<real, S> v;
107 for (unsigned int k = 0; k < S; ++k)
108 v[k] = extremes[k][0]
109 + (qrand_weyl(i, alpha[k])
110 * (extremes[k][1] - extremes[k][0]));
111
112 sum_y += f(v);
113 }
114
115 real volume = 1;
116 for (unsigned int i = 0; i < S; ++i)
117 volume *= (extremes[i][1] - extremes[i][0]);
118
119 return volume * sum_y / static_cast<real>(N);
120 }
121
122
130 template<unsigned int S>
131 inline real integral_quasi_crude(
132 real(*f)(vec<real, S>), vec<vec2, S> extremes,
133 unsigned int N = 1000, real alpha = 0) {
134
135 if(alpha == 0) {
136
137 // Find the only positive root of the polynomial
138 // x^s+1 - x - 1 = 0
139
140 alpha = 1.0 / root_bisect(
141 [](real x) {
142 return pow(x, S + 1) - x - 1;
143 }, 0, 2);
144 }
145
146 vec<real, S> v;
147 for (unsigned int i = 0; i < S; ++i)
148 v[i] = pow(alpha, i + 1);
149
150 return integral_quasi_crude(f, extremes, N, v);
151 }
152
153
162 inline real integral_hom(
163 real_function f,
164 real a, real b, real c, real d,
165 PRNG& g, unsigned int N = 1000) {
166
167 int N_inside = 0;
168
169 for (unsigned int i = 0; i < N; ++i) {
170
171 const real x_n = rand_uniform(a, b, g);
172 const real y_n = rand_uniform(c, d, g);
173
174 const real f_x = f(x_n);
175
176 if(f_x >= 0) {
177 if(f_x >= y_n)
178 N_inside++;
179 } else {
180 if(f_x < y_n)
181 N_inside--;
182 }
183 }
184
185 return (N_inside / static_cast<real>(N)) * (b - a) * (d - c);
186 }
187
188
200 inline real integral_hom(
201 real_function f,
202 real a, real b, real f_max,
203 PRNG& g, unsigned int N = 1000) {
204
205 unsigned int N_inside = 0;
206
207 for (unsigned int i = 0; i < N; ++i) {
208
209 const real x_n = rand_uniform(a, b, g);
210 const real y_n = rand_uniform(0, f_max, g);
211
212 if(f(x_n) > y_n)
213 N_inside++;
214 }
215
216 return (N_inside / static_cast<real>(N)) * (b - a) * f_max;
217 }
218
219
228 inline real integral_quasi_hom(
229 real_function f,
230 real a, real b, real c, real d,
231 unsigned int N = 1000) {
232
233 int N_inside = 0;
234
235 for (unsigned int i = 0; i < N; ++i) {
236
237 vec2 v = qrand_weyl2(i);
238
239 const real x_n = a + (b - a) * v[0];
240 const real y_n = c + (d - c) * v[1];
241
242 const real f_x = f(x_n);
243
244 if(f_x >= 0) {
245 if(f_x >= y_n)
246 N_inside++;
247 } else {
248 if(f_x < y_n)
249 N_inside--;
250 }
251 }
252
253 return (N_inside / static_cast<real>(N)) * (b - a) * (d - c);
254 }
255
256
265 inline real integral_quasi_hom(
266 real_function f,
267 real a, real b, real f_max,
268 unsigned int N = 1000) {
269
270 unsigned int N_inside = 0;
271
272 for (unsigned int i = 0; i < N; ++i) {
273
274 vec2 v = qrand_weyl2(i);
275
276 const real x_n = a + (b - a) * v[0];
277 const real y_n = v[1] * f_max;
278
279 if(f(x_n) > y_n)
280 N_inside++;
281 }
282
283 return (N_inside / static_cast<real>(N)) * (b - a) * f_max;
284 }
285
286
297 inline real integral_hom_2d(
298 real(*f)(real, real),
299 real a, real b,
300 real c, real d,
301 real f_max, PRNG& g,
302 unsigned int N = 1000) {
303
304 unsigned int N_inside = 0;
305
306 for (unsigned int i = 0; i < N; ++i) {
307
308 real x = rand_uniform(a, b, g);
309 real y = rand_uniform(c, d, g);
310 real z = rand_uniform(0, f_max, g);
311
312 if(f(x, y) > z)
313 N_inside++;
314 }
315
316 return (N_inside / static_cast<real>(N)) * (b - a) * (d - c) * f_max;
317 }
318
319
327 inline real integral_impsamp(
328 real_function f, real_function g, real_function Ginv,
329 PRNG& gen, unsigned int N = 1000) {
330
331 real sum_y = 0;
332
333 for (unsigned int i = 0; i < N; ++i) {
334 const real z = Ginv(rand_uniform(0, 1, gen));
335 sum_y += f(z) / g(z);
336 }
337
338 return sum_y / static_cast<real>(N);
339 }
340
341
349 inline real integral_quasi_impsamp(
350 real_function f, real_function g, real_function Ginv,
351 unsigned int N = 1000) {
352
353 real sum_y = 0;
354
355 for (unsigned int i = 0; i < N; ++i) {
356 const real z = Ginv(qrand_weyl(i + 1));
357 sum_y += f(z) / g(z);
358 }
359
360 return sum_y / static_cast<real>(N);
361 }
362
363
372 template <
373 typename Vector = std::vector<real>,
374 typename Function = std::function<real(vec<real>)>
375 >
376 Vector sample_mc(Function f, std::vector<pdf_sampler>& rv, unsigned int N) {
377
378 Vector sample;
379 sample.resize(N);
380
381 vec<real> v;
382 v.resize(rv.size());
383
384 for (unsigned int i = 0; i < N; ++i) {
385
386 for (unsigned int j = 0; j < rv.size(); ++j)
387 v[j] = rv[j]();
388
389 sample[i] = f(v);
390 }
391
392 return sample;
393 }
394
395}
396
397
398#endif
theoretica::PRNG
A pseudorandom number generator.
Definition prng.h:18
theoretica::vec::resize
void resize(size_t n) const
Compatibility function to allow for allocation or resizing of dynamic vectors.
Definition vec.h:416
theoretica
Main namespace of the library which contains all functions and objects.
Definition algebra.h:27
theoretica::real
double real
A real number, defined as a floating point type.
Definition constants.h:198
theoretica::integral_crude
real integral_crude(real_function f, real a, real b, PRNG &g, unsigned int N=1000)
Approximate an integral by using Crude Monte Carlo integration.
Definition montecarlo.h:24
theoretica::integral_quasi_impsamp
real integral_quasi_impsamp(real_function f, real_function g, real_function Ginv, unsigned int N=1000)
Approximate an integral by using Crude Quasi-Monte Carlo integration with importance sampling,...
Definition montecarlo.h:349
theoretica::integral_quasi_crude
real integral_quasi_crude(real_function f, real a, real b, unsigned int N=1000)
Approximate an integral by using Crude Quasi-Monte Carlo integration by sampling from the Weyl sequen...
Definition montecarlo.h:76
theoretica::integral_hom_2d
real integral_hom_2d(real(*f)(real, real), real a, real b, real c, real d, real f_max, PRNG &g, unsigned int N=1000)
Use the Hit-or-Miss Monte Carlo method to approximate a double integral.
Definition montecarlo.h:297
theoretica::vec2
vec< real, 2 > vec2
A 2-dimensional vector with real elements.
Definition algebra_types.h:39
theoretica::integral_quasi_hom
real integral_quasi_hom(real_function f, real a, real b, real c, real d, unsigned int N=1000)
Approximate an integral by using Hit-or-miss Quasi-Monte Carlo integration, sampling points from the ...
Definition montecarlo.h:228
theoretica::vector_element_t
std::remove_reference_t< decltype(std::declval< Structure >()[0])> vector_element_t
Extract the type of a vector (or any indexable container) from its operator[].
Definition core_traits.h:134
theoretica::integral_hom
real integral_hom(real_function f, real a, real b, real c, real d, PRNG &g, unsigned int N=1000)
Approximate an integral by using Hit-or-miss Monte Carlo integration.
Definition montecarlo.h:162
theoretica::integral_impsamp
real integral_impsamp(real_function f, real_function g, real_function Ginv, PRNG &gen, unsigned int N=1000)
Approximate an integral by using Crude Monte Carlo integration with importance sampling.
Definition montecarlo.h:327
theoretica::rand_uniform
real rand_uniform(real a, real b, PRNG &g, uint64_t prec=STATISTICS_RAND_PREC)
Generate a pseudorandom real number in [a, b] using a preexisting generator.
Definition sampling.h:33
theoretica::sample_mc
Vector sample_mc(Function f, std::vector< pdf_sampler > &rv, unsigned int N)
Generate a Monte Carlo sample of values of a given function of arbitrary variables following the give...
Definition montecarlo.h:376
theoretica::root_bisect
real root_bisect(RealFunction f, real a, real b, real tol=OPTIMIZATION_TOL, unsigned int max_iter=OPTIMIZATION_BISECTION_ITER)
Find the root of a univariate real function using bisection inside a compact interval [a,...
Definition roots.h:58
theoretica::qrand_weyl
real qrand_weyl(unsigned int n, real alpha=INVPHI)
Weyl quasi-random sequence.
Definition quasirandom.h:24
theoretica::qrand_weyl2
vec2 qrand_weyl2(unsigned int n, real alpha=0.7548776662466927)
Weyl quasi-random sequence in 2 dimensions.
Definition quasirandom.h:70
theoretica::real_function
std::function< real(real)> real_function
Function pointer to a real function of real variable.
Definition function.h:20
theoretica::pow
dual2 pow(dual2 x, int n)
Compute the n-th power of a second order dual number.
Definition dual2_functions.h:41
prng.h
Pseudorandom number generator class.
quasirandom.h
Quasi-random sequences.
sampling.h
Sampling from probability distributions.