How to find n average number of trials before criteria are met with differing probabilities per outcome?

Question

I've spent a few days trying to figure this out and looking up tutorials, but everything I've found so far seems like it's close to what I need but don't give the results I need.

I have a device that produces a single letter, A-F. For simplicity's sake, you can think of it like a die with letters. It will always produce one and only one letter each time it is used. However it has one major difference: each letter can have a differing known probability of being picked:

A: 25%
B: 5%
C: 20%
D: 15%
E: 20%
F: 15%

These probabilities remain constant throughout all attempts.

Additionally, I have a specific combination I must accrue before I am "successful":

As needed: 1
Bs needed: 3
Cs needed: 0
Ds needed: 1
Es needed: 2
Fs needed: 3

I need find the average number of letter picks (i.e. rolls/trials/attempts) that have to happen for this combination of letters to be accrued. It's completely fine for any individual outcome to have more than the required number of letters, but success is only counted when each letter has been chosen at least its minimum amount of times.

I've looked at plenty of tutorials for multinomial probability distribution and similar things, but I haven't found anything that explains how to find average number of trials for a scenario like this. Please kindly explain answers clearly as I'm not a wiz with statistics.

Answer 1

In addition to Severin's answer that logically looks good to me but might be costly to evaluate (i.e. infinite sum of factorials). Let me provide some intuition that should give a good approximation.

Considering each category at a time. Refer this math stackexchange question/ answer. Expected number of tosses in which you would get the k number of successes for each category (i) can be calculated as k(i)/ P(i):

Given,
p(A): 25% ; Expected number of tosses to get 1 A = 1/ 0.25 = 4
p(B): 5% ; Expected number of tosses to get 3 B's = 3/ 0.05 = 60
p(C): 20% ; Expected number of tosses to get 0 C = 0/ 0.20 = 0
p(D): 15% ; Expected number of tosses to get 1 D = 1/ 0.15 = 6.67 ~ 7
p(E): 20% ; Expected number of tosses to get 2 E's = 2/ 0.20 = 10
p(F): 15% ; Expected number of tosses to get 3 F's = 3/ 0.15 = 20

you get an idea that getting 3 B's is your bottleneck, you can expect on average 60 tosses for your scenario to play out.

Answer 2

Well, minimum number of throws is 10. Average would be infinite sum

A=10•P(done in 10)+11•P(done in 11)+12•P(done in 12) + ...

For P(done in 10) we could use multinomial

P(10)=Pm(1,3,0,1,2,3|probs), where probs=[.25, .05, .20, .15, .20, .15]

For P(11) you have one more throw which you could distribute like this

P(11)=Pm(2,3,0,1,2,3|probs)+Pm(1,4,0,1,2,3|probs)+Pm(1,3,0,2,2,3|probs)+ Pm(1,3,0,1,3,3|probs)+Pm(1,3,0,1,2,4|probs)

For P(12) you have to distribute 2 more throws. Note, that there are combinations of throws which are impossible to get, like Pm(2,3,0,2,2,3|probs), because you have to stop earlier

And so on and so forth

Answer 3

Your process can be described as a Markov chain with a finite number of states, and an absorbing state.

The number of steps before reaching the absorbing state is called the hitting time. The expected hitting time can be calculated easily from the transition matrix of the Markov chain.

• Enumerate all possible states (a, b, c, d, e, f). Consider only a finite number of states, because "b >= 3" is effectively the same as "b = 3", etc. The total number of states is (1+1)*(3+1)*(0+1)*(2+1)*(3+1) = 192.

• Make sure that in your enumeration, starting state (0, 0, 0, 0, 0, 0) comes first, with index 0, and absorbing state (1, 3, 0, 1, 2, 3) comes last.

• Build the transition matrix P. It's a square matrix with one row and column per state. Entry P[i, j] in the matrix gives the probability of going from state i to state j when rolling a die. There should be at most 6 non-zero entries per row.

• For example, if i is the index of state (1, 0, 0, 1, 2, 2) and j the index of state (1, 1, 0, 1, 2, 2), then P[i, j] = probability of rolling face B = 0.05. Another example: if i is the index of state (1,3,0,0,0,0), then P[i,i] = probability of rolling A, B or C = 0.25+0.05+0.2 = 0.5.

