scipy Gamma累积分布函数的正确拟合

xmq68pz9 于 2023-03-02 发布在其他

关注(0)|答案(1)|浏览(196)

我有两个Numpy数组（两者总共210个条目）的降雨量值，一个是观测值，另一个是预测值。我的目标是创建一个最适合的伽马CDF（我第一次深入研究gamma CDF），并确定随后提供的值将落入的相关百分位数。下图提供了我尝试使用这两个数组创建的gamma CDF的更简单的图形参考。需要注意的是，y轴引用直方图中每个值的百分位数，因此范围从第1位到第99位：

这些数组如下所示：

guess = [0.02 0.03 0.02 0.04 0.01 0.01 0.04 0.01 0.   0.   0.01 0.03 0.03 0.04
 0.05 0.03 0.   0.02 0.03 0.03 0.04 0.03 0.04 0.04 0.04 0.04 0.01 0.01
 0.01 0.03 0.04 0.03 0.02 0.05 0.03 0.   0.   0.04 0.05 0.03 0.05 0.03
 0.03 0.   0.01 0.02 0.01 0.05 0.01 0.05 0.05 0.04 0.04 0.02 0.02 0.04
 0.04 0.04 0.02 0.04 0.02 0.03 0.04 0.04 0.   0.15 0.07 0.08 0.15 0.08
 0.13 0.14 0.07 0.13 0.13 0.08 0.14 0.1  0.08 0.12 0.14 0.11 0.15 0.14
 0.14 0.16 0.15 0.15 0.06 0.1  0.1  0.09 0.09 0.11 0.07 0.12 0.11 0.15
 0.06 0.11 0.09 0.09 0.08 0.09 0.12 0.07 0.07 0.09 0.12 0.16 0.13 0.11
 0.1  0.08 0.13 0.06 0.09 0.13 0.16 0.12 0.23 0.35 0.33 0.28 0.24 0.33
 0.25 0.25 0.24 0.25 0.28 0.28 0.34 0.24 0.33 0.17 0.25 0.24 0.35 0.24
 0.24 0.22 0.29 0.23 0.2  0.32 0.25 0.25 0.33 0.21 0.18 0.22 0.27 0.18
 0.25 0.22 0.29 0.27 0.33 0.2  0.31 0.29 0.17 0.17 0.29 0.39 0.65 0.84
 0.71 0.64 0.52 0.91 0.82 0.36 0.37 0.95 0.87 0.73 0.67 0.73 0.8  0.91
 0.63 0.58 0.6  0.75 0.53 0.88 0.84 0.98 1.2  1.2  1.02 1.02 1.17 1.14
 1.02 1.13 1.15 1.25 1.03 1.04 1.25 1.12 1.02 1.26 1.44 1.33 1.33 1.49]

actual = [0.04 0.03 0.03 0.02 0.04 0.01 0.03 0.02 0.01 0.01 0.04 0.01 0.   0.05
 0.03 0.03 0.05 0.04 0.02 0.04 0.02 0.01 0.05 0.   0.01 0.05 0.01 0.02
 0.04 0.   0.01 0.01 0.04 0.04 0.03 0.01 0.03 0.04 0.   0.03 0.03 0.05
 0.05 0.01 0.05 0.05 0.03 0.02 0.02 0.05 0.04 0.05 0.04 0.04 0.01 0.03
 0.02 0.01 0.01 0.   0.03 0.02 0.05 0.03 0.04 0.13 0.06 0.07 0.14 0.11
 0.1  0.15 0.14 0.15 0.07 0.13 0.08 0.07 0.07 0.1  0.15 0.1  0.11 0.08
 0.09 0.06 0.15 0.12 0.1  0.12 0.14 0.16 0.16 0.11 0.07 0.06 0.15 0.1
 0.15 0.14 0.14 0.09 0.13 0.13 0.15 0.09 0.11 0.11 0.13 0.15 0.14 0.12
 0.12 0.06 0.08 0.13 0.07 0.16 0.09 0.1  0.21 0.17 0.27 0.24 0.33 0.24
 0.28 0.28 0.19 0.17 0.29 0.27 0.22 0.35 0.19 0.28 0.3  0.33 0.29 0.31
 0.17 0.27 0.34 0.26 0.22 0.3  0.22 0.22 0.32 0.34 0.21 0.21 0.3  0.19
 0.27 0.22 0.19 0.23 0.26 0.33 0.23 0.31 0.18 0.34 0.35 0.55 0.76 0.37
 0.92 0.86 0.72 0.78 0.54 0.7  0.4  0.45 0.37 1.   0.48 0.92 0.45 0.57
 0.55 0.56 0.75 0.5  0.41 0.71 0.82 0.73 1.04 1.17 1.17 1.09 1.06 1.04
 1.14 1.18 1.09 1.03 1.08 1.16 1.09 1.12 1.22 1.32 1.38 1.39 1.37 1.37]

