laziness, impatience, and hubris PerlMonks

### Re^3: Monte Carlo - Coin Toss

by roboticus (Chancellor)
 on Mar 12, 2011 at 16:11 UTC ( #892833=note: print w/replies, xml ) Need Help??

in reply to Re^2: Monte Carlo - Coin Toss
in thread Monte Carlo - Coin Toss

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...

Replies are listed 'Best First'.
Re^4: Monte Carlo - Coin Toss
by repellent (Priest) on Mar 13, 2011 at 04:31 UTC
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%        --

Sweet! So of course I had to take another look:

```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 robo_2 {
my \$row = shift;
my @cols = (1);
++\$row;
\$cols[\$_] = \$cols[\$_-1] * (\$row-\$_)/\$_ for 1 .. \$row-1;
return @cols;
}

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[1..\$#triangle-1];
}

use Benchmark qw(cmpthese);

print "robo_1: ", join(" ",triangle(8)), "\n";
print "repel1: ", join(" ",pascal_tri_row(8)), "\n";
print "robo_2: ", join(" ",robo_2(8)), "\n";

cmpthese -1, {
robo_tri  => sub { triangle(32) },
repel_tri => sub { pascal_tri_row(32) },
robo_2    => sub { robo_2(32) },
};

First I put in a trace so I could be sure I wasn't generating useless values. Next, I couldn't bear the zeroes you had to add to your function to match my original return value. They were just sentinel values to simplify the calculation, anyway. So I fixed the original to return only the values of interest. Finally, I had to squeeze a little more speed out of it:

```\$ perl 892898.pl
robo_1: 1 8 28 56 70 56 28 8 1
repel1: 1 8 28 56 70 56 28 8 1
robo_2: 1 8 28 56 70 56 28 8 1
Rate  robo_tri repel_tri    robo_2
robo_tri   2196/s        --      -42%      -93%
repel_tri  3775/s       72%        --      -88%
robo_2    31946/s     1355%      746%        --

You should've read a bit further down in the wikipedia article you linked. It had a much better algorithm for calculating the coefficients of any row. It saved a nested loop.

;^)

...roboticus

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

Wonderful! That is the one. Thanks, roboticus.

Create A New User
Node Status?
node history
Node Type: note [id://892833]
help
Chatterbox?

How do I use this? | Other CB clients
Other Users?
Others examining the Monastery: (9)
As of 2017-10-17 17:18 GMT
Sections?
Information?
Find Nodes?
Leftovers?
Voting Booth?
My fridge is mostly full of:

Results (235 votes). Check out past polls.

Notices?