pythonnumpyscipycdfgamma-distribution

Properly Fitting a Gamma Cumulative Distribution Function

I have two Numpy arrays (both 210 entries in total) of rainfall values, one observed and the other forecast. My goal is to create a best-fit gamma CDF (my first time diving into gamma CDFs) to both of these arrays and determine the relevant percentile that values then provided would fall into. The image below provides a simpler graphical reference of the gamma CDF I'm attempting to create with these two arrays. An important note is that the y-axis references the percentile of each value in the histogram, so ranging from 1st to 99th:

Gamma CDFs

These arrays are as follows:

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]

I've created a histogram for both of these arrays, binned in increments of 0.05 for a total of 30 bins. The code snippet to achieve this from the data supplied above is as follows:

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)])

I've also graphed both of these histograms, which looks as follows:

Dual Histogram Bar Chart

I am unsure where to proceed from here in order to create the best-fit gamma CDFs for both arrays, though I've initially found a stats.gamma function in scipy. Any help on how to proceed would be appreciated.

Solution

  • Use the built-ins made for this purpose in scipy. Also, a histogram is less illustrative than an ECPDF, which shows every data point and is more easily comparable to the fit 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()

    graphs