求助：如何绘制双变量的Joint Return Period图（附R语言示例代码）

我来帮你搞定这个联合重现期图的绘制！要实现你想要的效果，我们可以结合极值统计方法和等值线可视化，下面我会给出完整的R代码，再一步步拆解关键步骤，保证你能轻松复现。

先看完整的可运行代码：

# 加载需要的工具包
library(extRemes)
library(ggplot2)
library(dplyr)

# 生成你提供的示例数据
set.seed(101)
dates = seq(as.Date("1951-12-31"), as.Date("2000-12-31"), by="day")
var1 = rnorm(length(dates), 0.5)
var2 = rnorm(length(dates),5)

# 1. 拟合两个变量的边缘极值分布（这里选Gumbel，可根据数据调整）
fit_var1 = fevd(var1, type = "Gumbel")
fit_var2 = fevd(var2, type = "Gumbel")

# 2. 计算不同重现期对应的等值线数据
# 定义要展示的重现期（单位：年）
target_rps = c(1,5,10,50,100)
# 总年数（用于计算重现期）
total_years = length(dates)/365

# 创建覆盖变量范围的网格点
var1_grid = seq(min(var1), max(var1), length.out = 100)
var2_grid = seq(min(var2), max(var2), length.out = 100)
grid_df = expand.grid(var1 = var1_grid, var2 = var2_grid)

# 计算每个网格点的联合超越概率（这里先假设变量独立，后面会说如何处理相关情况）
grid_df$exceed_prob = (1 - pgev(grid_df$var1, loc = fit_var1$results$par[1], scale = fit_var1$results$par[2])) *
                      (1 - pgev(grid_df$var2, loc = fit_var2$results$par[1], scale = fit_var2$results$par[2]))

# 计算联合重现期：T = 总年数 / 联合超越概率
grid_df$return_period = total_years / grid_df$exceed_prob

# 3. 绘制最终图形
ggplot() +
  # 先画原始数据散点
  geom_point(aes(x = var1, y = var2), size = 0.5, alpha = 0.3) +
  # 画重现期等值线（蓝色）
  geom_contour(aes(x = var1, y = var2, z = return_period), 
               breaks = target_rps, color = "blue", linewidth = 1) +
  # 给等值线添加重现期标签
  geom_contour_label(aes(x = var1, y = var2, z = return_period), 
                     breaks = target_rps, color = "blue", fontface = "bold") +
  # 调整标签和主题
  labs(x = "变量1 (var1)", y = "变量2 (var2)", title = "双变量联合重现期图") +
  theme_minimal()

关键步骤详解 #

1. 拟合边缘极值分布 #

我们用extRemes包的fevd()函数给每个变量拟合极值分布（这里选Gumbel分布，适合多数极值场景；如果你的数据有明显下界/上界，可以换成Weibull或广义极值分布）。这一步是为了准确计算单个变量的超越概率（即某个值被超过的概率）。

2. 计算联合重现期 #

联合重现期的核心逻辑是：

联合重现期 ( T = \frac{总年数}{P(X > x, Y > y)} )
其中 ( P(X > x, Y > y) ) 是var1超过x 且 var2超过y的联合概率。

上面的代码用了变量独立假设，如果你的两个变量存在相关性（比如降水和径流），建议用copula包建模依赖关系，这样结果更准确：

# 示例：用高斯Copula处理变量相关性
library(copula)
# 先把原始数据转换为均匀分布样本（基于经验CDF）
u = ecdf(var1)(var1)
v = ecdf(var2)(var2)
# 拟合高斯Copula
cop_fit = fitCopula(gaussianCopula(), data = cbind(u, v))
# 计算网格点的联合超越概率
grid_df$u = ecdf(var1)(grid_df$var1)
grid_df$v = ecdf(var2)(grid_df$var2)
grid_df$exceed_prob = 1 - pCopula(cbind(grid_df$u, grid_df$v), cop_fit@copula)
# 重新计算重现期
grid_df$return_period = total_years / grid_df$exceed_prob

3. 可视化优化 #

如果你习惯用基础绘图系统，也可以这样实现：

# 基础绘图版本
plot(var1, var2, cex=0.1, xlab="变量1", ylab="变量2", main="双变量联合重现期图")
# 添加蓝色等值线
contour(var1_grid, var2_grid, matrix(grid_df$return_period, nrow=100), 
        levels=target_rps, col="blue", add=TRUE)
# 给等值线添加标签
clabel(contour(var1_grid, var2_grid, matrix(grid_df$return_period, nrow=100), 
               levels=target_rps, drawlabels=FALSE), 
       col="blue", font=2)

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