Beefy Boxes and Bandwidth Generously Provided by pair Networks
Clear questions and runnable code
get the best and fastest answer
 
PerlMonks  

Re: Monte Carlo - Coin Toss

by roboticus (Canon)
on Mar 12, 2011 at 13:53 UTC ( #892824=note: print w/ replies, xml ) Need Help??


in reply to Monte Carlo - Coin Toss

James_H:

Obviously your "modest laptop" is a bit zippier than my Atom-based desktop, as it took 3:10 to complete. But when I looked at your code, I was wondering why you did a monte-carlo simulation. Since there are only 20 coins, there are a little over a million combinations, so why not compute it exactly?

So I whipped this up:

$ cat 892293.pl use strict; use warnings; my $numTosses = 20; # Coin tosses per experiment (20 in this example) my $runs = (1<<$numTosses)-1; my @tailCnt; for (my $collection=0; $collection<1<<$numTosses; ++$collection) { my $tails=0; $tails += ($collection & 1<<$_) ? 1 : 0 for 0 .. $numTosses; $tailCnt[$tails]++; } print <<EOHDR; Tails Count % ----- -------- ------ EOHDR for (my $i = 0; $i < $numTosses+1; $i++) { printf "% 4u % 8u %5.2f\n", $i, $tailCnt[$i], 100*$tailCnt[$i]/$runs; } $ time perl 892293.pl 0x00100000, 1048576 Tails Count % ----- -------- ------ 0 1 0.00 1 20 0.00 2 190 0.02 3 1140 0.11 4 4845 0.46 5 15504 1.48 6 38760 3.70 7 77520 7.39 8 125970 12.01 9 167960 16.02 10 184756 17.62 11 167960 16.02 12 125970 12.01 13 77520 7.39 14 38760 3.70 15 15504 1.48 16 4845 0.46 17 1140 0.11 18 190 0.02 19 20 0.00 20 1 0.00 real 0m20.372s user 0m20.273s sys 0m0.016s

It runs quite a bit faster on my machine, and it tried each combination once. Of course, I'm limited by the number of bits available, so large experiments aren't going to work well. You can use Bit::Vector, but I suspect that would be a good deal slower. (I wouldn't know, I didn't try it.)

...roboticus

When your only tool is a hammer, all problems look like your thumb.


Comment on Re: Monte Carlo - Coin Toss
Download Code
Re^2: Monte Carlo - Coin Toss
by BrowserUk (Pope) on Mar 12, 2011 at 14:38 UTC
    You can use Bit::Vector, but I suspect that would be a good deal slower.

    There wouldn't be any need to do that. With 32-bit ints, a full run would take around 7 hours. With 64-bit ints, it would take around 3 million years.

    If you switched to using unpack to count the bits:

    $tailCnt[ unpack '%32b*', pack 'N', $collection ]++;

    Then the full 32-bit comes down to a very reasonable 40 minutes, but that still leaves the full 64-bit requiring several hundred thousand years.


    Examine what is said, not who speaks -- Silence betokens consent -- Love the truth but pardon error.
    "Science is about questioning the status quo. Questioning authority".
    In the absence of evidence, opinion is indistinguishable from prejudice.

      BrowserUk:

      Heh, yep, that would be a good deal slower, all right! Of course, if speed gets to be an issue, you can always skip counting all the iterations, and directly compute the binomial expansion:

      use strict; use warnings; use List::Util qw(reduce); my $numTosses = 20; #my $runs = (1<<$numTosses)-1; my @triangle = (0, 1, 0); for (1 .. $numTosses) { my @newTriangle=(0); push @newTriangle, $triangle[$_]+$triangle[$_+1] for 0 .. $#triang +le-1; push @newTriangle, 0; @triangle = @newTriangle; } print <<EOHDR; Tails Count % ----- ---------- ------ EOHDR my $runs = reduce { $a + $b } @triangle; for (my $i = 0; $i < $numTosses+1; $i++) { printf "% 4u % 10u %5.2f\n", $i, $triangle[$i+1], 100*$triangle[$i+1]/$runs; }

      It runs a good bit faster:

      $ time perl 892293.pl Name "main::a" used only once: possible typo at 892293.pl line 21. Name "main::b" used only once: possible typo at 892293.pl line 21. Tails Count % ----- ---------- ------ 0 1 0.00 1 20 0.00 2 190 0.02 3 1140 0.11 4 4845 0.46 5 15504 1.48 6 38760 3.70 7 77520 7.39 8 125970 12.01 9 167960 16.02 10 184756 17.62 11 167960 16.02 12 125970 12.01 13 77520 7.39 14 38760 3.70 15 15504 1.48 16 4845 0.46 17 1140 0.11 18 190 0.02 19 20 0.00 20 1 0.00 real 0m0.034s user 0m0.024s sys 0m0.012s

      Even if you use 32 bits:

      $ time perl 892293.pl Name "main::a" used only once: possible typo at 892293.pl line 21. Name "main::b" used only once: possible typo at 892293.pl line 21. Tails Count % ----- ---------- ------ 0 1 0.00 1 32 0.00 2 496 0.00 3 4960 0.00 4 35960 0.00 5 201376 0.00 6 906192 0.02 7 3365856 0.08 8 10518300 0.24 9 28048800 0.65 10 64512240 1.50 11 129024480 3.00 12 225792840 5.26 13 347373600 8.09 14 471435600 10.98 15 565722720 13.17 16 601080390 13.99 17 565722720 13.17 18 471435600 10.98 19 347373600 8.09 20 225792840 5.26 21 129024480 3.00 22 64512240 1.50 23 28048800 0.65 24 10518300 0.24 25 3365856 0.08 26 906192 0.02 27 201376 0.00 28 35960 0.00 29 4960 0.00 30 496 0.00 31 32 0.00 32 1 0.00 real 0m0.034s user 0m0.028s sys 0m0.004s

      ...roboticus

      When your only tool is a hammer, all problems look like your thumb.

      Update: When I use bignum;, it prints the PDF for 64 tosses per run in just over a second, but I haven't got it formatted nicely...

      use strict; use warnings; use List::Util qw(reduce); use bignum; my $numTosses = 64; #my $runs = (1<<$numTosses)-1; my @triangle = (0, 1, 0); for (1 .. $numTosses) { my @newTriangle=(0); push @newTriangle, $triangle[$_]+$triangle[$_+1] for 0 .. $#triang +le-1; push @newTriangle, 0; @triangle = @newTriangle; } print <<EOHDR; Tails Count % ----- ---------- ------ EOHDR my $runs = reduce { $a + $b } @triangle; for (my $i = 0; $i < $numTosses+1; $i++) { print "$i\t$triangle[$i+1]\t",100*$triangle[$i+1]/$runs, "\n"; }

      I guess I'll have to figure out how to make bugnum and printf cooperate better...

        Here's a slightly speedier alternative to get Pascal's triangle row for the number of tosses:
        sub pascal_tri_row { my $r = shift; return () if $r < 0; my @row = (1) x ($r + 1); for my $i (1 .. $r - 1) { $row[$_] += $row[$_ - 1] for reverse 1 .. $i; } return @row; } sub triangle { my $numTosses = shift; my @triangle = (0, 1, 0); for (1 .. $numTosses) { my @newTriangle=(0); push @newTriangle, $triangle[$_]+$triangle[$_+1] for 0 .. $#tr +iangle-1; push @newTriangle, 0; @triangle = @newTriangle; } return @triangle; } use Benchmark qw(cmpthese); cmpthese -1, { robo_tri => sub { triangle(32) }, repel_tri => sub { 0, pascal_tri_row(32), 0 }, # 0's to match tria +ngle()'s return }; __END__ Rate robo_tri repel_tri robo_tri 1128/s -- -38% repel_tri 1811/s 60% --
Re^2: Monte Carlo - Coin Toss
by eye (Chaplain) on Mar 12, 2011 at 16:41 UTC
    ...find the probability density function of distribution (pdf) of tossing 20 coins at a time.

    Let me add a pedantic point. When asked to find a PDF (e.g., in a statistical theory class), the requester is likely to mean that the function should be found analytically. That said, it is equally valid to compute it exactly by enumeration (as roboticus did).

    Someone looking for a Monte Carlo approach will typically ask you to estimate the PDF.

Log In?
Username:
Password:

What's my password?
Create A New User
Node Status?
node history
Node Type: note [id://892824]
help
Chatterbox?
and the web crawler heard nothing...

How do I use this? | Other CB clients
Other Users?
Others romping around the Monastery: (6)
As of 2014-12-22 04:37 GMT
Sections?
Information?
Find Nodes?
Leftovers?
    Voting Booth?

    Is guessing a good strategy for surviving in the IT business?





    Results (110 votes), past polls