4
$\begingroup$

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]]]]
$\endgroup$

1 Answer 1

5
$\begingroup$

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.

$\endgroup$
3
  • $\begingroup$ Why the multiplication factor of $3$? $\endgroup$
    – Amzoti
    Dec 15, 2014 at 5:54
  • 2
    $\begingroup$ (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]. $\endgroup$ Dec 15, 2014 at 7:13
  • 4
    $\begingroup$ 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. ) $\endgroup$
    – xzczd
    Dec 15, 2014 at 8:44

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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