|
Re: Monte Carlo approximation of PI
by robartes (Priest) on Apr 04, 2003 at 11:40 UTC
|
Very nice, but shouldn't it be rand()**2 + rand()**2 <= 1? Given the fact that the range of numbers returned from a random generator is discrete, given enough iterations, doing just < 1 would underestimate pi ever so slightly. Then again, nobody would probably run this algorithm long enough to see an effect from this :).
CU
Robartes- | [reply]
[d/l]
[select] |
|
|
It is probable that imperfections in Perl's rand will result in a consistent slight theoretical bias. It would also take an insane amount of time to have any chance of detecting this fact. Remember that to get the error down to 1 part in n, the number of iterations that you need is proportional to n*n. (The appropriate constant depends on what probability of being outside of your interval you are willing to accept.) A theoretical error on the order of one part in 50 million takes 2,500 trillion steps to find. Assuming that you can invoke the necessary code 100 million times per section (probably rather optimistic on today's hardware), it would take on the order of a year or so to find it.
If the error is smaller than that, it would take longer...
(I don't have sufficient interest to find out what Perl's rand algorithm is to calculate the order of magnitude of the exact theoretical bias. It will depend on how close together the grid of possible answers are that rand can return.)
| [reply] |
|
|
But if you want to factor in the discreteness, with the
same argumenation it follows that using <= 1
will lead to an overestimated pi.
Abigail
| [reply]
[d/l] |
|
|
Indeed, well spotted. So to make the algorithm even slower, we might actually have to run it twice, once with < 1 and once with <= 1, and then take the average of the two. Oh well :).
CU
Robartes-
| [reply]
[d/l]
[select] |
|
|
|
|
(jeffa) Re: Monte Carlo approximation of PI
by jeffa (Bishop) on Apr 04, 2003 at 20:08 UTC
|
Very nice Abigail-II, but not slow enough! ;) If you will
allow me to lift a few lines of your code, we can make it
slower:
use strict;
use warnings;
use POE;
my $in = 0;
my $count = shift || 1000;
POE::Session->create(
inline_states => {
_start => sub {
rand () ** 2 + rand () ** 2 < 1 && $in ++
},
},
) for 1 .. $count;
POE::Kernel->run();
printf "Pi equals %f\n", 4 * $in / $count;
UPDATE: yes, this is very much an attempt at humor. :O) No
disrespect meant to you or my new favorite
app framework, POE. ;)
jeffa
L-LL-L--L-LL-L--L-LL-L--
-R--R-RR-R--R-RR-R--R-RR
B--B--B--B--B--B--B--B--
H---H---H---H---H---H---
(the triplet paradiddle with high-hat)
| [reply]
[d/l] |
|
|
| [reply] |
|
|
| [reply] |
Re: Monte Carlo approximation of PI
by YuckFoo (Abbot) on Apr 04, 2003 at 18:00 UTC
|
With a little modification, a 1.0 approximator:
OneFoo
#!/usr/bin/perl
# Two Monte Carlo approximations of 1. Slow!
use strict;
use warnings;
my $in = 0;
my $count = @ARGV ? shift : 1000;
rand () + rand () < 1 && $in ++ for 1 .. $count;
printf "One equals %f\n", 2 * $in / $count;
$in = 0;
rand () < .5 && rand () < .5 && $in ++ for 1 .. $count;
printf "Or maybe %f\n", 4 * $in / $count;
| [reply]
[d/l] |