原文链接：/?p=15301
结合POT模型的洪水风险评估能够从有限的实测资料中获取更多的洪水风险信息,得到更贴近事实的风险评估结果,能为决策者提供更多的依据,从而使决策结果更加可靠实用。
对于这些同样面临挑战的人，我希望这个博客将有助于简化工作。
案例POT序列在47年的记录期内提供了高于74 m3/ s阈值的47个峰值。
我们的目标是将概率模型拟合到这些数据并估算洪水分位数。
我从获取了每次洪水的日期，并将其包含在文件中。有趣的是，最早的洪水流量是1943年，而最后一次是1985年，是43年的记录，而不是47年。这是因为1939年至1943年的洪水都小于74 m3/ s的阈值。
首先计算这些数据点的绘制位置。
T给定排放超标之间的平均间隔（年）R是POT系列中的流量等级（最大流量是等级1）n是数据的年数。
请注意，这是记录的年数，而不是峰值数。
同样，重要的是要认识到，方程式1对POT系列的作用与年度系列不同。让我们看一个显示这种差异的示例。考虑以下情况：我们根据47年的数据分析了POT系列的94个峰。在这种情况下，最小的峰的等级为94。重复间隔为：
这大约是半年或6个月，这似乎是合理的（47年中有94个高峰，因此平均每年有2个高峰，平均相隔约6个月）。
将绘图位置解释为年度超出概率将得出以下结果：
也就是说，概率大于1，这没有意义。因此，我们不能使用绘图位置公式来计算阈值峰值序列中的数据的AEP。取而代之的是，方程式1的逆可以解释为EY，即每年的预期超出次数。
ARR示例将指数分布拟合为概率模型。
为了计算L2，我们使用QJ Wang（Wang，1996）的公式
L2 <- function(q){q <- sort(q)n <- length(q)0.5*(1/choose(n,2))*sum((0:(n-1) - (n-1):0)*q)}
qi从最小到最大的顺序是流量（POT）n是流的数量
L2 = 79.12
指数分布的参数可以用L矩表示。我们使用的是广义帕累托（GP）公式。
对于指数分布：
这些参数估计值的置信区间可以使用bootstrapping计算得出。
Beta的95％置信区间是（37.4，89.4）和（120.6，244.7）。参数之间的相关性约为-0.5。参数的不确定性如图1所示。
param_errors_df %>% ggplot(aes(x = V1, y = V2)) +geom_point(size = 0.1) +scale_x_continuous(name = 'beta') +scale_y_continuous(name = bquote('q'['*'])) +stat_ellipse(colour = 'red') + # 95% 置信区间theme_gray(base_size = 7)
图1：参数的不确定性。椭圆显示置信限度为95％
指数分布将超出概率与流的大小相关。在这种情况下，在任何POT事件中，峰值流量超过某个值的概率为：
这是针对超额概率的。在水文学中，我们通常使用超出概率（洪水大于特定值的概率），因此所需方程式为一个减去所示方程式。
通过将每年超过阈值的洪峰平均数乘以POT概率，我们可以将POT概率转换为每年的预期超标次数。
74 m3/ s阈值，POT系列中有47个值，并且有47年的数据，因此每年的平均峰值数为1。因此，EY可以表示为：
其中，q是每年POT洪水的平均数量，
w也可以与EY以下内容相关：
我们还可以EY通过以下公式与年度超出概率相关：
因此，通过以下等式，洪水幅度可以与AEP相关，而AEP可以与洪水幅度相关。
这些方程式可用于估计标准EY值的分位数。使用bootstrap自举法估计了置信区间（95％）（表1）。
res
现在绘制数据和拟合度（图2）。x值是根据等式1的逆计算的EY；y值是流量。拟合基于等式6。使用bootstrap自举法计算分位数的置信区间。
图2：河流的部分序列显示契合度和置信区间
我个人更希望该图向右增加，这通常是洪水频率曲线的绘制方式。这仅涉及使用ARI作为纵坐标（图3）。
图3：河流部分序列显示契合度和置信区间
参考文献
1.R语言基于ARMA-GARCH-VaR模型拟合和预测实证研究
2.R语言时变参数VAR随机模型
3.R语言时变参数VAR随机模型
4.R语言基于ARMA-GARCH过程的VAR拟合和预测
5.GARCH（1,1），MA以及历史模拟法的VaR比较
6.R语言时变参数VAR随机模型
7.R语言实现向量自动回归VAR模型
8.R语言随机搜索变量选择SSVS估计贝叶斯向量自回归（BVAR）模型
9.R语言VAR模型的不同类型的脉冲响应分析
