1
$\begingroup$

I want to make an application that gives you the probability of reaching a certain threshold with a number of rolls of special dice. Just to clarify, I am not asking for help on programming, I need a formula to program from.

The variables are

  • d: the number of dice, all with the sides [1, 2, 2, 3, 3, 4]

  • x: the number the roll needs to reach

I need to work out the probability of rolling d dice (of the special type noted above) and getting a total equal to or greater than x

$\endgroup$
  • 2
    $\begingroup$ I'd try simulation. Make a virtual die and roll $d$ of them $10000$ times. Then see how many times you get to that total.The alternative is just solving a generic equation which I'm sure you can do. $\endgroup$ – stuart stevenson Jan 1 '18 at 23:52
1
$\begingroup$

Let $p(t,d)$ be the probability that when you roll $d$ dice, the sum of their numbers is exactly $t.$

Then $p(t,d) = 0$ whenever $t < 0.$ Also $p(t,0) = 0$ when $t > 0,$ but $p(0,0) = 1.$

For $d \geq 0,$ $$p(t,d+1) = \frac16 p(t-1,d) + \frac13 p(t-2,d) + \frac13 p(t-3,d) + \frac16 p(t-4,d).$$

You can use that formula to find any value of $p(t,d)$ recursively. An efficient way to do this is to compute $p(t,1)$ for every relevant value of $t,$ then $p(t,2)$ for every relevant $t,$ and so forth until you have computed the necessary values of $p(t,d)$ in order to find the probability of rolling at least $x$; that probability is $$ \sum_{t\geq x} p(t,d).$$

| cite | improve this answer | |
$\endgroup$
  • $\begingroup$ Note that for $t \in \{d,d+1,\cdots, 4d\}$, $p(t,d)$ is the coefficient of $x^t$ in the polynomial $\phi^d(x)$ below. Exactly the same recursion appears, unsurprisingly, it is essentially a convolution. Do note that $p(t,d)$ is nonzero iff for $t \in \{d,d+1,\cdots, 4d\}$. $\endgroup$ – copper.hat Jan 3 '18 at 5:15
1
$\begingroup$

For small values of $d$ you could use the following generating function approach:

Let $\phi(x) = {x+ 2 x^2 + 2 x^3 + x^4 \over 6}$, then the probability that $d$ rolls gives exactly $m$ is the coefficient of $x^m$ in $\phi^d$, where $m=d,d+1,...,4d$.

To illustrate: If $d=7$ then $\phi^d(x) = \sum_{k=0}^{28} a_k x^k$ where $a={1 \over 6^7}(0 , 0 , 0 , 0 , 0 , 0 , 0 , 1 , 14 , 98 , 455 , 1568 , 4256 , 9429 , 17446 , 27370 , 36771 , 42560 , 42560 , 36771 , 27370 , 17446 , 9429 , 4256 , 1568 , 455 , 98 , 14 , 1)$.

For example, the probability of getting a sum of exactly $26$ is the coefficient of $x^{26}$ which is ${98 \over 6^7}$.

The probability of the sum being $\ge x$ will be the sum of the coefficients of $x^k$ where $k \ge x$.

There are various ways of computing the polynomial coefficients depending on what you are trying to achieve.

If you have access to Octave/Matlab, you could try the following:

# poly of phi*6, constant term last
f = [1 2 2 1 0];
# 7 rolls
d = 7;
# compute polynomial of (phi*6)^d
p = 1;
for i = 1:d
    p = conv(p, f);
endfor
# scale result so the coefficients are probabilities
# constant term last
p = p/6^d;

If not, you could using something like Neville's algorithm to compute the interpolating polynomial of the points $\phi^d(0),\cdots, \phi^d(4d)$.

Another way would be to compute the coefficients by multiplying the polynomials $\phi^2 = \phi \cdot \phi$, $\phi^{k+1} = \phi^k \cdot \phi$. This is essentially the convolution of coefficients of the polynomials being multiplied.

Aside: My original Octave/Matlab code looked something like

f = [1 2 2 1 0];
r = roots(f);
d = 7;
v = poly(repmat(r',[1, d]));
v = v/6^d;

however, for $d \in \{5,6,7\}$, the poly call resulted in some small imaginary noise. The conv call avoids this entirely.

| cite | improve this answer | |
$\endgroup$
  • $\begingroup$ How small? d will always be less than seven. $\endgroup$ – TheOnlyMrCat Jan 2 '18 at 8:46
  • $\begingroup$ The call to poly produces an error with both Octave and MATLAB because the matrix argument is not square. $\endgroup$ – Fabio Somenzi Jan 2 '18 at 21:54
  • $\begingroup$ @FabioSomenzi: My apologies, I will fix momentarily. $\endgroup$ – copper.hat Jan 2 '18 at 22:02
  • $\begingroup$ @FabioSomenzi: Thanks for catching that, I gave up on my previous approach as it failed for $d=7$. There were a number of errors in my previous code, I should know not to wing it at this stage! $\endgroup$ – copper.hat Jan 2 '18 at 23:05
  • $\begingroup$ Glad to be of some help! $\endgroup$ – Fabio Somenzi Jan 2 '18 at 23:13

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.