如何在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 and inf for0和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

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