python - 来自广义极值 (GEV) 最大似然拟合数据的奇怪 pdf

Question

我正在进行一些数据分析，涉及将数据集拟合为广义极值 (GEV) 分布，但我得到了一些奇怪的结果。这是我正在做的:

from scipy.stats import genextreme as gev
import numpy
data = [1.47, 0.02, 0.3, 0.01, 0.01, 0.02, 0.02, 0.12, 0.38, 0.02, 0.15, 0.01, 0.3, 0.24, 0.01, 0.05, 0.01, 0.0, 0.06, 0.01, 0.01, 0.0, 0.05, 0.0, 0.09, 0.03, 0.22, 0.0, 0.1, 0.0]
x = numpy.linspace(0, 2, 20)
pdf = gev.pdf(x, *gev.fit(data))
print(pdf)

输出:

array([  5.64759709e+05,   2.41090345e+00,   1.16591714e+00,
         7.60085002e-01,   5.60415578e-01,   4.42145248e-01,
         3.64144425e-01,   3.08947114e-01,   2.67889183e-01,
         2.36190826e-01,   2.11002185e-01,   1.90520108e-01,
         1.73548832e-01,   1.59264573e-01,   1.47081601e-01,
         1.36572220e-01,   1.27416958e-01,   1.19372442e-01,
         1.12250072e-01,   1.05901466e-01,   1.00208313e-01,
         9.50751375e-02,   9.04240603e-02,   8.61909342e-02,
         8.23224528e-02,   7.87739599e-02,   7.55077677e-02,
         7.24918532e-02,   6.96988348e-02,   6.71051638e-02,
         6.46904782e-02,   6.24370827e-02,   6.03295277e-02,
         5.83542648e-02,   5.64993643e-02,   5.47542808e-02,
         5.31096590e-02,   5.15571710e-02,   5.00893793e-02,
         4.86996213e-02,   4.73819114e-02,   4.61308575e-02,
         4.49415891e-02,   4.38096962e-02,   4.27311763e-02,
         4.17023886e-02,   4.07200140e-02,   3.97810205e-02,
         3.88826331e-02,   3.80223072e-02])

问题是第一个值很大，完全扭曲了所有结果，它在图中显示得很清楚:

我试验过其他数据和随机样本，在某些情况下它是有效的。我的数据集中的第一个值明显高于其余值，但它是一个有效值，所以我不能就这样放弃它。

有人知道为什么会这样吗？

更新

这是另一个更清楚地说明问题的例子:

In [1]: from scipy.stats import genextreme as gev, kstest

In [2]: data = [0.01, 0.0, 0.28, 0.0, 0.0, 0.0, 0.01, 0.0, 0.0, 0.13, 0.07, 0.03
, 0.01, 0.42, 0.11, 0.0, 0.0, 0.0, 0.0, 0.25, 0.0, 0.0, 0.26, 1.32, 0.06, 0.02,
1.57, 0.07, 1.56, 0.04]

In [3]: fit = gev.fit(data)

In [4]: kstest(data, 'genextreme', fit)
Out[4]: (0.48015007915450658, 6.966510064376763e-07)

In [5]: x = linspace(0, 2, 200)

In [6]: plot(x, gev.pdf(x, *fit))
Out[6]: [<matplotlib.lines.Line2D at 0x97590f0>]

In [7]: hist(data)

请特别注意，第 4 行显示的 p 值约为 7e-7，远低于通常认为可接受的值。这是制作的情节:

Answer 1

首先，我认为您可能希望将位置参数固定为 0。

其次，您的数据中有零，结果拟合可能有 +inf pdf 在 x=0 例如对于 GEV 拟合或 Weibull 拟合。因此，拟合实际上是正确的，但是当您绘制 pdf(包括 x=0)时，生成的图会失真。

第三，我真的认为 scipy 应该放弃对 x=0 的支持，以支持 Weibull 等一系列分布。对于 x=0，R 给出了一个很好的警告:error in fitdistr(data, "weibull") : Weibull values must be > 0，这很有帮助。

In [103]:

p=ss.genextreme.fit(data, floc=0)
ss.genextreme.fit(data, floc=0)
Out[103]:
(-1.372872096699608, 0, 0.011680600795499299)
In [104]:

plt.hist(data, bins=20, normed=True, alpha=0.7, label='Data')
plt.plot(np.linspace(5e-3, 1.6, 100),
         ss.genextreme.pdf(np.linspace(5e-3, 1.6, 100), p[0], p[1], p[2]), 'r--',
         label='GEV Fit')
plt.legend(loc='upper right')
plt.savefig('T1.png')

