用scipy拟合具有固定均值的伽马分布？

编程问答 2022-06-04

问题描述

scipy.stats.rv_continuous.fit，可让您在拟合分布时固定参数，但这取决于scipy的参数化选择。由于伽马分布使用k，theta（形状，比例）参数化，因此例如在保持theta恒定的情况下很容易拟合。我想适合我知道平均值的数据集，但是观察到的平均值可能会由于采样误差而有所不同。如果scipy使用使用mu = k * theta而不是theta的参数化，这将很容易。有办法让scipy做到这一点吗？如果没有，还有另一个库可以吗？

下面是一些带有数据集的示例代码，其观察到的平均值为9.952，但我知道基础分布的实际平均值为11：

from scipy.stats import gamma
observations = [17.6,24.9,3.9,17.6,11.8,10.4,4.1,11.7,5.7,1.6,8.6,12.9,8.0,7.4,1.2,11.3,1.0,1.9,6.0,9.3,13.3,5.4,9.1,4.0,12.8,11.1,23.1,4.2,7.9,10.0,3.4,27.8,7.2,14.9,2.9,5.5,7.0,12.3,10.6,22.1,5.0,21.3,15.9,34.5,8.1,19.6,10.8,13.4,22.8,27.6,6.8,5.9,9.0,7.1,21.2,14.6,16.9,6.5,14.1,15.2,7.8,4.9,2.1,9.5,5.6,7.7,18.3,3.8,11.0,12.5,8.4,3.2,2.0,24.7,24.6,4.3,7.6,8.3,14.5,14.0,9.0]
shape,_,scale = gamma.fit(observations,floc = 0)
print(shape*scale)

这给了

9.952

但是我想要适合shape*scale = 11.0

解决方法

SciPy分布的fit方法提供了参数的最大似然估计。您是正确的，它仅适合形状，位置和比例。 （实际上，您说的是形状和比例，但SciPy还包含一个位置参数。有时称为三参数伽马分布。）

对于SciPy中的大多数分布，fit方法使用数值优化器来最小化nnlf方法中定义的负对数似然率。除了使用fit方法之外，您还可以只用几行代码来完成此操作。这样，您就可以仅使用一个参数（例如形状k）创建目标函数，然后在该函数中设置theta = mean/k，其中mean是期​​望的均值，然后调用{{1 }来评估负对数可能性。这是您可以做到的一种方法：

gamma.nnlf

此脚本打印

import numpy as np
from scipy.stats import gamma
from scipy.optimize import fmin


def nll(k,mean,x):
    return gamma.nnlf(np.array([k[0],mean/k[0]]),x)


observations = [17.6,24.9,3.9,17.6,11.8,10.4,4.1,11.7,5.7,1.6,8.6,12.9,8.0,7.4,1.2,11.3,1.0,1.9,6.0,9.3,13.3,5.4,9.1,4.0,12.8,11.1,23.1,4.2,7.9,10.0,3.4,27.8,7.2,14.9,2.9,5.5,7.0,12.3,10.6,22.1,5.0,21.3,15.9,34.5,8.1,19.6,10.8,13.4,22.8,27.6,6.8,5.9,9.0,7.1,21.2,14.6,16.9,6.5,14.1,15.2,7.8,4.9,2.1,9.5,5.6,7.7,18.3,3.8,11.0,12.5,8.4,3.2,2.0,24.7,24.6,4.3,7.6,8.3,14.5,14.0,9.0]


# This is the desired mean of the distribution.
mean = 11

# Initial guess for the shape parameter.
k0 = 3.0
opt = fmin(nll,k0,args=(mean,np.array(observations)),xtol=1e-11,disp=False)

k_opt = opt[0]
theta_opt = mean / k_opt
print(f"k_opt:     {k_opt:9.7f}")
print(f"theta_opt: {theta_opt:9.7f}")

或者，可以修改wikipedia中所示的对数似然极值的一阶条件，以便只有一个参数k_opt: 1.9712604 theta_opt: 5.5801861 。然后可以将极值的条件实现为一个标量方程，其根可以用k找到。以下是上述使用此技术的脚本的变体。

scipy.optimize.fsolve

输出：

import numpy as np
from scipy.special import digamma
from scipy.optimize import fsolve


def first_order_eq(k,x):
    mean_logx = np.mean(np.log(x))
    return (np.log(k) - digamma(k) + mean_logx - np.mean(x)/mean
            - np.log(mean) + 1)


observations = [17.6,9.0]


# This is the desired mean of the distribution.
mean = 11

# Initial guess for the shape parameter.
k0 = 3
sol = fsolve(first_order_eq,observations),xtol=1e-11)

k_opt = sol[0]
theta_opt = mean / k_opt
print(f"k_opt:     {k_opt:9.7f}")
print(f"theta_opt: {theta_opt:9.7f}")
pythonscipyscipyscipy.stats