Theoretica
A C++ numerical and automatic mathematical library
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 
15 namespace 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_bisection(
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
A statically allocated N-dimensional vector with elements of the given type.
Definition: vec.h:88
theoretica::vec::resize
void resize(size_t n) const
Compatibility function to allow for allocation or resizing of dynamic vectors.
Definition: vec.h:412
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:188
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::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::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::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
theoretica::root_bisection
real root_bisection(RealFunction f, real a, real b, real tolerance=OPTIMIZATION_TOL)
Approximate a root of an arbitrary function using bisection inside a compact interval [a,...
Definition: roots.h:56
prng.h
Pseudorandom number generator class.
quasirandom.h
Quasi-random sequences.
sampling.h
Sampling from probability distributions.