scipy 中的负二项式参数化通过均值和标准差

parameterization of the negative binomial in scipy via mean and std

我正在尝试使用scipy中的包使我的数据符合负二项分布 ]Python。但是,我的验证似乎失败了。

这些是我的步骤:

  1. 我有一些需求数据,是用统计描述的:
mu = 1.4
std = 1.59
print(mu, std)
  1. 我使用下面的参数化函数,取自这个 post 来计算两个 NB 参数。
def convert_params(mu, theta):
    """
    Convert mean/dispersion parameterization of a negative binomial to the ones scipy supports

    See https://en.wikipedia.org/wiki/Negative_binomial_distribution#Alternative_formulations
    """
    r = theta
    var = mu + 1 / r * mu ** 2
    p = (var - mu) / var
    return r, 1 - p

我通过了(希望正确...)我的两个统计数据 - 不同来源之间的命名约定在这一点上相当混乱 prk

firstParam, secondParam = convert_params(mu, std)
  1. 然后我将使用这两个参数来拟合分布:
from scipy.stats import nbinom

rv = nbinom(firstParam, secondParam)

然后我用 百分比函数 .ppf(0.95) 计算出一个值 R。我的问题上下文中的值 R 是再订购点。

R = rv.ppf(0.95)
  1. 现在我希望验证前面的步骤,但我没能检索我的原始统计数据 mustd 分别对应 meanmath.sqrt(var)
import math

mean, var = nbinom.stats(firstParam, secondParam, moments='mv')
print(mean, math.sqrt(var))

我错过了什么?关于 Scipy 中实施的参数化的任何反馈?

您似乎使用了不同的转换方式。引用 wikipedia section 处的最后一个项目符号给出了如下所示的公式。使用这些公式,您可以得到完全相同的 mustd:

import numpy as np
from scipy.stats import nbinom

def convert_mu_std_to_r_p(mu, std):
    r = mu ** 2 / (std ** 2 - mu)
    p = 1 - mu / std ** 2
    return r, 1 - p

mu = 1.4
std = 1.59
print("mu, std:", mu, std)
firstParam, secondParam = convert_mu_std_to_r_p(mu, std)
mean, var = nbinom.stats(firstParam, secondParam, moments='mv')
print("mean, sqrt(var):", mean, np.sqrt(var))

rv = nbinom(firstParam, secondParam)
print("reorder point:", rv.ppf(0.95))

输出:

mu, std: 1.4 1.59
mean, sqrt(var): 1.4 1.59
reorder point: 5.0

转换代码是错误的,我相信,SciPy 不是使用 Wiki 约定,而是 Mathematica 约定

#%%
import numpy as np
from scipy.stats import nbinom

def convert_params(mean, std):
    """
    Convert mean/dispersion parameterization of a negative binomial to the ones scipy supports

    See https://mathworld.wolfram.com/NegativeBinomialDistribution.html
    """
    p = mean/std**2
    n = mean*p/(1.0 - p)
    return n, p

mean = 1.4
std  = 1.59

n, p = convert_params(mean, std)

print((n, p))

#%%

m, v = nbinom.stats(n, p, moments='mv')
print(m, np.sqrt(v))

代码打印回 1.4、1.59 对

再订货点计算为

rv = nbinom(n, p)
print("reorder point:", rv.ppf(0.95))

输出5