9
$\begingroup$

I want to do numerical integration of some functions using the Monte Carlo method. The default setting for the Monte Carlo method is to use a uniform distribution as far as I know.

How can I change the random number generation to Gaussian when doing numerical integration?

$\endgroup$
  • 2
    $\begingroup$ You can look up the examples here. $\endgroup$ – b.gatessucks Apr 1 '14 at 9:00
  • 1
    $\begingroup$ If you can figure out the PointGenerator option for MonteCarloRule, I think that's a path. Zero documentation of it (at least as far as I've seen) other than the mention in the advanced numerical integration tutorial docs. $\endgroup$ – ciao Apr 1 '14 at 10:08
  • $\begingroup$ Attempting to use the mentioned but undocumented sub-option "PointGenerator" of option Method of NIntegrate can be dangerous to your kernel. At least, my attempts to set this sub-option anything other than its default resulted in a kernel crash. $\endgroup$ – m_goldberg Apr 1 '14 at 11:08
  • $\begingroup$ why don't you implement it yourself? There is a mathematica example on the wiki page by the way en.wikipedia.org/wiki/Monte_Carlo_integration. (I didn't study it but it looks like a good start ) $\endgroup$ – george2079 Apr 1 '14 at 12:48
  • 4
    $\begingroup$ @m_goldberg These work "PointGenerator" -> {Random, "Sobol", "Niederreiter", "Lattice"} $\endgroup$ – Dr. belisarius Feb 22 '16 at 8:55
12
$\begingroup$

Because of the comment by @ciao I think it is a good idea to give a solution/answer using MonteCarloRule and PointGenerator.

PointGenerator is an object that expects a sub-value function definition for "MakePoint" with the signature:

"MakePoint"[dim_, totalNumberOfPoints_, i_, wprec_]

where the arguments are:

  1. dim - integration dimension,
  2. totalNumberOfPoints - desired total number of points per rule application,
  3. i - point index,
  4. wprec - working precision.

Note that "MakePoint" provides very fine granularity, it makes only one point. For optimization purposes custom definitions (like the one given below) have to provide some sort of optimization for generating a bulk of points. The consequence of such an approach is that if AdaptiveMonteCarlo is used then the sampling points of each application of MonteCarloRule are going to be rescaled versions of the same set of points.

Here are definitions that make PointGenerator work with a general distribution argument.