我已经为这两个数组创建了一个直方图，以0.05为增量进行分箱，总共有30个分箱，从上面提供的数据实现这一点的代码片段如下：

rngst = 0.00
rngend = 1.50
gushist = np.histogram(guess, bins = [round(x, 2) for x in np.arange(rngst,(rngend + 0.05),0.05)])
acthist = np.histogram(actual, bins = [round(x, 2) for x in np.arange(rngst,(rngend + 0.05),0.05)])

我还绘制了这两个直方图，如下所示：

我不确定从这里开始要为两个数组创建最适合的gamma CDF，尽管我最初在scipy中找到了一个stats.gamma函数。

scipy

来源：https://stackoverflow.com/questions/75598344/properly-fitting-a-gamma-cumulative-distribution-function

1条答案

按热度按时间

tv6aics11#

使用scipy中为此目的而设计的内置函数。此外，直方图不如ECPDF直观，ECPDF显示每个数据点，并且更容易与拟合CDF进行比较：

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

guess = (
    0.02, 0.03, 0.02, 0.04, 0.01, 0.01, 0.04, 0.01, 0.00, 0.00, 0.01, 0.03, 0.03, 0.04,
    0.05, 0.03, 0.00, 0.02, 0.03, 0.03, 0.04, 0.03, 0.04, 0.04, 0.04, 0.04, 0.01, 0.01,
    0.01, 0.03, 0.04, 0.03, 0.02, 0.05, 0.03, 0.00, 0.00, 0.04, 0.05, 0.03, 0.05, 0.03,
    0.03, 0.00, 0.01, 0.02, 0.01, 0.05, 0.01, 0.05, 0.05, 0.04, 0.04, 0.02, 0.02, 0.04,
    0.04, 0.04, 0.02, 0.04, 0.02, 0.03, 0.04, 0.04, 0.00, 0.15, 0.07, 0.08, 0.15, 0.08,
    0.13, 0.14, 0.07, 0.13, 0.13, 0.08, 0.14, 0.10, 0.08, 0.12, 0.14, 0.11, 0.15, 0.14,
    0.14, 0.16, 0.15, 0.15, 0.06, 0.10, 0.10, 0.09, 0.09, 0.11, 0.07, 0.12, 0.11, 0.15,
    0.06, 0.11, 0.09, 0.09, 0.08, 0.09, 0.12, 0.07, 0.07, 0.09, 0.12, 0.16, 0.13, 0.11,
    0.10, 0.08, 0.13, 0.06, 0.09, 0.13, 0.16, 0.12, 0.23, 0.35, 0.33, 0.28, 0.24, 0.33,
    0.25, 0.25, 0.24, 0.25, 0.28, 0.28, 0.34, 0.24, 0.33, 0.17, 0.25, 0.24, 0.35, 0.24,
    0.24, 0.22, 0.29, 0.23, 0.20, 0.32, 0.25, 0.25, 0.33, 0.21, 0.18, 0.22, 0.27, 0.18,
    0.25, 0.22, 0.29, 0.27, 0.33, 0.20, 0.31, 0.29, 0.17, 0.17, 0.29, 0.39, 0.65, 0.84,
    0.71, 0.64, 0.52, 0.91, 0.82, 0.36, 0.37, 0.95, 0.87, 0.73, 0.67, 0.73, 0.80, 0.91,
    0.63, 0.58, 0.60, 0.75, 0.53, 0.88, 0.84, 0.98, 1.20, 1.20, 1.02, 1.02, 1.17, 1.14,
    1.02, 1.13, 1.15, 1.25, 1.03, 1.04, 1.25, 1.12, 1.02, 1.26, 1.44, 1.33, 1.33, 1.49,
)

