Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

Consider the following example:

4.*Mean@Table[ Boole[Random[]^2 + Random[]^2 < 1], {10^6}] // AbsoluteTiming
(*{0.408023, 3.14105}*)

With[{n = 10^6},
  4. Mean@Boole@Thread[(RandomReal[1, n]^2 + RandomReal[1, n]^2) < 1]] // AbsoluteTiming
(*{1.460084, 3.14276}*)

I expected the second method to be more efficient since it used vectorized operations, but it works slowly, and I don't know why. Is it possible to make it faster? I used Mathematica 9 on Windows 7.

share|improve this question

2 Answers 2

up vote 7 down vote accepted

First, the old Random function uses an inferior PRN generator and should not be used other than for legacy compatibility, if I recall correctly.

On my machine (version 7, Windows 7) the second method is faster than the first, though not by a great amount:

4.*Mean@Table[Boole[RandomReal[]^2 + RandomReal[]^2 < 1], {10^6}] // AbsoluteTiming

With[{n = 10^6}, 
  4. Mean@Boole@Thread[(RandomReal[1, n]^2 + RandomReal[1, n]^2) < 1]] // AbsoluteTiming

{1.2840735, 3.14154}

{0.4570262, 3.1401}

Numeric methods will be much faster:

With[{n = 10^6},
 With[{r = RandomReal[1, n]^2 + RandomReal[1, n]^2},
  4. * Mean[1 - UnitStep[r - 1]]
 ]
] // AbsoluteTiming

{0.0510029, 3.14232}

This is partly because the array is not unpacked, as it is in your code as can be seen by first running On["Packing"]:

On["Packing"]

With[{n = 10^6}, 
  4. Mean@Boole@Thread[(RandomReal[1, n]^2 + RandomReal[1, n]^2) < 1]] // AbsoluteTiming

Developer`FromPackedArray::unpack: Unpacking array in call to Less. >>

Developer`FromPackedArray::punpack1: Unpacking array with dimensions {1000000} to level 1. >>

share|improve this answer
1  
"First, the old Random function uses an inferior PRN generator and should not be used other than for legacy compatibility, if I recall correctly." - right, the old Marsaglia-Zaman method (which is now also available with SeedRandom[seed, Method -> "Legacy"]) has been shown to be statistically inferior to the current default method. The OP might also consider switching to a different PRNG if need be, using also SeedRandom[]. –  J. M. May 18 '13 at 13:47
    
Thank you. Is this right? With[{n = 10^6}, With[{r = RandomReal[1, n]^2 + RandomReal[1, n]^2}, 4.*Mean[UnitStep[1 - r]]]] –  expression May 19 '13 at 2:55
    
@explorer I believe you are asking: is UnitStep[1 - r] equivalent to 1 - UnitStep[r - 1]? Not precisely, because of the case of 1: {UnitStep[1 - r], 1 - UnitStep[r - 1]} /. r -> 1. It is unlikely to make a difference in your application however, IMHO, but to be rigorous I used the longer form. –  Mr.Wizard May 19 '13 at 5:49
    
@Mr.Wizard Thanks again for being so helpful. –  expression May 19 '13 at 8:14

Can't compete with Mr Wizard, but here's a vectorized version that's not too shabby:

n = 10^6;
-Mean[Sign@Clip[Total[RandomReal[1, {2, n}]^2] - 1, {-100, 0}]] // AbsoluteTiming

{0.133305, 3.14322}

share|improve this answer

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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