Monte Carlo integration with random numbers generated from a Gaussian distribution

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?

-
You can look up the examples here. – b.gatessucks Apr 1 '14 at 9:00
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. – ciao Apr 1 '14 at 10:08
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. – m_goldberg Apr 1 '14 at 11:08
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 ) – george2079 Apr 1 '14 at 12:48
@m_goldberg These work "PointGenerator" -> {Random, "Sobol", "Niederreiter", "Lattice"} – Dr. belisarius Feb 22 at 8:55

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["IntegrationNIntegrateUtilities"]

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}}}]


-
+1 I love to learn about these under documented methods from folks with inside knowledge! – Jason is no longer a postdoc Feb 23 at 17:57
@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. – Anton Antonov Feb 23 at 18:04
Wow! If only stuff like this customization was in the docs… – J. M. Feb 23 at 18:26
@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.) – Anton Antonov Feb 24 at 2:01
@AntonAntonov thank you for this wonderful answer and your blog "mathematicaforprediction" :) – ubpdqn Feb 24 at 4:14

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


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


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


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


-
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] σ)]. – J. M. Feb 23 at 14:24
And today I learned about the 2-argument version of Erf, thank you very much – Jason is no longer a postdoc Feb 23 at 14:28
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}]]. – J. M. Feb 23 at 14:46
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 – Jason is no longer a postdoc Feb 23 at 14:53
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. – J. M. Feb 23 at 14:57