3
$\begingroup$

I have the following functions:

T = Compile[{n}, Boole[Sum[RandomReal[]^2, {n}] <= 1]];
Monte[d_, n_] := Sum[T[d], {n}]/n*2^d

T is a function that generates a random n-dimensional point in the n-dimensional unit hypercube and determines if it is within the n-dimensional unit hypersphere. Monte is a function that samples $n$ such points in $d$-dimenional space. For "small" numbers of points, Monte works fine:

Monte[2, 10^4]*1.0
(* 3.1244 *)

Please note that every time Monte is run the output will be different when Monte is working correctly, so the above output is just one of a number of possible outputs. Anyway: Monte breaks when I try inputting a larger value for $n$:

Monte[2, 10^7]*1.0
(* 4.*10^-7 Sum[1, {10000000}] *)

Why does this occur?

Edit: I tested it when T is declared as a regular function and not a compiled function:

T = Boole[Sum[RandomReal[]^2, {#}] <= 1] &;
Monte[d_, n_] := Sum[T[d], {n}]/n*2^d;
Monte[2, 10^7]*1.0
(* 4.*10^-7 Sum[1, {10000000}] *)
$\endgroup$
7
$\begingroup$

The help page for Sum (in v11.0 anyway) states "This example gives an unexpected result above the threshold value of 10^6:" but doesn't mention that threshold anywhere else on the page! Why is it THE threshold?

Anyway, their explanation is that the first argument is symbolically evaluated once the second argument is above that threshold. So that's why you see just the 1 in 4.*10^-7 Sum[1, {10000000}].

Wolfram's solution (workaround?) worked for me:

T = Compile[{n}, Boole[Sum[RandomReal[]^2, {n}] <= 1]];
Monte[d_, n_] := Sum[T[d], {n}, Method -> "Procedural"]/n*2^d
Monte[2, 10^7]*1.0
(* 3.14127 *)

and they give an alternative solution, which is to use _Integer after the variable definition to prevent symbolic evaluation, although that didn't work for me here. The reason is likely that T is not defined as a true function above.

| improve this answer | |
$\endgroup$
2
$\begingroup$

In the documentation for Sum there's a section about Possible Issues, where you'll find an example of a sum that gives unexpected results above the threshold value $10^6$.

Using your code you'll see that the strange behaviour, that, according to the documentation is due to symbolic evaluation of the first argument. Indeed the strange behaviour starts at the threshold value stated above:

Monte[2, 10^6] // N
Out[]:= 3.14004
Monte[2, 10^6+1]
Out[]:= (4 Sum[1, {1000001}])/1000001

The solution suggested in the documentation is to specify the method that sum should use as Method -> "Procedural", i.e.

Monte[d_, n_] := Sum[T[d], {n}, Method -> "Procedural"]/n*2^d;
Monte[2, 10^6 + 1]//N
Out[]:= 3.14013

Why exactly the symbolic evaluation results in what it does I don't know. Maybe someone else can elaborate on that.

About your approach: My first hunch was to use Total and Table instead of your Sum which is about 15 times faster than your approach. I'm sure there are even faster methods but Sum seems to me to be a bad choice.

Cheers

| improve this answer | |
$\endgroup$

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