2
$\begingroup$

I was given this problem by a colleague, and I cannot figure it out.

Given nine cups and an unlimited supply of six-sided dice rolled one at a time, what is the probability that you will, at some point, obtain a value of exactly seven in each of the nine cups, if the cup into which any given dice lands is completely random?

$\endgroup$

closed as off-topic by boboquack, Wen1now, JMP, user58, Glorfindel Oct 23 '17 at 11:28

This question appears to be off-topic. The users who voted to close gave this specific reason:

  • "This question is off-topic as it appears to be a mathematics problem, as opposed to a mathematical puzzle. For more info, see "Are math-textbook-style problems on topic?" on meta." – boboquack, JMP, user58, Glorfindel
If this question can be reworded to fit the rules in the help center, please edit the question.

  • $\begingroup$ Welcome to Puzzling! Why don't you take the tour and earn your first badge? Just a heads-up, if this puzzle is purely mathematical, it may be closed, unless there is a more lateral solution. $\endgroup$ – boboquack Oct 21 '17 at 21:35
  • $\begingroup$ This would make a nice Euler Project question. Maybe is already is. If it is, I'll remove the detailed answer, and just show the idea. Let me know if anyone knows. $\endgroup$ – Dr Xorile Oct 22 '17 at 21:20
  • 1
    $\begingroup$ I'm voting to close this question as off-topic because it appears to be an exhausting bash with no nice solving method $\endgroup$ – Wen1now Oct 23 '17 at 1:01
  • $\begingroup$ Can this be transferred to the math se? It's an interesting probability problem, and I'd be interested to see what they do? $\endgroup$ – Dr Xorile Oct 29 '17 at 18:50
  • $\begingroup$ Math.SE would want you to provide context and a little more clarity to what the question is looking for. If you add that here I'll migrate it. (@DrXorile) $\endgroup$ – Rubio Nov 1 '17 at 22:36
5
$\begingroup$

So the probability is vanishingly small. It is (I believe) $1.274\times10^{-10}$. Or, to be exact:

1580941466784207610825588390565400459910497074701370825267293858713643202540644069402625/  
12408574807343622185064460000854792067100180333248751160777882097868222836057400196426369022820352

I used recursive programming with python to get this. I'll throw the code below for people to audit. But remember, to understand recursive programming, you must first understand recursive programming.

The idea of the code is to create a state of the total throw in the 9 cups. I keep the state sorted to reduce the number of possible states to a manageable number (11,440).

I then have a function (prob) that figures out the probability of eventually getting (7,7,7,7,7,7,7,7,7) from a given starting point. It works in terms of itself (which is the recursive part).

