使用Scipy拟合Weibull分布

Question

我试图重新创建最大似然分布拟合,我已经可以在Matlab和R中做到这一点,但现在我想使用scipy.特别是,我想估计我的数据集的Weibull分布参数.

我试过这个:

import scipy.stats as s
import numpy as np
import matplotlib.pyplot as plt

def weib(x,n,a):
    return (a / n) * (x / n)**(a - 1) * np.exp(-(x / n)**a)

data = np.loadtxt("stack_data.csv")

(loc, scale) = s.exponweib.fit_loc_scale(data, 1, 1)
print loc, scale

x = np.linspace(data.min(), data.max(), 1000)
plt.plot(x, weib(x, loc, scale))
plt.hist(data, data.max(), normed=True)
plt.show()

得到这个:

(2.5827280639441961, 3.4955032285727947)

并且看起来像这样的分布:

我一直exponweib在阅读http://www.johndcook.com/distributions_scipy.html.我也尝试了scipy中的其他Weibull函数(以防万一!).

在Matlab(使用分布拟合工具 - 参见屏幕截图)和R(使用MASS库函数fitdistr和GAMLSS包)中,我得到(loc)和b(比例)参数更像1.58463497 5.93030013.我相信所有三种方法都使用最大似然法进行分布拟合.

如果你想去,我已经在这里发布了我的数据!为了完整起见,我使用的是Python 2.7.5,Scipy 0.12.0,R 2.15.2和Matlab 2012b.

为什么我会得到不同的结果!？

Answer 1

My guess is that you want to estimate the shape parameter and the scale of the Weibull distribution while keeping the location fixed. Fixing loc assumes that the values of your data and of the distribution are positive with lower bound at zero.

floc=0 keeps the location fixed at zero, f0=1 keeps the first shape parameter of the exponential weibull fixed at one.

>>> stats.exponweib.fit(data, floc=0, f0=1)
[1, 1.8553346917584836, 0, 6.8820748596850905]
>>> stats.weibull_min.fit(data, floc=0)
[1.8553346917584836, 0, 6.8820748596850549]

The fit compared to the histogram looks ok, but not very good. The parameter estimates are a bit higher than the ones you mention are from R and matlab.

Update

The closest I can get to the plot that is now available is with unrestricted fit, but using starting values. The plot is still less peaked. Note values in fit that don't have an f in front are used as starting values.

>>> from scipy import stats
>>> import matplotlib.pyplot as plt
>>> plt.plot(data, stats.exponweib.pdf(data, *stats.exponweib.fit(data, 1, 1, scale=02, loc=0)))
>>> _ = plt.hist(data, bins=np.linspace(0, 16, 33), normed=True, alpha=0.5);
>>> plt.show()

Answer 2

很容易验证哪个结果是真正的MLE,只需要一个简单的函数来计算对数似然:

>>> def wb2LL(p, x): #log-likelihood
    return sum(log(stats.weibull_min.pdf(x, p[1], 0., p[0])))
>>> adata=loadtxt('/home/user/stack_data.csv')
>>> wb2LL(array([6.8820748596850905, 1.8553346917584836]), adata)
-8290.1227946678173
>>> wb2LL(array([5.93030013, 1.57463497]), adata)
-8410.3327470347667

来自和R (@Warren)的fit方法的结果更好并且具有更高的对数似然性.它更可能是真正的MLE.GAMLSS的结果不同并不奇怪.它是一个完全不同的统计模型:广义加法模型.exponweibfitdistr

还是不相信？我们可以围绕MLE绘制2D置信限制图,详见Meeker和Escobar的书.

同样,这验证array([6.8820748596850905, 1.8553346917584836])了正确的答案,因为对数似然低于参数空间中的任何其他点.注意:

>>> log(array([6.8820748596850905, 1.8553346917584836]))
array([ 1.92892018,  0.61806511])

BTW1,MLE拟合可能看起来不符合分布直方图.考虑MLE的一种简单方法是,MLE是给定观测数据最可能的参数估计.它不需要在视觉上很好地拟合直方图,这将是最小化均方误差的东西.

