# Monte Carlo Volume Calculation and Speed Up

I am trying to do a Monte Carlo area calculation for an irregular area defined by:

$$0 \le x \le 1, 10 \le y \le 13, y \ge 12 \cos(x), y \ge 10 + x^3$$

I used the code with modifications (so I can actually read it) from: Finding the volume of a sphere using the Monte Carlo algorithm

  vol[num_] := Module[{hit, miss, index, x, y}, hit = 0; miss = 0;
For[index = 1, index <= num, ++index,
x = RandomReal[{0, 1}]; y = RandomReal[{10, 13}];
If[y >= 12 Cos[x] && y >= 10 + x^3, ++hit]];
hit/num]

Print["time and value...... :", Timing[N[vol[1000]]]]


It appears to be giving incorrect results (off by a factor of $\approx 3$):

$$\text{time and value...... :}\{0.015600,0.663\}$$

Any ideas why off by a factor of $\approx 3$ (the book claims $2.000346869$).

Can we easily print the hit and miss ratios?

Also, what is the best way to speed this up to do a several million trials? Certainly the code works for that value of $n$, but it is likely slow compared to what is possible.

Update

Curiously, if you look at: Cheney and Kincaid: Ice Cream Cone Problem, a similar program is giving different results. Here is the program in Fortran (http://www.ma.utexas.edu/CNA/cheney-kincaid/f90code/CHP13/cone.f90)

  vol3[num_] := Module[{hit, miss, index, x, y, z}, hit = 0; miss = 0;
For[index = 1,
index <= num, ++index, {x, y} = RandomReal[{-1, 1}, 2];
z = RandomReal[{0, 2}];
If[x^2 + y^2 <= z^2 && x^2 + y^2 <= z (2 - z), ++hit]];
8 hit/num]

Print["time and value...... :", Timing[N[vol3[10000]]]]


For loop is always slow.

You may try this:

f = Compile[{{x, _Real}, {y, _Real}}, If[y >= 12. Cos[x] && y >= 10 + x^3, 1, 0]];
vol[n_Integer /; n <= 10^6] := 3.* Total[f@@@Transpose@{RandomReal[{0, 1}, n], RandomReal[{10, 13}, n]}]/n;


The calculation of 1000000 samples takes 1.1 s on my i5-3210M.

• Why the multiplication factor of $3$? – Amzoti Dec 15 '14 at 5:54
• (hit/total) represents the probability that a random point falls onto the region. To obtain the area of the region, the probability has to be multiplied by the area of the domain [0, 1]*[10, 13]. – njpipeorgan Dec 15 '14 at 7:13
• If you add RuntimeAttributes -> Listable inside Compile, then you can use f[RandomReal[{0, 1}, n], RandomReal[{10, 13}, n]] instead of f@@@Transpose@…… and get a 7X speed up. (Parallelization -> True may also help a little. ) – xzczd Dec 15 '14 at 8:44