如何在满足每个结果的不同概率的标准之前找到 n 次平均试验次数？
How to find n average number of trials before criteria are met with differing probabilities per outcome?
我花了几天时间试图解决这个问题并查找教程，但到目前为止我发现的一切似乎都接近我需要的，但没有给出我需要的结果。
我有一个可以生成单个字母 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
我需要找到要累积的这种字母组合必须发生的平均字母选择次数（即 rolls/trials/attempts）。个别结果的字母数超过要求的字母数，但只有当每个字母至少被选择了最少次数时才算成功。
我看过很多关于多项式概率分布和类似内容的教程，但我还没有找到任何内容来解释如何为这种情况找到平均试验次数。请清楚地解释答案，因为我不是统计专家。
嗯，最小投掷次数是 10。平均数将是无限总和
A=10•P（10 日完成）+11•P（11 日完成）+12•P（12 日完成）+ ...
对于 P（在 10 内完成）我们可以使用多项式
P(10)=Pm(1,3,0,1,2,3|概率)，其中概率=[.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,3,0,1,2,4|概率)
对于 P(12)，您必须再分配 2 次掷球。请注意，有些组合是不可能获得的，例如 Pm(2,3,0,2,2,3|probs)，因为您必须提前停止
以此类推
除了 Severin 的 答案在我看来在逻辑上看起来不错，但评估起来可能代价高昂（即阶乘的无穷和）。
让我提供一些直觉，应该给出一个很好的近似值。
一次考虑每个类别。参考这个math stackexchange question/ answer。对于每个类别 (i)，您将获得 k 成功次数的预期投掷次数可以计算为 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 次 让你的场景发挥作用。
您的过程可以描述为具有有限数量状态和吸收状态的马尔可夫链。
达到吸收状态之前的步数称为击球时间。预期命中时间可以很容易地从马尔可夫链的转移矩阵中计算出来。
枚举所有可能的状态（a、b、c、d、e、f）。仅考虑有限数量的状态，因为“b >= 3”实际上与“b = 3”等同。状态总数为(1+1)*(3+1)*(0+1)*(2+1)*(3+1) = 192.
确保在您的枚举中，起始状态 (0, 0, 0, 0, 0, 0) 首先出现，索引为 0，吸收状态 (1, 3, 0, 1 , 2, 3) 最后。
构建转换矩阵P。它是一个方阵，每个状态一行一行。矩阵中的条目 P[i, j] 给出了掷骰子时从状态 i 进入状态 j 的概率。每行最多应有 6 non-zero 个条目。
例如，如果i是状态(1, 0, 0, 1, 2, 2)的索引，j是状态(1, 1, 0, 1, 2, 2)的索引，那么P[i, j] = 滚动面 B 的概率 = 0.05。又如：如果i是状态(1,3,0,0,0,0)的索引，那么P[i,i]=滚动A、B或C的概率=0.25+0.05+0.2=0.5.
调用Q去掉P的最后一行最后一列得到的方阵。
调用I与Q相同维度的单位矩阵。
计算矩阵M = (I - Q)^-1，其中^-1是矩阵求逆。
在矩阵 M 中，条目 M[i, j] 是从状态开始时在吸收状态之前达到状态 j 的预期次数i.
由于我们的实验从状态 0 开始，我们对矩阵 M 的第 0 行特别感兴趣。
矩阵M第0行之和为吸收状态前预计达到的状态总数。这正是我们寻求的答案：达到吸收状态的步数。
要了解其工作原理，您应该阅读有关马尔可夫链的课程！也许这个：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.
下面是 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。
- 然后我们计算矩阵
Q 和 M 正如我在上面指出的那样。
- 我们return第一行所有单元格的总和
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)，以便您可以单独查看变量。
一个Monte-Carlo模拟：
- 实际掷骰子直到我们达到要求；
- 数一数我们做了多少卷；
- 重复实验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
我花了几天时间试图解决这个问题并查找教程，但到目前为止我发现的一切似乎都接近我需要的，但没有给出我需要的结果。
我有一个可以生成单个字母 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
我需要找到要累积的这种字母组合必须发生的平均字母选择次数（即 rolls/trials/attempts）。个别结果的字母数超过要求的字母数，但只有当每个字母至少被选择了最少次数时才算成功。
我看过很多关于多项式概率分布和类似内容的教程，但我还没有找到任何内容来解释如何为这种情况找到平均试验次数。请清楚地解释答案，因为我不是统计专家。
嗯，最小投掷次数是 10。平均数将是无限总和
A=10•P（10 日完成）+11•P（11 日完成）+12•P（12 日完成）+ ...
对于 P（在 10 内完成）我们可以使用多项式
P(10)=Pm(1,3,0,1,2,3|概率)，其中概率=[.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,3,0,1,2,4|概率)
对于 P(12)，您必须再分配 2 次掷球。请注意，有些组合是不可能获得的，例如 Pm(2,3,0,2,2,3|probs)，因为您必须提前停止
以此类推
除了 Severin 的 答案在我看来在逻辑上看起来不错，但评估起来可能代价高昂（即阶乘的无穷和）。 让我提供一些直觉，应该给出一个很好的近似值。
一次考虑每个类别。参考这个math stackexchange question/ answer。对于每个类别 (i)，您将获得 k 成功次数的预期投掷次数可以计算为 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 次 让你的场景发挥作用。
您的过程可以描述为具有有限数量状态和吸收状态的马尔可夫链。
达到吸收状态之前的步数称为击球时间。预期命中时间可以很容易地从马尔可夫链的转移矩阵中计算出来。
枚举所有可能的状态（a、b、c、d、e、f）。仅考虑有限数量的状态，因为“b >= 3”实际上与“b = 3”等同。状态总数为
(1+1)*(3+1)*(0+1)*(2+1)*(3+1) = 192.确保在您的枚举中，起始状态 (0, 0, 0, 0, 0, 0) 首先出现，索引为 0，吸收状态 (1, 3, 0, 1 , 2, 3) 最后。
构建转换矩阵
P。它是一个方阵，每个状态一行一行。矩阵中的条目P[i, j]给出了掷骰子时从状态i进入状态j的概率。每行最多应有 6 non-zero 个条目。例如，如果
i是状态(1, 0, 0, 1, 2, 2)的索引，j是状态(1, 1, 0, 1, 2, 2)的索引，那么P[i, j]= 滚动面 B 的概率 = 0.05。又如：如果i是状态(1,3,0,0,0,0)的索引，那么P[i,i]=滚动A、B或C的概率=0.25+0.05+0.2=0.5.调用
Q去掉P的最后一行最后一列得到的方阵。调用
I与Q相同维度的单位矩阵。计算矩阵
M = (I - Q)^-1，其中^-1是矩阵求逆。在矩阵
M中，条目M[i, j]是从状态开始时在吸收状态之前达到状态j的预期次数i.由于我们的实验从状态 0 开始，我们对矩阵 M 的第 0 行特别感兴趣。
矩阵M第0行之和为吸收状态前预计达到的状态总数。这正是我们寻求的答案：达到吸收状态的步数。
要了解其工作原理，您应该阅读有关马尔可夫链的课程！也许这个：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.
下面是 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。 - 然后我们计算矩阵
Q和M正如我在上面指出的那样。 - 我们return第一行所有单元格的总和
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)，以便您可以单独查看变量。
一个Monte-Carlo模拟：
- 实际掷骰子直到我们达到要求；
- 数一数我们做了多少卷；
- 重复实验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