ClearAll[MCAbscGenerator, MCDistributionPoints];
MCDistributionPoints[distr_, dim_, totalNumberOfPoints_] :=
  MCDistributionPoints[distr, dim, totalNumberOfPoints] =
   Block[{points = RandomVariate[distr, {totalNumberOfPoints, dim}]},
    Transpose[
     Map[Rescale[#, {Min[#], Max[#]}, {0, 1}] &, 
      Transpose[points]]]];
MCAbscGenerator["MakePoint"[dim_, totalNumberOfPoints_, n_, wprec_]] :=
    MCDistributionPoints[NormalDistribution[], dim, 
    totalNumberOfPoints][[n]];
MCAbscGenerator[distr_][
  "MakePoint"[dim_, totalNumberOfPoints_, i_, wprec_]] :=
 MCDistributionPoints[distr, dim, totalNumberOfPoints][[i]];

(The definitions are probably not very didactic, but I think they provide good general properties: signatures and funcitonality.)

Here is how we apply the custom point generator with MonteCarloRule:

NIntegrate[x^2 + y^3, {x, 0, 1}, {y, 0, 1}, PrecisionGoal -> 2,
   Method -> {"MonteCarlo", 
   Method -> {"MonteCarloRule", 
     "PointGenerator" -> 
      MCAbscGenerator[TruncatedDistribution[{-1, 1}, NormalDistribution[]]]}}]

(* 0.561876 *)

The real value is 7/12:

Integrate[x^2 + y^3, {x, 0, 1}, {y, 0, 1}] // N

(* 0.583333 *)

Let us visualize the sampling points for different truncated versions of the Normal Distribution and different Monte-Carlo methods. Notice that the integral gets more underestimated for the wider truncated distribution ranges.

Needs["Integration`NIntegrateUtilities`"]

ds = {TruncatedDistribution[{-1, 1}, NormalDistribution[]], 
   TruncatedDistribution[{-2, 2}, NormalDistribution[]], 
   TruncatedDistribution[{-3, 3}, NormalDistribution[]], NormalDistribution[]};
Grid[Riffle[Table[{Style[distr, Blue], SpanFromLeft}, {distr, ds}], 
  Table[(est =
     NIntegrate[x^2 + y^3, {x, 0, 1}, {y, 0, 1}, PrecisionGoal -> 2, 
      Method -> {meth, 
        Method -> {"MonteCarloRule", 
          "PointGenerator" -> MCAbscGenerator[distr]}}];
    gr = NIntegrateSamplingPoints@
      NIntegrate[x^2 + y^3, {x, 0, 1}, {y, 0, 1}, PrecisionGoal -> 2, 
       Method -> {meth, 
         Method -> {"MonteCarloRule", 
           "PointGenerator" -> MCAbscGenerator[distr]}}];
    Labeled[Append[gr, Frame -> True], 
     Grid[{{"method:", meth}, {"estimate:", est}}]]
    ), {distr, ds}, {meth, {"MonteCarlo", "AdaptiveMonteCarlo"}}]], 
 Dividers -> {None, {{True, False}}}]

enter image description here

$\endgroup$
  • 3
    $\begingroup$ +1 I love to learn about these under documented methods from folks with inside knowledge! $\endgroup$ – Jason B. Feb 23 '16 at 17:57
  • $\begingroup$ @JasonB Thanks. MonteCarloRule and PointGenerator were designed/developed before the functions RandomVariate and TruncatedDistribution were introduced. So, PointGenerator was designed with Sobol and Hammersley sequence generators in mind. I also considered Importance Sampling Integration to be easily implemented through PointGenerator. $\endgroup$ – Anton Antonov Feb 23 '16 at 18:04
  • 3
    $\begingroup$ Wow! If only stuff like this customization was in the docs… $\endgroup$ – J. M. is away Feb 23 '16 at 18:26
  • $\begingroup$ @J.M. Thanks. I did prepare a few documents for method top level programming within NIntegrate's framework, but other things became more important. Specifically, for PointGenerator I did not document it because I was not sure how good its design is. At that point I was only sure that it is a good fit for the pseudo random algorithms of Sobol, Niederreiter, Hammersley, and high dimensional lattices. (And of course RandomReal.) $\endgroup$ – Anton Antonov Feb 24 '16 at 2:01
  • 1
    $\begingroup$ @Jordan I will try to answer in the next few days. $\endgroup$ – Anton Antonov Apr 17 '17 at 1:36
5
$\begingroup$

So I took the code that is on the Wikipedia page and corrected one problem (why it had 1.1*x - 0.1 instead of x is beyond me, but it throws the answer way off), and then rewrote it to be more general and (to me) easier to read.

Here is the function,

monteCarloGaussian[func_, x1_, x2_, n_, mean_: 0, σ_: 1] := 
 1/2 (Erf[(mean - x2)/(Sqrt[2] σ), (mean - x1)/(
    Sqrt[2] σ)]) 
    Mean[(func[#]/(E^(-((-mean + #)^2/(2 σ^2)))/(
       Sqrt[2 π] σ)) ) & /@ 
    RandomVariate[
     TruncatedDistribution[{x1, x2}, 
      NormalDistribution[mean, σ]], n]]

The arguments are a pure function to be integrated, the limits of integration, the number of points to sample, and optionally the mean and standard deviation of the underlying normal distribution. The prefactor is the integral of the normal distribution over the given range, and the integral is the sum of the function evaluated at the random points divided by their probability density at that point.

You can evaluate it for some functions and see how large a sample is needed for convergence,

