我花了几天时间试图弄清楚这一点，并查找教程，但到目前为止我发现的一切似乎都接近我需要的，但没有给出我需要的结果。
我有一个能产生一个字母a-F的装置。为了简单起见，你可以把它想象成一个有字母的骰子。每次使用时，它总是会产生一个且只有一个字母。然而，它有一个主要区别：每个字母被选中的已知概率不同：
A: 25%
B: 5%
C: 20%
D: 15%
E: 20%
F: 15%
这些概率在所有尝试中都保持不变。
此外，我有一个特定的组合，我必须积累，我是"；"成功"：
As needed: 1
Bs needed: 3
Cs needed: 0
Ds needed: 1
Es needed: 2
Fs needed: 3
我需要找到要累积的字母组合必须发生的字母拾取(即滚动/试用/尝试)的平均数量任何一个单独的结果都可以有超过所需数量的字母，但只有当每个字母至少被选择了最少的次数时，成功才算在内。
我看了很多关于多项式概率分布和类似内容的教程，但我没有找到任何解释如何在这样的场景中找到平均试验次数的东西。请把答案解释清楚，因为我不是统计学专家。
除了Severin的答案在逻辑上看起来不错，但评估(即阶乘的无限和)可能会很昂贵。让我提供一些直觉，应该会给出一个很好的近似值。
一次考虑每个类别。请参阅此数学堆叠交换问题/答案。预期投掷次数，其中您将获得k每个类别的成功次数(i)可以计算为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
如果你认为获得3个B是你的瓶颈，那么你的场景平均会有60次投掷。
好吧，最小投掷次数是10次。平均值将是的无穷和
A=10•p(10年完成)+11•p(11年完成)+12•p(12年完成)+。。。
对于p(在10中完成)，我们可以使用多项式
p(10)=Pm(1,3,0,1,2,3| probs)，其中probs=[0.25，.05，.20，.15，.20，.15]
对于p(11)，你还有一次投掷，你可以像这个一样分发
p(11)=Pm(2,3,0,1,2,3|问题)+Pm(1,4,0,1,2.3|问题)+Pm(1,3,0,2,2,3 |问题)+Pm(1,3,0,1,3,3|问题)+Pm(1,30,0,1,2,4|问题)
对于p(12)，你必须再投2次球。注意，有些投掷组合是不可能得到的，比如Pm(2,3,0,2,3|probs)，因为你必须提前停止
等等
您的过程可以被描述为具有有限状态和吸收状态的马尔可夫链。
达到吸收状态之前的步骤数称为撞击时间。根据马尔可夫链的转移矩阵可以很容易地计算出预期的命中时间。
-
枚举所有可能的状态(a、b、c、d、e、f)。只考虑有限数量的状态，因为"；b>=3〃；实际上与"；b＝3"；，状态的总数是CCD_ 1。
-
请确保在枚举中，起始状态(0，0，0、0、0，0)排在第一位，索引为0，吸收状态(1，3，0，1，2，3)排在最后。
-
建立转换矩阵
P。它是一个正方形矩阵，每个状态有一行和一列。矩阵中的条目P[i, j]给出了当轧制模具时从状态i进入状态j的概率。每行最多应有6个非零条目。 -
例如，如果
i是状态(1, 0, 0, 1, 2, 2)的索引，而j是状态(1, 1, 0, 1, 2, 2)的索引，则(1+1)*(3+1)*(0+1)*(2+1)*(3+1) = 1920＝滚动面的概率B＝0.05。另一个例子：如果i是状态(1,3,0,0,0,0)的索引，则P[i,i]=滚动A、B或C的概率=0.25+0.05+0.2=0.5。 -
将去掉
P的最后一行和最后一列得到的平方矩阵称为Q。 -
将
I称为与Q具有相同维度的恒等矩阵。 -
计算矩阵
M = (I - Q)^-1，其中^-1是矩阵求逆。 -
在矩阵
M中，条目M[i, j]是从状态i开始时在吸收状态之前达到状态j的预期次数。 -
由于我们的实验从状态0开始，我们对矩阵M的第0行特别感兴趣。
-
矩阵M的第0行之和是在吸收状态之前达到的预期状态总数。这正是我们寻求的答案：达到吸收状态的步骤数。
要理解为什么这样做，你应该阅读一门关于马尔可夫链的课程！也许是这样的：詹姆斯·诺里斯关于马尔可夫链的课程笔记。关于"；击球次数"；(这是达到目标状态之前的步骤数的名称)是第1.3章。
下面是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
代码说明：
- 首先，我们使用笛卡尔乘积
itertools.product构建状态列表 - 对于给定的状态
i和模具面f，我们需要计算j=当加上f时从i达到的状态。我有两种计算方法，要么作为j = i + base[f] * (s[f] < targets[f])，要么作为j = lookup[s[:f] + (min(s[f]+1, targets[f]),) + s[f+1:]]。因为我有妄想症，所以我用两种方法计算，并检查两种方法是否给出了相同的结果。但你只需要一种方式。如果需要，可以删除j1 = ...到assert(j == j1)行 - 矩阵
P从零开始填充，我们用P[i, j] += p每行最多填充六个单元格，其中p是滚动面f的概率 - 然后我们计算矩阵CCD_ 36和CCD_
- 我们返回
M第一行上所有单元格的总和
为了帮助您更好地了解发生了什么，我鼓励您检查所有变量的值。例如，您可以将return M[0, :].sum()替换为return states, base, lookup, P, Q, I, M，然后在python交互式shell中编写states, base, lookup, P, Q, I, M = get_expected_n_trials(targets, dice_weights)，这样您就可以单独查看变量。
蒙特卡罗模拟：
- 实际滚动模具，直到我们达到要求
- 数数我们做了多少卷
- 重复实验1000次，得到经验平均值
在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