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

1 Answer 1

Reset to default
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.

Improve this answer
$\endgroup$
3
  • $\begingroup$ Why the multiplication factor of $3$? $\endgroup$ Commented 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$ Commented 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$ Commented Dec 15, 2014 at 8:44

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.