exact = NIntegrate[1/(1 + Sinh[2 x] Log[x]^2), {x, 0, 2}]; 
Table[{exact, 
    monteCarloGaussian[Function[x, 1/(1 + Sinh[2 x] Log[x]^2)], 0, 
     2.0, n]}, {n, 10^Range[1, 6]}] // Transpose // ListLinePlot

enter image description here

When the limits of integration are far from zero, it is necessary to change the mean and standard deviation of the distribution:

exact = NIntegrate[1/(1 + Sinh[2 x/4] Log[x/2]^2), {x, 5, 10}]; 
Table[{exact, 
    monteCarloGaussian[Function[x, 1/(1 + Sinh[2 x/4] Log[x/2]^2)], 5,
      10.0, n]}, {n, 10^Range[1, 6]}] // Transpose // ListLinePlot

enter image description here

that looks bad, but this

exact = NIntegrate[1/(1 + Sinh[2 x/4] Log[x/2]^2), {x, 5, 10}]; 
Table[{exact, 
    monteCarloGaussian[Function[x, 1/(1 + Sinh[2 x/4] Log[x/2]^2)], 5,
      10.0, n, 7.5, 5]}, {n, 10^Range[1, 6]}] // 
  Transpose // ListLinePlot

enter image description here

is much better.

It should be straightforward to extend this for functions of multiple variables, and I can do that if there is interest.

The above function is optimized to use a normal distribution, but we could modify it to use any distribution,

monteCarloDist[func_, x1_, x2_, n_, dist_: NormalDistribution[0, 1]] :=
 Module[{int, vars, probs},
  int = Integrate[PDF[dist, xx], {xx, x1, x2}];
  vars = RandomVariate[
    TruncatedDistribution[{x1, x2},
     dist], n];
  probs = PDF[dist, vars];
   int  Mean[(func[#1]/#2 ) & @@@ Transpose[{vars, probs}]]
  ]

With this you can use any of the built-in or user-defined distribution functions. It is a little slower than the one above for a NormalDistribution since it has to work out the integral on its own and it has to get the probabilities on its own, but not too much slower. The more complicated you make your distribution, the longer it takes,

exact = NIntegrate[1/(1 + Sinh[2 x] Log[x]^2), {x, 0, 2}];
Table[{exact, 
    monteCarloDist[Function[x, 1/(1 + Sinh[2 x] Log[x]^2)], 0, 2.0, n,
      ProbabilityDistribution[Sqrt[2]/Pi/(1 + x^4), {x, -4, 4}]]}, {n,
     10^Range[1, 6]}] // Transpose // ListLinePlot

enter image description here

$\endgroup$
  • 1
    $\begingroup$ It is more compact and more numerically stable to use Erf[(mean - x2)/(Sqrt[2] σ), (mean - x1)/(Sqrt[2] σ)] instead of Erf[(mean - x1)/(Sqrt[2] σ)] - Erf[(mean - x2)/(Sqrt[2] σ)]. $\endgroup$ – J. M. is away Feb 23 '16 at 14:24
  • $\begingroup$ And today I learned about the 2-argument version of Erf, thank you very much $\endgroup$ – Jason B. Feb 23 '16 at 14:28
  • $\begingroup$ With respect to your more general monteCarloDist: why not int = CDF[dist, x2] - CDF[dist, x1]? Also, the final result can be done as int Mean[(func[#1]/#2) & @@@ Transpose[{vars, pdfs}]]. $\endgroup$ – J. M. is away Feb 23 '16 at 14:46
  • $\begingroup$ It would be easier if there were a sandbox on here I could just show these to you in first before I post. :-) So in the case of the normal distribution, the CDF is a couple of order of magnitude faster, but for complicated user-defined definitions it is slower $\endgroup$ – Jason B. Feb 23 '16 at 14:53
  • $\begingroup$ I'd think that since CDF[] does table lookup (at least for the known distributions) instead of direct integration, it'd perform much faster... but I guess no advantage is accrued for arbitrary ProbabilityDistribution[] objects. $\endgroup$ – J. M. is away Feb 23 '16 at 14:57

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.