python - 如何在 python 中正确地适应 beta 分布？

Question

我正在尝试获得拟合 Beta 分布的正确方法。这不是现实世界中的问题，我只是在测试几种不同方法的效果，而在这样做的过程中，有些事情让我感到困惑。

这是我正在处理的 Python 代码，我在其中测试了 3 种不同的方法:1>:使用矩拟合(样本均值和方差)。2>:通过最小化负对数似然来拟合(通过使用 scipy.optimize.fmin())。3>: 只需调用 scipy.stats.beta.fit()

from scipy.optimize import fmin
from scipy.stats import beta
from scipy.special import gamma as gammaf
import matplotlib.pyplot as plt
import numpy


def betaNLL(param,*args):
    '''Negative log likelihood function for beta
    <param>: list for parameters to be fitted.
    <args>: 1-element array containing the sample data.

    Return <nll>: negative log-likelihood to be minimized.
    '''

    a,b=param
    data=args[0]
    pdf=beta.pdf(data,a,b,loc=0,scale=1)
    lg=numpy.log(pdf)
    #-----Replace -inf with 0s------
    lg=numpy.where(lg==-numpy.inf,0,lg)
    nll=-1*numpy.sum(lg)
    return nll

#-------------------Sample data-------------------
data=beta.rvs(5,2,loc=0,scale=1,size=500)

#----------------Normalize to [0,1]----------------
#data=(data-numpy.min(data))/(numpy.max(data)-numpy.min(data))

#----------------Fit using moments----------------
mean=numpy.mean(data)
var=numpy.var(data,ddof=1)
alpha1=mean**2*(1-mean)/var-mean
beta1=alpha1*(1-mean)/mean

#------------------Fit using mle------------------
result=fmin(betaNLL,[1,1],args=(data,))
alpha2,beta2=result

#----------------Fit using beta.fit----------------
alpha3,beta3,xx,yy=beta.fit(data)

print '\n# alpha,beta from moments:',alpha1,beta1
print '# alpha,beta from mle:',alpha2,beta2
print '# alpha,beta from beta.fit:',alpha3,beta3

#-----------------------Plot-----------------------
plt.hist(data,bins=30,normed=True)
fitted=lambda x,a,b:gammaf(a+b)/gammaf(a)/gammaf(b)*x**(a-1)*(1-x)**(b-1) #pdf of beta

xx=numpy.linspace(0,max(data),len(data))
plt.plot(xx,fitted(xx,alpha1,beta1),'g')
plt.plot(xx,fitted(xx,alpha2,beta2),'b')
plt.plot(xx,fitted(xx,alpha3,beta3),'r')

plt.show()

我遇到的问题是关于规范化过程 (z=(x-a)/(b-a))，其中 a 和 b 是分别为样本的最小值和最大值。

当我不做归一化时，一切正常，不同的拟合方法之间存在细微差异，还算不错。

但是当我进行归一化时，这是我得到的结果图。

只有 moment 方法(绿线)看起来没问题。

无论我使用什么参数生成随机数，scipy.stats.beta.fit() 方法(红线)始终是统一的。

MLE(蓝线)失败。

因此，规范化似乎造成了这些问题。但我认为在 beta 发行版中使用 x=0 和 x=1 是合法的。如果给定一个现实世界的问题，难道不是将样本观察值标准化以使其介于 [0,1] 之间的第一步吗？那么，我该如何拟合曲线呢？

Answer 1

问题是 beta.pdf() 有时返回 0 and inf for 0 和 1。例如:

>>> from scipy.stats import beta
>>> beta.pdf(1,1.05,0.95)
/usr/lib64/python2.6/site-packages/scipy/stats/distributions.py:1165: RuntimeWarning: divide by zero encountered in power
  Px = (1.0-x)**(b-1.0) * x**(a-1.0)
inf
>>> beta.pdf(0,1.05,0.95)
0.0

您保证在规范化过程中您将在 0 和 1 处获得一个数据样本。尽管您“更正”了 pdf 为 0 的值，但您并未更正那些返回 inf 的值。为了解决这个问题，您可以删除所有不是有限的值:

def betaNLL(param,*args):
    """
    Negative log likelihood function for beta
    <param>: list for parameters to be fitted.
    <args>: 1-element array containing the sample data.

    Return <nll>: negative log-likelihood to be minimized.
    """

    a, b = param
    data = args[0]
    pdf = beta.pdf(data,a,b,loc=0,scale=1)
    lg = np.log(pdf)
    mask = np.isfinite(lg)
    nll = -lg[mask].sum()
    return nll

不过，实际上您不应该像这样进行归一化，因为您实际上是在将两个数据点不合身。