离散数据拟合：负二项式、泊松分布、几何分布

在 scipy 中，不支持使用数据拟合离散分布。我知道有很多关于这个的话题。

例如，如果我有一个如下所示的数组：

x = [2,3,4,5,6,7,0,1,1,0,1,8,10,9,1,1,1,0,0]

我无法申请这个数组：

from scipy.stats import nbinom
param = nbinom.fit(x)

但我想问您最新的情况，有没有什么方法可以拟合这三个离散分布，然后选择最适合离散数据集的方法？

您可以使用矩量法以适应任何特定的分布。

基本思想：得到经验一阶矩、二阶矩等，然后从这些矩导出分布参数。

因此，在所有这些情况下，我们只需要两个时刻。让我们得到它们：

import pandas as pd
# for other distributions, you'll need to implement PMF
from scipy.stats import nbinom, poisson, geom

x = pd.Series(x)
mean = x.mean()
var = x.var()
likelihoods = {}  # we'll use it later

注意：我使用 pandas 而不是 numpy。那是因为 numpy 的var() and std()不申请贝塞尔修正，而熊猫则这样做。如果您有 100 多个样本，应该不会有太大差异，但对于较小的样本，这可能很重要。

现在，让我们获取这些分布的参数。负二项式有两个参数：p，r。让我们估计它们并计算数据集的可能性：

# From the wikipedia page, we have:
# mean = pr / (1-p)
# var = pr / (1-p)**2
# without wiki, you could use MGF to get moments; too long to explain here
# Solving for p and r, we get:

p = 1 - mean / var  # TODO: check for zero variance and limit p by [0, 1]
r = (1-p) * mean / p

UPD:维基百科和 scipy 使用不同的 p 定义，一种将其视为成功的概率，另一种将其视为失败的概率。因此，为了与 scipy 概念保持一致，请使用：

p = mean / var
r = p * mean / (1-p)

UPD 结束

UPD2:

我建议改用@thilak 的代码日志可能性。它可以避免精度损失，这对于大样本尤其重要。

UPD2 结束

计算可能性：

likelihoods['nbinom'] = x.map(lambda val: nbinom.pmf(val, r, p)).prod()

同样适用于Poisson，只有一个参数：

# from Wikipedia,
# mean = variance = lambda. Nothing to solve here
lambda_ = mean
likelihoods['poisson'] = x.map(lambda val: poisson.pmf(val, lambda_)).prod()

同样适用于几何分布:

# mean = 1 / p  # this form fits the scipy definition
p = 1 / mean

likelihoods['geometric'] = x.map(lambda val: geom.pmf(val, p)).prod()

最后，让我们获得最佳拟合：

best_fit = max(likelihoods, key=lambda x: likelihoods[x])
print("Best fit:", best_fit)
print("Likelihood:", likelihoods[best_fit])

如果您有任何疑问，请告诉我