通过指定位置和尺度参数的协变相关指数函数，缩放降水数据的GEV分布

Question

背景：我正在进行一项气候归因研究，该研究应用了Philip等人发表的协议中提出的方法。(2020)。我正在使用的数据是来自气象站的大约100年的每日降水量数据。特别关注的是1998年发生的一次强降水事件。目标是找出这一事件发生的可能性和规模是否在多年内发生了变化"对统计模型的拟合结果提供了对感兴趣事件的回报期的估计，以及是否存在超出自然变异性预期偏差范围的趋势"根据降水数据计算区块最大值(每年)。该极值数据集随后被拟合到广义极值(GEV)分布，并用平滑全局平均温度(GMST)的指数函数对其进行缩放。因此，整个分布利用以下函数中表示的位置和比例参数进行缩放：µ=µ0 exp(αT’/µ0)和σ=σ0 exp"协变相关函数可以反转，并对给定年份的分布进行评估，例如过去的一年(T’=T0)或当前年份(T’=T1)。这给出了至少与这2年中观测到的事件一样极端的事件的概率，p0和p1，或表示为重现期τ0=1/p0和τ1=1/p1；

问题：我的问题在于指数标度函数的应用。我拟合了平稳GEV模型，并从中获得了位置、尺度和形状参数。我的理解是，非平稳拟合(以GMST值为协变量)的位置和尺度参数由分别求解T0和T1的函数组成。我还没有弄清楚如何在python中实现这一点来计算概率p0和p1。

如果有人能帮我解决这个问题，我将不胜感激！

数据集：

数据集1：重新采样以毫米为单位的日降水量值，以显示仅春季月份(3、4、5)的三天总和(毫米)不幸的是，我不被允许发布我的数据集，它的值范围从0-45mm

数据集2：以摄氏度为单位的年度GMST(温度)异常值年度数据集可在此处下载：https://www.ncei.noaa.gov/access/monitoring/global-temperature-anomalies/anomalies

我尝试了什么：

• 我在模块"1"的帮助下计算了数据集1的极值；pyextremes">

from pyextremes import EVA
# transform df into 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

• 我导入了软件包climtextremes(最初为R:climexRemes开发)，以使用其允许协变的fit_gev函数

import climextremes

# dataframe conversion for R usage
import pandas as pd
import rpy2.robjects as ro
from rpy2.robjects.packages import importr
from rpy2.robjects import pandas2ri
from rpy2.robjects.conversion import localconverter
with localconverter(ro.default_converter + pandas2ri.converter):
r_from_df_extremes = ro.conversion.py2rpy(df_extremes.threeday_prcp)
r_from_df_extremes
with localconverter(ro.default_converter + pandas2ri.converter):
r_from_GMST = ro.conversion.py2rpy(GMST)
r_from_GMST

# stationary fit
result_stat = climextremes.fit_gev(r_from_df_extremes, getParams=True, bootSE=False)
result_stat
# Output:
# mle = maximum likelihood estimation
# {'mle_names': array(['location', 'scale', 'shape'], dtype='<U8'),
#  'mle': array([17.57752926,  6.09155533,  0.14391791]),
#  'se_mle_names': array(['location', 'scale', 'shape'], dtype='<U8'),
#  'se_mle': array([0.71883326, 0.56333485, 0.09658412]),
#  'nllh': array([340.05145515]),
#  'info': {'convergence': array([0], dtype=int32),
#   'counts': 'Compute Error',
#   'message': None,
#   'failure': array([0], dtype=int32)}}

# Parameters of stationary fit:
loc = 17.57752926
scale = 6.09155533
shape = 0.14391791
alpha = 0.9 --> assumed to be 0.9 

• 根据上面给出的指数公式，我指定了以下位置和比例函数：

# Defining functions for non-stationary fit
locfun = loc**(alpha*r_from_GMST/loc)
scalefun = scale**(alpha*r_from_GMST/loc)
# Defining functions for non-stationary fit - year 1913
locfun_t0 = loc**(alpha*-0.35/loc)
scalefun_t0 = scale**(alpha*-0.35/loc)

对于1913年，我插入了GMST异常值，该值为-0.35

• 对于T0(1913)的非平稳拟合，我应用以下公式，得到RRuntimeError

# Non-stationary fit
result_nonstat = climextremes.fit_gev(r_from_df_extremes, r_from_GMST, locationFun=locfun_t0, scaleFun=scalefun_t0, getParams=True, bootSE=False)
result_nonstat

File "C:UsersXanaconda3envsclimRtestlibsite-packagesrpy2rinterface.py", line 810, in __call__
raise embedded.RRuntimeError(_rinterface._geterrmessage())
RRuntimeError: Error in parseParamInput(locationFun, names(x), .allowNoInt) : 
parseParamInput: expecting integer-valued indices in locationFun.

我想要实现的是：

"[…]使用与GMST成比例的GEV分布，并对1900年和2016年进行评估。根据许可协议，转载自van der Wiel等人(2017)https://creativecommons.org/licenses/by/4.0/(最后访问时间：2018年10月31日)；有关更多详细信息，请参阅原始文件">

我愿意就如何到达上面显示的图形的其他包裹或方式提出建议

提前感谢您的任何建议：)

Answer 1

我是气候极限的开发者。当别人向我提到这篇文章时，我才姗姗来迟地注意到这篇文章

不幸的是，你不能在极端气候中使用Philip等人的参数化。特别是参数化共享参数跨越地点和规模，这是不可能的在极端气候中，它提供了协变的线性函数分别针对每个位置、比例、形状。

你可以联系世界天气归因的人，看看他们是否有可以共享的代码。更好的是，也许他们已经在某个地方公开发布了。

Answer 2

在联系了climextremes的作者后，我发现，到目前为止(2023年4月)，不可能像你试图做的那样指定指数依赖关系。只能通过提供locationFun=1来指定线性位置依赖关系。可以通过指定locationFun=1, scaleFun=1来另外添加比例参数的对数依赖性。如果你致力于指定位置和尺度参数的指数依赖性，你可以尝试在R中结合rpy2和fevd函数。然而，如果我理解正确，你也可以使用线性依赖性，而不会对你的科学造成任何重大不利影响。

顺便说一句，我不太确定你是如何读取数据的，但你可以很容易地将numpy.array提供给climextremes.fit_gev函数的x和y，而无需将它们转移到R对象中。