您可以使用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

注意：我用panda代替numpy。这是因为numpy的var()和std()不适用于Bessel's correction，而panda的适用于Bessel's correction。如果你有100多个样本，应该不会有太大的差异，但在较小的样本上，这可能很重要。
现在，让我们来获取这些分布的参数。Negative binomial有两个参数：让我们估计它们并计算数据集的似然性：

# 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结束
更新版本2：

我建议使用@thilak的代码log likelihood来代替，它可以避免精度损失，这在大样本情况下尤其重要。

UPD 2结束

计算可能性：

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

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

马拉特的回答很棒。
除了马拉特的帖子，我肯定会推荐取概率密度函数的对数。一些关于为什么对数似然比似然更受欢迎的信息-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
• 用和代替积

为了找出最佳的分布-

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