如何在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分别是样本的最小值和最大值。

当我不做规范化时，一切正常，不同的拟合方法之间有细微的差别，相当好。

但是当我做标准化时，这里是我得到的结果图。

只有矩法（绿线）看起来可以。

无论使用什么参数来生成随机数，scipy.stats.beta.fit（）方法（红线）总是一致的。

MLE（蓝线）失败了。

因此，正常化似乎造成了这些问题。但我认为β分布中有x=0和x=1是合法的。如果给定一个真实世界的问题，那么将样本观测值规范化使其介于[0,1]之间不是第一步吗？在这种情况下，我应该如何拟合曲线？

Answer 1

问题是beta.pdf()有时返回0和1的inf。例如：

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

实际上，您不应该像这样进行规范化，因为您实际上是在抛出两个不合适的数据点。

Answer 2

如果没有用于beta.fit的docstring，查找有点困难，但是如果您知道要对beta.fit强制的上限和下限，则可以使用kwargs floc和fscale。

我只使用beta.fit方法运行了您的代码，但是使用和不使用floc和fscale kwargs。另外，我用int和float形式的参数检查了它，以确保这不会影响您的答案。在这次测试中没有。我不能说它是否永远不会。）

>>> from scipy.stats import beta
>>> 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

>>> data=beta.rvs(5,2,loc=0,scale=1,size=500)
>>> beta.fit(data)
(5.696963536654355, 2.0005252702837009, -0.060443307228404922, 1.0580278414086459)
>>> beta.fit(data,floc=0,fscale=1)
(5.0952451826831462, 1.9546341057106007, 0, 1)
>>> beta.fit(data,floc=0.,fscale=1.)
(5.0952451826831462, 1.9546341057106007, 0.0, 1.0)

总之，这似乎不会改变您的数据（通过规范化）或抛出数据。我只是觉得应该注意的是，使用这个的时候要小心。在您的例子中，您知道限制是0和1，因为您从定义的分布中获取的数据介于0和1之间。在其他情况下，限制可能是已知的，但如果它们是未知的，beta.fit将提供它们。在这种情况下，在不指定0和1的限制的情况下，beta.fit将它们计算为loc=-0.06和scale=1.058。