基于GMST协变量的GEV分布降水数据缩放实现求助

背景 #

开展气候归因研究，采用Philip et al. (2020)提出的方法处理某气象站近百年逐日降水数据，重点分析1998年极端强降水事件，目标是探究该事件发生概率和强度随年份的变化。通过块最大值法（逐年）提取降水极值数据集，拟拟合到广义极值分布（GEV），并采用平滑后的全球平均温度距平（GMST）的指数函数对分布进行缩放，位置和尺度参数公式为：
µ = µ₀ exp(αT'/µ₀)
σ = σ₀ exp(αT'/µ₀)
需要计算过去年份（T'=T₀）和当前年份（T'=T₁）中极端事件的发生概率p₀、p₁及重现期。

数据集 #

• 数据集1：逐日降水数据，已重采样为春季（3-5月）的三日累计降水量，范围0-45mm
• 数据集2：年度GMST距平数据（可从指定渠道获取）

已尝试操作 #

1. 用pyextremes模块提取降水极值：

from pyextremes import EVA
# 将DataFrame转换为Series
threeday_evaseries = pd.Series(threeday_springseries['threeday_prcp'])
model = EVA(threeday_evaseries)
model.get_extremes(method="BM", errors="coerce")
print(model.extremes.head())
extremes = model.extremes

2. 用climextremes包拟合平稳GEV模型，得到参数：loc=17.57752926，scale=6.09155533，shape=0.14391791，假设α=0.9
3. 定义非平稳拟合的位置和尺度函数，针对1913年（GMST距平-0.35）设置后调用fit_gev时出现报错：

RRuntimeError: Error in parseParamInput(locationFun, names(x), .allowNoInt) : 
  parseParamInput: expecting integer-valued indices in locationFun.

需求 #

实现基于GMST协变量的非平稳GEV拟合，计算不同年份的极端事件概率，并绘制类似van der Wiel et al. (2017)展示的随GMST变化的GEV分布图表，接受其他Python包的方案建议。

方案1：纯Python自定义拟合（用scipy.stats） #

避开climextremes的R底层参数解析问题，用纯Python实现更灵活：

1. 对齐数据：确保逐年提取的极值序列与对应年份的GMST距平序列长度一致
2. 定义非平稳GEV对数似然函数：

import numpy as np
from scipy.stats import genextreme
from scipy.optimize import minimize

def nonstationary_gev_loglike(params, data, gmst_anomaly):
    mu0, sigma0, shape, alpha = params
    # 计算随GMST变化的位置和尺度参数
    mu = mu0 * np.exp(alpha * gmst_anomaly / mu0)
    sigma = sigma0 * np.exp(alpha * gmst_anomaly / mu0)
    # 返回负对数似然（供scipy极小化）
    return -np.sum(genextreme.logpdf(data, c=shape, loc=mu, scale=sigma))

3. 拟合模型：用平稳GEV参数作为初始值

# 初始参数：平稳GEV的loc、scale、shape，加上假设的alpha
initial_params = [17.5775, 6.0916, 0.1439, 0.9]
# 传入极值数据和对应GMST距平，设置参数边界避免不合理值
result = minimize(nonstationary_gev_loglike, initial_params, 
                  args=(extremes.values, gmst_anomaly_values),
                  bounds=[(10, 25), (3, 10), (-0.5, 0.5), (0.5, 1.5)])

# 提取拟合后的参数
mu0_fit, sigma0_fit, shape_fit, alpha_fit = result.x

4. 计算极端事件概率与重现期：

def calculate_probability(threshold, gmst, mu0, sigma0, shape, alpha):
    mu = mu0 * np.exp(alpha * gmst / mu0)
    sigma = sigma0 * np.exp(alpha * gmst / mu0)
    # GEV生存函数（1-CDF）即超过阈值的概率
    return genextreme.sf(threshold, c=shape, loc=mu, scale=sigma)

# 示例：计算30mm阈值在过去（gmst=-0.35）和当前（gmst=0.6）的概率
p0 = calculate_probability(30, -0.35, mu0_fit, sigma0_fit, shape_fit, alpha_fit)
p1 = calculate_probability(30, 0.6, mu0_fit, sigma0_fit, shape_fit, alpha_fit)

# 重现期（单位：年）
return_period0 = 1 / p0
return_period1 = 1 / p1

方案2：用pyextremes内置协变量支持 #

pyextremes原生支持非平稳模型，无需手动写似然函数：

from pyextremes import EVA
import pandas as pd

# 确保极值序列带年度时间索引，GMST数据与极值时间对齐
model = EVA(extremes)
# 添加GMST距平作为协变量
model.add_covariate("gmst", gmst_data)
# 定义非平稳参数公式，用Patsy语法结合自定义变换
model.fit_model(
    model="GEV",
    distribution_kwargs={
        "loc": "mu0 * np.exp(alpha * gmst / mu0)",
        "scale": "sigma0 * np.exp(alpha * gmst / mu0)",
        "shape": "shape",
    },
    initial_params={
        "mu0": 17.5775,
        "sigma0": 6.0916,
        "shape": 0.1439,
        "alpha": 0.9,
    }
)

# 查看拟合结果
print(model.summary)

绘制GMST相关的GEV分布图表 #

用matplotlib实现类似van der Wiel et al. (2017)的可视化：

import matplotlib.pyplot as plt
import seaborn as sns

# 生成GMST距平范围和降水阈值范围
gmst_range = np.linspace(-0.5, 1.0, 100)
prcp_range = np.linspace(0, 45, 200)

plt.figure(figsize=(10,6))
# 每隔10个GMST值绘制一条分布曲线
for gmst in gmst_range[::10]:
    mu = mu0_fit * np.exp(alpha_fit * gmst / mu0_fit)
    sigma = sigma0_fit * np.exp(alpha_fit * gmst / mu0_fit)
    pdf = genextreme.pdf(prcp_range, c=shape_fit, loc=mu, scale=sigma)
    plt.plot(prcp_range, pdf, label=f"GMST距平 = {gmst:.2f}°C")

plt.xlabel("春季三日累计降水量 (mm)")
plt.ylabel("概率密度")
plt.title("GEV分布随GMST距平的变化")
plt.legend(bbox_to_anchor=(1.05, 1), loc='upper left')
plt.grid(alpha=0.3)
plt.show()

修复climextremes的报错（若坚持使用） #

报错原因是climextremes要求位置/尺度函数用协变量的整数索引，而非直接代入数值。需将GMST距平作为协变量列传入：

from climextremes import fit_gev

# 假设extremes是包含'prcp'和'gmst_anomaly'列的DataFrame
locationFun = "mu0 * exp(alpha * x[,2]/mu0)"
scaleFun = "sigma0 * exp(alpha * x[,2]/mu0)"
shapeFun = "shape"

fit_result = fit_gev(
    x=extremes[['prcp', 'gmst_anomaly']].values,
    locationFun=locationFun,
    scaleFun=scaleFun,
    shapeFun=shapeFun,
    startParams=[17.5775, 6.0916, 0.1439, 0.9]
)

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