BTW2,您的数据似乎是leptokurtic和左倾斜,这意味着Weibull分布可能不适合您的数据.尝试,例如Gompertz-Logistic,它将对数可能性提高了大约100. 干杯!

Answer 3

我知道这是一个老帖子,但我刚遇到了类似的问题,这个帖子帮我解决了.认为我的解决方案可能对像我这样的人有帮助:

# Fit Weibull function, some explanation below
params = stats.exponweib.fit(data, floc=0, f0=1)
shape = params[1]
scale = params[3]
print 'shape:',shape
print 'scale:',scale

#### Plotting
# Histogram first
values,bins,hist = plt.hist(data,bins=51,range=(0,25),normed=True)
center = (bins[:-1] + bins[1:]) / 2.

# Using all params and the stats function
plt.plot(center,stats.exponweib.pdf(center,*params),lw=4,label='scipy')

# Using my own Weibull function as a check
def weibull(u,shape,scale):
    '''Weibull distribution for wind speed u with shape parameter k and scale parameter A'''
    return (shape / scale) * (u / scale)**(shape-1) * np.exp(-(u/scale)**shape)

plt.plot(center,weibull(center,shape,scale),label='Wind analysis',lw=2)
plt.legend()

一些额外的信息,帮助我理解:

Scipy Weibull函数可以采用四个输入参数:(a,c),loc和scale.你想修复loc和第一个形状参数(a),这是用floc = 0,f0 = 1完成的.然后拟合将给出参数c和比例,其中c对应于双参数威布尔分布的形状参数(通常用于风数据分析),并且比例对应于其比例因子.

来自docs:

exponweib.pdf(x, a, c) =
    a * c * (1-exp(-x**c))**(a-1) * exp(-x**c)*x**(c-1)

如果a是1,那么

exponweib.pdf(x, a, c) =
    c * (1-exp(-x**c))**(0) * exp(-x**c)*x**(c-1)
  = c * (1) * exp(-x**c)*x**(c-1)
  = c * x **(c-1) * exp(-x**c)

由此,与"风分析"威布尔函数的关系应该更加明确

Answer 4

我对你的问题感到好奇,尽管这不是一个答案,但它将Matlab结果与你的结果和使用的结果进行了比较,结果leastsq显示与给定数据的最佳相关性:

代码如下:

import scipy.stats as s
import numpy as np
import matplotlib.pyplot as plt
import numpy.random as mtrand
from scipy.integrate import quad
from scipy.optimize import leastsq

## my distribution (Inverse Normal with shape parameter mu=1.0)
def weib(x,n,a):
    return (a / n) * (x / n)**(a-1) * np.exp(-(x/n)**a)

def residuals(p,x,y):
    integral = quad( weib, 0, 16, args=(p[0],p[1]) )[0]
    penalization = abs(1.-integral)*100000
    return y - weib(x, p[0],p[1]) + penalization

#
data = np.loadtxt("stack_data.csv")


x = np.linspace(data.min(), data.max(), 100)
n, bins, patches = plt.hist(data,bins=x, normed=True)
binsm = (bins[1:]+bins[:-1])/2

popt, pcov = leastsq(func=residuals, x0=(1.,1.), args=(binsm,n))

loc, scale = 1.58463497, 5.93030013
plt.plot(binsm,n)
plt.plot(x, weib(x, loc, scale),
         label='weib matlab, loc=%1.3f, scale=%1.3f' % (loc, scale), lw=4.)
loc, scale = s.exponweib.fit_loc_scale(data, 1, 1)
plt.plot(x, weib(x, loc, scale),
         label='weib stack, loc=%1.3f, scale=%1.3f' % (loc, scale), lw=4.)
plt.plot(x, weib(x,*popt),
         label='weib leastsq, loc=%1.3f, scale=%1.3f' % tuple(popt), lw=4.)

plt.legend(loc='upper right')
plt.show()