13
$\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$
6
  • 2
    $\begingroup$ You can look up the examples here. $\endgroup$ Commented Apr 1, 2014 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
    Commented Apr 1, 2014 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
    Commented Apr 1, 2014 at 11:08
  • 4
    $\begingroup$ @m_goldberg These work "PointGenerator" -> {Random, "Sobol", "Niederreiter", "Lattice"} $\endgroup$ Commented Feb 22, 2016 at 8:55
  • 1
    $\begingroup$ This is a very helpful discussion on Monte Carlo integration. I'm wondering if there are any recent updates on these methods. In particular, Antonov's answer is attractive given its leveraging of the NIntegrate function, while the clarity of Jason's answer is very simple to understand. Have any advances in Mathematica made these implements with NIntegrate simpler? $\endgroup$
    – Kenric
    Commented Mar 11, 2021 at 11:59

2 Answers 2

15
$\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$
10
  • 3
    $\begingroup$ +1 I love to learn about these under documented methods from folks with inside knowledge! $\endgroup$
    – Jason B.
    Commented Feb 23, 2016 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$ Commented Feb 23, 2016 at 18:04
  • 3
    $\begingroup$ Wow! If only stuff like this customization was in the docs… $\endgroup$ Commented Feb 23, 2016 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$ Commented Feb 24, 2016 at 2:01
  • 1
    $\begingroup$ @Jordan I will try to answer in the next few days. $\endgroup$ Commented Apr 17, 2017 at 1:36
7
$\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$
6
  • 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$ Commented Feb 23, 2016 at 14:24
  • $\begingroup$ And today I learned about the 2-argument version of Erf, thank you very much $\endgroup$
    – Jason B.
    Commented Feb 23, 2016 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$ Commented Feb 23, 2016 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.
    Commented Feb 23, 2016 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$ Commented Feb 23, 2016 at 14:57

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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