• Call Q the square matrix obtained by removing the last row and last column of P.

• Call I the identity matrix of the same dimensions as Q.

• Compute matrix M = (I - Q)^-1, where ^-1 is matrix inversion.

• In matrix M, the entry M[i, j] is the expected number of times that state j will be reached before the absorbing state, when starting from state i.

• Since our experiment starts in state 0, we're particularly interested in row 0 of matrix M.

• The sum of row 0 of matrix M is the expected total number of states reached before the absorbing state. That is exactly the answer we seek: the number of steps to reach the absorbing state.

To understand why this works, you should read a course on Markov chains! Perhaps this one: James Norris' course notes on Markov chains. The chapter about "hitting times" (which is the name for the number of steps before reaching target state) is chapter 1.3.

Below, an implementation in python.

from itertools import product, accumulate
from operator import mul
from math import prod
import numpy as np

dice_weights = [0.25, 0.05, 0.2, 0.15, 0.2, 0.15]
targets = [1, 3, 0, 1, 2, 3]

def get_expected_n_trials(targets, dice_weights):
    states = list(product(*(range(n+1) for n in targets)))
    base = list(accumulate([n+1 for n in targets[:0:-1]], mul, initial=1))[::-1]
    lookup = dict(map(reversed, enumerate(states)))
    P = np.zeros((len(states), len(states)))
    for i, s in enumerate(states):
        a,b,c,d,e,f = s
        for f, p in enumerate(dice_weights):
            #j = index of state reached from state i when rolling face f
            j = i + base[f] * (s[f] < targets[f]) 
            j1 = lookup[s[:f] + (min(s[f]+1, targets[f]),) + s[f+1:]]
            if (j != j1):
                print(i, s, f, ' --> ' , j, j1)
            assert(j == j1)
            P[i,j] += p
    Q = P[:-1, :-1]
    I = np.identity(len(states)-1) 
    M = np.linalg.inv(I - Q)
    return M[0,:].sum()

print(get_expected_n_trials(targets, dice_weights))
# 61.28361802372382

Explanations of code:

• First we build the list of states using Cartesian product itertools.product
• For a given state i and die face f, we need to calculate j = state reached from i when adding f. I have two ways of calculating that, either as j = i + base[f] * (s[f] < targets[f]) or as j = lookup[s[:f] + (min(s[f]+1, targets[f]),) + s[f+1:]]. Because I'm paranoid, I calculated it both ways and checked that the two ways gave the same result. But you only need one way. You can remove lines j1 = ... to assert(j == j1) if you want.
• Matrix P begins filled with zeroes, and we fill up to six cells per row with P[i, j] += p where p is probability of rolling face f.
• Then we compute matrices Q and M as I indicated above.
• We return the sum of all the cells on the first row of M.

To help you better understand what is going on, I encourage you to examine the values of all variables. For instance you could replace return M[0, :].sum() with return states, base, lookup, P, Q, I, M and then write states, base, lookup, P, Q, I, M = get_expected_n_trials(targets, dice_weights) in the python interactive shell, so that you can look at the variables individually.

Answer 4

A Monte-Carlo simulation:

• Actually roll the die until we hit the requirements;
• Count how many rolls we did;
• Repeat experiment 1000 times to get the empirical average value.

Implementation in python:

from collections import Counter
from random import choices
from itertools import accumulate
from statistics import mean, stdev

dice_weights = [0.25, 0.05, 0.2, 0.15, 0.2, 0.15]
targets = [1, 3, 0, 1, 2, 3]

def avg_n_trials(targets, dice_weights, n_experiments=1000):
    dice_faces = range(len(targets))
    target_state = Counter(dict(enumerate(targets)))
    cum_weights = list(accumulate(dice_weights))
    results = []
    for _ in range(n_experiments):
        state = Counter()
        while not state >= target_state:
            f = choices(dice_faces, cum_weights=cum_weights)[0]
            state[f] += 1
        results.append(state.total()) # python<3.10: sum(state.values())
    m = mean(results)
    s = stdev(results, xbar=m)
    return m, s

m, s = avg_n_trials(targets, dice_weights, n_experiments=10000)
print(m)
# 61.4044