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

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:

# 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

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