Calculating odds distribution with 6-sided dice

I'm trying to calculate the odds distribution of a changing number of 6-sided die rolls. For example, 3d6 ranges from 3 to 18 as follows:

3:1, 4:3, 5:6, 6:10, 7:15, 8:21, 9:25, 10:27, 11:27, 12:25, 13:21, 14:15, 15:10, 16:6, 17:3, 18:1

I wrote this php program to calculate it:

function distributionCalc($numberDice,$sides=6) {
for ( $i=0; $i<pow($sides,$numberDice); $i++)
    {
    $sum=0;
    for  ($j=0; $j<$numberDice; $j++)
        { $sum+=(1+(floor($i/pow($sides,$j))) % $sides); }
    $distribution[$sum]++;
    }
return $distribution;
}

The inner $j for-loop uses the magic of the floor and modulus functions to create a base-6 counting sequence with the number of digits being the number of dice, so 3d6 would count as:

111,112,113,114,115,116,121,122,123,124,125,126,131,etc.

The function takes the sum of each, so it would read as: 3,4,5,6,7,8,4,5,6,7,8,9,5,etc. It plows through all 6^3 possible results and adds 1 to the corresponding slot in the $distribution array between 3 and 18. Pretty straightforward. However, it only works until about 8d6, afterward i get server time-outs because it's now doing billions of calculations.

But I don't think it's necessary because die probability follows a sweet bell-curve distribution. I'm wondering if there's a way to skip the number crunching and go straight to the curve itself. Is there a way to do this with, for example, 80d6 (range: 80-480)? Can the distribution be projected without doing 6^80 calculations?

I'm not a professional coder and probability is still new to me, so thanks for 开发者_运维知识库all the help!

Stephen

In PERL:

#!
my( $DieType, $NumDice, $Loaded ) = @ARGV;

my $subname = "D" . $DieType . ( ( $Loaded eq "Loaded" ) ? "Loaded" : "Normal" );
my $Prob = \&$subname;

my $width = 12;
my $precision = $width - 2;

printf "%5s  %-${width}s \n", "Pip:", "Frequency:"; 
for ( my $j = $NumDice; $j <= $DieType * $NumDice ; $j++ ) {
  printf "%5d  %${width}.${precision}f \n", $j, Frequency( $DieType, $NumDice, $j );
}

sub D6Normal {
  my $retval = 1/6;
}

sub D6Loaded {
  my $retval = 1/6;

  CASE: for ($_[0]) {
    /1/    && do { $retval -= 0.02/6;   last CASE; }; 
    /2..5/ && do { $retval += 0.0025/6; last CASE; }; 
    /6/    && do { $retval += 0.01/6;   last CASE; }; 
  }
  return $retval;
}

sub D8Normal {
  my $retval = 1/8;
}

sub D10Normal {
  my $retval = 1/10;
}

sub D10Loaded {
  my $retval = 1/10;

  CASE: for ($_[0]) {
    /1..8/ && do { last CASE; }; 
    /9/    && do { $retval -= 0.01/10;  last CASE; }; 
    /10/   && do { $retval += 0.01/10;  last CASE; }; 
  }
  return $retval;
}

sub D12Normal {
  my $retval = 1/12;
}

sub D20Normal {
  my $retval = 1/20;
}

sub D32Normal {
  my $retval = 1/32;
}

sub D100Normal {
  my $retval = 1/100;
}

sub Frequency {
  my( $DieType, $NumberofDice, $PipCount ) = @_;

  if ( ( $PipCount > ($DieType * $NumberofDice) ) || ( $PipCount < $NumberofDice ) ) { 
    return 0; 
  }

  if ( ! exists $Freq{$NumberofDice}{$PipCount} ) {
    if ( $NumberofDice > 1 ) {
      for ( my $i = max( 1, $PipCount - $DieType ); $i <= min( $DieType * ($NumberofDice - 1), $PipCount - 1 ); $i++ ) {
        $Freq{$NumberofDice}{$PipCount} += &$Prob( $PipCount - $i ) * Frequency( $DieType, $NumberofDice - 1, $i );
      }
    } else {
      $Freq{$NumberofDice}{$PipCount} = &$Prob( $PipCount );
    }
  }
  return $Freq{$NumberofDice}{$PipCount}; 
}

sub max {
  my $max = shift(@_);
  foreach my $arg (@_) {
    $max = $arg if $max < $arg;
  }
  return $max;
}

sub min {
  my $min = shift(@_);
  foreach my $arg (@_) {
    $min = $arg if $min > $arg;
  }
  return $min;
}

You are looking for a Binomial Distribution

Ok, so let's start with rolling just one die. We know that the mean is 3.5. We can also calculate the variance,

sum(p(x) * (x - M)^2), where M is the mean, x is a dice result, and p is the probability of that dice result.

Using this formula, the variance of a single dice roll is 35/12 = 1/6*((-2.5)^2 + (-1.5)^2 + (-0.5)^2 + 0.5^2 + 1.5^2 + 2.5^2)

It's also a fact that for multiple independent samples from the same distribution, their variances add. So, if you roll N dice, you should get a new distribution with mean 3.5*N and variance 35*N/12.

So, if you generate a normal distribution with mean 3.5*N and variance 35*N/12, it will be a pretty good fit, assuming you're rolling a decent number of dice.

I'm wondering if there's a way to skip the number crunching and go straight to the curve itself. Is there a way to do this with, for example, 80d6 (range: 80-480)? Can the distribution be projected without doing 6^80 calculations?

Yes. The probability function of the sum of independent variables is the convolution of the probability function of each variable.

The convolution in this case is just a special summation. (More generally, the convolution is an integral.) Let p and q be two discrete probability functions. Convolution is conventionally indicated by an asterisk.

(p * q)[i] = sum_{j=1}^(n_p) p[j] q[i - j + 1]

where i ranges from 1 to (n_p + n_q - 1) with n_p is the number of elements of p and n_q the number of elements of q. If (i - j + 1) is less than 1 or greater than n_q, then let q[i - j + 1] be zero (so those terms just disappear from the summation).

In the case at hand, you have p = q = [1/6, 1/6, 1/6, 1/6, 1/6, 1/6], n_p = n_q = 6. The distribution of the sum of 3 rolls is (p * p * p). The distribution of the sum of 80 rolls is (p * p * p * ... (76 more p's) ... * p).

I don't know PHP so I wrote a little program in Maxima.

discrete_conv (p, q) := makelist (discrete_conv1 (p, q, i), i, 1, length (p) + length (q) - 1);
discrete_conv1 (p, q, i) := sum (p [j] * foo (q, i - j + 1), j, 1, length (p));
foo (a, i) := if 1 <= i and i <= length (a) then a [i] else 0;
r : [1/6, 1/6, 1/6, 1/6, 1/6, 1/6];
discrete_conv (r, discrete_conv (r, r));
 => [1/216,1/72,1/36,5/108,5/72,7/72,25/216,1/8,1/8,25/216,7/72,
     5/72,5/108,1/36,1/72,1/216]

If you continue repeating discrete_conv, you should find the numbers become more and more like a normal distribution. This is an illustration of the central limit theorem.

It's entirely possible that I've made some mistake with indexing so you'll want to check that. Hope this sheds some light on the problem.