如何为降雨数据拟合最优Gamma累积分布函数（CDF）

核心实现步骤 #

Scipy的scipy.stats.gamma内置了拟合、CDF（累积分布函数）和PPF（百分位数函数）工具，直接用这些方法就能完成需求，具体步骤如下：

1. 导入依赖库

import numpy as np
from scipy.stats import gamma

2. 拟合Gamma分布参数

Gamma分布的核心参数为形状参数（shape/a）、位置参数（loc）和尺度参数（scale）。考虑到降雨数据的非负特性，建议固定loc=0（符合物理意义，同时减少拟合参数提升稳定性），用最大似然估计（MLE）拟合两组数据：

# 拟合预报值数组guess的Gamma分布
guess_shape, guess_loc, guess_scale = gamma.fit(guess, floc=0)

# 拟合观测值数组actual的Gamma分布
actual_shape, actual_loc, actual_scale = gamma.fit(actual, floc=0)

注：如果数据中存在0值，Gamma分布的定义域是x>0，可以给所有数据加极小值（如1e-8）避免拟合报错：guess = guess + 1e-8

3. 计算给定降雨值对应的百分位数

用gamma.cdf()计算目标值在拟合分布中的累积概率，乘以100就是对应的百分位数：

# 示例：计算降雨值5.0在两组分布中的百分位数
target_rainfall = 5.0

# 预报值分布的百分位数
guess_percentile = gamma.cdf(target_rainfall, a=guess_shape, loc=guess_loc, scale=guess_scale) * 100
# 观测值分布的百分位数
actual_percentile = gamma.cdf(target_rainfall, a=actual_shape, loc=actual_loc, scale=actual_scale) * 100

print(f"预报值{target_rainfall}对应的百分位数: {guess_percentile:.2f}%")
print(f"观测值{target_rainfall}对应的百分位数: {actual_percentile:.2f}%")

4. 反向计算：给定百分位数求对应降雨值

如果需要根据百分位数（如90%分位数）反推降雨值，用gamma.ppf()（分位数函数）：

# 示例：计算90百分位数对应的降雨值
target_percentile = 90

guess_90th = gamma.ppf(target_percentile / 100, a=guess_shape, loc=guess_loc, scale=guess_scale)
actual_90th = gamma.ppf(target_percentile / 100, a=actual_shape, loc=actual_loc, scale=actual_scale)

print(f"预报分布的90百分位数: {guess_90th:.2f}")
print(f"观测分布的90百分位数: {actual_90th:.2f}")

5. 验证拟合效果（可选）

可以绘制经验CDF和拟合CDF对比，确认拟合质量：

import matplotlib.pyplot as plt

# 生成覆盖数据范围的x轴序列
x_vals = np.linspace(min(min(guess), min(actual)), max(max(guess), max(actual)), 1000)

# 计算拟合的CDF曲线
guess_fitted_cdf = gamma.cdf(x_vals, a=guess_shape, loc=guess_loc, scale=guess_scale)
actual_fitted_cdf = gamma.cdf(x_vals, a=actual_shape, loc=actual_loc, scale=actual_scale)

# 绘制经验CDF与拟合CDF
plt.figure(figsize=(10, 6))
# 经验CDF（排序后的数据+累积概率）
plt.plot(np.sort(guess), np.linspace(0, 1, len(guess)), label='Guess Empirical CDF')
plt.plot(np.sort(actual), np.linspace(0, 1, len(actual)), label='Actual Empirical CDF')
# 拟合CDF（虚线）
plt.plot(x_vals, guess_fitted_cdf, label='Guess Fitted Gamma CDF', linestyle='--')
plt.plot(x_vals, actual_fitted_cdf, label='Actual Fitted Gamma CDF', linestyle='--')

plt.xlabel('Rainfall')
plt.ylabel('Cumulative Probability')
plt.legend()
plt.title('Empirical vs Fitted Gamma CDF')
plt.show()

内容的提问来源于stack exchange，提问作者TornadoEric