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$
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$
  • $\begingroup$ Why the multiplication factor of $3$? $\endgroup$ – Amzoti Dec 15 '14 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$ – njpipeorgan Dec 15 '14 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 '14 at 8:44

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.