In [105]:

p=ss.expon.fit(data, floc=0)
ss.expon.fit(data, floc=0)
Out[105]:
(0, 0.14838807003769505)
In [106]:

plt.hist(data, bins=20, normed=True, alpha=0.7, label='Data')
plt.plot(np.linspace(0, 1.6, 100),
         ss.expon.pdf(np.linspace(0, 1.6, 100), p[0], p[1]), 'r--',
         label='Expon. Fit')
plt.legend(loc='upper right')
plt.savefig('T2.png')

In [107]:

p=ss.weibull_min.fit(data[data!=0], floc=0)
ss.weibull_min.fit(data[data!=0], floc=0)
Out[107]:
(0.67366030738733995, 0, 0.10535422201164378)
In [108]:

plt.hist(data[data!=0], bins=20, normed=True, alpha=0.7, label='Data')
plt.plot(np.linspace(5e-3, 1.6, 100),
         ss.weibull_min.pdf(np.linspace(5e-3, 1.6, 100), p[0], p[1], p[2]), 'r--',
         label='Weibull_Min Fit')
plt.legend(loc='upper right')
plt.savefig('T3.png')

---

编辑

您的第二个数据(包含更多 0)是一个很好的例子，当涉及位置参数的 MLE 拟合变得非常具有挑战性时，尤其是可能涉及大量浮点溢出/下溢:

In [122]:
#fit with location parameter fixed, scanning loc parameter from 1e-8 to 1e1
L=[] #stores the Log-likelihood
P=[] #stores the p value of K-S test
for LC in np.linspace(-8, 1, 200):
    fit = gev.fit(data, floc=10**LC)
    L.append(np.log(gev.pdf(data, *fit)).sum())
    P.append(kstest(data, 'genextreme', fit)[1])
L=np.array(L)
P=np.array(P)
In [123]:
#plot log likelihood, a lot of overflow/underflow issues! (see the zigzag line?)
plt.plot(np.linspace(-8, 1, 200), L,'-')
plt.ylabel('Log-Likelihood')
plt.xlabel('$log_{10}($'+'location parameter'+'$)$')

In [124]:
#plot p-value
plt.plot(np.linspace(-8, 1, 200), np.log10(P),'-')
plt.ylabel('$log_{10}($'+'K-S test P value'+'$)$')
plt.xlabel('$log_{10}($'+'location parameter'+'$)$')
Out[124]:
<matplotlib.text.Text at 0x107e3050>

In [128]:
#The best fit between with location paramter between 1e-8 to 1e1 has the loglikelihood of 515.18
np.linspace(-8, 1, 200)[L.argmax()]
fit = gev.fit(data, floc=10**(np.linspace(-8, 1, 200)[L.argmax()]))
np.log(gev.pdf(data, *fit)).sum()
Out[128]:
515.17663678368604
In [129]:
#The simple MLE fit is clearly bad, loglikelihood is -inf
fit0 = gev.fit(data)
np.log(gev.pdf(data, *fit0)).sum()
Out[129]:
-inf

In [133]:
#plot the fit
x = np.linspace(0.005, 2, 200)
plt.plot(x, gev.pdf(x, *fit))
plt.hist(data,normed=True, alpha=0.6, bins=20)
Out[133]:
(array([ 8.91719745,  0.8492569 ,  0.        ,  1.27388535,  0.        ,
        0.42462845,  0.        ,  0.        ,  0.        ,  0.        ,
        0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
        0.        ,  0.42462845,  0.        ,  0.        ,  0.8492569 ]),
 array([ 0.    ,  0.0785,  0.157 ,  0.2355,  0.314 ,  0.3925,  0.471 ,
        0.5495,  0.628 ,  0.7065,  0.785 ,  0.8635,  0.942 ,  1.0205,
        1.099 ,  1.1775,  1.256 ,  1.3345,  1.413 ,  1.4915,  1.57  ]),
 <a list of 20 Patch objects>)

---

编辑，GEV 的拟合优度检验

关于 KS 测试的旁注。您正在测试 GEV 的拟合优度，其参数 ESTIMATED FROM THE DATA 本身。在这种情况下，p 值无效，请参阅:itl.nist.gov/div898/handbook/eda/section3/eda35g.htm

似乎有很多关于 GEV 拟合优度检验主题的研究，我还没有找到任何可用的实现。

http://onlinelibrary.wiley.com/doi/10.1029/98WR02364/abstract http://onlinelibrary.wiley.com/doi/10.1029/91WR00077/abstract http://www.idrologia.polito.it/~laio/articoli/16-WRR%20EDFtest.pdf