R语言极值分析：GEV与GPD模型与MCMC的海洋观测数据极值模拟可视化研究

全文链接：https://tecdat.cn/?p=37007

原文出处：拓端数据部落公众号

在海洋科学领域，极端天气和海洋事件如极端海浪、风暴潮和海啸等，对沿海社区、基础设施及生态环境构成了重大威胁。准确预测和评估这些极端事件的强度和频率，对于制定有效的防灾减灾策略至关重要。极值分析作为统计学的一个重要分支，专门用于处理和分析极端值的出现规律，近年来在海洋观测数据的处理中得到了广泛应用。

GEV分布是一个包含三种基本极值分布（Gumbel、Fréchet和Weibull）的通用框架，能够适用于广泛的极值数据场景。而GPD分布则在阈值超量建模中表现出色，特别适用于超过某一特定阈值的极值数据分析。结合这两种模型，研究人员可以更全面地理解和模拟海洋观测数据中的极值行为。

此外，马尔可夫链蒙特卡洛（Markov Chain Monte Carlo, MCMC）方法作为一种强大的统计模拟技术，在复杂分布的参数估计和预测中发挥着重要作用。将MCMC与GEV或GPD模型相结合，不仅可以提高参数估计的精度，还能够生成符合特定分布规律的模拟数据，为海洋极值事件的情景分析和风险评估提供有力支持。

本文旨在帮助客户探讨R语言环境下，实现GEV与GPD模型结合MCMC方法在海洋观测数据极值模拟可视化研究中的应用课题研究。通过具体案例，展示如何构建这些模型，进行参数估计，并利用MCMC方法进行模拟数据生成和可视化分析。本文的研究不仅有助于深入理解海洋极值事件的统计规律，还为海洋灾害的预警和应对提供了科学依据和技术支持。

数据清理与统计分析

在本文中，我们首先对一组海洋观测数据进行了预处理，随后进行了基本的统计分析和可视化展示。这些数据包含了不同时间点的海浪高度（Hs）、最大海浪高度（Hmax）、水温（SST）以及其他海洋环境参数，如潮汐（Tz）、波峰周期（Tp）和波峰方向（Dir_Tp.TRUE）。数据的时间跨度为1996年11月21日的多个时间点，每行数据代表一个观测点。

在这里插入图片描述

数据预处理

为了确保数据分析的准确性和可靠性，我们首先对数据进行了清理，移除了任何可能存在的缺失值。在R语言中，使用na.omit()函数可以轻松实现这一步骤，该函数能够自动删除包含NA值的行。

 
  data <- na.omit(data)  

  
  summary_stats <- c("Mean", "Median", "Sd")
  data_summary <- sapply(data[, c("Hs")], function(x) summary(x)[summary_stats])
  data_summary

可视化分析

为了进一步直观展示海浪高度的分布情况，我们采用了箱线图和直方图两种可视化方法。

箱线图

箱线图是一种用于展示数据分布情况的图表，它能够清晰地显示数据的四分位数、中位数以及异常值。在本例中，我们绘制了Hs的箱线图，并通过设置outlier.shape = NA来隐藏异常值的标记，使图表更加简洁。在这里插入图片描述

直方图

直方图则是另一种常用的数据分布展示方式，它通过条形图的形式展示了数据的频数分布。在本例中，我们设置了30个分组来绘制Hs的直方图，并通过设置颜色和填充样式来增强图表的可读性。

在这里插入图片描述

通过上述分析和可视化展示，我们可以对海浪高度（Hs）的分布情况有一个较为全面的了解，为后续的研究和预测提供了有力的数据支持。

GPD-MCMC模型在极端海浪高度分析中的应用

在本文中，我们利用广义帕累托分布（Generalized Pareto Distribution, GPD）模型对极端海浪高度数据（Hs）进行了详细的分析与拟合。GPD模型因其能够灵活捕捉极端值分布的特性，在极端事件统计分析中得到了广泛应用。

模型拟合

首先，我们使用fevd函数从extRemes或类似R包中拟合GPD模型至海浪高度数据（Hs）。在拟合过程中，我们指定了类型为"GP"（即广义帕累托分布），并设置了阈值（threshold）为6，以筛选出极端海浪高度数据进行专门分析。

模型诊断与可视化

为了评估GPD模型的拟合效果，我们进行了多个诊断步骤和可视化展示。通过distill函数，我们提取了模型的关键统计信息，以便进一步分析。

在这里插入图片描述

接着，我们使用mrlplot函数绘制了超出量均值剩余寿命（Mean Residual Life, MRL）图，以检查数据是否满足GPD模型的假设。通过设定xlim参数，我们专注于海浪高度在0到12米范围内的分析。

在这里插入图片描述

置信区间与极值预测

为了评估模型参数的不确定性，我们计算了参数的置信区间（CI）。这有助于我们理解模型估计的稳健性。
在这里插入图片描述
进一步地，我们利用ci函数预测了不同重现期（return period）下的极端海浪高度及其置信区间。这对于海洋工程设计和风险评估尤为重要。

阈值敏感性分析

为了探究阈值选择对模型结果的影响，我们进行了阈值范围分析。通过threshrange.plot函数，我们绘制了不同阈值下的模型拟合效果，并特别关注了阈值在-75.5到-68.5范围内（注意：这里假设的阈值范围可能基于特定数据集或分析需求，实际应用时需根据数据特性调整）的变化情况，采用“PP”（Pickands Plot）类型进行展示。

在这里插入图片描述

MCMC 采样

我们开发了MCMC采样函数，通过提议分布（如正态分布）生成新的参数候选值，并根据Metropolis-Hastings算法接受或拒绝这些候选值。对数似然函数用于计算候选值和当前值的接受概率。

我们使用R语言实现了GPD的对数似然函数、提议分布和MCMC采样函数。具体实现细节已在前文给出。

# 定义 GPD 的概率密度函数
  gpd_pdf <- function(x, shape, scale) {
    if (shape == 0) {
      return(exp(-x / scale) / scale)
    } else {
      return((1 + shape * x / scale)^(-1 - 1 / shape) * shape / scale)
    }
  }
  
  # 定义提议分布（这里使用正态分布）
  proposal_dist <- function(theta, var) {
    return(rnorm(length(theta), theta, var))
  }
  
  # MCMC 采样函数
  mcmc_sample <- function(log_likelihood, initial_theta, data, n_iter, var) {
    theta <- initial_theta