actual = (
    0.04, 0.03, 0.03, 0.02, 0.04, 0.01, 0.03, 0.02, 0.01, 0.01, 0.04, 0.01, 0.00, 0.05,
    0.03, 0.03, 0.05, 0.04, 0.02, 0.04, 0.02, 0.01, 0.05, 0.00, 0.01, 0.05, 0.01, 0.02,
    0.04, 0.00, 0.01, 0.01, 0.04, 0.04, 0.03, 0.01, 0.03, 0.04, 0.00, 0.03, 0.03, 0.05,
    0.05, 0.01, 0.05, 0.05, 0.03, 0.02, 0.02, 0.05, 0.04, 0.05, 0.04, 0.04, 0.01, 0.03,
    0.02, 0.01, 0.01, 0.00, 0.03, 0.02, 0.05, 0.03, 0.04, 0.13, 0.06, 0.07, 0.14, 0.11,
    0.10, 0.15, 0.14, 0.15, 0.07, 0.13, 0.08, 0.07, 0.07, 0.10, 0.15, 0.10, 0.11, 0.08,
    0.09, 0.06, 0.15, 0.12, 0.10, 0.12, 0.14, 0.16, 0.16, 0.11, 0.07, 0.06, 0.15, 0.1,
    0.15, 0.14, 0.14, 0.09, 0.13, 0.13, 0.15, 0.09, 0.11, 0.11, 0.13, 0.15, 0.14, 0.12,
    0.12, 0.06, 0.08, 0.13, 0.07, 0.16, 0.09, 0.10, 0.21, 0.17, 0.27, 0.24, 0.33, 0.24,
    0.28, 0.28, 0.19, 0.17, 0.29, 0.27, 0.22, 0.35, 0.19, 0.28, 0.30, 0.33, 0.29, 0.31,
    0.17, 0.27, 0.34, 0.26, 0.22, 0.30, 0.22, 0.22, 0.32, 0.34, 0.21, 0.21, 0.30, 0.19,
    0.27, 0.22, 0.19, 0.23, 0.26, 0.33, 0.23, 0.31, 0.18, 0.34, 0.35, 0.55, 0.76, 0.37,
    0.92, 0.86, 0.72, 0.78, 0.54, 0.70, 0.40, 0.45, 0.37, 1.00, 0.48, 0.92, 0.45, 0.57,
    0.55, 0.56, 0.75, 0.50, 0.41, 0.71, 0.82, 0.73, 1.04, 1.17, 1.17, 1.09, 1.06, 1.04,
    1.14, 1.18, 1.09, 1.03, 1.08, 1.16, 1.09, 1.12, 1.22, 1.32, 1.38, 1.39, 1.37, 1.37,
)

fig, ax = plt.subplots()

for label, rv in (('guess', guess), ('actual', actual)):
    x = np.sort(rv)
    ecpdf = np.linspace(0, 1, len(x), endpoint=False)
    ax.step(x, ecpdf, label=f'{label}, ecpdf')

    param = scipy.stats.gamma.fit(rv)
    x = np.linspace(0, 1.5, 500)
    cdf = scipy.stats.gamma.cdf(x, *param)
    ax.plot(x, cdf, label=f'{label}, gamma cdf')

ax.set_title('Rainfall, 11 Aug 2011')
ax.set_xlabel('Rainfall (furlongs per fortnight)')
ax.legend()
plt.show()

赞(0）回复(0）举报 2023-03-02

我来回答

scipy Gamma累积分布函数的正确拟合

1条答案

相关问题

热门标签

最新问答