Почему функция Scipy's fit () для гамма-распределения создает совершенно другое распределение? - PullRequest
Новая
       58

Почему функция Scipy's fit () для гамма-распределения создает совершенно другое распределение?

0 голосов
/ 23 мая 2019

Я пытаюсь использовать функцию подгонки scipy для подгонки гамма-распределения к наблюдаемым данным.Гистограмма для наблюдаемых данных ниже:

enter image description here

Среднее значение и дисперсия 

print("mean",np.mean(observed_data)) # mean 0.427611176580073
print("Var",np.var(observed_data)) # Var 0.6898193689790143

Однако, если я использую scipy.stats.gamma.fit(), чтобы согласовать эти наблюдаемые данные с гамма-распределением, а затем снова произвести выборку из этого распределения, среднее значение и дисперсия совершенно разные:

enter image description here

Я знаючто scipy подходит на основе MLE, но я не понимаю интуицию, объясняющую, почему эти ключевые статистические данные настолько отклонены - среднее значение и дисперсия совершенно разные.На самом деле, я могу получить гораздо лучшие результаты, просто запустив это через свой собственный решатель: 

from scipy.optimize import fsolve
from typing import List
def fit_gamma_distribution(data: List[float]):
    mean = np.mean(data)
    variance = np.var(data)

    def equations(p):
        k, theta = p
        return (k * theta - mean, k * theta **2 - variance)
    solved_k, solved_theta = fsolve(equations, (1,1))

    if np.isclose(np.array([solved_k * solved_theta]), np.array([mean]), rtol=0.01):
        return fsolve(equations, (1,1))

k, theta = fit_gamma_distribution(observed_data)

new_dist = np.random.gamma(shape=k, scale=
theta, size=len(observed_data))
plt.hist(new_dist, alpha= 0.5, bins=40)
plt.hist(observed_data, alpha=0.2, bins=40)
plt.xlim(0,5)
plt.title(f"New sampled distribution: μ = {round(np.mean(new_dist),2)} observed μ = {round(np.mean(observed_data), 2)}")

enter image description here

Это гораздо лучше, чем scipy один на мой взгляд.Почему это так?

...