I am working on doing some monte carlo simulations. I want to do a particular manipulation n times, but I want to constrain what I do based on three parameters, x, y, and z, which are probability distributions coded as arrays. For example, if I want to run this simulation 1000 times, then 24 should be xayaza, 24 xayazb, 72 xaybzc and so on. My code seems to work right when everything is integers, but not when some of the numbers are non-integers reals (always positive). I have tried a couple different strategies of rounding, ceiling, and floor, but it seems to always be off by a few. The correct solution should output $total = $trials;
I will offer 100 yootles to the first person who finds a solution for me. (I am including the perl code I have been playing around with, but my problem lies in the algorithm, not in the syntax.)
#!/usr/bin/perl -w
use strict;
use POSIX qw(floor ceil);
my $trials=795;
$trials = shift;
my @x = (.4,.4,.2);
my @y = (.3,.5,.2);
my @z = (.2,.2,.6);
my $trial =0;
my $total =1;
my $a=0;
while ($a<scalar @x && $trial <= round($trials*$x[$a])) {
my $b=0;
while ($b< scalar @y && $trial <= round($trials*$x[$a]*$y[$b])) {
my $c=0;
while ($c< scalar @z && $trial <= round($trials*$x[$a]*$y[$b]*$z[$c])) {
$total++;
$trial++;
if ($trial >= round($trials*$x[$a]*$y[$b]*$z[$c])) {
$c++;
$trial=0;
}
}
$b++;
$trial=0;
}
$a++;
$trial=0;
}
print "trials = $trials, total = $total\n";
sub round {
my($number) = shift;
return int($number + .5 * ($number <=> 0));
}
Update
My friend Danny Reeves, along with some help from David Yang solved my problem in a completely different way. Here is Danny’s solution:
#!/usr/bin/perl
# Rob's monstronsity that is surely the solution to the wrong problem.
# But for 100 yootles, we'll just do as we're told.
# This would be much nicer in a properly functional-style language!
my $trials = 795;
my @x = (.4,.4,.2);
my @y = (.3,.5,.2);
my @z = (.2,.2,.6);
@tuples = cross(\@x,\@y,\@z);
# compute idealized tuple counts:
@counts = map { $trials*prod(@$_) } @tuples;
@ic = map { int($_) } @counts; # floors of counts.
@fc = map { $_-int($_) } @counts; # fractional parts.
$f = sum(@fc); # sum of fractional parts, to redistribute.
# redistribute...
for($i=1; $i<=$f; $i++) {
$ic[posmax(deltas(\@counts,\@ic))]++;
}
$total = 0;
for($i=0; $i<scalar(@ic); $i++) {
($x, $y, $z) = @{$tuples[$i]};
for($j=0; $j<$ic[$i]; $j++) {
print "do something with ($x, $y, $z)\n";
$total++;
}
}
print "trials = $trials, total = $total\n";
# Return a cross product from its arguments. Arguments are array refs.
# Result is a list of array refs. [found this on the web; damn slick]
# (note that this returns the tuples in not quite canonical order)
sub cross {
my @r = [];
@r = map {my $s = $_; map {[@$_ => $s]} @r} @$_ for @_;
@r
}
# Return sum of arguments. Ie, reduce(+, args, 0).
sub sum { my $x = 0; for(@_) { $x += $_; } $x }
# Return product of arguments. Ie, reduce(*, args, 1).
sub prod { my $x = 1; for(@_) { $x *= $_; } $x }
# Return a list of differences between 2 lists, passed as refs.
# (assumes the lists have the same length)
sub deltas {
my @ans;
for(my $i=0; $i<scalar(@{$_[0]}); $i++) {
push(@ans, $_[0]->[$i] - $_[1]->[$i]);
}
@ans
}
# Takes list, return the position of the largest element.
sub posmax {
if (scalar(@_)==0) { return -1; }
my $x = 0; # index of best so far.
for(my $i=0; $i<scalar(@_); $i++) {
if($_[$i] > $_[$x]) { $x = $i; }
}
$x
}
And because Danny is a huge fan of Mathematica, he also included a Mathematica version
trials = 795;
x = {.4, .4, .2};
y = {.3, .5, .2};
z = {.2, .2, .6};
(* index of the largest element; i'm sure there's an adorable one-liner
for this if i thought hard enough *)
posmax[{}] = -1;
posmax[l_] := Module[{i, x = 1}, (* x is index of best so far *)
For[i = 1, i < Length[l], i++, If[l[[i]] > l[[x]], x = i]];
x]
tuples = Tuples[{x, y, z}];
counts = (trials*Times @@ # &) /@ tuples;
ic = IntegerPart /@ counts;
fc = FractionalPart /@ counts;
f = Total[fc];
For[i=1, i<=f, i++, ic[[posmax[counts - ic]]]++];
total = 0;
MapThread[Do[{"do something with ", #1}; total++, {#2}]&, {tuples,ic}];
total
