# Number of realizations in NExpectation, Method Monte Carlo

I have been looking in NExpectation about how large the sample size is when Mathematica performs Monte Carlo simulation. For example, if we write

NExpectation[x Sin[x] + 1, x \[Distributed] NormalDistribution[], Method -> "MonteCarlo"]


which is the number of realizations of $x$ that Mathematica makes?

Also, I do not understand well how you can specify the number of realizations you want. For instance, if I want $1000$ realizations of $x$, what should I write?

• If you have looked at the documentation, you already know that NExpectation does integration of expression if the distribution is continuous, or discrete sum if this is a discrete one. In your case, it will do integration of $(x \sin(x)+1)$ with the NormalDistribution – José Antonio Díaz Navas Dec 19 '17 at 8:58
• @JoséAntonioDíazNavas But Monte Carlo is based on sampling: http://mathworld.wolfram.com/MonteCarloIntegration.html. I would like to know which $N$ Mathematica takes. – user54284 Dec 19 '17 at 9:01
• maybe you will be interested in this – José Antonio Díaz Navas Dec 19 '17 at 12:33

The number of samples taken by the "MonteCarlo" Method of NExpectation[] is controlled through an undocumented option. To see what happens, let's use a function with Sow[]:

f[x_?NumericQ] := (Sow[x]; x Sin[x] + 1)


and then evaluate

tst = Reap[NExpectation[f[x], x \[Distributed] NormalDistribution[],
Method -> "MonteCarlo"]];


We then find that

Length[tst[[-1, 1]]]
1012001


Now, let's adjust that undocumented option to use $10^7$ samples:

tst2 = Reap[NExpectation[f[x], x \[Distributed] NormalDistribution[],
Method -> {"MonteCarlo", "SamplingIncrement" -> 1*^7}]];


and we find that

Length[tst2[[-1, 1]]]
10012001


so, $10^7$ samples plus a little extra was used.

To see all the other adjustable options:

 Options[StatisticsLibraryNExpectationMonteCarloMethod]
{PrecisionGoal -> 2, AccuracyGoal -> 3, MaxIterations -> 50000,
"RandomSeed" -> Automatic, EvaluationMonitor -> None,
ConfidenceLevel -> 19/20, WorkingPrecision -> MachinePrecision,
"ReportingMethod" -> Automatic, "SamplingIncrement" -> 1000000}


You can use NIntegrate and specify for the method "MonteCarlo" the maximum number of points:

NIntegrate[(x*Sin[x] + 1) PDF[NormalDistribution[],
x], {x, -Infinity, Infinity},
Method -> {"MonteCarlo", "MaxPoints" -> 999}]

(* 1.63975 *)


This shows that 1000 points were used:

res =
Reap@NIntegrate[(x*Sin[x] + 1) PDF[NormalDistribution[],
x], {x, -Infinity, Infinity},
Method -> {"MonteCarlo", "MaxPoints" -> 999},
EvaluationMonitor :> Sow[x]];
Length@res[[2, 1]]

(* 1000 *)


Also, you might want to take a look into this answer of "Monte Carlo integration with random numbers generated from a Gaussian distribution".

Looking at official doc and NExpectation options:

Options[NExpectation]

{AccuracyGoal -> \[Infinity], Compiled -> Automatic,
Method -> Automatic, PrecisionGoal -> Automatic,
WorkingPrecision -> MachinePrecision, TargetUnits -> Automatic}


it is true that we do not find how to have a fine control (number of realizations etc.).

Concerning the number of realizations, here 100, you can use this trick:

NExpectation[Block[{$IterationLimit = 100}, x*Sin[x] + 1], x \[Distributed] NormalDistribution[], Method -> "MonteCarlo"]  1.60641 to be compared with the more accurate result: 1.60624 • Since NExpectation does not hold its arguments, Block[{$IterationLimit = 100}, x*Sin[x] + 1] evaluates to x*Sin[x] + 1 before it is passed to NExpectation, and the change to $IterationLimit is lost before it can have any effect. (I don't think $IterationLimit works in the way implied here. I'd be surprised if it is used in any other way than what is described in this tutorial.) – Michael E2 Mar 19 '18 at 12:36