Marbles arranged in a circle probability question

Question

You have a bag with an infinite amount of marbles, 36% of them are red and 64% of them are green.

You draw 46 marbles out of the bag and place them in a circle in random order.

What are the chances of having 3, 4, or 5 red marbles next to each other?

I tried solving this graphically, but I just can't work through all the possible combinations. I can do it if the circle is small, but I don't know how to turn it into an equation for larger numbers.

Answer 1

This can be modeled as a finite-state absorbing Markov chain. If we are looking for $k$ or more reds in a row, then the states are the $m=2^k$ sequences or red and green marbles of length $k$ and the absorbing state is the one with $k$ red marbles. For programming, it is more convenient to make the states the numbers from $0$ through $m-1$ represented as $k$-bit binary strings ($0$-padded on the left) so that the absorbing state is $m-1$.

If the probability of a $1$ is $p$ then the probability that the chain starts in state $S$ is $p^s(1-p)^{(k-s)}$ where $s$ is the number of $1$'s in the binary representation of $S$. This gives us the initial vector $V$. It's easy to get the transition matrix $P$. From a non-absorbing state $S$ the system can only go to two states, each starting with the last $k-1$ bits of $S$ and followed by a $0$ or $1$ with probability $1-p$ or $p$ respectively.

If there are $n$ marbles in the circle, we have to compute $$X=VP^{n-k}$$ since the initial state already has $k$ marbles. The final entry of $X$ gives the probability that we have encountered at least $k$ red marbles in a row.

So far, we have not accounted for the possibility that there are $k$ red marbles that "close the ring." That is, that the chain starts with $0<j<k$ red marbles, and ends with at least $k-j$ marbles. To compute this, we have to consider each individual state $S$ that starts with $1\leq j<k$ $1's$, compute the probability of ending in a non-absorbing state that ends in at least $k-j$ $1$'s, given that we start in state $S$, and multiply by the probability of starting in state $S$. The conditional probability is computed in a manner analogous the the calculation described above.

I wrote a python script that does these calculations. The output from the script is

3 (0.7745555576353794, 0.012089246740243075)
4 (0.3882662366474866, 0.017389132719839696)
5 (0.1547312088582764, 0.011606489722908718)

The first number in each line is $k$. The second number is the computed probability of ending in the absorbing state, and the third number is the additional probability from closing the ring. I broke it out this way to try to judge whether the additional computation to compute the probability of closing the ring is worthwhile. Since the number of cases we must consider increases exponentially with $k$, this is worth thinking about. Looking at the $k=5$ case, it appears that the answer is "yes", and leak should probably be changed to return the sum of the two probabilities.

Here's my script:

import numpy as np
'''
n marbles are arranged in a ring
Each marble is red (1) with probability p and green (0) with
probability 1-p
What is the probability that at least k consecutive marbles are red

We view this as an absorbing Markov chain with 2^k states, consisting
of the possible bit strings of length k.  The state 11...11 is absorbing.
Once we have the transition matrix P, and initial vector X, we can essentially
solve the problem by computing P^nX and checking the probability of
ending in the absorbing state.

Then we have to add in the probability of "closing the ring", that is 
the probability that the chain starts in a state with 0<j<k 1's and
ends in a non-absorbing state that ends in at least k-j 1's.  
'''
def prob(state, p, k):
    '''Probability of starting in state'''
    fmt = '0%db'%k
    s = format(state, fmt).count('1')
    return p**s *(1-p)**(k-s)

def initial(p,k):
    m = 2**k
    a = np.zeros(m)
    for n in range(m):
        a[n] = prob(n, p, k)
    return a

def leak(n, p, k):
    m= 2**k
    P = np.zeros((m,m))
    V = initial(p, k) 
    for i in range(m-1):
        j = (2*i) % m
        P[i, j] = 1-p
        P[i, j+1] = p
    # all red is absorbing
    P[m-1,m-1] = 1
    N = np.linalg.matrix_power(P,n-k)
    V = V@N
    # Now compute probabilities of k reds to close the ring
    # State must start with red, but not be absorbing
    base = V[-1]
    fmt = '0%db'%k
    addl = 0
    for s in range(m//2, m-1):
        pr = prob(s, p, k)
        V = np.zeros(m) 
        V[s] = 1
        V = V@N
        f = format(s,fmt)
        f0 = f.index('0')  # s starts with f0 1's
        need = k-f0       # need this many 1's at the end
        md = 2**need
        addl += pr*sum(V[md-1:-1:md])
    return base, addl

for k in range(3,6):
    print(k, leak(46, .36, k))