在Python中实现广义生日悖论

Question

我的问题是关于实现概率函数时遇到的数值问题，而不是与其背后的概率/数学有关。我也知道下面的代码可能没有很好地优化（例如，如果在“comb”中使用“exact=False”，则可以对第一个函数进行矢量化）。所以我愿意接受优化建议，但这并不是我目前的主要关注点。

我试图在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])

输出的图像： 感谢任何想法/建议/建议。

- combinatoricky

Answer 1

您可以使用内置的decimal模块来提高精度。

from decimal import *

getcontext().prec = 10000

def factorial(n):
    res = Decimal(1)
    for i in range(int(n)):
        res = res * Decimal(i + 1)
    return res

def binomial_coefficient(n, k):
    return factorial(n) / factorial(k) / factorial(n - k)

def p_unique_birthdays(m, k, n):
    m = Decimal(m)
    k = Decimal(k)
    n = Decimal(n)
    total = Decimal(0)
    for i in range(int(m) + 1):
        total += Decimal((-1) ** i) * binomial_coefficient(m, i) * binomial_coefficient(k, m) * ((m - i) / k) ** n
    return total

print(p_unique_birthdays(49, 365, 50))

上述代码输出0.11484925，与http://www.wolframalpha.com/input/?i=sum+combination(49,x)combination(365,49)++(((49-x)%2F365)%5E50)+*+(-1)%5Ex,+x%3D0+to+49相同。

- shaun shia

Answer 2

你是正确的，这是一些奇怪的数字原因。

将这行代码更改为： total += (-1)**i * comb(m, i, exact=True) * ((m-i)**n)/(k**n) 由于某种原因，如果强制使用不同的操作顺序，事情就会变得很好。

你可能需要花更多时间来找出如何修改你的“log'd”版本，但是考虑到上面的更改可以解决问题，你可能只想放弃“log'd”版本。

希望有所帮助！

- Felipe D.

---

谢谢您的建议！它很有效。实际上，在这种特殊情况下，似乎 total += (-1)**i * comb(m, i, exact=True) / k**n * (m-i)**n 的效果更好。不幸的是，无论哪个版本都无法处理稍大的 n 和 k 值，所以我猜这个数值计算确实很难克服！ - combinatoricky