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

Posted 2023-03-12

【中文标题】拟合离散数据：负二项式、泊松、几何分布

【英文标题】：Fitting For Discrete Data: Negative Binomial, Poisson, Geometric Distribution

【发布时间】：2020-04-06 01:44:17

【问题描述】：

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

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

【问题讨论】：

什么意思，没有支持？ docs.scipy.org/doc/scipy/reference/generated/…呢？

我编辑了我的问题@mkrieger1

x 应该是什么意思？你期望nbinom.fit(x) 做什么？ scipy.stats.nbinom 没有 fit 方法。

我知道“不合适的方法”。我想知道有什么方法可以拟合这些离散分布并获取其参数...@mkrieger1

没有适合任意离散分布的通用方法，因为它们的数量是无限的，并且可能具有无限的参数。但是，有一些方法可以适应特定的分布，例如矩量法。如果你只需要这三个我可以展示如何使用它

【参考方案1】：

您可以使用Method of Moments 来适应任何特定的分布。

基本思想：获取经验的第一、第二等矩，然后从这些矩中推导出分布参数。

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

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() 不适用 Bessel's correction，而 pandas 则适用。如果您有 100 多个样本，则应该不会有太大差异，但对于较小的样本，这可能很重要。

现在，让我们获取这些分布的参数。 Negative binomial 有两个参数：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： Wikipedia 和 scipy 使用 p 的不同定义，一种将其视为成功的概率，另一种将其视为失败的概率。因此，为了与 scipy 概念保持一致，请使用：

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

UPD 结束

计算可能性：

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

Geometric distribution 也一样：

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

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

【讨论】：

非常感谢。我还有一个问题，如果你能回答我将不胜感激。我知道我的数据集是离散的，但假设我想看看它是否符合正态分布。可能吗？有没有什么方法可以做到这一点，比如矩量法？

另外，我们可以将这种方法应用于混合模型吗？像二项式混合？

@Salih 是的，它适用于像高斯这样的连续分布。然而，在某些情况下，很难甚至不可能估计所有参数。例如，对于两个二项式的混合，您需要三个参数，因此需要三个矩；解决起来已经很不爽了。如果在混合物中添加更多成分，情况会变得更糟。

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

@Salih 对不起，但我不会。它已经偏离了原始问题的范围，因此最好为此发布一个单独的问题。此外，它实际上是五个参数，而不是三个，正如我所提到的，以封闭的形式解决它真的很不愉快。在实践中，混合模型通常使用 EM 算法和高斯来估计。