параметризация отрицательного бинома в scipy через среднее и стандартное - PullRequest
Новая
       149

параметризация отрицательного бинома в scipy через среднее и стандартное

1 голос
/ 18 июня 2020

Я пытаюсь подогнать свои данные к отрицательному биномиальному распределению с пакетом scipy в Python. Однако моя проверка, похоже, не удалась.

Это мои шаги:

  1. У меня есть данные о спросе, которые описываются статистикой: 
mu = 1.4
std = 1.59
print(mu, std)
Я использую функцию параметризации ниже, взятую из этого сообщения для вычисления двух параметров 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

Я передаю (надеюсь, правильно ...) две мои статистические данные - соглашение об именах между разными источниками на данном этапе довольно сбивает с толку p, r, k

firstParam, secondParam = convert_params(mu, std)
Затем я бы использовал эти два параметра для соответствия распределению: 
from scipy.stats import nbinom

rv = nbinom(firstParam, secondParam)

Затем я вычисляю значение R с помощью функции процентной точки .ppf(0.95). Значение R в контексте моей проблемы является точкой изменения порядка. 

R = rv.ppf(0.95)
Сейчас я ожидаю проверки предыдущих шагов , но мне не удается получить мою исходную статистику mu и std с mean и math.sqrt(var) соответственно. 
import math

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

Что мне не хватает? Любые отзывы о параметризации, реализованной в Scipy?

Ответы [ 2 ]

1 голос
/ 19 июня 2020

Код преобразования неверен, я считаю, 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

1 голос
/ 19 июня 2020

Похоже, вы используете другое преобразование. Последний пункт в процитированном разделе википедии дает формулы, показанные ниже. С помощью этих формул вы получите точно такие же mu и std: 

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
...