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

Question

在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)

但我想问一下您最新的情况，是否有办法适用于这三个离散分布，并选择最适合离散数据集的拟合方式？

- Salih

---

你的意思是说没有支持？那 https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.curve_fit.html 呢？ - mkrieger1

我知道"no fit method"。我想学习是否有任何方法可以拟合这些离散分布并获取其参数... @mkrieger1 - Salih

FYI：有关泊松分布，请参见https://dev59.com/i5ffa4cB1Zd3GeqP84Ed#37500643。 - Warren Weckesser

谢谢，但我正在尝试为所有离散分布寻找合适的拟合，并选择最佳拟合方法 @Weckesser。 - Salih

没有通用的方法来拟合任意离散分布，因为存在无限多个这样的分布，其参数可能是无限的。不过，有一些方法可以拟合特定的分布，例如矩法。如果你只需要这三个分布，我可以展示如何使用它们。 - Marat

Answer 1

您可以使用矩法来拟合任何特定的分布。

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

因此，在所有这些情况下，我们只需要两个矩。让我们来获取它们：

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()和std()没有应用贝塞尔矫正，而pandas的应用了。如果你有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

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

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

更新结束

更新2：

我建议使用@thilak的代码来计算对数似然，它可以避免在大样本中丢失精度，这尤其重要。

更新2结束

计算似然：

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

对于泊松分布也是一样，只有一个参数：

# 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])

如果您有任何问题，请告知我

- Marat

---

非常感谢。我还有一个问题，如果您能回答就感激不尽了。我知道我的数据集是离散的，但假设我想看看它是否符合正态分布。这可能吗？是否有像矩估计法这样的方法可以做到这一点？ - Salih

@Salih 是的，它适用于像高斯分布这样的连续分布。然而，在某些情况下，估计所有参数是困难甚至不可能的。例如，对于两个二项式的混合物，您将需要三个参数和三个矩; 这已经很难解决了。一旦您将更多组分添加到混合物中，情况就会变得更糟。 - Marat

@Maral，如果我想为两个二项混合物做这件事，我该怎么做？你能给我指条路吗？ - Salih

@Salih 抱歉，我不会这样做。它已经偏离了原始问题的范围，所以最好为此发布一个单独的问题。另外，实际上有五个参数，而不是三个，正如我所提到的，在封闭形式中解决它变得非常不愉快。在实践中，混合模型通常使用EM算法和高斯函数进行估计。 - Marat

@Maral，你确定负二项式的代码是正确的吗？当我生成随机数据时，它并没有给出该分布的正确结果。 - Salih

Answer 2

Marat的回答很好。

除了Marat的帖子之外，我肯定会建议将概率质量函数取对数。关于为什么优先考虑对数似然而不是似然的一些信息-https://math.stackexchange.com/questions/892832/why-we-consider-log-likelihood-instead-of-likelihood-in-gaussian-distribution

我会重写负二项式的代码为-

log_likelihoods = {}
log_likelihoods['nbinom'] = x.map(lambda val: nbinom.logpmf(val, r, p)).sum()

请注意，我使用了以下内容进行翻译 -

• 使用logpmf代替pmf
• 使用sum代替product

为了找到最佳分布 -

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

- thilak g