So, for a simplified example: \begin{eqnarray} prob((1,0,0)) &=& \frac{1}{18}(prob((2,0,0))+prob((3,0,0))+\cdots+prob((7,0,0)))\\ &+& \frac{2}{18}(prob((1,1,0))+prob((2,1,0))+\cdots+prob((6,1,0)) \end{eqnarray} where the denominator of the fractions is the number of states (3 in the simple example, 9 in the actual example) multiplied by the number of throws (6); and the numberators is the number of cups effected.

By using the fractions module we keep everything exact. And the memoize code just means that when I call prob more than once, it saves the answer from the previous occasion.

Code is below:

from collections import defaultdict
import fractions

class memoize:
  def __init__(self, function):
    self.function = function
    self.memoized = {}

  def __call__(self, *args):
    try:
      return self.memoized[args]
    except KeyError:
      self.memoized[args] = self.function(*args)
      return self.memoized[args]

def add_dice(d,tot,add):
    """d is a dictionary, tot is the total in the cup(s)
    d[tot] is the number of cups with tot in them, and add
    is the new throw added to one of them.
    Returns state= (n1,n2,...,n9). n1>=n2>=...>=n9"""
    d[tot]+=-1
    d[tot+add]+=1
    s = []
    for k in d:
        s+=[k]*d[k]
    s.sort(reverse=True)
    return tuple(s)

@memoize
def prob(state):
    """state = (n1,n2,...,n9). n1>=n2>=...>=n9
    is the score in the cups.
    Returns p(success)"""
    if state[0]==7 and state[-1]==7:
        return fractions.Fraction(1,1)
    #Convert to dictionary 
    a = defaultdict(int)
    for s in state:
        a[s]+=1
    probability=fractions.Fraction(0)
    for k in a:
        for d in range(1,min([7,8-k])):
            suc = prob(add_dice(a.copy(),k,d))
            probability += suc*fractions.Fraction(a[k],len(state)*6)
    return probability

print(prob((0,0,0,0,0,0,0,0,0)))

It is easy enough to confirm this in simpler instances. For example: prob((0,)) gives you 70993/279936, which can be confirmed by multiplying out a Markov chain matrix, or just looking at the tree.

Here's the other ways of working out the 1 cup scenario, which is already a challenging puzzle. In this scenario, we are assuming just 1 cup and asking how often it ends up on 7.

We can, in this scenario, consider 7 throws. There are $6^7=279936$ possible throws. We then have to look at how many of them add to 7 when we are allowed to throw away as many as we need to. So (6,1,x,x,...) and (1,6,x,x,...) would work, and would each have $6^5=7776$ items. You would need to continue in this way. There might be a clever way of counting them all, but I can't see it right now. The generalization to 9 cups would be horrible.

You might wonder why we can't just look at $7\times9=63$ throws in the 9 cups case (I say “just” but obviously this is already a lot). The problem is that you wouldn't necessarily be limited to just 7 throws in a single cup, so the number of possibilities grows dramatically from there. Hence the huge denominator above.

The alternative is to look at a Markov chain. One might be able to generalize a Markov chain more easily, but it would involve a matrix that was $11,441\times11441$ being multiplied to itself a number of times. The matrix is relatively sparse, but this is probably still a bit of work.

Anyway, it's a nice example in the one cup case, and it confirms that the above code works at least in this simple(r) case.

So to do a Markov Chain, we need to envision 9 states:

  1. Score of 0
  2. Score of 1
  3. Score of 2
  4. Score of 3
  5. Score of 4
  6. Score of 5
  7. Score of 6
  8. Score of 7
  9. Bust.

We then need a matrix to hold the probabilities of transitioning from one state to the next with an additional throw:

\begin{equation}M=\begin{bmatrix} 0 & \frac{1}{6} & \frac{1}{6} & \frac{1}{6} & \frac{1}{6} & \frac{1}{6} & \frac{1}{6} & 0 & 0 \\ 0 & 0 & \frac{1}{6} & \frac{1}{6} & \frac{1}{6} & \frac{1}{6} & \frac{1}{6} & \frac{1}{6} & 0 \\ 0 & 0 & 0 & \frac{1}{6} & \frac{1}{6} & \frac{1}{6} & \frac{1}{6} & \frac{1}{6} & \frac{1}{6} \\ 0 & 0 & 0 & 0 & \frac{1}{6} & \frac{1}{6} & \frac{1}{6} & \frac{1}{6} & \frac{2}{6} \\ 0 & 0 & 0 & 0 & 0 & \frac{1}{6} & \frac{1}{6} & \frac{1}{6} & \frac{3}{6} \\ 0 & 0 & 0 & 0 & 0 & 0 & \frac{1}{6} & \frac{1}{6} & \frac{4}{6} \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & \frac{1}{6} & \frac{5}{6} \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 \\ \end{bmatrix}\end{equation}

Where the value in row $i$ and column $j$ shows the probability of moving from the $i$th state to the $j$th state. Note that if we get to a score of 7 or if we get to bust we stay there.

Our starting position is 0, which we represent as \begin{equation}S=\begin{bmatrix} 1 \\ 0 \\ \vdots \\ 0 \end{bmatrix}\end{equation}

We know that we will be done within 7 throws, one way or another, so we need to work out $S.M^7$. Actually $S.M^8$ will also be fine, and is easier computationally because we can just keep squaring our result. But computers can handle these little calculations easily.

\begin{equation} S.M^n=\begin{bmatrix}0 \\ 0 \\ 0 \\ 0 \\ 0 \\ 0 \\ 0 \\ \frac{70993}{279936} \\ \frac{208943}{279936} \end{bmatrix} \end{equation} for any $n\geq7$. So the probability of ending on 7 with one cup is $\frac{70993}{279936}$ in agreement with our running prob((0,)). The code to do this is as follows:

import numpy as np
import fractions

S = np.matrix([1,0,0,0,0,0,0,0,0])
f = fractions.Fraction(1,6)
M = np.matrix([[0,f,f,f,f,f,f,0,0],
          [0,0,f,f,f,f,f,f,0],
          [0,0,0,f,f,f,f,f,f],
          [0,0,0,0,f,f,f,f,2*f],
          [0,0,0,0,0,f,f,f,3*f],
          [0,0,0,0,0,0,f,f,4*f],
          [0,0,0,0,0,0,0,f,5*f],
          [0,0,0,0,0,0,0,1,0],
          [0,0,0,0,0,0,0,0,1]])

print(S*M**7)
$\endgroup$

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