1
$\begingroup$

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?

$\endgroup$
  • $\begingroup$ 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 $\endgroup$ – José Antonio Díaz Navas Dec 19 '17 at 8:58
  • $\begingroup$ @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. $\endgroup$ – user54284 Dec 19 '17 at 9:01
  • 1
    $\begingroup$ maybe you will be interested in this $\endgroup$ – José Antonio Díaz Navas Dec 19 '17 at 12:33
1
$\begingroup$

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[Statistics`Library`NExpectationMonteCarloMethod]
   {PrecisionGoal -> 2, AccuracyGoal -> 3, MaxIterations -> 50000, 
    "RandomSeed" -> Automatic, EvaluationMonitor -> None, 
    ConfidenceLevel -> 19/20, WorkingPrecision -> MachinePrecision, 
    "ReportingMethod" -> Automatic, "SamplingIncrement" -> 1000000}
$\endgroup$
1
$\begingroup$

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".

$\endgroup$
0
$\begingroup$

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

$\endgroup$
  • 1
    $\begingroup$ 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.) $\endgroup$ – Michael E2 Mar 19 '18 at 12:36

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.