我的问题是关于实现概率函数时遇到的数值问题,而不是与其背后的概率/数学有关。我也知道下面的代码可能没有很好地优化(例如,如果在“comb”中使用“exact=False”,则可以对第一个函数进行矢量化)。所以我愿意接受优化建议,但这并不是我目前的主要关注点。
我试图在Python 3.6.5中验证给出的公式here,用于“在选择n次时从[0,k)获得m个唯一值的概率”。
为此,我使用
对于较小的k和n值,模拟结果与公式非常吻合,因此我相当满意它基本上是正确的。但是,当k和n稍微大一些时,我从公式中获得负值。我怀疑这是因为它包含了微小分数和非常大的阶乘的乘积,因此在某些阶段可能会失去精度。
为了尝试解决这个问题,我实现了相同的公式,但是在可能的情况下使用对数,最后再进行指数运算。令人烦恼的是,它并没有真正帮助,可以在下面给出的代码输出中看到。
因此,我的问题是,是否有人能提出建议,以便我可以继续实现更大的n和k值的这个公式?我是否正确地认为这是由于大数和小数的乘积引入的数字异常?
我的代码:
输出的图像:
感谢任何想法/建议/建议。
我试图在Python 3.6.5中验证给出的公式here,用于“在选择n次时从[0,k)获得m个唯一值的概率”。
为此,我使用
numpy.ramdom.choice(k, n, replace=True)
在Python 3.6.5中获取一个多重集合,然后计算多重集合中唯一值的数量,并保存该数字。然后重复。对于较小的k和n值,模拟结果与公式非常吻合,因此我相当满意它基本上是正确的。但是,当k和n稍微大一些时,我从公式中获得负值。我怀疑这是因为它包含了微小分数和非常大的阶乘的乘积,因此在某些阶段可能会失去精度。
为了尝试解决这个问题,我实现了相同的公式,但是在可能的情况下使用对数,最后再进行指数运算。令人烦恼的是,它并没有真正帮助,可以在下面给出的代码输出中看到。
因此,我的问题是,是否有人能提出建议,以便我可以继续实现更大的n和k值的这个公式?我是否正确地认为这是由于大数和小数的乘积引入的数字异常?
我的代码:
import numpy as np
import numpy.random as npr
from scipy.special import comb, gammaln
import matplotlib.pyplot as plt
def p_unique_birthdays(m, k, n):
"""PMF for obtaining m unique elements when selecting from [0,k) n times.
I wanted to use exact=True to see if that helped, hence why this is not
vectorised.
"""
total = 0
for i in range(m):
total += (-1)**i * comb(m, i, exact=True) * ((m-i)/k)**n
return comb(k, m, exact=True) * total
def p_unique_birthdays_logs(m, k, n):
"""PMF for obtaining m unique elements when selecting from [0,k) n times.
I use logs to try and deal with some of the numerical craziness that seems
to arise.
"""
total = 0
for i in range(m):
log_mCi = gammaln(m+1) - gammaln(i+1) - gammaln(m-i+1)
log_exp_bit = n * (np.log(m-i) - np.log(k))
total += (-1)**i * np.exp(log_mCi + log_exp_bit)
return comb(k, m, exact=True) * total
def do_stuff(k, n, pmf):
n_samples = 50000
p_ms = np.zeros(n)
for i in range(n):
temp_p = pmf(i+1, k, n)
p_ms[i] = temp_p
print("Sum of probabilities:", p_ms.sum())
samples = np.zeros(n_samples)
for i in range(n_samples):
samples[i] = np.unique(npr.choice(k, n, replace=True)).size
# So that the histogram is centered on the correct integers.
d = np.diff(np.unique(samples)).min()
left_of_first_bin = samples.min() - float(d)/2
right_of_last_bin = samples.max() + float(d)/2
fig = plt.figure(figsize=(8,5))
ax = fig.add_subplot(111)
ax.grid()
ax.bar(range(1, n+1), p_ms, color="C0",
label=labels[j])
ax.hist(samples, np.arange(left_of_first_bin, right_of_last_bin + d, d),
alpha=0.5, color="C1", density=True, label="Samples")
ax.legend()
ax.set_xlabel("Unique birthdays")
ax.set_ylabel("Normalised frequency")
ax.set_title(f"k = {k}, n = {n}")
#fig.savefig(f"k{k}_n{n}_{labels[j]}.png")
plt.show()
random_seed = 1234
npr.seed(random_seed)
labels = ["PMF", "PMF (logs)"]
pmfs = [p_unique_birthdays, p_unique_birthdays_logs]
for j in range(2):
for k, n in [(30, 20), (60, 40)]:
do_stuff(k, n, pmfs[j])
输出的图像:
![output](https://istack.dev59.com/FzGqc.webp)
![output](https://istack.dev59.com/4Gjro.webp)
![output](https://istack.dev59.com/WQE5K.webp)
![output](https://istack.dev59.com/8SI2O.webp)
total += (-1)**i * comb(m, i, exact=True) / k**n * (m-i)**n
的效果更好。不幸的是,无论哪个版本都无法处理稍大的 n 和 k 值,所以我猜这个数值计算确实很难克服! - combinatoricky