如何在Python中正确拟合Beta分布？

Question

我正在尝试找到正确的拟合beta分布的方法。这不是一个真实的问题，我只是在测试几种不同方法的效果，在做这件事时有些困惑。

这是我正在处理的Python代码，我测试了三种不同的途径： 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（蓝线）失败了。

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

- Jason

---

科学家是否会在操作符之间使用空格来格式化他们的代码...或者他们只是太忙了 :) - Erik Kaplun

@Ffisegydd 感谢您的帮助。 - Jason

Answer 1

问题在于beta.pdf()有时会返回0和inf，对于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。虽然您“纠正”了概率密度函数为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

但你不应该像这样进行规范化，因为这样做实际上是在逐步淘汰拟合中的两个数据点。

- user545424

---

谢谢您的回答，很有道理。但是还应该使用哪些其他规范化呢？ - Jason

@Jason 我有同样的问题。 - Chi

Answer 2

没有为 beta.fit 编写 docstring，所以有点难以找到，但如果您知道要强制施加的上限和下限，就可以使用 kwargs floc 和 fscale。我只使用了 beta.fit 方法运行了您的代码，但是使用和不使用 floc 和 fscale kwargs 进行了检查。此外，我使用整数和浮点数作为参数进行了检查，以确保这不会影响您的答案（在这个测试中是不会的，但不能保证以后都不会）。

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

- jdj081

---

通用的*.fit()文档：https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.rv_continuous.fit.html#scipy.stats.rv_continuous.fit - mathisfun

Answer 3

我使用了doi:10.1080/00949657808810232中提出的方法来拟合beta参数：

from scipy.special import psi
from scipy.special import polygamma
from scipy.optimize import root_scalar
from numpy.random import beta
import numpy as np

def ipsi(y):
    if y >= -2.22:
        x = np.exp(y) + 0.5
    else:
        x = - 1/ (y + psi(1))  
    for i in range(5):
        x = x - (psi(x) - y)/(polygamma(1,x))
    return x
        
#%%
# q satisface
# psi(q) - psi(ipsi(lng1 - lng2 + psi(q)) + q) -lng2 = 0 
# O sea, busco raíz de 
# f(q) = psi(q) - psi(ipsi(lng1 - lng2 + psi(q)) + q) -lng2
# luego:
# p = ipsi(lng1 - lng2 + psi(q))
def f(q,lng1,lng2):
    return psi(q) - psi(ipsi(lng1 - lng2 + psi(q)) + q) -lng2

#%%
def ml_beta_pq(sample):
    lng1 = np.log(sample).mean()
    lng2 = np.log(1-sample).mean()
    def g(q):
        return f(q,lng1,lng2)
    q=root_scalar(g,x0=1,x1=1.1).root
    p = ipsi(lng1 - lng2 + psi(q))
    return p, q

#%%
p = 2
q = 5
n = 1500
sample = beta(p,q,n)
ps,qs = ml_beta_pq(sample) #s de sombrero

print(f'Estimación de parámetros de una beta({p}, {q}) \na partir de una muestra de tamaño n = {n}')
print(f'\nn ={n:5d} |   p   |   q')
print(f'---------+-------+------')
print(f'original | {p:2.3f} | {q:2.3f}')
print(f'estimado | {ps:2.3f} | {qs:2.3f}')

- rgrimson

---

你能解释一下这段代码吗？它应该是有信息量的，而不需要阅读任何文章。此外，你在代码中使用了西班牙语，这对某些人来说可能非常困惑，请考虑将其翻译成英语，因为这是英文 SO :) - Ruli