zehua

Zehua

广义线性模型实现降水统计降尺度

Downscaling
15
0
2025-05-13

简介

目标

侧重于PP方法下广义线性模型(GLM)的降水统计降尺度，重点要实现区域PCA以及交叉验证

方法

模型: 降水被建模为伯努利-伽马（Bernoulli-Gamma）分布，包括三个参数：p, α (shape), 和 β (rate)。

  • 第一阶段：二项 GLM（逻辑连接函数）拟合降水发生（概率 p ，即日降水量 >1mm 的概率）。

  • 第二阶段：伽马 GLM（对数连接函数）拟合降水强度（湿日降水量，即 \alpha \beta ）。并且强度模型假设 rate 参数 \beta 在给定预测因子条件下是恒定的，并通过残差估计。

Predictor:

  • 训练数据: ERA-Interim 再分析数据（1979-2008）的 5 个变量在 4 个气压层上的数据: 1000、850、700、500 hPa 四个气压层上的温度、位势高度、北向风、东向风和比湿。

  • 预测数据:

    • 3.1 使用交叉验证,

    • 3.2 使用EC-Earth GCM的历史时期数据(1979–2008),

    • 3.3 使用未来RCP8.5情景数据(2071–2100)。

数据预处理(仅针对交叉验证):

  • 对每个目标站点，预测因子来源于其所属 PRUDENCE 区域内解释了 95% 方差的主成分（PCs）

  • GLM需要标准化预测变量(均值为0，标准差为1)。

Predictand:

ECA&D的83个站点的日降水量(原始86个，由于缺失数据清理后为83个)。

数据结构

st86_precip.rda


组件

类型

描述

$Data

数值矩阵 [10958 × 86]

表示 1979–2008 年间每日降水数据，共 86 个站点、10958 天

$Dates$start

字符向量 [10958]

与降水数据时间轴一一对应

$xyCoords

data.frame [86 × 2]

每个站点的经纬度（x 为经度，y 为纬度）

$Metadata$station_id 等

字符/数值向量 [86]

站点 ID、名称、海拔等元信息

era_interim.rda

组件

类型

描述

$Data

数组 [20, 1, 10958, 20, 30]

[变量数, 成员数, 时间步, 纬度, 经度]

20 个变量（例如 z@500、q@850 等），10958 天，20×30 个格点

$Dates

列表 [20]

每个变量对应一个日期列表，每个列表含 start/end，时间格式为 GMT

$xyCoords$x

向量 [30]

经度 [-14, ..., 44]，步长 2°

$xyCoords$y

向量 [20]

纬度 [36, ..., 74]，步长 2°

实现

核心实现1：区域PCA -- prepareData

对每个目标站点，不使用全局预测因子，而是使用该站点所属的 PRUDENCE 区域内的预测因子数据进行 PCA（主成分分析），保留解释 95% 方差的主成分作为该站点的 GLM 输入，得到该区域最具有代表性的主成分PCs。

这里细节比较多，下面详细解释一下。

ERA-Interim是一个全球气候再分析数据集，对应于我们的 x，即 predictors 数据，在我们的情况下，它是一个具有2°×2°空间分辨率的网格，覆盖了整个欧洲地区不同压力层(1000, 850, 700, 500 hPa)上的气候要素信息，包括温度、位势高度、北向风、东向风和比湿，共计20个变量(Predictors)。

我们的目标是建立这些再分析变量 x （低分辨率）与观测降水 y （高分辨率）之间的统计关系，基于83个地面观测站点逐站建模。由于ERA-Interim数据的空间分辨率较粗(2°×2°)且覆盖面积大(整个欧洲)，如果直接对整个欧洲的预测因子进行降维（全局PCA），存在以下明显问题：

  • 忽略局部差异: 全局PCA提取的是整个欧洲最显著的主成分。例如，北欧天气可能对整个欧洲气候有很大影响，权重很大，但对于我们感兴趣的巴黎站点，我们只关心巴黎周围的气候信息。北欧的气候特征可能和巴黎没关系，甚至是干扰。

  • 计算成本极高: 使用整个欧洲区域的所有20个变量的网格点进行降维会导致特征空间非常大，PCA的计算成本非常高。

为了解决全局PCA存在的问题，人们引入了 PRUDENCE 区域的概念（例如，法国对应 FR 区域）。针对每一个站点，只选取其所属 PRUDENCE 区域内的 ERA-Interim 网格数据（共20个变量）进行PCA，从而生成该特定区域的主成分。

首先确定站点所属的 PRUDENCE 区域（例如 巴黎 站属于 FR 区域），提取该 PRUDENCE 区域内 ERA-Interim 的 20 个 predictors 数据，将这些变量按时间展开，这个展开形式取决于参数 which.combine ，一般情况下用TRUE，那么就是对特定 PRUDENCE 区域内的 20 个 predictors 进行组合 PCA，比如说，对于第一个 predictor ，把其所有格点数据展开成一个长向量，其长度是 N1(取决于ERA-Interim 数据的空间分辨率，分辨率越高， prudence 区域内对应 predictors 格点就越多，最后得到的向量长度就大)，同理对于第二个 predictor 也可以把其所有格点数据展开成一个长向量N2，然后一直到N20，最后把这20个长向量 concatenate 起来。因此输入到PCA的矩阵维度就是，行数是时间步长 （10958天），列数是N1+N2+...+N20（变量1的区域内格点数 + 变量2的区域内格点数 + ... + 变量20的区域内格点数）

因为我们的 ERA-Interim 数据是规则网格，即所有20个 Predictors 都在相同的空间网格上，在区域内有相同的网格点。我们可以通过代码验证ERA-Interim数据是在规则网格上。

# 加载数据
load("era_interim.rda")
​
# 查看其空间坐标
coords <- getCoordinates(era_interim)
lon <- coords$x
lat <- coords$y
​
# 判断是否是规则网格：是否间距一致。
is_regular_lon <- all(diff(lon) == diff(lon)[1])
is_regular_lat <- all(diff(lat) == diff(lat)[1])
​
if (is_regular_lon && is_regular_lat) {
  cat("ERA-Interim是规则网格\n")
  cat("经度范围：", range(lon), "，步长为：", diff(lon)[1], "\n")
  cat("纬度范围：", range(lat), "，步长为：", diff(lat)[1], "\n")
} else {
  cat("ERA-Interim不是规则网格\n")
}

对应结果:

> ERA-Interim是规则网格
> 经度范围： -14 44 ，步长为： 2
> 纬度范围： 36 74 ，步长为： 2

代码实现1：prepareData

原始的 prepareData 函数默认使用全局 PCA，通过参数spatial.predictors ，并调用transformeR::prinComp 对整个 ERA-Interim 数据（空间范围覆盖整个欧洲，经度 -14 到 44，纬度 36 到 74）进行 PCA，对应的$xyCoords$x长度为 30、$xyCoords$y 长度为 20，共计 600 个空间格点。每个格点包含 20 个预测因子变量（包括温度、位势高度、风速、比湿等），因此在进行全局 PCA 时，输入矩阵的特征维度为 600 × 20 = 12,000，计算负担比较大

因此为了实现区域PCA，我们对prepareData进行了以下修改：

1. PRUDENCE区域定义

setwd("/Users/zehua/Downloads/data_precipitation")
​
# 加载原始数据
load("st86_precip.rda")
​
# 获取站点坐标和ID
station_coords <- as.data.frame(st86_precip$xyCoords)
colnames(station_coords) <- c("lon", "lat")
station_coords$station_id <- st86_precip$Metadata$station_id
​
# 定义PRUDENCE区域
prudence_regions <- list(
  BI = list(lon = c(-10, 2), lat = c(50, 59)),
  IP = list(lon = c(-10, 3), lat = c(36, 44)),
  FR = list(lon = c(-5, 5), lat = c(44, 50)),
  ME = list(lon = c(2, 16), lat = c(48, 55)),
  SC = list(lon = c(5, 30), lat = c(55, 70)),
  AL = list(lon = c(5, 16), lat = c(44, 48)),
  MD = list(lon = c(3, 25), lat = c(36, 44)),
  EA = list(lon = c(16, 30), lat = c(44, 55))
)

2. 将每个站点分配到其对应区域

assign_station_to_region <- function(lon, lat, regions) {
  for (region in names(regions)) {
    bbox <- regions[[region]]
    if (lon >= bbox$lon[1] && lon <= bbox$lon[2] &&
        lat >= bbox$lat[1] && lat <= bbox$lat[2]) {
      return(region)
    }
  }
  return(NA)
}
station_coords$region <- mapply(assign_station_to_region,
                                 station_coords$lon,
                                 station_coords$lat,
                                 MoreArgs = list(regions = prudence_regions))

4. 查找每个当前站点对应的PRUDENCE区域

辅助函数map_station_to_prudence_region()，该函数将每个站点的经纬度坐标与其对应的区域匹配。如果站点不完全在任何区域边界内，则使用最近邻方法选择最近的PRUDENCE区域。

# ==================================================================================
# 辅助函数：区域查找
# ==================================================================================
​
map_station_to_prudence_region <- function(lon, lat, regions_definition) {
  region_name_inside <- NA
  for (region_code in names(regions_definition)) {
    bounds <- regions_definition[[region_code]]
    if (lon >= bounds$lonLim[1] && lon < bounds$lonLim[2] &&
        lat >= bounds$latLim[1] && lat < bounds$latLim[2]) {
      region_name_inside <- region_code
      return(region_name_inside)
    }
  }
  min_dist <- Inf; closest_region <- NA; point_coord <- c(lon, lat)
  for (region_code in names(regions_definition)) {
    bounds <- regions_definition[[region_code]]
    center_lon <- mean(bounds$lonLim); center_lat <- mean(bounds$latLim)
    region_center_coord <- c(center_lon, center_lat)
    current_dist <- geosphere::distHaversine(point_coord, region_center_coord)
    if (current_dist < min_dist) {
      min_dist <- current_dist
      closest_region <- region_code
    }
  }
  return(closest_region)
}

后续调用关键行代码：

station_coords <- getCoordinates(y)
station_regions <- sapply(1:nrow(station_coords), function(i) {
  map_station_to_prudence_region(station_coords$x[i], station_coords$y[i], prudence_regions)
})

5. 按PRUDENCE区域提取ERA子网格并执行PCA

对spatial.predictors$grid的使用进行了修改。不对整个x网格使用 prinComp进行PCA，而是先遍历每个PRUDENCE区域，然后使用subsetGrid()函数从全欧洲 Predictors中 提取只属于该 PRUDENCE 区域的 Predictors，然后对该区域进行独立的PCA。

示例：

# --- 提取区域子集 ---
x_regional <- subsetGrid(x, lonLim = region_bounds$lonLim, latLim = region_bounds$latLim)
​
# 准备此区域的参数
current_spatial_predictors <- spatial.predictors 
current_spatial_predictors[["grid"]] <- x_regional # 使用区域网格
​
message("... 为区域调用prinComp:", region_code)
full.pca <- do.call("prinComp", current_spatial_predictors)

6. 主成分提取和保存

对于每个区域的full.pca输出，我们提取主成分矩阵(如full.pca$COMBINED[[1]]$PCs)并构建两个主要结构：

  • regional_pca_models[[region_code]]：存储PCA模型对象(用于新数据投影)

  • regional_pcs_matrices[[region_code]]：存储主成分得分矩阵 EOFs(用于GLM建模输入)

最终以下列形式保存：

predictor.list <- list("y" = y, "x.global" = x, "x.local" = nn, "pca" = full.pca)

主函数中的调用示例：

############################################################################
## 4) 使用PCA维度降低(prepareData)
############################################################################
data.pca <- prepareData(
  x = era_interim,  # ERA-Interim再分析数据用作大尺度预测因子。
  y = st86_precip,  # 观测站点数据用作本地小尺度预测量。
  spatial.predictors = list(
    v.exp = 0.95,                      # 保留解释方差≥95%的主成分
    which.combine = getVarNames(era_interim)  # 使用所有预测因子进行联合PCA
  ),
  combined.only = TRUE   # 使用联合PCA而不是单独降低每个预测因子的维度
)

系统对 8 个区域（AL, EA, SC, FR, MD, ME, BI, IP）依次进行 PCA 分析，进行如下所示:

> 处理区域：FR
> ... 为区域 FR 调用prinComp
> [timestamp] 对20个变量及其组合执行PC分析...
> [timestamp] 完成。

7. 检查data.pca的输出

执行结束后，保存区域PCA结果，我们详细查看这个 data.pca 结果，首先他是一个有四个元素的列表list

y：观测数据(包括83个站点×10958天)

x.global：按PRUDENCE区域获得的主成分矩阵列表

  • x.global$FR → [10958 × 20] 法国区域，10958天 x 20个PCs

  • x.global$SC → [10958 × 41]

  • x.global$MD → [10958 × 49]

  • 列数即该区域保留的主成分数量（满足解释 95% 方差）。每个主成分自动命名为 "PC_1", "PC_2", ...

pca：对应区域的 prinComp 模型对象，便于在 prepareNewData 中将新数据投影到已训练好的主成分空间中。例如 data.pca$pca$FR 包含法国区域的 PCA 模型细节：

  • 预测变量层级信息（level: 500/700/850/1000 hPa）

  • 区域网格坐标（xCoords, yCoords）

  • 每个变量对应的主成分权重与解释方差

  • 日期与季节属性等

x.local 为 NULL ,这个没用，因为关于local predictor的代码我没改动，在我们的工作中不需要使用这个，可以忽略

TO DO (prepareData)

1.data.pca中还有一个对象属性，由于我没改动对应的逻辑，目前为止 nature = "spatial" ，这里我后续会新加一个判断逻辑，改成nature = "region"

2.regional_pca_info信息目前为NULL，因为在prepareData中，我之前没弄这部分，导致regional_pca_info信息没有传过来（见下面），回来再次训练的时候要生成带有这个info信息的结果，这个结果可以送给prepareNewData，因为在prepareNewData中，data.structure$pca 和 data.structure$x.global 都是列表（包含不同区域），因此 prepareNewData 需要先知道在哪个prudence区域里进行转换，然后在data.structure$pca 列表中提取对应区域的 PCA 模型对象，用该PCA模型对象来转换mapping对应区域的 newdata 子集

因此，对应的改动是，不再使用 is.null(data.structure$pca)，而是检查 spatialPred.pars 属性，添加循环遍历data.structure$pca 中的区域，对于每个循环内部，首先我们应该获取该区域的PCA模型对象(regional_pca_model)，然后获得这个区域对应的 newdata 子集(newdata_regional)，而这个区域边界信息就需要从prepareData中通过 regional_pca_info 传过来 ✅

核心实现2：区域PCA -- prepareNewData

prepareNewData 的作用主要是将新的预测数据 newdata 按照 data.structure（由 prepareData）进行标准化、降维，以便用于后续的降尺度预测环节，确保新数据（例如3.2 using historical period of EC-Earth GCM (1979-2008), 3.3 using future RCP8.5 scenario data (2071-2100)）能够与训练时的特征空间一致

大体流程介绍为, 首先 local.predictors一点不改动，对于 global.Predictors，原代码分两种情况处理，第一是如果原 prepareData没用使用PCA，那么直接提取原始变量数据etc...，这部分与我们的工作无关，不改动，忽略。第二，如果 prepareData 使用了PCA，那么判断是全局PCA（同样不改动）还是区域PCA，如果是区域PCA，这就是我们的工作了，则对新数据 newdata在对应 prudence 区域内格点提取 PCs，并然后根据训练时保存的进行center 与 scale 进行标准化，用主成分矩阵（EOFs）进行投影mapping，得到新的 PCs 表示。

下面我们更详细的解释代码每部分流程思路:

首先进行初始检查，检查变量名（训练/预测使用的预测因子变量名一致）和空间坐标匹配（确保 newdata 的原始格点与训练时 prepareData 使用的格点具有相同的空间坐标系），

然后针对 local.predictors ，它从 data.structure$x.local 中获取局地邻近点的索引 local.index.list，然后调用 predictor.nn.values 函数从 newdata 中提取这些局地预测因子的值，它独立于global.predictors ，因此即使这里删除也毫无影响，我这里直接跳过了，没改动过

然后针对 global.predictors ，进行了区域PCA的功能添加，不要在意这个名称，后面可以改region.predictors ，我只是现在原代码基础上改了，还没优化代码。前面说了，这里有很多种情况，同理我们只修改了区域PCA的情况，其他未改动。对应区域PCA入口为

  # 检查数据结构中是否包含区域PCA信息
  regional_info <- attr(data.structure, "regional_pca_info")
  
  # --- 区域PCA处理路径 ---
  if (!is.null(regional_info)) {

代码通过检查 data.structure 是否拥有一个名为 regional_pca_info 的属性来判断 prepareData 阶段是否执行了区域PCA

我们还需要导入PRUDENCE区域定义信息

    # 获取区域定义列表（包含每个区域的经纬度边界信息）
    regions_definition_list <- regional_info$prudence_definitions_used

所有信息都在 regional_pca_info 这里面，所以之前说过，这个regional_pca_info 很重要。

    # 获取需要处理的区域代码列表
    regions_to_process <- names(data.structure$pca)

在我们的数据里，区域代码列表就是 AL, EA, SC, FR, MD, ME, BI, IP

由于我们按区域来进行投影，因此然后得到的结果列表是按区域来投影得到的PC序列，因此我们必须有一个结果列表，元素是每个区域名称，再每个区域里有对应投影后得到的PCs。这个结果列表叫做 newdata.global.list.regional <- list()

然后我们就可以开始循环了，循环处理每个区域 for (region_code in regions_to_process) ，然后我们要获得当前区域的 PCA 模型，即从 data.structure$pca 列表中提取当前 region_code 对应的PCA模型对象（这是之前 prinComp 函数在该区域上的输出）

      # 获取当前区域的PCA模型
      regional_pca_model <- data.structure$pca[[region_code]]

随后我们要获取当前区域训练时生成的PC矩阵，主要用来获取原始的 predictorNames 或 colnames

      # 获取当前区域训练时的PC输出矩阵（用于提取列名和属性）
      regional_pcs_template <- data.structure$x.global[[region_code]]

获取区域边界:

      # 获取当前区域的边界信息
      region_bounds <- regions_definition_list[[region_code]]

提取newdata的区域子集:

      # 提取当前区域 predictors 的空间子集
      # 区域PCA的关键步骤：从全欧洲数据中提取特定区域的网格数据
      message("... 提取 newdata 的区域子集 ", region_code)
      newdata_regional <- subsetGrid(newdata, lonLim = region_bounds$lonLim, latLim = region_bounds$latLim)

最关键的一步，将输入的 newdata 按照当前 region_code 的地理边界进行裁剪，得到只包含该区域数据的 newdata_regional。后续的PCA投影都将在这个区域子集上进行。

后面我们有两种情况，一种是联合PCA，一种是非联合PCA

      # 检查当前区域是使用联合PCA还是非联合PCA
      if (is.null(regional_pca_model$COMBINED)) {

检查当前区域的PCA模型 regional_pca_model 是否包含 COMBINED 组件。

对于非联合PCA路径 (每个变量单独PCA)，只简单介绍解释

首先得到需要单独进行PCA的变量名，对于每个变量，从 newdata_regional 中提取该变量的格点数据 grid_var_regional，然后使用 transformeR::grid2PCs 函数将其投影到该区域、该变量的EOFs上。最后将不同变量投影得到的PC序列（可能包含多个成员 member）通过 cbind 按成员合并起来。

对于联合PCA路径 (变量合并后进行PCA)

我们需要先得到在该区域 联合PCA中所涉及的变量名和联合EOFs信息

        # 获取联合PCA使用的变量列表和特征向量矩阵
        combined_vars <- attr(regional_pca_model$COMBINED, "combined_variables")
        regional_combined_EOFs <- regional_pca_model$COMBINED[[1]]$EOFs

然后根据获得的 regional_combined_EOFs 信息计算 newdata_regional的投影。但是！这里有一个区别，即我是根据原始代码逻辑修改的，原始代码在这里不再调用 grid2PCs函数，而是使用手动投影。因此我们也采用手动投影过程，总的来说，逻辑处理流程：标准化->合并标准化变量->手动矩阵乘法投影

对于联合PCA中的每个变量 var_name 和 newdata_regional 中的每个成员 j_mem_idx，从 regional_pca_model[[var_name]][[1]]$orig 中获取该变量在训练时的标准化参数 (center 和 scale)，然后提取 newdata_regional 中对应变量和成员的数据，对其进行标准化

然后将同一成员内所有变量标准化后的数据通过 cbind 合并成一个大的标准化矩阵。

          # 合并每个成员的所有标准化变量
          combined_std_mat_list_by_member <- lapply(1:n.mem, function(j_mem_idx) {

将上一步得到的每个成员的组合标准化矩阵与该区域的联合EOFs (regional_combined_EOFs) 通过矩阵乘法 %*% 进行投影，得到最终的PC序列。

          # PCA投影：通过矩阵乘法将标准化数据投影到主成分空间
          projected_pcs_list_by_member <- lapply(1:n.mem, function(j_mem_idx) {
              combined_std_mat_list_by_member[[j_mem_idx]] %*% regional_combined_EOFs,

最后为当前区域的投影结果设置适当的列名和属性，尝试从之前获取的 regional_pcs_template (即 data.structure$x.global[[region_code]]) 中提取 predictorNames 属性或 colnames。这是为了确保新投影的PCs与训练时的PCs具有一致的名称。将获取到的名称作为属性赋给当前区域的PC列表:

        # 首先，为整个列表设置predictorNames属性
        attr(projected_pcs_list_by_member, "predictorNames") <- original_predictor_names

进一步为列表中每个成员的PC矩阵（member_pc_matrix）设置列名。

        # 然后，为每个成员矩阵设置列名
        projected_pcs_list_by_member <- lapply(projected_pcs_list_by_member, function(member_pc_matrix) {

最后储存结果，将处理完毕的当前区域的PC数据（projected_pcs_list_by_member）存入以区域代码 region_code 为键的总列表 newdata.global.list.regional 中

      # 将当前区域的处理结果存储到总列表中
      newdata.global.list.regional[[region_code]] <- projected_pcs_list_by_member
      message("... 区域 '", region_code, "' 的投影和赋名处理完成。")
    }

当所有区域都处理完毕后，这个包含了所有区域投影结果的 newdata.global.list.regional 列表被赋给 newdata.global.list，它是函数返回结果之一。

    # 区域循环结束，将区域结果列表赋给最终输出变量
    newdata.global.list <- newdata.global.list.regional

这样区域PCA流程就结束了。

这里有一个关键点值得一提:

新数据 (newdata) 投影计算得到的PCs需要与训练期间的PCs保持完全相同的 predictorNames (例如 "PC_1", "PC_2", ..., "PC_15")

当 prepareData 函数处理完训练数据后，对于某个特定的区域 region_code，我们得到了一组主成分（PCs），例如 k 个PCs。这些PCs通常被命名为 "PC_1", "PC_2", ..., "PC_k"。它们存储在 data.structure$x.global[[region_code]] 中。这些PCs随后被用作输入，去训练一个统计模型，比如广义线性模型 (GLM)，用以降尺度预测特定变量（例如降水）。模型的简化形式可能如下：

g\left( \mathbb{E}[Y \mid \text{PCs}] \right) = \beta_0 + \beta_1 \cdot \text{PC}_1 + \beta_2 \cdot \text{PC}_2 + \dots + \beta_k \cdot \text{PC}_k

β_j ​ (对于 j=1,…,k) 是模型为第 j 个主成分 PC_{j} 学习到的系数， β_2 严格对应 PC_2

当有新的输入数据 newdata 进入 prepareNewData 函数时，对于同一个 region_code，我们会使用在 prepareData阶段存储的该区域的EOF（经验正交函数/加载向量）和标准化参数（来自 data.structure$pca[[region_code]]），将 newdata 投影到与训练阶段完全相同的 k 维主成分空间。这个投影过程会为 newdata 生成一组新的PC时间序列值，我们称之为\text{PC}_1',\text{PC}_2', ...,\text{PC}_k'。这些值代表了 newdata 在那原始 k 个主成分方向上的表现。这个新的PC矩阵 member_pc_matrix的维度也是 [时间序列 × k个PCs]。

接下来，这些新计算得到的PC时间序列值 ( \text{PC}_1', ...,\text{PC}_k') 将被输入到先前训练好的GLM模型中，以获得对Y 的预测值：

\hat{Y} = g^{-1}\left(\beta_0 + \beta_1 \cdot \text{PC}_1' + \beta_2 \cdot \text{PC}_2' + \dots + \beta_k \cdot \text{PC}_k'\right)

因此输入的\text{PC}_j' 必须准确无误地对应于训练时模型学习\beta_j 系数时所使用的原始\text{PC}_j

因此，prepareNewData 中为新投影得到的PC矩阵（member_pc_matrix）的列赋予与训练阶段完全相同的名称（ "PC_1", "PC_2", ..., "PC_k"）

对于非区域PCA部分代码，我也没改动，这就是原本处理全局PCA的部分。

最后返回结果，newdata.refdates 用于获取 newdata 的起止日期

  # 处理日期信息
  newdata.refdates <- list(start = getRefDates(newdata, "start"), end = getRefDates(newdata, "end"))

nature 属性没改，后续可以改

  # 复制'nature'属性（指示数据的性质，如spatial、temporal等）
  attr(newdata.global.list, "nature") <- attr(data.structure, "nature")

将 regional_info 附加到最终输出的 newdata.global.list 的属性中

    attr(newdata.global.list, "regional_pca_info") <- regional_info

最后结果如下

  message("prepareNewData 函数处理完成")
  # 返回的 x.global:
  # - 如果是区域 PCA: 是一个列表，元素为区域代码，值为该区域的 (成员列表 -> PC矩阵)
  # - 如果是全局 PCA 或 原始预测器: 是一个列表，元素为成员名，值为该成员的 (PC矩阵 或 原始特征矩阵)
  return(list("x.global" = newdata.global.list,
              "x.local" = newdata.local.list, # 忽略
              "Dates" = newdata.refdates)
  )
}

  • 如果是区域PCA：x.global 是一个列表，其名称是区域代码 (FR, BI, ...)。每个区域代码对应的值是另一个列表，这个列表的名称是成员名 (member_1, member_2, ...)，每个成员对应的值是该区域、该成员投影后得到的PC矩阵 (时间序列 x PCs数量)。

  • 如果是非区域PCA (全局PCA或原始数据)：x.global 的结构与原始 prepareNewData 结果一样，以成员名 (member_1, ...) 为名称的列表，每个成员对应一个PC矩阵或原始特征矩阵。

核心实现3：GLM模型进行训练和预测 -- glimpr.R

该函数用于实现我们之前说的双模型降水结构:

  • 二项分布GLM (Binomial GLM)：使用 logit 连接函数，预测降水发生与否（湿日/干日）。

  • 伽马分布GLM (Gamma GLM)：使用 log 连接函数，预测湿日的降水量。

这个代码原本已经实现了双模型结构，及训练/预测，按道理不需要更改，但是由于我们最后需要评估指标，其中包含AUC COR，这两个指标就必须需要降水发生概率和湿日条件下的降水强度，因此我们只针对这点进行改动，把代码中中间变量概率和强度也保存到结果中

首先我们解释一下整体流程:

处理来自 modelPars 的 pred.mat (训练期 PCs) 和 sim.mat (预测/模拟期 PCs)

  # ==================================================================================
  #  数据预处理 - 解析 predictors 和模拟矩阵
  # ==================================================================================
  # 处理预报因子矩阵，兼容单矩阵和区域列表两种输入格式
  if (is.matrix(modelPars$pred.mat)) { # 如果是单矩阵格式
    predMat <- modelPars$pred.mat
  } else if (is.list(modelPars$pred.mat) && is.matrix(modelPars$pred.mat[[1]])) {# 如果是区域列表，使用第一个区域的数据
    predMat <- modelPars$pred.mat[[1]]
  } else {
    stop("modelPars$pred.mat must be a matrix or list of matrices")
  }
  
  # 如果未指定主成分个数，使用全部主成分
  if (is.null(n.pcs)) n.pcs <- ncol(predMat)

将观测数据 y$Data 转换成统一的二维矩阵 ymat [时间, 站点数]

  # ==================================================================================
  #  数据预处理 - 准备观测数据矩阵
  # ==================================================================================
  # 根据数据类型转换为标准矩阵格式
  ymat <- if (isTRUE(modelPars$stations)) { y$Data } # 站点数据直接使用
              
          else {    # 格点数据需要转换
            
                  if (is.null(dim(y$Data))) {as.matrix(y$Data) }  # 一维数据转为矩阵
    
                  else {transformeR::array3Dto2Dmat(y$Data)} # 三维数组转为二维矩阵
                    
               }

由于降水发生模型是预测0/1二值事件，因此根据阈值 wet.threshold 把 ymat 转换成二值化形式 ymat.bin。

  # 创建二值化矩阵，用于降水发生概率建模
  ymat.bin <- ymat
  ymat.bin[ymat <  wet.threshold] <- 0L  # 低于阈值记为0（无降水）
  ymat.bin[ymat >= wet.threshold] <- 1L  # 超过阈值记为1（有降水）

随后再计算每个站点的历史湿润日频率。

  # 计算各站点/格点的降水发生频率
  wet.prob <- apply(ymat.bin, 2, function(x) sum(x, na.rm = TRUE) / length(stats::na.exclude(x)))

确定在观测期内有效数据的站点的索引，这里对于我们就是83个，无效数据站点之前早删除了

  # 剔除全部缺测的站点/格点
  rm.ind <- which(is.na(wet.prob))
  cases  <- if (length(rm.ind) > 0) (1:ncol(ymat))[-rm.ind] else 1:ncol(ymat)

然后创建新的列表来额外存储降水发生概率prob.list 和条件强度cond_int.list ，而不仅仅是最终降水预测pred.list。从这里开始进行了改动。

  # ==================================================================================
  #  初始化结果保存矩阵
  # ==================================================================================
  err        <- numeric()              # 记录模型拟合错误
  prob.list  <- vector("list", length(cases))   # 降水概率列表
  cond.list  <- vector("list", length(cases))   # 条件降水强度列表
  pred.list  <- vector("list", length(cases))   # 最终预报值列表
  mod.bin    <- vector("list", length(cases))   # 二值GLM模型列表
  mod.gamma  <- vector("list", length(cases))   # Gamma GLM模型列表

然后遍历每个有效站点/格点进行模型拟合和预测

  # ==================================================================================
  #  模型拟合主循环 - 逐站点建模
  # ==================================================================================
  message("[", Sys.time(), "] Fitting models ...")
  suppressWarnings(
    for (k in seq_along(cases)) {  # 对 cases 中的每个有效站点/格点索引 k 进行循环

拟合二项GLM (降水发生模型)，判断某天是否降水 （未改动），即为当前正在处理的站点 cases[k] 建立并拟合一个逻辑回归GLM模型。

      #  第一步：拟合 伯努利 二值GLM模型 - 预报降水发生概率
      # ------------------------------------------------------------------------------
      # 使用逻辑回归模型，以 区域PCs 为预报因子
      mod.bin.k <- stats::glm(ymat.bin[, cases[k]] ~ ., # 响应变量: 当前站点 cases[k] 在训练期的0/1发生序列
                              data   = as.data.frame(predMat)[, 1:n.pcs], # 预测变量: 训练PCs predMat 的前 n.pcs 列
                              family = stats::binomial(link = "logit"))

为当前正在处理的站点建立并拟合一个广义线性模型，更具体的话，是逻辑回归GLM模型

ymat.bin 是一个二维矩阵 [时间, 站点]，其中存储的是观测期间每日的降水发生情况（1代表湿日，0代表干日），cases[x] 是当前站点的索引，所以，ymat.bin[ ,cases[x]] 提取了当前站点在整个观测期（训练期）的每日降水发生序列（一列0和1的值）

predMat[, 1:n.pcs] 提供了用于训练的预测因子（即训练期的PCA主成分的前 n.pcs 个）。

最后stats::glm 函数进行参数估计，计算出模型中各个预测因子（PCs）的系数，以获得统计关系，将拟合得到的整个模型对象（包含系数、拟合优度信息等）存储在变量 mod.bin.k（当前站点 k 的二项分布模型），这个 mod.bin.k 对象后续可以被 stats::predict() 函数用来对新的预测因子数据进行预测，得到降水发生的概率。

注意:

这里有一个问题，就是当前modelPars$pred.mat 是一个单一的、全局的矩阵，但是在区域PCA中，理论上应该变成一个按区域组织的list，但是我没修改这里，我修改了主函数的调用方法，主函数中针对每个当前区域，来调用 glimpr ，这样就使得:

  • glimpr 函数接收到的 predMat (源自 modelPars$pred.mat) 是为该特定区域准备的训练期PCs (单一矩阵)。

  • glimpr 函数接收到的 modelPars$sim.mat 是一个列表(即 simMat.list)，其元素（这里只有一个元素）也是一个单一矩阵(sim_mat_regional)，这个矩阵是为该特定区域准备的预测期PCs。

  • glimpr 函数接收到的 y (y_subset_for_region) 只包含该区域内的站点数据。因此，glimpr 内部的 for (x in 1:length(cases)) 循环实际上只会遍历当前区域的站点，并使用当前区域的预测因子（pred_mat_regional 和 sim_mat_regional）为这些站点建模和预测。

预测降水发生概率并生成0/1发生序列

使用训练好的发生模型 mod.bin.x 对新数据（modelPars$sim.mat）进行预测，得到原始发生概率，并将其存储。同时，根据这些概率生成0/1的发生序列。

      # 计算第一个成员的降水概率（用于输出）
      probs.first <- stats::predict(mod.bin.k,
                                    newdata = as.data.frame(simMat.first)[, 1:n.pcs],
                                    type = "response")

使用 mod.bin.k 对 simMat.first (预测/模拟期的第一个成员的PCs) 进行预测，得到原始的发生概率 \hat{π}。这个结果将用于最终输出。

      # 对所有集合成员进行降水发生预报
      sims.bin <- sapply(seq_along(simMat.list), function(i) {
        # 计算降水概率
        p <- stats::predict(mod.bin.k,
                            newdata = as.data.frame(simMat.list[[i]])[, 1:n.pcs],
                            type = "response")
        
        # 根据模拟方式来决定如何处理概率
        if (simulate != "no") 
          {
          # 随机模拟：基于概率随机生成0/1
          rnd <- stats::runif(length(p))
          p[p > rnd] <- 1L; p[p <= rnd] <- 0L
          } 
        
        else 
          {
          # 确定性模式：基于历史频率阈值判断
          p[p >= stats::quantile(mod.bin.k$fitted.values,
                                 1 - wet.prob[cases[k]],
                                 na.rm = TRUE)] <- 1L
          p <- as.integer(p)
          }
        p
      })

对于 simMat.list 中的每个成员 i，首先计算原始发生概率 p，然后，根据 simulate 方法将概率 p 转换为0/1的发生序列:

  • 随机模式 (simulate != "no"): 通过与一个均匀分布的随机数比较来决定发生与否。

  • 确定性模式 (simulate == "no"): 将概率与一个基于训练期拟合值和观测湿日频率计算得到的阈值 threshold_prob 进行比较来决定发生与否。

最终，sims.bin 存储了每个预测成员在每个时间点的0/1降水发生情况。

拟合伽马GLM -- 为当前站点训练一个预测湿日降水量的模型

伽马模型只在观测到的湿日 (subset = wet.idx) 上进行训练，使用 predMat (训练期的PCs) 作为预测因子。

      # ------------------------------------------------------------------------------
      #  第二步：拟合Gamma GLM模型 - 预报降水强度（仅针对降水日）
      # ------------------------------------------------------------------------------
      # 筛选有降水的样本
      wet.idx <- which(ymat[, cases[k]] >= wet.threshold)
      
      # 逐步减少主成分数直到模型收敛
      Npcs    <- n.pcs + 1
      mod.g.k <- NULL
      while (is.null(mod.g.k) && Npcs > 1) {
        Npcs <- Npcs - 1
        mod.g.k <- tryCatch(
          stats::glm(ymat[, cases[k]] ~ .,
                     data   = as.data.frame(predMat)[, 1:Npcs],
                     family = stats::Gamma(link = "log"),
                     subset = wet.idx),
          error = function(e) NULL)
      }
      
      # 记录是否使用了全部主成分
      if (length(mod.g.k$model) - 1 < n.pcs) err[k] <- 1 else err[k] <- 0

先找到训练期中当前站点 cases[k] 的湿日索引wet.idx，对于while 循环，这是原本函数里的鲁棒性设计，从 Npcs_gamma_initial（通常等于 n.pcs）个主成分开始尝试拟合Gamma GLM，如果出现拟合失败的情况，比如说湿润日太少，或者预测因子问题导致模型不收敛，那么tryCatch 会捕获错误并使 mod.g.k 保持/变为 NULL ，然后循环会减少 Npcs 的数量再试一次，直到模型成功拟合或者 Npcs 减到1。 此逻辑在原代码基础上未做改动，只做适配

  • stats::glm(...) subset = wet.idx: 核心的Gamma模型拟合，只使用了湿润日的观测数据。

  • err[k]: 用于记录Gamma模型是否成功使用最初数量的PCs (Npcs_gamma_initial) 拟合。如果PC数量被减少或模型最终未能拟合，err[k] 设为1。

预测降水量并存储条件强度

使用训练好的伽马模型 mod.g.k 对新数据 (simMat.list) 进行预测，得到原始条件强度，并计算最终的降水预测值。

      # ------------------------------------------------------------------------------
      #  第三步：生成最终预报
      # ------------------------------------------------------------------------------
      if (is.null(mod.g.k)) {
        # 模型拟合失败，返回缺测值
        cond.first <- rep(NA_real_, nrow(simMat.first))
        preds.all  <- sapply(seq_along(simMat.list), function(i)
          rep(NA_real_, nrow(simMat.list[[i]])))
      } else {
        # 计算条件降水强度（第一个成员）
        cond.first <- stats::predict(mod.g.k,
                                     newdata = as.data.frame(simMat.first)[, 1:n.pcs],
                                     type = "response")
        
        # 生成所有成员的降水预报
        preds.all <- sapply(seq_along(simMat.list), function(i) {
          # 预报降水强度
          pg <- stats::predict(mod.g.k,
                               newdata = as.data.frame(simMat.list[[i]])[, 1:n.pcs],
                               type = "response")
          
          # 根据模拟模式生成最终降水量
          if (simulate == "yes") {
            # 随机模拟：从Gamma分布采样
            rgamma(length(pg),
                   shape = 1 / summary(mod.g.k)$dispersion,
                   scale = summary(mod.g.k)$dispersion * pg) * sims.bin[, i]
          } else {
            # 确定性预报：强度乘以发生概率
            pg * sims.bin[, i]
          }
        })
      }

如果伽马模型 mod.g.k 未能成功拟合，则该站点的条件强度 cond.first 和所有成员的最终预测 preds.all 都将被设为 NA。

如果伽马模型 mod.g.k 拟合成功，使用 mod.g.k 对 simMat.first (第一个预测成员的PCs，注意这里使用 Npcs，即Gamma模型实际使用的PC数) 进行预测，得到条件平均降水量 \hat{μ​}

preds.all: 这是一个矩阵，维度为 [预测时间步数, 成员数]。具体计算步骤: 对于 simMat.list 中的每个成员 i ，首先计算其条件平均降水量 pg ，然后根据 simulate 方式计算最终降水量:

  • simulate == "yes" (完全随机模拟): 从一个形状参数为 1 / summary(mod.g.k)$dispersion、尺度参数为 summary(mod.g.k)$dispersion * pg 的伽马分布中随机抽样得到降水量，然后乘以该成员的0/1发生指示 sims.bin[, i]。

  • simulate == "no" 或 "occurrence" (确定性强度): 直接将预测的条件平均降水量 pg乘以该成员的0/1发生指示 sims.bin[, i]。

存储当前站点的所有类型结果

      # ------------------------------------------------------------------------------
      #  保存当前站点/格点的结果
      # ------------------------------------------------------------------------------
      prob.list[[k]] <- probs.first         # 降水概率
      cond.list[[k]] <- cond.first          # 条件强度
      pred.list[[k]] <- preds.all[, 1]      # 第一个成员的最终预报
      
      mod.bin[[k]]   <- mod.bin.k           # 二值模型
      mod.gamma[[k]] <- mod.g.k             # 强度模型
    }
  )
  
  message("[", Sys.time(), "] Done.")

补充知识: 第一个成员的意思就是 member ，在我们的数据中 member=1 ，这个member的意思涉及到气候模型的运行模式。简要来说，气候模型为了捕捉预报的不确定性，所以不可能只运行一次模型，通常会运行好几次，每次运行的时候，模型的初始条件都会略微不同，或者说模型参数会有所调整，那么就导致了得到多个预报结果，这多个预报结果就组成了一个集合，即 member列表。

后处理：填充NA并组合输出数组

确保所有输出列表（pred.list, prob.list, cond_int.list）都包含与输入y同样数量的站点/格点，对于那些因数据不足而未进行建模的站点，其对应的元素填充为NA矩阵。

  # ==================================================================================
  #  结果整理 - 输出矩阵
  # ==================================================================================
  # 创建完整的站点矩阵
  n.test <- length(prob.list[[1]])
  probs.mat <- cond.mat <- preds.mat <- matrix(NA_real_, nrow = n.test, ncol = ncol(ymat))
  
  # 填充有效站点的结果
  probs.mat[, cases] <- do.call(cbind, prob.list)
  cond.mat [, cases] <- do.call(cbind, cond.list)
  preds.mat[, cases] <- do.call(cbind, pred.list)

n.test 是预测期的时间长度。初始化三个与原始 ymat 列数相同的 NA 矩阵：probs.mat，cond.mat，preds.mat。

使用 do.call(cbind, ...) 将之前为每个有效站点 (cases) 存储在 prob.list、cond.list 和 pred.list 中的向量结果按列合并，并填充到这三个大矩阵的相应 cases 列中。那些无效站点（不在 cases 中）的列将保持为 NA。

5. 调整并返回最终输出

将所有处理好的结果（最终预测、发生概率、条件强度以及可选的模型对象）打包成一个列表返回。

  # ==================================================================================
  #  构建返回对象
  # ==================================================================================
  # 对于单站点情况，返回向量；多站点返回矩阵
  out <- list(
    probabilities    = probs.mat[, 1],        # 降水发生概率
    cond_intensities = cond.mat [, 1],        # 条件降水强度（给定降水发生）
    predictions      = preds.mat[, 1],        # 最终降水预报值
    models           = list(binary = mod.bin, gamma = mod.gamma)  # 拟合的模型对象
  )
  
  return(out)
}

这里都取第一列的原因其实之前解释过，主要原因是glimpr函数的调用方式，在主函数中，我们对站点内部循环中调用glimpr 函数，也就是每个站点都构建一个glm模型，每次调用 glimpr 时，传递给它的观测数据 y (即主脚本中的 y_subset_for_station) 只包含单个站点的数据，因此，在 glimpr 函数内部，ymat 变量（由 y$Data 转换而来）实际上只会有一个列 (ncol(ymat) 等于1)，cases （有效站点的索引）也只会包含一个元素，即 1。内部的循环 for (k in seq_along(cases)) 也只会执行一次 (k=1)。

所以最后最终生成的 probs.mat、cond.mat 和 preds.mat 这些矩阵在这种单站点调用情景下，它们都只会有一列。

因为我原本打算在glimpr里实现所有站点循环构建模型，后面想通了在主函数里挨个调用算了，就没改动这部分代码。

主函数调用

六月三日补充，这部分主函数没有交叉验证，且第二部分指标计算方面有错误，后续要改，目前放在这里做备份

1. 环境设置与脚本加载

# ############################################################################
# ## Première partie: PCA régionale sans validation croisée -
# ## Traitement des données + Obtention des paramètres de distribution + Simulation aléatoire paramétrée et moyenne
# ############################################################################

# ############################################################################
# ## 1) Définir le répertoire de travail & Charger les packages
# ############################################################################
setwd("/Users/zehua/Downloads/PCA_regional")

library(transformeR)
library(downscaleR)
library(magrittr)
library(abind)
library(pROC)
library(ggplot2)
library(dplyr)
library(tidyr)
library(maps)
library(fields)
library(RColorBrewer)
library(geosphere)  # La fonction map_station_to_prudence_region est nécessaire de ça

# ############################################################################
# ## 2) Charger le script (prepareData de la version régionale et glimpr modifié)
# ############################################################################
source("glimpr.R")          # Modifié, peut renvoyer prob, cond_int, modèles
source("prepareData.R")     # Zone PCA version prepareData
source("prepareNewData.R")  # Zone PCA version prepareNewData
source("helpers.R")

2. PRUDENCE 区域定义与站点映射函数

定义地理区域的边界，并创建一个函数将具体的气象站点映射到这些预定义的区域中

# Les fonctions auxiliaires predictor.nn.indices et predictor.nn.values sont définies dans prepareData

# --- Définir la fonction de mappage des régions et des sites PRUDENCE ---
prudence_regions <- list(
  BI = list(name = "British Isles", lonLim = c(-10, 2), latLim = c(50, 59)),
  IP = list(name = "Iberian Peninsula", lonLim = c(-10, 3), latLim = c(36, 44)),
  FR = list(name = "France", lonLim = c(-5, 5), latLim = c(44, 50)),
  ME = list(name = "Mid-Europe", lonLim = c(2, 16), latLim = c(48, 55)),
  SC = list(name = "Scandinavia", lonLim = c(5, 30), latLim = c(55, 70)),
  AL = list(name = "Alps", lonLim = c(5, 16), latLim = c(44, 48)),
  MD = list(name = "Mediterranean", lonLim = c(3, 25), latLim = c(36, 44)),
  EA = list(name = "Eastern Europe", lonLim = c(16, 30), latLim = c(44, 55))
)

map_station_to_prudence_region <- function(lon, lat, regions_definition) {
  region_name_inside <- NA
  for (region_code in names(regions_definition)){
    bounds <- regions_definition[[region_code]]
    if (lon >= bounds$lonLim[1] && lon < bounds$lonLim[2] &&
        lat >= bounds$latLim[1] && lat < bounds$latLim[2] ) {
      region_name_inside <- region_code
      return(region_name_inside) # 找到即返回
    }
  }
  # Si le site ne se trouve dans aucune zone rectangulaire prédéfinie, il est attribué à la zone la plus proche du point central
  min_dist <- Inf; closest_region <- NA; point_coord <- c(lon, lat)
  for (region_code in names(regions_definition)){
    bounds <- regions_definition[[region_code]]; center_lon <- mean(bounds$lonLim); center_lat <- mean(bounds$latLim)
    region_center_coord <- c(center_lon, center_lat)
    current_dist <- geosphere::distHaversine(point_coord, region_center_coord)
    if (current_dist < min_dist){ min_dist <- current_dist; closest_region <- region_code }
  }
  if (is.na(closest_region) && length(names(regions_definition)) > 0) {
    closest_region <- names(regions_definition)[1]
    warning("Site (", lon, ",", lat, ") ne affecté pas à une zone spécifique et incapable de déterminer la zone la plus proche, affecté à la première zone : ", closest_region)
  }
  return(closest_region)
}

这部分在 prepareData 里也有，这里为了图方面直接搬过来，后面优化的时候需要改

3. 数据加载与站点数据清洗


# ############################################################################
# ## 3) Charger les données & Nettoyer les sites non valides
# ############################################################################
load("st86_precip.rda")
load("era_interim.rda")

remove_ids <- c("000017", "000274", "000275")
st86_precip_original_coords <- st86_precip$xyCoords
st86_precip_original_ids <- st86_precip$Metadata$station_id

idx_remove <- which(st86_precip$Metadata$station_id %in% remove_ids)

if (length(idx_remove) > 0) {
  st86_precip$Data <- st86_precip$Data[, -idx_remove, drop = FALSE]
  st86_precip$xyCoords <- st86_precip$xyCoords[-idx_remove, , drop = FALSE]
  st86_precip$Metadata$station_id <- st86_precip$Metadata$station_id[-idx_remove]
  st86_precip$Metadata$name <- st86_precip$Metadata$name[-idx_remove]
  st86_precip$Metadata$altitude <- st86_precip$Metadata$altitude[-idx_remove]
  st86_precip$Metadata$source <- st86_precip$Metadata$source[-idx_remove]
}

if (is.null(attr(st86_precip$Data, "dimensions"))) {
  attr(st86_precip$Data, "dimensions") <- c("time", "loc")
}
message("Le chargement et le nettoyage des données sont terminés. Nombre de sites après nettoyage:", ncol(st86_precip$Data))

# --- Obtenir les coordonnées et l'identifiant du site final utilisé (après nettoyage) ---
final_station_coords <- as.matrix(st86_precip$xyCoords) # Assure que c'est une matrice
final_station_ids <- st86_precip$Metadata$station_id
n_final_stations <- nrow(final_station_coords)

4. 使用区域PCA进行降维 (调用 prepareData)

使用支持区域PCA的 prepareData 函数，对大尺度预测因子 era_interim 进行降维，为每个PRUDENCE区域生成主成分（PCs）。

# ############################################################################
# ## 4) Utilisation de l'ACP régionale pour la réduction de dimensionnalité
# ############################################################################
message("En train d'utiliser PCA de zone pour la réduction de dimensionnalité (prepareData)...")

data.pca <- prepareData(
  x = era_interim,
  y = st86_precip, # Données du site nettoyées
  spatial.predictors = list(
    v.exp = 0.95,
    which.combine = getVarNames(era_interim)
  ),
  combined.only = TRUE
)
message("La réduction de dimension PCA de la zone a été effectuée")
# Résultat de sortie 1 : data.pca$x.global est une liste avec la clé comme le code de la région et la valeur comme la matrice PC de cette région (ou une liste de membres->matrice PC)
# Résultat de sortie 2 : data.pca$pca est également une liste de structure similaire, contenant des objets de modèle PCA pour chaque région.
# Résultat de sortie 3: attr(data.pca, "regional_pca_info") contient des informations de traitement régional

5. 将新数据映射到区域PC空间

: 使用支持区域PCA的 prepareNewData 函数，将（通常是预测期的）大尺度预测因子数据投影到由 prepareData 建立的各个区域的PC空间上。

这里对新数据的选择有问题，仍然选择 era_interim 因为没有新数据，拿这个顶一下

# ############################################################################
# ## 5) Mapper les nouvelles données (era_interim) à l'espace PC régional
# ############################################################################
message("En train de mapper era_interim (réanalyse) à l'espace PC régional (prepareNewData)...")
# data.structure est maintenant data.pca, qui contient des informations de PCA régionalisées.
newdata.pca <- prepareNewData(era_interim, data.pca)
message("La cartographie des données est terminée.")
# newdata.pca$x.global est également une liste organisée par région, avec une structure similaire à celle de data.pca$x.global

if(!is.list(newdata.pca$x.global) || !identical(names(newdata.pca$x.global), names(data.pca$x.global))){
  stop("La structure de newdata.pca$x.global retournée par prepareNewData ne correspond pas à la région.")
}

newdata.pca$x.global: 也是一个按区域组织的列表。每个区域代码下是一个列表，该列表的元素是对应区域、对应成员（目前是单成员）的预测期的PC矩阵 [时间, 该区域PC数量]

6.A 按区域调用 glimpr 获取降水分布参数

对每个区域内的每个站点，使用该区域的PCs作为预测因子，调用 glimpr 函数来拟合GLM模型并获取降水分布的关键参数（发生概率 \pi 和条件强度 \mu），并收集训练好的伽马模型。

# ############################################################################
# ## 6) Obtenir les paramètres de distribution de chaque site, puis effectuer une simulation aléatoire paramétrique
# ############################################################################
message("En train d'obtenir les paramètres de distribution des précipitations pour chaque site (appel régional à glimpr)...")

# --- A. Résumé des paramètres de tous les sites ---
# Initialiser le conteneur stockant tous les paramètres du site
n_total_days <- nrow(data.pca$y$Data)
prob_estimates_all_stations <- matrix(NA_real_, nrow = n_total_days, ncol = n_final_stations)
cond_mean_intensity_all_stations <- matrix(NA_real_, nrow = n_total_days, ncol = n_final_stations)
gamma_models_all_stations <- vector("list", n_final_stations)

# Déterminez à quelle région chaque site appartient
station_region_mapping <- sapply(1:n_final_stations, function(i) {
  map_station_to_prudence_region(final_station_coords[i, 1], final_station_coords[i, 2], prudence_regions)
})
names(station_region_mapping) <- final_station_ids

# Appeler glimpr région par région
processed_regions <- names(data.pca$x.global) # Obtenir la zone réellement traitée à partir du résultat de prepareData

for (current_region_code in processed_regions) {
  message("  En cours de traitement de la région: ", current_region_code)

  # Trouver l'index du site appartenant à la zone actuelle
  station_indices_in_region <- which(station_region_mapping == current_region_code)

  if (length(station_indices_in_region) == 0) {
    message("    Région ", current_region_code, " a pas de site, sauter")
    next
  }

  # Préparez y (observation) de la région actuelle
  y_regional <- data.pca$y

  # Préparez pred.mat et sim.mat de la région actuelle
  pred_mat_regional <- data.pca$x.global[[current_region_code]]
  # Membre unique, donc prendre le premier élément de la liste
  sim_mat_regional  <- newdata.pca$x.global[[current_region_code]][[1]]

  # # Assure que pred_mat_regional et sim_mat_regional sont des matrices
  # if (is.null(pred_mat_regional) || !is.matrix(pred_mat_regional)) { # 额外检查 NULL
  #   stop("区域 ", current_region_code, " 的 pred_mat_regional 不是一个有效的矩阵。可能该区域在 prepareData 中处理失败。")
  # }
  # if (is.null(sim_mat_regional) || !is.matrix(sim_mat_regional)) { # 额外检查 NULL
  #   stop("区域 ", current_region_code, " 的 sim_mat_regional (成员1) 不是一个有效的矩阵。可能该区域在 prepareNewData 中处理失败。")
  # }

  # Obtenir le nombre actuel de PC dans la région
  n_pcs_regional <- ncol(pred_mat_regional)
  if (is.null(n_pcs_regional)) stop("Impossible de déterminer le nombre de PCs de la région ", current_region_code, " ")


  # Appeler glimpr pour le site de la région actuelle
  #
  # 重要: glimpr 内部的站点循环需要能正确对应到 y_regional 中属于当前区域的站点
  # 实际上，glimpr 是对 y_regional$Data 的所有列进行操作。
  # 如果 y_regional$Data 包含所有站点，而 pred_mat_regional 只对应一个区域，这是不匹配的。

  # Préparer y_subset pour glimpr contenant uniquement le site actuel de la région
  y_subset_for_region <- list(
    Data = data.pca$y$Data[, station_indices_in_region, drop = FALSE],
    xyCoords = data.pca$y$xyCoords[station_indices_in_region, , drop = FALSE],
    Dates = data.pca$y$Dates,
    Metadata = list(
      station_id = data.pca$y$Metadata$station_id[station_indices_in_region]

    )
  )
  # Définir les attributs de dimension pour y_subset_for_region (erreur si non ça)
  if (is.null(attr(y_subset_for_region$Data, "dimensions")) && !is.null(data.pca$y$Data)) {
    attr(y_subset_for_region$Data, "dimensions") <- attr(data.pca$y$Data, "dimensions")
  }


  message("    Pour la région '", current_region_code, "' la site", length(station_indices_in_region), " appel à glimpr...")

  # --- Code de diagnostic: Cette partie est inutile, utilisée pour le débogage ---
  message("---- Région: ", current_region_code, " ----")

  if (!is.null(pred_mat_regional) && is.matrix(pred_mat_regional)) {
    message("Les noms de colonnes de pred_mat_regional (colnames):")
    print(colnames(pred_mat_regional))
    message("La propriété 'predictorNames' de pred_mat_regional:")
    print(attributes(pred_mat_regional)$predictorNames)
  } else {
    message("pred_mat_regional n'est pas une matrice ou est NULL. La structure est la suivante:")
    str(pred_mat_regional)
  }

  if (!is.null(sim_mat_regional) && is.matrix(sim_mat_regional)) {
    message("Les noms de colonnes de sim_mat_regional (colnames):")
    print(colnames(sim_mat_regional))
    message("La propriété 'predictorNames'de sim_mat_regional:")
    # sim_mat_regional est une matrice de membres, dont l'attribut predictorNames est basé sur newdata.pca$x.global[[current_region_code]]
    print(attributes(newdata.pca$x.global[[current_region_code]])$predictorNames)
  } else {
    message("sim_mat_regional n'est pas une matrice ou NULL. La structure est la suivante:")
    str(sim_mat_regional)
  }
  message("---- Fin des informations de diagnostic ----")



  glimpr_region_result <- glimpr(
    y               = y_subset_for_region, #y contient uniquement des sites de la région actuelle
    modelPars       = list(
      stations      = TRUE,
      multi.member  = FALSE,
      pred.mat      = pred_mat_regional,  # PCs d'entraînement de la région actuelle
      sim.mat       = list(sim_mat_regional), # Prévisions des PC pour la région actuelle
      sim.dates     = newdata.pca$Dates    # Objet de date global
    ),
    simulate        = "no",
    return.models   = TRUE,
    wet.threshold   = 1,
    n.pcs           = n_pcs_regional
  )

  # Enregistrez les résultats de la région actuelle dans la matrice/liste de paramètres globale
  # Le résultat de glimpr_region_result est pour les sites dans y_subset_for_region,Donc, le nombre de colonnes est égal à la longueur (station_indices_in_region)
  prob_estimates_all_stations[, station_indices_in_region] <- glimpr_region_result$probabilities
  cond_mean_intensity_all_stations[, station_indices_in_region] <- glimpr_region_result$cond_intensities
  # gamma_models_list: Il est nécessaire de stocker selon l'index du site d'origine
  for(k in 1:length(station_indices_in_region)){
    original_station_idx <- station_indices_in_region[k]
    gamma_models_all_stations[[original_station_idx]] <- glimpr_region_result$models$gamma[[k]]
  }
} # Terminer la boucle de la zone

message("Les paramètres glimpr de toutes les régions ont été récupérés")

6.B 参数化随机模拟

利用上一步从 glimpr获得的每个站点的降水发生概率 \pi 和条件强度 \mu以及伽马模型的离散度参数，进行500次随机模拟，以生成降水的集合预报或更好地表征不确定性。

# --- B. Effectuer 500 simulations paramétriques ---
message("Commencer la simulation aléatoire paramétrique 500 fois...")
n_simulations <- 500

all_simulations_array <- array(NA_real_, dim = c(n_total_days, n_final_stations, n_simulations))

for (s_idx in 1:n_simulations) {
  current_simulation_run_matrix <- matrix(0, nrow = n_total_days, ncol = n_final_stations)
  for (i_station in 1:n_final_stations) {
    pi_daily   <- prob_estimates_all_stations[, i_station]
    mu_daily   <- cond_mean_intensity_all_stations[, i_station]
    gamma_model <- gamma_models_all_stations[[i_station]]

    if (is.null(gamma_model) || any(is.na(pi_daily)) || any(is.na(mu_daily))) {
      current_simulation_run_matrix[, i_station] <- NA; next
    }
    dispersion_i <- summary(gamma_model)$dispersion
    if (is.na(dispersion_i) || dispersion_i <= 0) {
      occurrence_sim_daily <- rbinom(n_total_days, 1, pi_daily)
      wet_days_idx_fallback <- which(occurrence_sim_daily == 1 & !is.na(mu_daily) & mu_daily > 0)
      current_simulation_run_matrix[wet_days_idx_fallback, i_station] <- mu_daily[wet_days_idx_fallback]
      next
    }
    shape_i <- 1 / dispersion_i
    occurrence_sim_daily <- rbinom(n_total_days, 1, pi_daily)
    wet_days_indices <- which(occurrence_sim_daily == 1)
    if (length(wet_days_indices) > 0) {
      for (d_idx in wet_days_indices) {
        if (!is.na(mu_daily[d_idx]) && mu_daily[d_idx] > 0) {
          scale_daily_d <- dispersion_i * mu_daily[d_idx]
          if (shape_i > 0 && scale_daily_d > 0 && is.finite(shape_i) && is.finite(scale_daily_d)) {
            current_simulation_run_matrix[d_idx, i_station] <- rgamma(1, shape = shape_i, scale = scale_daily_d)
          } else {
            current_simulation_run_matrix[d_idx, i_station] <- 0 # 或 NA
          }
        } else { current_simulation_run_matrix[d_idx, i_station] <- 0 }
      }
    }
  } # Terminer la boucle de la station
  all_simulations_array[,, s_idx] <- current_simulation_run_matrix
  if (s_idx %% 50 == 0) message("  Terminer le nombre de simulations paramétriques: ", s_idx, "/", n_simulations)
} # Terminer la boucle de simulation
message("La simulation aléatoire paramétrée est terminée")

message("Calcul de la moyenne des résultats de la simulation paramétrique en cours...")
pred_mean <- apply(all_simulations_array, MARGIN = c(1, 2), FUN = mean, na.rm = TRUE)
message("Les résultats moyens de la simulation paramétrique ont été calculés.")

rm(all_simulations_array, current_simulation_run_matrix)
gc()

7. 保存结果

# ############################################################################
# ## 7) Les données à sauvegarder
# ############################################################################
obs  <- data.pca$y$Data # Utilisez data.pca$y$Data comme observations finalement alignées

prob_for_auc      <- prob_estimates_all_stations
cond_int_for_cor  <- cond_mean_intensity_all_stations
expected_value    <- prob_for_auc * cond_int_for_cor
expected_value[prob_for_auc == 0 & !is.na(prob_for_auc)] <- 0
expected_value[is.na(cond_int_for_cor) & !is.na(prob_for_auc)] <- NA


save(obs, pred_mean, prob_for_auc, cond_int_for_cor, expected_value, final_station_coords, # 保存 final_station_coords
     file = "regional_pca_glm_param_sim_averaged_results.rda") 

message("Terminé : exécution réussie de la simulation paramétrique et enregistrement des résultats moyens de prédiction, des estimations de paramètres et des valeurs attendues dans 'regional_pca_glm_param_sim_averaged_results.rda' ")

将关键的中间和最终结果保存到 .rda 文件

第二部分：计算所有评估指标

1. 加载之前保存的结果

# ############################################################################
# ## Deuxième partie: Calculer tous les indicateurs d'évaluation
# ############################################################################
message("\n--- Commencer à calculer les indicateurs d'évaluation ---")
# --- Chargement ---
load("regional_pca_glm_param_sim_averaged_results.rda")
# Il y a actuellement: obs, pred_mean, prob_for_auc, cond_int_for_cor, expected_value, final_station_coords

2. 数据有效性检查

确保所有必需的数据都已成功加载，并且维度匹配，以避免在后续计算中出错

# Renommez final_station_coords en station_coords pour correspondre au code ultérieur
station_coords <- final_station_coords 

wet_threshold <- 1

# Vérification de la validité des données
if (!exists("station_coords") || !is.matrix(station_coords) || ncol(station_coords) != 2) {
  stop("Impossible de charger 'station_coords' ou son format est incorrect")
}
if (!exists("obs") || !exists("pred_mean") || !exists("prob_for_auc") || !exists("expected_value")) {
  stop("Impossible de charger les données 'obs', 'pred_mean', 'prob_for_auc' ou 'expected_value'.")
}
if (nrow(station_coords) != ncol(obs) || ncol(obs) != ncol(pred_mean) || ncol(obs) != ncol(prob_for_auc) || ncol(obs) != ncol(expected_value)) {
  stop("Le nombre de coordonnées du site ne correspond pas au nombre de sites dans les données d'observation/prédiction/probabilité/valeur attendue.")
}
message("Les données d'évaluation ont été chargées.")

3. 定义评估指标计算函数

# calculate_R01, SDII, P98, DW, WW, calculate_AUC, calculate_COR_paper
calculate_R01 <- function(precip_series, threshold) {
  valid_days <- sum(!is.na(precip_series)); if (valid_days == 0) return(NA)
  wet_days <- sum(precip_series >= threshold, na.rm = TRUE); return(wet_days / valid_days)
}
calculate_SDII <- function(precip_series, threshold) {
  wet_days_precip <- precip_series[precip_series >= threshold & !is.na(precip_series)]
  if (length(wet_days_precip) == 0) return(NA); return(mean(wet_days_precip))
}
calculate_P98 <- function(precip_series, threshold) {
  wet_days_precip <- precip_series[precip_series >= threshold & !is.na(precip_series)]
  if (length(wet_days_precip) < 2) return(NA)
  return(quantile(wet_days_precip, probs = 0.98, na.rm = TRUE, type = 8))
}
create_transition_data <- function(precip_series, threshold) {
  is_wet <- precip_series >= threshold & !is.na(precip_series); if (length(is_wet) < 2) return(NULL)
  data.frame(current = is_wet[-1], previous = is_wet[-length(is_wet)])
}
calculate_DW <- function(precip_series, threshold) {
  trans_data <- create_transition_data(precip_series, threshold); if (is.null(trans_data)) return(NA)
  n_dry_days <- sum(!trans_data$previous); if (n_dry_days == 0) return(NA)
  n_dry_to_wet <- sum(!trans_data$previous & trans_data$current); return(n_dry_to_wet / n_dry_days)
}
calculate_WW <- function(precip_series, threshold) {
  trans_data <- create_transition_data(precip_series, threshold); if (is.null(trans_data)) return(NA)
  n_wet_days <- sum(trans_data$previous); if (n_wet_days == 0) return(NA)
  n_wet_to_wet <- sum(trans_data$previous & trans_data$current); return(n_wet_to_wet / n_wet_days)
}
calculate_AUC <- function(probabilities, observed_binary) {
  valid_idx <- !is.na(probabilities) & !is.na(observed_binary)
  if (sum(valid_idx) < 10 || length(unique(observed_binary[valid_idx])) < 2) { return(NA) }
  roc_obj <- tryCatch(roc(response = observed_binary[valid_idx], predictor = probabilities[valid_idx], quiet = TRUE), error = function(e) NULL)
  if (is.null(roc_obj)) return(NA); return(as.numeric(auc(roc_obj)))
}
calculate_COR_paper <- function(observed_precip, predicted_expected_value) {
  cor_val <- tryCatch(cor(observed_precip, predicted_expected_value, method = "spearman", use = "pairwise.complete.obs"), error = function(e) NA)
  return(cor_val)
}

4. 为每个站点计算评估指标

遍历所有站点，使用上面定义的函数分别为观测数据和模型预测数据（参数化模拟的均值 pred_mean和期望值 expected_value）计算各项评估指标


message("En train de calculer tous les indicateurs d'évaluation...")
n_stations_eval <- ncol(obs)
metrics_results_all <- list()
obs_binary <- ifelse(obs >= wet_threshold & !is.na(obs), 1L, 0L)
obs_binary <- apply(obs_binary, 2, as.integer)

for (i in 1:n_stations_eval) {
  obs_series <- obs[, i]
  pred_series_for_distrib <- pred_mean[, i]
  prob_series_for_auc <- prob_for_auc[, i]
  exp_series_for_cor <- expected_value[, i]
  obs_bin_series <- obs_binary[, i]
  valid_obs_count <- sum(!is.na(obs_series)); valid_pred_dist_count <- sum(!is.na(pred_series_for_distrib))
  valid_prob_auc_count <- sum(!is.na(prob_series_for_auc)); valid_exp_cor_count <- sum(!is.na(exp_series_for_cor))
  if (valid_obs_count < 30 || valid_pred_dist_count < 30 || valid_prob_auc_count < 30 || valid_exp_cor_count < 30) {
    metrics_results_all[[i]] <- list(station_id = i, lon = station_coords[i, 1], lat = station_coords[i, 2], R01_obs = NA, SDII_obs = NA, P98_obs = NA, DW_obs = NA, WW_obs = NA, R01_pred = NA, SDII_pred = NA, P98_pred = NA, DW_pred = NA, WW_pred = NA, AUC = NA, COR_paper = NA); next
  }
  metrics_results_all[[i]] <- list(
    station_id = i, lon = station_coords[i, 1], lat = station_coords[i, 2],
    R01_obs = calculate_R01(obs_series, wet_threshold), R01_pred = calculate_R01(pred_series_for_distrib, wet_threshold),
    SDII_obs = calculate_SDII(obs_series, wet_threshold), SDII_pred = calculate_SDII(pred_series_for_distrib, wet_threshold),
    P98_obs = calculate_P98(obs_series, wet_threshold), P98_pred = calculate_P98(pred_series_for_distrib, wet_threshold),
    DW_obs = calculate_DW(obs_series, wet_threshold), DW_pred = calculate_DW(pred_series_for_distrib, wet_threshold),
    WW_obs = calculate_WW(obs_series, wet_threshold), WW_pred = calculate_WW(pred_series_for_distrib, wet_threshold),
    AUC = calculate_AUC(prob_series_for_auc, obs_bin_series), COR_paper = calculate_COR_paper(obs_series, exp_series_for_cor)
  )
}
metrics_df_all_raw <- bind_rows(metrics_results_all)
message("Le calcul de l'indicateur est terminé")

5. 计算相对偏差并整理最终评估数据框

calculate_relative_bias <- function(pred_val, obs_val) {
  bias <- ifelse(is.na(obs_val) | is.na(pred_val) | obs_val == 0, NA, 100 * (pred_val - obs_val) / obs_val); return(bias)
}
metrics_df_eval <- metrics_df_all_raw %>%
  mutate(R01_Bias = calculate_relative_bias(R01_pred, R01_obs), SDII_Bias = calculate_relative_bias(SDII_pred, SDII_obs), P98_Bias = calculate_relative_bias(P98_pred, P98_obs), DW_Bias = calculate_relative_bias(DW_pred, DW_obs), WW_Bias = calculate_relative_bias(WW_pred, WW_obs)) %>%
  select(station_id, lon, lat, AUC, COR_paper, R01_Bias, SDII_Bias, P98_Bias, DW_Bias, WW_Bias)

6. 绘制评估指标的空间分布图

使用自定义的绘图函数 plotMetricMap_fixed，将 metrics_df_eval 中的各项评估指标值绘制在地图上，直观展示模型性能的空间差异，并将所有图保存到一个PDF文件中。

plotMetricMap_fixed <- function(metric.values, station_coords, title, colorscale = "PuBu", reverse = FALSE,
                                xlim = c(-14, 44), ylim = c(36, 74), zlim = NULL, legend_text = "Value") {
  valid_idx <- !is.na(metric.values); if (sum(valid_idx) == 0) { warning(paste("Il n'y a pas de valeurs non-NA valides pour dessiner la carte:", title)); plot(1, type="n", xlim=xlim, ylim=ylim, xlab="Longitude", ylab="Latitude", main=paste(title, "\n(No valid data)")); map("world", add=TRUE, fill=TRUE, col="gray95", border="gray70", lwd=0.2); return() }
  metric.values <- metric.values[valid_idx]; station_coords_plot <- station_coords[valid_idx, , drop = FALSE]
  palette <- brewer.pal(9, colorscale); if (reverse) palette <- rev(palette)
  if (is.null(zlim)) { zlim <- range(metric.values, na.rm = TRUE) }
  if (diff(zlim) == 0) { offset <- if (zlim[1] != 0) abs(zlim[1] * 0.1) else 0.1; zlim <- c(zlim[1] - offset, zlim[2] + offset); if(diff(zlim)==0) zlim <- c(zlim[1]-1, zlim[2]+1) }
  par(mar = c(4, 4, 3, 5.5)); plot(station_coords_plot[, 1], station_coords_plot[, 2], xlim = xlim, ylim = ylim, type = "n", main = title, xlab = "Longitude", ylab = "Latitude")
  map("world", add = TRUE, fill = TRUE, col = "gray95", border = "gray70", lwd = 0.2)
  color_fun <- colorRampPalette(palette); n_colors_smooth <- 200; color_bins <- seq(zlim[1], zlim[2], length.out = n_colors_smooth + 1)
  color_map <- color_fun(n_colors_smooth); metric.values_clamped <- pmax(zlim[1], pmin(zlim[2], metric.values))
  color_idx <- findInterval(metric.values_clamped, color_bins, rightmost.closed = TRUE, all.inside = TRUE); color_idx[color_idx == 0] <- 1
  point_colors <- color_map[color_idx]
  points(station_coords_plot[, 1], station_coords_plot[, 2], pch = 21, bg = point_colors, col = "black", cex = 1.2)
  image.plot(legend.only = TRUE, zlim = zlim, col = color_map, axis.args = list(cex.axis = 0.8), legend.args = list(text = legend_text, side = 4, line = 3, cex = 0.9), legend.shrink = 0.8, legend.width = 1.2)
}

# Exporter toutes les cartes d'indicateurs au format PDF
message("En train de générer la carte des indicateurs d'évaluation...")

pdf("regional_pca_all_metrics_maps_separate_pages.pdf", width = 8, height = 7)

coords_matrix <- as.matrix(metrics_df_eval[, c("lon", "lat")])

# --- Tracer le graphique AUC  ---
auc_colors <- "YlGnBu"
auc_zlim <- c(min(metrics_df_eval$AUC, 0.5, na.rm = TRUE), max(metrics_df_eval$AUC, 1.0, na.rm = TRUE))
# 处理 zlim 范围为0的情况
if(diff(auc_zlim) < 1e-6 && !is.na(diff(auc_zlim))) auc_zlim <- auc_zlim + c(-0.05, 0.05) # 确保不为NA
plotMetricMap_fixed(metrics_df_eval$AUC, coords_matrix, "GLM RegionalPCA - AUC (Param. Sim.)",
                    colorscale = auc_colors, reverse = FALSE, zlim = auc_zlim, legend_text = "AUC")

# --- Tracer le graphique COR  ---
cor_colors <- "RdYlBu"
cor_values_no_na <- metrics_df_eval$COR_paper[!is.na(metrics_df_eval$COR_paper)] 
if (length(cor_values_no_na) > 0) {
  cor_limit <- max(abs(range(cor_values_no_na)))
  cor_limit <- min(1, ceiling(cor_limit * 10) / 10)
  cor_zlim <- c(-cor_limit, cor_limit)
  if(diff(cor_zlim) < 1e-6 && !is.na(diff(cor_zlim))) cor_zlim <- cor_zlim + c(-0.1, 0.1)
} else {
  cor_zlim <- c(-1, 1) 
}
plotMetricMap_fixed(metrics_df_eval$COR_paper, coords_matrix, "GLM RegionalPCA - Spearman COR (Param. Sim.)",
                    colorscale = cor_colors, reverse = FALSE, zlim = cor_zlim, legend_text = "Spearman Corr.")

# --- Créer un graphique des écarts relatifs ---
bias_colors <- "BrBG"
bias_bias_values_no_na <- unlist(metrics_df_eval[, grep("_Bias", names(metrics_df_eval))])
bias_bias_values_no_na <- bias_bias_values_no_na[!is.na(bias_bias_values_no_na)]
if (length(bias_bias_values_no_na) > 0) {
  bias_limit_abs <- quantile(abs(bias_bias_values_no_na), 0.98, na.rm=TRUE)
  bias_limit <- ceiling(bias_limit_abs / 10) * 10
  bias_limit <- max(bias_limit, 10, na.rm = TRUE)  
  bias_zlim <- c(-bias_limit, bias_limit)
  if(diff(bias_zlim) < 1e-6 && !is.na(diff(bias_zlim))) bias_zlim <- bias_zlim + c(-5, 5) # 处理范围为0
} else {
  bias_zlim <- c(-50, 50)  
}

# R01 Bias
plotMetricMap_fixed(metrics_df_eval$R01_Bias, coords_matrix, "GLM RegionalPCA - R01 Rel. Bias (%)",
                    colorscale = bias_colors, reverse = FALSE, zlim = bias_zlim, legend_text = "Relative Bias (%)")
# SDII Bias 
plotMetricMap_fixed(metrics_df_eval$SDII_Bias, coords_matrix, "GLM RegionalPCA - SDII Rel. Bias (%)",
                    colorscale = bias_colors, reverse = FALSE, zlim = bias_zlim, legend_text = "Relative Bias (%)")
# P98 Bias 
plotMetricMap_fixed(metrics_df_eval$P98_Bias, coords_matrix, "GLM RegionalPCA - P98 Rel. Bias (%)",
                    colorscale = bias_colors, reverse = FALSE, zlim = bias_zlim, legend_text = "Relative Bias (%)")
# WW Bias
plotMetricMap_fixed(metrics_df_eval$WW_Bias, coords_matrix, "GLM RegionalPCA - WW Rel. Bias (%)",
                    colorscale = bias_colors, reverse = FALSE, zlim = bias_zlim, legend_text = "Relative Bias (%)")
# DW Bias 
plotMetricMap_fixed(metrics_df_eval$DW_Bias, coords_matrix, "GLM RegionalPCA - DW Rel. Bias (%)",
                    colorscale = bias_colors, reverse = FALSE, zlim = bias_zlim, legend_text = "Relative Bias (%)")


dev.off()
message("Toutes les cartes d'indicateurs ont été enregistrées dans: regional_pca_all_metrics_maps_separate_pages.pdf")

7. 输出指标的汇总统计信息

将各项关键评估指标的描述性统计（如最小值、中位数、均值、最大值等）输出


cat("\n=== Summary RegionalPCA (Parameterized Simulations): AUC ===\n"); summary(metrics_df_eval$AUC)
cat("\n=== Summary RegionalPCA (Parameterized Simulations): Spearman COR (Paper Definition) ===\n"); summary(metrics_df_eval$COR_paper)
cat("\n=== Summary RegionalPCA (Parameterized Simulations): R01 Relative Bias (%) ===\n"); summary(metrics_df_eval$R01_Bias)
cat("\n=== Summary RegionalPCA (Parameterized Simulations): SDII Relative Bias (%) ===\n"); summary(metrics_df_eval$SDII_Bias)
cat("\n=== Summary RegionalPCA (Parameterized Simulations): P98 Relative Bias (%) ===\n"); summary(metrics_df_eval$P98_Bias)
cat("\n=== Summary RegionalPCA (Parameterized Simulations): WW Relative Bias (%) ===\n"); summary(metrics_df_eval$WW_Bias)
cat("\n=== Summary RegionalPCA (Parameterized Simulations): DW Relative Bias (%) ===\n"); summary(metrics_df_eval$DW_Bias)

message("Fini")

交叉验证

在上述代码流程基础上，加入交叉验证

# ############################################################################
# ## Version régionale de l'ACP avec validation croisée : Reproduction de la section 3.1 de l'article
# ## Traitement des données + Obtention des paramètres de distribution + Simulation aléatoire des paramètres et moyennage
# ############################################################################

# ==================================================================================
# ## 1) Définition du répertoire de travail & Chargement des paquets
# ==================================================================================

# Définir le répertoire de travail, veuillez modifier selon votre chemin réel
setwd("~/Downloads/GCM_regoion_PCA_complet/")

# Chargement des paquets R requis
library(transformeR)
library(downscaleR)
library(magrittr)
library(abind)
library(pROC)
library(ggplot2)
library(dplyr)
library(tidyr)
library(maps)
library(fields)
library(RColorBrewer)
library(geosphere)

# ==================================================================================
# ## 2) Chargement des scripts (glimpr modifiée et traitement des données modifiées)
# ==================================================================================

# Chargement des fichiers de script R
source("glimpr.R")       # Chargement de la version modifiée de la fonction glimpr, potentiellement pour le downscaling des modèles linéaires généralisés, peut retourner probabilités, intensités conditionnelles et modèle.         s
source("prepareData.R")  # Chargement de la version ACP régionale de la fonction prepareData, pour préparer les données nécessaires à l'analyse ACP

# Note : Les fonctions auxiliaires predictor.nn.indices et predictor.nn.values sont définies dans prepareData.R

source("prepareNewData.R")# Chargement de la version ACP régionale de la fonction prepareNewData, pour projeter de nouvelles données dans l'espace ACP existant
source("helpers.R")       # Chargement d'autres fonctions


cv_folds <- 5

# --- Définition de la fonction de mappage des régions PRUDENCE ---
# Liste des étendues des régions européennes définies par le projet PRUDENCE
# Chaque région contient un nom (name), une plage de longitudes (lonLim) et une plage de latitudes (latLim)
prudence_regions <- list(
  BI = list(name = "British Isles", lonLim = c(-10, 2), latLim = c(50, 59)),
  IP = list(name = "Iberian Peninsula", lonLim = c(-10, 3), latLim = c(36, 44)),
  FR = list(name = "France", lonLim = c(-5, 5), latLim = c(44, 50)),
  ME = list(name = "Mid-Europe", lonLim = c(2, 16), latLim = c(48, 55)),
  SC = list(name = "Scandinavia", lonLim = c(5, 30), latLim = c(55, 70)),
  AL = list(name = "Alps", lonLim = c(5, 16), latLim = c(44, 48)),
  MD = list(name = "Mediterranean", lonLim = c(3, 25), latLim = c(36, 44)),
  EA = list(name = "Eastern Europe", lonLim = c(16, 30), latLim = c(44, 55))
)

# Fonction : Mapper les coordonnées d'une station à une région PRUDENCE
# Entrée : lon (longitude), lat (latitude), regions_definition (liste des définitions de régions)
# Sortie : Code de la région à laquelle appartient la station (ex: "BI", "IP", etc.)
map_station_to_prudence_region <- function(lon, lat, regions_definition) {
  region_name_inside <- NA # Initialisation du nom de la région où se trouve la station
  # Parcourir toutes les régions prédéfinies
  for (region_code in names(regions_definition)){
    bounds <- regions_definition[[region_code]] # Obtention des limites de la région actuelle
    # Vérifier si la station se trouve dans l'emprise rectangulaire de la région actuelle
    if (lon >= bounds$lonLim[1] && lon < bounds$lonLim[2] &&
        lat >= bounds$latLim[1] && lat < bounds$latLim[2] ) {
      region_name_inside <- region_code # Si oui, enregistrer le code de la région
      return(region_name_inside) # Retourner dès que trouvé, une station n'appartient qu'à une seule région rectangulaire prédéfinie
    }
  }

  # Si la station n'est dans aucune région rectangulaire prédéfinie, l'assigner au centre de la région la plus proche
  min_dist <- Inf; closest_region <- NA; point_coord <- c(lon, lat) # Initialisation de la distance minimale et de la région la plus proche
  # Parcourir toutes les régions prédéfinies
  for (region_code in names(regions_definition)){
    bounds <- regions_definition[[region_code]]; # Obtention des limites de la région actuelle
    center_lon <- mean(bounds$lonLim); center_lat <- mean(bounds$latLim) # Calcul de la longitude et latitude du centre de la région
    region_center_coord <- c(center_lon, center_lat) # Coordonnées du centre de la région
    # Utilisation de la fonction distHaversine du paquet geosphere pour calculer la distance sphérique (en mètres) entre la station et le centre de la région
    current_dist <- geosphere::distHaversine(point_coord, region_center_coord)
    # Si la distance actuelle est inférieure à la distance minimale enregistrée
    if (current_dist < min_dist){
      min_dist <- current_dist # Mettre à jour la distance minimale
      closest_region <- region_code # Mettre à jour le code de la région la plus proche
    }
  }
  return(closest_region) # Retourner le code de la région la plus proche
}

# ==================================================================================
# ## 3) Chargement des données & Nettoyage des stations invalides
# ==================================================================================

# Chargement des données de précipitations observées (données de stations)
load("st86_precip.rda")
# Chargement des données de réanalyse ERA-Interim (données grillées, utilisées comme prédicteurs)
load("era_interim.rda")

orig_era_dates <- era_interim$Dates
# ───────────────────────────────────────────────────────────────────────────────

# --- Alignement temporel manuel (basé sur les objets Date, en ignorant les différences d'heure, minute, seconde) ---
message("Alignement temporel manuel de st86_precip et era_interim (basé sur la partie date)...")

# 1. Extraction des chaînes de dates originales de st86_precip ($start / $end)
if (!is.list(st86_precip$Dates) || !("start" %in% names(st86_precip$Dates))) {
  stop("La structure de st86_precip$Dates n'est pas conforme aux attentes, ce n'est pas une liste contenant $start.")
}
st86_orig_char_dates_start <- st86_precip$Dates$start
st86_orig_char_dates_end   <- st86_precip$Dates$end

# 2. Extraction des chaînes de dates originales de era_interim (du $start du premier membre)
if (!is.list(era_interim$Dates) ||
    length(era_interim$Dates) == 0 ||
    !is.list(era_interim$Dates[[1]]) ||
    !("start" %in% names(era_interim$Dates[[1]]))) {
  message("Structure de era_interim$Dates :")
  print(str(era_interim$Dates))
  stop("La structure de era_interim$Dates n'est pas conforme aux attentes, ou le premier élément ne contient pas $start.")
}
era_orig_char_dates_start <- era_interim$Dates[[1]]$start
era_orig_char_dates_end   <- era_interim$Dates[[1]]$end   # Conservation facultative

# 3. Conversion en Date (ignorer heures, minutes, secondes)
st86_Date_obj_for_match  <- as.Date(st86_orig_char_dates_start)
era_Date_obj_for_match   <- as.Date(era_orig_char_dates_start)

message("Longueur de st86_Date_obj_for_match après conversion : ", length(st86_Date_obj_for_match))
message("Longueur de era_Date_obj_for_match après conversion : ", length(era_Date_obj_for_match))

# 4. Trouver les objets Date communs
common_days_Date_obj <- intersect(st86_Date_obj_for_match,
                                  era_Date_obj_for_match)
if (length(common_days_Date_obj) == 0) {
  stop("Échec de l'alignement des dates basé sur les objets Date, aucune date commune ! Veuillez vérifier la conversion as.Date.")
}
message("Nombre d'objets Date communs trouvés : ", length(common_days_Date_obj))

# Forcer la reconversion en Date du vecteur potentiellement numérique issu de intersect
common_days_Date_obj <- as.Date(common_days_Date_obj, origin = "1970-01-01")
# Puis trier et dédupliquer
common_days_Date_obj <- sort(unique(common_days_Date_obj))
message("Classe de common_days_Date_obj après tri et déduplication : ", class(common_days_Date_obj))
message("5 premiers common_days_Date_obj : "); print(head(common_days_Date_obj, 5))

# 5. Faire correspondre aux séquences originales en utilisant les objets Date communs
st86_indices <- match(common_days_Date_obj, st86_Date_obj_for_match)
era_indices  <- match(common_days_Date_obj, era_Date_obj_for_match)

# 6. Vérifier la sortie de match()
message("--- st86_indices (après match) ---")
message("Longueur : ", length(st86_indices))
message("Présence de NA : ", anyNA(st86_indices))
if (anyNA(st86_indices)) {
  message("Nombre de NA : ", sum(is.na(st86_indices)))
  message("Valeur de common_days_Date_obj correspondant au premier NA : ")
  print(common_days_Date_obj[which(is.na(st86_indices))[1]])
}
message("5 premières valeurs d'indice :"); print(head(st86_indices, 5))
message("5 dernières valeurs d'indice :"); print(tail(st86_indices, 5))

message("--- era_indices (après match) ---")
message("Longueur : ", length(era_indices))
message("Présence de NA : ", anyNA(era_indices))
if (anyNA(era_indices)) {
  message("Nombre de NA : ", sum(is.na(era_indices)))
  message("Valeur de common_days_Date_obj correspondant au premier NA : ")
  print(common_days_Date_obj[which(is.na(era_indices))[1]])
}
message("5 premières valeurs d'indice :"); print(head(era_indices, 5))
message("5 dernières valeurs d'indice :"); print(tail(era_indices, 5))

if (length(st86_indices) != length(common_days_Date_obj) ||
    length(era_indices)  != length(common_days_Date_obj)) {
  stop("Problème lors de la mise en correspondance des dates communes avec les séquences originales, longueurs d'indices incohérentes.")
}

# 7. Filtrage et remplacement de st86_precip
st86_precip$Data <- st86_precip$Data[st86_indices, , drop = FALSE]
unified_start_strings <- format(common_days_Date_obj, "%Y-%m-%d 00:00:00")
unified_end_strings_st86 <- format(
  as.Date(st86_orig_char_dates_end[st86_indices]),
  "%Y-%m-%d 00:00:00"
)
st86_precip$Dates$start <- unified_start_strings
st86_precip$Dates$end   <- unified_end_strings_st86


# --- 8. Sous-ensemble manuel de era_interim$Data (identification dynamique de la dimension temporelle) ---
original_era_dims  <- dim(era_interim$Data)
era_dims_names     <- attr(era_interim$Data, "dimensions")
time_len           <- length(era_Date_obj_for_match)

# 1) Trouver quelle dimension est le temps
time_dim <- which(original_era_dims == time_len)
if (length(time_dim) != 1) {
  stop("Impossible de localiser de manière unique la dimension temporelle : dim = ", paste(original_era_dims, collapse="×"))
}

# 2) Construire la liste d'indices : utiliser era_indices uniquement pour la dimension temporelle, ":" pour les autres dimensions
subset_args <- lapply(seq_along(original_era_dims), function(d) {
  if (d == time_dim) era_indices else substitute()
})

# 3) Réalisation effective du sous-ensemble (drop=FALSE pour conserver toutes les dimensions)
era_interim$Data <- do.call(
  `[`,
  c(list(era_interim$Data), subset_args, list(drop = FALSE))
)

# 4) Restaurer l'attribut dimensions
attr(era_interim$Data, "dimensions") <- era_dims_names

# 5) Mettre à jour Dates pour qu'elles correspondent exactement à celles de st86_precip
unified_end_strings_era <- format(
  as.Date(era_orig_char_dates_end[era_indices]),
  "%Y-%m-%d 00:00:00"
)
era_interim$Dates <- list(
  start = unified_start_strings,
  end   = unified_end_strings_era
)



# 9. Vérification finale
if (!identical(st86_precip$Dates$start, era_interim$Dates$start)) {
  message("st86_precip$Dates$start (5 premiers) :"); print(head(st86_precip$Dates$start))
  message("era_interim$Dates$start (5 premiers) :"); print(head(era_interim$Dates$start))
  stop("Après l'alignement temporel manuel, les chaînes $Dates$start ne correspondent toujours pas !")
}

message("Alignement temporel réussi — Nombre de dates valides communes : ", length(common_days_Date_obj), " jours")
# —————— Recalcul des jours et des indices nécessaires pour la validation croisée ——————
# Utiliser la longueur de st86_precip$Dates$start juste après l'alignement
total_days <- length(st86_precip$Dates$start)
fold_size  <- floor(total_days / cv_folds)

cv_indices <- vector("list", cv_folds)
for (i in seq_len(cv_folds)) {
  start_idx <- (i - 1) * fold_size + 1
  end_idx   <- i * fold_size
  # (Si vous voulez absolument garantir de ne pas dépasser, vous pouvez aussi écrire end_idx <- min(i*fold_size, total_days))
  cv_indices[[i]] <- seq(start_idx, end_idx)
}
# —————— Fin de alignement temporel ——————

message("Longueur de st86_precip$Dates$start après alignement : ", length(st86_precip$Dates$start))
message("Longueur de era_interim$Dates$start après alignement : ", length(era_interim$Dates$start))

# --- Nettoyage des stations invalides ---
# Définition de la liste des ID de stations à supprimer
remove_ids <- c("000017", "000274", "000275")
# Sauvegarde des coordonnées et ID originaux des stations (nul)
st86_precip_original_coords <- st86_precip$xyCoords
st86_precip_original_ids <- st86_precip$Metadata$station_id

# Trouver les indices des stations à supprimer dans les données
idx_remove <- which(st86_precip$Metadata$station_id %in% remove_ids)

# Si des stations à supprimer ont été trouvées
if (length(idx_remove) > 0) {
  # Supprimer les colonnes de données de ces stations de la matrice de données, drop=FALSE garantit que le résultat reste une matrice même s'il ne reste qu'une colonne
  st86_precip$Data <- st86_precip$Data[, -idx_remove, drop = FALSE]
  # Supprimer les lignes de coordonnées de ces stations de la matrice de coordonnées, drop=FALSE garantit que le résultat reste une matrice même s'il ne reste qu'une ligne
  st86_precip$xyCoords <- st86_precip$xyCoords[-idx_remove, , drop = FALSE]
  # Supprimer les ID de ces stations des métadonnées
  st86_precip$Metadata$station_id <- st86_precip$Metadata$station_id[-idx_remove]
  # Supprimer les noms de ces stations des métadonnées
  st86_precip$Metadata$name <- st86_precip$Metadata$name[-idx_remove]
  # Supprimer les altitudes de ces stations des métadonnées
  st86_precip$Metadata$altitude <- st86_precip$Metadata$altitude[-idx_remove]
  # Supprimer les informations de source de ces stations des métadonnées
  st86_precip$Metadata$source <- st86_precip$Metadata$source[-idx_remove]
}

# S'assurer que la matrice de données st86_precip$Data possède l'attribut "dimensions"
# Ceci est requis pour certaines fonctions des paquets transformeR et downscaleR
if (is.null(attr(st86_precip$Data, "dimensions"))) {
  attr(st86_precip$Data, "dimensions") <- c("time", "loc") # "time" représente la dimension temporelle, "loc" la dimension station/localisation
}
message("Chargement et nettoyage des données terminés. Nombre de stations après nettoyage :", ncol(st86_precip$Data)) # Affichage du nombre de stations après nettoyage

# --- Obtention des coordonnées et ID des stations finales utilisées (après nettoyage) ---
final_station_coords <- as.matrix(st86_precip$xyCoords) # S'assurer que les coordonnées des stations sont sous forme de matrice
final_station_ids <- st86_precip$Metadata$station_id   # Obtention des ID des stations finales utilisées
n_final_stations <- nrow(final_station_coords)         # Obtention du nombre de stations finales utilisées

# ==================================================================================
# ## 4) Définition des paramètres de validation croisée
# ==================================================================================

cv_folds <- 5  # Définition du nombre de plis pour la validation croisée (K-fold cross-validation)
set.seed(123)  # Définition de la graine aléatoire pour assurer la reproductibilité des résultats

# Obtention de la plage temporelle des données d'observation (c.-à-d. la séquence de dates)
# Note : Si st86_precip$Dates est une liste (contenant $start et $end),
# observation_dates sera cette liste elle-même, et non un vecteur de dates.
# length(observation_dates) sera alors la longueur de la liste (généralement 2).
observation_dates <- st86_precip$Dates


# Liste pour stocker tous les résultats de la validation croisée
cv_results <- list()

这部分没有什么好说的，设置工作环境并加载包，PRUDENCE区域定义与站点映射函数，数据加载与无效站点清理，特别是时间对齐，这里卡了好几天，原本使用函数对齐的方式没做出来，一直报错，所以用的手动对齐方法。大体流程如下

首先提取 st86_precip 的原始日期字符串，检查 st86_precip$Dates 是否是一个列表并且包含名为 start 的元素，因为我们的数据结构就是这样，然后，将start（开始日期字符串向量）和end（结束日期字符串向量）分别赋值给 st86_orig_char_dates_start 和 st86_orig_char_dates_end。

同样提取 era_interim 的原始日期字符串，检查 era_interim$Dates[[1]]$start 是否存在，然后提取其开始和结束日期字符串。

将日期字符串转换为 Date 对象

st86_Date_obj_for_match  <- as.Date(st86_orig_char_dates_start) # 将观测日期的开始字符串转为Date对象
era_Date_obj_for_match   <- as.Date(era_orig_char_dates_start)  # 将再分析日期的开始字符串转为Date对象
message("Longueur de st86_Date_obj_for_match après conversion : ", length(st86_Date_obj_for_match)) # "转换后 st86_Date_obj_for_match 的长度："
message("Longueur de era_Date_obj_for_match après conversion : ", length(era_Date_obj_for_match)) # "转换后 era_Date_obj_for_match 的长度："

as.Date() 会自动解析多种常见的日期格式，并只保留年月日信息，忽略时间部分（时、分、秒），之所以忽略时、分、秒，是因为两个数据时间小时级别不对齐，x 是00:00开始 y 是12:00开始

转换后的 Date 对象存储在 st86_Date_obj_for_match 和 era_Date_obj_for_match 中。

然后找到两个数据集中共有的日期

intersect() 函数找出在 st86_Date_obj_for_match 和 era_Date_obj_for_match 两个 Date 对象向量中都存在的日期，结果存入 common_days_Date_obj。如果没有任何共同日期，脚本自动停止报错，最后要确保 common_days_Date_obj 是 Date 类型，并对其进行排序和去重，以得到一个干净、有序的共同日期序列。

获取共同日期在原始日期序列中的索引

st86_indices <- match(common_days_Date_obj, st86_Date_obj_for_match) # 获取共同日期在观测数据原始Date序列中的索引
era_indices  <- match(common_days_Date_obj, era_Date_obj_for_match)  # 获取共同日期在再分析数据原始Date序列中的索引

match() 函数返回 common_days_Date_obj 中的每个日期在 st86_Date_obj_for_match (或 era_Date_obj_for_match) 中首次出现的位置（索引）。这些索引 (st86_indices 和 era_indices) 将用于从原始数据中精确提取与共同日期对应的数据行/时间片。

然后调试信息不用多说，检查 st86_indices 和 era_indices 的长度是否与 common_days_Date_obj 的长度一致，以及是否包含 NA 值（理论上不应该有，因为是基于 intersect 的结果来匹配的）

根据索引过滤 st86_precip (站点观测数据)

st86_precip$Data (观测数据矩阵) 通过行索引 st86_indices 进行子集筛选，只保留与共同日期对应的时间点。drop = FALSE 确保结果仍为矩阵。

原始的结束日期字符串 st86_orig_char_dates_end 也根据 st86_indices 筛选，转换为 Date 对象后再格式化为 unified_end_strings_st86。

最后，st86_precip$Dates$start 和 st86_precip$Dates$end 被更新为这些统一和筛选后的日期字符串。

st86_precip$Dates$start <- unified_start_strings # 更新观测数据的开始日期信息
st86_precip$Dates$end   <- unified_end_strings_st86 # 更新观测数据的结束日期信息

根据索引过滤 era_interim$Data (再分析数据，多维数组):

这部分比较复杂，因为 era_interim$Data 是一个多维数组，时间维度不固定。大体来说，代码先确定 era_interim$Data 中哪一维是时间维，然后，构造一个参数列表 subset_args。这个列表用于 do.call 函数，以便在正确的时间维度上使用 era_indices 进行索引，而其他维度则保持不变（通过 substitute() 实现，它在 do.call 上下文中表示一个空参数，代表取该维度的全部）。子集操作后，恢复可能丢失的 dimensions 属性，并像处理 st86_precip 一样，更新 era_interim$Dates使其与共同日期对齐。

最终的日期一致性检查:再次严格比较对齐后两个数据对象的start 字符串是否完全一致。如果不一致则报错。最后展示输出时间对齐成功的消息和共同有效日期的数量。

下面是交叉验证核心流程部分

外层循环：遍历每一折 (fold)，初始化当前折的结果存储容器

# ==================================================================================
#  交叉验证主循环 - 执行k折交叉验证评估模型性能
# ==================================================================================
message("开始执行交叉验证，", cv_folds, " 折交叉验证...")

for (fold in seq_len(cv_folds)) {
  message("\n处理第 ", fold, "/", cv_folds, " 折交叉验证")

  # 初始化当前折的模型存储容器，每个站点对应一个Gamma模型
  gamma_models_all_stations <- vector("list", n_final_stations)

  # 创建当前折的结果容器
  cv_results[[fold]] <- list()

为当前这一折 fold 初始化一个名为 gamma_models_all_stations 的列表。这个列表的长度等于最终有效站点的数量 n_final_stations。它的目的是在后续逐站点调用 glimpr 函数后，存储每个站点拟合得到的Gamma模型对象，这些模型将在蒙特卡洛模拟步骤中用到。

在之前创建的总结果列表 cv_results 中，为当前折 fold 创建一个新的空子列表。这个子列表 cv_results[[fold]] 将用于存储与当前这一折相关的所有计算结果和元数据。

数据划分 - 将数据集划分为训练集和测试集:

  # ==================================================================================
  #  数据划分 - 将数据集划分为训练集和测试集
  # ==================================================================================
  # 获取当前折的测试集时间索引
  test_idx  <- cv_indices[[fold]]
  # 训练集索引为除测试集外的所有时间点
  train_idx <- setdiff(seq_along(common_days_Date_obj), test_idx)

test_idx <- cv_indices[[fold]]: 从之前在数据预处理阶段（时间对齐后）创建的 cv_indices 列表中，提取出第 fold 个元素。这个元素本身是一个整数向量，包含了属于当前这一折 测试集 的所有时间点在完整对齐后数据序列 (common_days_Date_obj) 中的索引。

train_idx <- setdiff(seq_along(common_days_Date_obj), test_idx): 确定当前这一折的 训练集 索引。seq_along(common_days_Date_obj) 生成一个从1到总天数的完整索引序列。setdiff() 函数从这个完整序列中移除掉 test_idx（测试集索引），剩下的索引就构成了训练集的时间索引 train_idx。

然后根据 train_idx 和 test_idx 将站点观测数据分割为训练集 (train_obs) 和测试集 (test_obs)

  # ------------------------------------------------------------------------------
  #  划分站点观测数据
  # ------------------------------------------------------------------------------
  # 复制完整数据结构
  train_obs <- st86_precip
  test_obs  <- st86_precip

  # 提取训练集时间段的观测数据
  train_obs$Data  <- st86_precip$Data[train_idx, , drop = FALSE]
  attr(train_obs$Data, "dimensions") <- c("time", "loc")

  # 提取测试集时间段的观测数据
  test_obs$Data   <- st86_precip$Data[test_idx,  , drop = FALSE]
  attr(test_obs$Data, "dimensions") <- c("time", "loc")

  # 更新日期信息
  train_obs$Dates <- list(
    start = unified_start_strings[train_idx],
    end   = unified_end_strings_st86[train_idx]
  )
  test_obs$Dates  <- list(
    start = unified_start_strings[test_idx],
    end   = unified_end_strings_st86[test_idx]
  )

首先，通过直接赋值创建 train_obs 和 test_obs 作为 st86_precip 的完整副本。这样做是为了保留原始对象的完整结构（如坐标、元数据等），然后只替换其数据和日期部分。然后，train_obs$Data 被赋予 st86_precip$Data 中由 train_idx 索引的行（时间维度）。drop = FALSE 确保即使只有一个站点或一个时间点，结果仍保持矩阵格式。

attr(train_obs$Data, "dimensions") <- c("time", "loc") 重新为子集化后的数据矩阵设置 "dimensions" 属性，指明行是时间，列是站点/位置。对 test_obs$Data 也这么做，不加就报错

划分预测因子数据 (ERA-Interim再分析数据):

同样，根据 train_idx 和 test_idx 将预测因子数据 (era_interim) 分割为训练集 (train_pred) 和测试集 (test_pred)。

  # ------------------------------------------------------------------------------
  #  划分预测因子数据（ERA-Interim再分析数据）
  # ------------------------------------------------------------------------------
  train_pred <- era_interim
  test_pred  <- era_interim

  # 提取对应时间段的ERA数据（注意维度顺序）
  train_pred$Data <- era_interim$Data[ , , train_idx, , , drop = FALSE]
  train_pred$Dates <- lapply(orig_era_dates, function(d) {
    list(
      start = d$start[train_idx],
      end   = d$end  [train_idx]
    )
  })
  attr(train_pred$Data, "dimensions") <- era_dims_names

  test_pred$Data  <- era_interim$Data[ , , test_idx,  , , drop = FALSE]
  test_pred$Dates  <- lapply(orig_era_dates, function(d) {
    list(
      start = d$start[test_idx],
      end   = d$end  [test_idx]
    )
  })
  attr(test_pred$Data, "dimensions") <- era_dims_names

  message("训练集 (", length(train_idx), " 天) 和测试集 (", length(test_idx), " 天) 已划分")

区域主成分分析 (Regional PCA) - 对训练数据进行PCA降维:

PCA仅在当前折的训练集预测因子 (train_pred) 上进行，然后将当前折的测试集预测因子 (test_pred) 投影到相同的PCA空间。

  # ==================================================================================
  #  区域主成分分析 - 对训练数据进行PCA降维
  # ==================================================================================
  # 使用prepareData函数进行区域PCA，保留95%的方差
  data.pca <- prepareData(
    x                = train_pred,      # 预测因子（ERA数据）
    y                = train_obs,       # 预报对象（站点观测）
    spatial.predictors = list(
      v.exp           = 0.95,           # 方差解释率阈值
      which.combine   = getVarNames(train_pred)  # 需要组合的变量
    ),
    combined.only    = TRUE             # 只使用组合后的主成分
  )

  # 将测试集数据投影到训练集的PCA空间
  message("将测试数据投影到PCA空间...")
  newdata.pca <- prepareNewData(test_pred, data.pca)

  # 保存当前折的时间范围信息
  cv_results[[fold]]$train_dates <- range(train_obs$Dates$start)
  cv_results[[fold]]$test_dates <- range(test_obs$Dates$start)

将当前折 fold 的训练集 (train_obs$Dates$start) 和测试集 (test_obs$Dates$start) 的日期范围（通过 range()函数获得，即最早和最晚日期）存储到总结果列表 cv_results 的当前折子列表cv_results[[fold]] 中，字段名为 train_dates 和 test_dates。

初始化用于存储当前折预测结果的矩阵

为当前折的测试期 (n_test_days 天) 和所有站点 (n_final_stations 个) 初始化三个矩阵，用 NA_real_ (实数型NA) 填充，用于后续存储GLM的预测输出。

  # ==================================================================================
  #  初始化预测结果存储矩阵
  # ==================================================================================
  n_test_days <- length(test_idx)
  # 降水概率估计
  fold_prob_estimates <- matrix(NA_real_, nrow = n_test_days, ncol = n_final_stations)
  # 条件降水强度（给定降水发生）
  fold_cond_mean_intensity <- matrix(NA_real_, nrow = n_test_days, ncol = n_final_stations)
  # 最终降水预测值
  fold_predictions <- matrix(NA_real_, nrow = n_test_days, ncol = n_final_stations)

n_test_days 是当前测试集包含的天数

  • fold_prob_estimates: 用于存储每个测试日、每个站点的降水发生概率。

  • fold_cond_mean_intensity: 用于存储每个测试日、每个站点的条件平均降水强度

  • fold_predictions: 用于存储每个测试日、每个站点的最终降水预测值（要么是 概率乘以条件强度，要么是经过随机模拟后的值）。

站点-区域映射

确定每个最终有效的气象站点属于哪个PRUDENCE气候区域。

  # ==================================================================================
  #  站点-区域映射 - 确定每个站点所属的PRUDENCE区域
  # ==================================================================================
  # 获取站点坐标矩阵
  station_coords <- as.matrix(train_obs$xyCoords)
  # 将每个站点映射到对应的PRUDENCE区域
  station_region_mapping <- sapply(1:n_final_stations, function(i) {
    map_station_to_prudence_region(station_coords[i, 1], station_coords[i, 2], prudence_regions)
  })
  names(station_region_mapping) <- train_obs$Metadata$station_id

  # 获取PCA处理后的区域列表
  processed_regions <- names(data.pca$x.global)

对 1 到 n_final_stations 的每个站点索引 i，调用之前定义的 map_station_to_prudence_region 函数，传入该站点的经度 (station_coords[i, 1])、纬度 (station_coords[i, 2]) 和 prudence_regions 定义。sapply 会将每次函数调用的返回结果（区域代码字符串）收集到向量 station_region_mapping 中。

data.pca$x.global 是一个列表，其每个元素的名称对应一个区域代码，元素内容是该区域训练集的主成分得分。最后processed_regions <- names(data.pca$x.global) 的意思是说从 prepareData 函数返回的 data.pca 对象中提取实际进行了PCA处理并生成了主成分的区域列表。

逐站点GLM建模与预测 - glimpr 函数:

  # ==================================================================================
  #  逐站点建模与预测 - 使用glimpr函数
  # ==================================================================================
  for(i in 1:n_final_stations) {
    message("     处理站点 ", i, " (第 ", fold, " 折)")

    # 获取当前站点所属区域
    station_region <- station_region_mapping[i]

    # 检查该区域是否在PCA处理范围内
    if(!(station_region %in% processed_regions)) {
      message("       站点 ", i, " 的区域 ", station_region, " 未被处理，跳过")
      next
    }

为调用 glimpr 函数准备符合其输入格式要求的观测数据 y

    # ------------------------------------------------------------------------------
    #  准备单个站点的训练数据
    # ------------------------------------------------------------------------------
    y_subset_for_station <- list(
      Data = train_obs$Data[, i, drop = FALSE],                    # 该站点的观测数据
      xyCoords = train_obs$xyCoords[i, , drop = FALSE],           # 站点坐标
      Dates = train_obs$Dates,                                     # 日期信息
      Metadata = list(
        station_id = train_obs$Metadata$station_id[i]             # 站点ID
      )
    )

    # 确保数据具有正确的维度属性
    if(is.null(attr(y_subset_for_station$Data, "dimensions"))) {
      attr(y_subset_for_station$Data, "dimensions") <- c("time", "loc")
    }

创建一个名为 y_subset_for_station 的列表，其结构与 glimpr 函数期望的输入参数 y 一致

获取当前站点所属区域的PCA结果

从 data.pca (训练集PCA结果) 和 newdata.pca (测试集在训练集PCA空间上的投影结果) 中提取当前站点 station_region 对应的PCA得分。

    # ------------------------------------------------------------------------------
    #  获取区域PCA结果
    # ------------------------------------------------------------------------------
    # 训练集的主成分得分矩阵
    pred_mat_regional <- data.pca$x.global[[station_region]]
    # 测试集的主成分得分矩阵
    sim_mat_regional <- newdata.pca$x.global[[station_region]][[1]]

    # 获取该区域的主成分数量
    n_pcs_regional <- ncol(pred_mat_regional)
    if(is.null(n_pcs_regional)) {
      warning("无法确定区域 ", station_region, " 的PC数量，跳过站点 ", i)
      next
    }

  • pred_mat_regional: 从训练PCA结果 data.pca$x.global 中提取当前站点所属区域 station_region 对应的元素。这是一个矩阵，行代表训练期的时间，列代表该区域的主成分，值是主成分得分。

  • sim_mat_regional: 从测试投影结果 newdata.pca$x.global 中提取 station_region 对应的元素，并取其第一个子元素 [[1]] ，这是测试集在该区域PCA空间中的主成分得分。

  • n_pcs_regional: 从 pred_mat_regional 的列数确定该区域实际使用的主成分数量。如果 pred_mat_regional 为 NULL (例如该区域PCA失败)，则 n_pcs_regional 也会是 NULL，此时会发出警告并跳过当前站点。

调用 glimpr 函数进行降尺度建模与预测:

    # ------------------------------------------------------------------------------
    #  调用glimpr函数进行降尺度建模
    # ------------------------------------------------------------------------------
    wet.threshold <- 1  # 降水阈值：1mm，用于区分有无降水

    glimpr_station_result <- glimpr(
      y = y_subset_for_station,         # 站点观测数据
      modelPars = list(
        stations = TRUE,                # 标记为站点数据
        multi.member = FALSE,           # 单成员预报
        pred.mat = pred_mat_regional,   # 训练集主成分
        sim.mat = list(sim_mat_regional), # 测试集主成分（列表格式）
        sim.dates = newdata.pca$Dates  # 预报日期
      ),
      simulate = "no",                  # 使用确定性预测模式
      return.models = TRUE,             # 返回拟合的模型对象
      wet.threshold = wet.threshold,    # 降水阈值
      n.pcs = n_pcs_regional           # 使用的主成分数
    )
    # 保存Gamma模型供后续模拟使用
    gamma_models_all_stations[[i]] <- glimpr_station_result$models$gamma[[1]]

    # 存储预测结果
    fold_prob_estimates[, i] <- glimpr_station_result$probabilities      # 降水概率
    fold_cond_mean_intensity[, i] <- glimpr_station_result$cond_intensities # 条件强度
    fold_predictions[, i] <- glimpr_station_result$predictions           # 最终预测值
  }

从 glimpr 函数返回的结果 glimpr_station_result 中，提取出拟合的Gamma模型 ($models$gamma)。因为 glimpr 内部的 mod.gamma 是一个列表（每个元素对应一个case，而这里case只有一个，即当前站点），所以取其第一个元素 [[1]]。这个Gamma模型被存储在 gamma_models_all_stations 列表的第 i 个位置，对应当前站点 i。这些保存的Gamma模型将在后续的蒙特卡洛模拟步骤中使用，用于从Gamma分布中抽样生成降水量。

将 glimpr 返回的预测概率、条件强度和最终预测值，分别存储到之前为当前折初始化的结果矩阵的第 i 列（对应当前站点）

蒙特卡洛模拟 - 生成多个降水场景

在对所有站点完成确定性预测后，接下来对每个站点进行蒙特卡洛模拟，以生成降水量的集合样本。

  # ==================================================================================
  #  蒙特卡洛模拟 - 生成500个降水场景
  # ==================================================================================
  message("\n开始进行500次模拟采样...")

  # 设置模拟次数
  n_simulations <- 500

  # 初始化三维数组存储所有模拟结果
  # 维度说明：[模拟次数, 测试天数, 站点数]
  fold_simulations <- array(NA_real_, dim = c(n_simulations, n_test_days, n_final_stations))

逐站点进行蒙特卡洛模拟 (内部循环)

  # ------------------------------------------------------------------------------
  #  逐站点进行蒙特卡洛模拟
  # ------------------------------------------------------------------------------
  for(i in 1:n_final_stations) {
    message("  站点 ", i, "/", n_final_stations, " - 开始500次采样")

    # 获取站点区域信息
    station_region <- station_region_mapping[i]

    # 检查区域有效性
    if(!(station_region %in% processed_regions)) {
      message("    站点 ", i, " 的区域 ", station_region, " 未被处理，跳过采样")
      next
    }

    # 获取测试集的区域主成分
    sim_mat_regional <- newdata.pca$x.global[[station_region]][[1]]
    n_pcs_regional <- ncol(sim_mat_regional)

获取当前站点所属区域 station_region 在测试集上的主成分得分 sim_mat_regional 和使用的主成分数量 n_pcs_regional。这些将作为模拟过程中模型的输入。

重新训练二元模型（专门用于模拟过程中的降水发生判断）

这里没用glimpr中的伯努利模型，主要原因是一开始写代码总报错，索性直接把glimpr里相关代码复制过来，后续改进

    # ------------------------------------------------------------------------------
    #  重新训练二元模型（用于模拟）
    # ------------------------------------------------------------------------------
    # 准备训练数据的二值化版本
    y_train_station <- train_obs$Data[, i]
    y_train_bin <- ifelse(y_train_station < wet.threshold, 0, 1)

    # 获取训练集主成分
    pred_mat_regional <- data.pca$x.global[[station_region]]

    # 训练逻辑回归模型预测降水发生
    mod_binary <- glm(y_train_bin ~ .,
                      data = as.data.frame(pred_mat_regional)[,1:n_pcs_regional],
                      family = binomial(link = "logit"))

获取之前保存的Gamma模型

    # 获取之前保存的Gamma模型
    mod_gamma <- gamma_models_all_stations[[i]]

    # # 计算历史降水频率和概率阈值
    # wet_prob_station <- sum(y_train_bin, na.rm = TRUE) / length(na.exclude(y_train_bin))
    # prob_threshold <- quantile(mod_binary$fitted.values, 1 - wet_prob_station, na.rm = TRUE)

取出之前为当前站点 i 调用 glimpr 时保存的Gamma模型对象，用于模拟下雨日的降水量。

对于历史降水频率和概率阈值的计算，后面没用到这里，可以删除的

执行n_simulations次 (500次) 随机模拟抽样

    # ------------------------------------------------------------------------------
    #  执行500次随机模拟
    # ------------------------------------------------------------------------------
    for(sim in 1:n_simulations) {
      # 步骤1：预测降水发生概率
      prob_rain <- predict(mod_binary,
                           newdata = as.data.frame(sim_mat_regional)[,1:n_pcs_regional],
                           type = "response")

      # 步骤2：基于概率随机生成降水发生（0或1）
      rain_occurrence <- ifelse(runif(n_test_days) < prob_rain, 1, 0)

      # 步骤3：为降水日生成降水量
      if(!is.null(mod_gamma)) {
        # 预测条件降水强度（给定降水发生）
        cond_intensity <- predict(mod_gamma,
                                  newdata = as.data.frame(sim_mat_regional)[,1:n_pcs_regional],
                                  type = "response")

        # 从Gamma GLM提取分布参数
        gamma_shape <- 1/summary(mod_gamma)$dispersion     # 形状参数
        gamma_scale <- summary(mod_gamma)$dispersion       # 尺度参数

        # 从Gamma分布中采样降水量
        rain_amount <- rgamma(n = n_test_days,
                              shape = gamma_shape,
                              rate = gamma_shape/cond_intensity)

        # 最终降水量 = 降水发生标记 × 采样的降水量
        fold_simulations[sim, , i] <- rain_occurrence * rain_amount

      } else {
        # Gamma模型拟合失败的情况
        fold_simulations[sim, , i] <- NA
      }
    }

步骤1：预测降水发生概率

使用刚刚为模拟训练的 mod_binary 模型和测试集的区域主成分 sim_mat_regional (的前 n_pcs_regional 列) 来预测测试集每一天的降水发生概率。结果 prob_rain 是一个长度为 n_test_days 的概率向量。

步骤2：基于概率随机生成降水发生（0或1）

为测试集的每一天生成一个 [0,1) 区间内的均匀分布随机数（通过 runif(n_test_days)）。如果当天预测的降水概率 prob_rain 大于该随机数，则 rain_occurrence 在该天被设为1（模拟为发生降水），否则设为0（模拟为不发生降水）。

步骤3：为模拟的降水日生成降水量

首先检查之前保存的 mod_gamma 是否有效，如果有效，使用 mod_gamma 和测试集的区域主成分 sim_mat_regional 预测条件平均降水强度 cond_intensity (即假设当天有雨，预期的雨量均值)。

然后从 mod_gamma 的 summary() 中提取离散度参数 summary(mod_gamma)$dispersion。Gamma分布的形状参数 gamma_shape 通常等于 1 / dispersion。使用 rgamma() 函数从参数化的Gamma分布中随机抽样生成 n_test_days 个降水量值 rain_amount。这里使用的参数是 shape = gamma_shape 和 rate = gamma_shape / cond_intensity (因为对于Gamma分布，均值 = shape / rate)。

最终，将该次模拟 (sim)、该站点 (i) 在测试集所有天数的模拟降水量存入 fold_simulations 数组。这个值是 rain_occurrence (0或1的发生指示) 乘以随机抽样的 rain_amount。如果当天rain_occurrence 为0，则模拟降水量也为0。

然后汇报进度，以及当前折所有站点模拟完成的汇报

    # 显示进度信息
    if(i %% 10 == 0) {
      message("    已完成 ", i, "/", n_final_stations, " 个站点的采样")
    }
  }

  message("500次模拟采样完成")

计算当前折500次模拟的集合平均值

  # 计算500次模拟的集合平均值作为期望降水量
  fold_expected_values <- apply(fold_simulations, c(2, 3), mean, na.rm = TRUE)

对三维数组 fold_simulations (维度为 [模拟次数, 测试天数, 站点数]) 进行操作。apply 函数的第二个参数 c(2, 3)表示在第一个维度（模拟次数维度）上进行操作，即对每个测试日、每个站点，计算其500次模拟值的平均值。结果 fold_expected_values 是一个二维矩阵 [测试天数, 站点数]，代表了当前交叉验证折的集合平均（期望）降水预测。

保存当前折 (fold) 的所有结果到 cv_results 列表:

  # ==================================================================================
  #  保存当前折的所有结果
  # ==================================================================================
  # 确定性预测结果
  cv_results[[fold]]$predictions <- fold_predictions
  # 基于500次采样的平均预测（期望值）
  cv_results[[fold]]$predictions_simulated <- fold_expected_values
  # 降水发生概率
  cv_results[[fold]]$prob_estimates <- fold_prob_estimates
  # 条件降水强度
  cv_results[[fold]]$cond_intensity <- fold_cond_mean_intensity
  # 期望值（与predictions_simulated相同）
  cv_results[[fold]]$expected_value <- fold_expected_values
  # 实际观测值
  cv_results[[fold]]$observations <- test_obs$Data
  # 保存所有500次模拟的完整结果
  cv_results[[fold]]$all_simulations <- fold_simulations

  # 保存元数据信息
  cv_results[[fold]]$station_ids <- final_station_ids             # 站点标识
  cv_results[[fold]]$station_coords <- final_station_coords       # 站点坐标
  cv_results[[fold]]$station_regions <- station_region_mapping    # 站点-区域映射
  cv_results[[fold]]$test_indices <- test_idx                    # 测试集索引
  cv_results[[fold]]$train_indices <- train_idx                  # 训练集索引

  message("第 ", fold, " 折交叉验证处理完成")
}

message("\n所有交叉验证折完成！")

结果融合与评估 - 整合所有交叉验证折的预测结果

# ==================================================================================
#  结果融合与评估 - 整合所有折的预测结果
# ==================================================================================
message("Fusion de tous les résultats de la validation croisée et évaluation...")

# ------------------------------------------------------------------------------
#  初始化完整时间序列的结果矩阵
# ------------------------------------------------------------------------------
# 所有矩阵维度：[总天数, 站点数]
cv_all_obs <- matrix(NA_real_, nrow = total_days, ncol = n_final_stations)       # 观测值
cv_all_pred <- matrix(NA_real_, nrow = total_days, ncol = n_final_stations)      # 确定性预测
cv_all_prob <- matrix(NA_real_, nrow = total_days, ncol = n_final_stations)      # 降水概率
cv_all_cond_int <- matrix(NA_real_, nrow = total_days, ncol = n_final_stations)  # 条件强度
cv_all_expected <- matrix(NA_real_, nrow = total_days, ncol = n_final_stations)  # 期望值
# 基于模拟的预测值
cv_all_pred_simulated <- matrix(NA_real_, nrow = total_days, ncol = n_final_stations)

# ------------------------------------------------------------------------------
#  按时间顺序填充结果矩阵
# ------------------------------------------------------------------------------
for(fold in 1:cv_folds) {
  # 获取当前折的测试集时间索引
  test_indices <- cv_indices[[fold]]

  # 将测试集结果填充到对应时间位置
  cv_all_obs[test_indices, ] <- cv_results[[fold]]$observations           # 观测值
  cv_all_pred[test_indices, ] <- cv_results[[fold]]$predictions          # 确定性预测
  cv_all_prob[test_indices, ] <- cv_results[[fold]]$prob_estimates       # 降水概率
  cv_all_cond_int[test_indices, ] <- cv_results[[fold]]$cond_intensity   # 条件强度
  cv_all_expected[test_indices, ] <- cv_results[[fold]]$expected_value   # 期望值
  cv_all_pred_simulated[test_indices, ] <- cv_results[[fold]]$predictions_simulated  # 模拟预测
}

# ------------------------------------------------------------------------------
#  重命名变量以便后续分析
# ------------------------------------------------------------------------------
obs <- cv_all_obs                          # 观测数据
pred_mean <- cv_all_pred                   # 确定性预测均值
prob_for_auc <- cv_all_prob               # 用于计算AUC的概率值
cond_int_for_cor <- cv_all_cond_int       # 用于相关性分析的条件强度
expected_value <- cv_all_expected          # 期望降水量
pred_mean_simulated <- cv_all_pred_simulated  # 基于模拟的预测均值

# 获取最终的站点坐标信息
final_station_coords <- as.matrix(st86_precip$xyCoords)

# ------------------------------------------------------------------------------
#  保存所有结果到RDA文件
# ------------------------------------------------------------------------------
save(obs,                    # 观测值
     pred_mean,              # 确定性预测
     pred_mean_simulated,    # 模拟预测
     prob_for_auc,           # 降水概率
     cond_int_for_cor,       # 条件强度
     expected_value,         # 期望值
     final_station_coords,   # 站点坐标
     cv_results,             # 详细的交叉验证结果
     cv_indices,             # 交叉验证索引
     file = "regional_pca_glm_cv_results.rda")

message("交叉验证已完成。结果已保存在'regional_pca_glm_cv_results.rda'中。")

按时间顺序填充（拼接）结果矩阵，遍历之前存储了每折结果的 cv_results 列表，将每折测试集的特定结果填充到上述全局结果矩阵的相应时间位置。

重命名变量只是为了后面用的时候更简洁直观

第二部分 评判指标计算

逐折合并模拟结果到完整时间序列的 sim_full 数组

# ==================================================================================
#  合并交叉验证结果 - 整合所有折的模拟数据
# ==================================================================================
# 定义模拟次数
n_simulations <- 500

# 逐折合并结果到完整时间序列
for (fold in seq_len(cv_folds)) {
  # 获取当前折的测试集时间索引
  test_idx <- cv_results[[fold]]$test_indices

  # 填充观测值到对应时间位置
  cv_all_obs[test_idx, ] <- cv_results[[fold]]$observations

  # 填充各类预测结果
  cv_all_pred[test_idx, ] <- cv_results[[fold]]$predictions              # 确定性预测
  cv_all_prob[test_idx, ] <- cv_results[[fold]]$prob_estimates          # 降水概率
  cv_all_cond_int[test_idx, ] <- cv_results[[fold]]$cond_intensity      # 条件强度
  cv_all_expected[test_idx, ] <- cv_results[[fold]]$expected_value      # 期望值
  cv_all_pred_simulated[test_idx, ] <- cv_results[[fold]]$predictions_simulated  # 模拟预测均值

  # ------------------------------------------------------------------------------
  #  创建完整的500次模拟三维数组
  # ------------------------------------------------------------------------------
  if (fold == 1) {
    # 第一折时初始化数组，维度：[模拟次数, 总天数, 站点数]
    sim_full <- array(NA_real_,
                      dim = c(n_simulations, total_days, n_final_stations))
  }
  # 将当前折的模拟结果填充到对应时间位置
  sim_full[, test_idx, ] <- cv_results[[fold]]$all_simulations
}

message("500次模拟结果合并完成")

和之前合并确定性预测和观测值的循环类似，但是它处理三维的模拟数据。

cv_results[[fold]]$all_simulations 是一个三维数组，存储了当前折 fold 的测试集上所有站点、所有模拟次数的降水值。其维度是 [n_simulations, length(test_idx), n_final_stations]。

sim_full[, test_idx, ] 表示对 sim_full 数组进行索引，sim_full 的维度被设为 [n_simulations, total_days, n_final_stations]，并用 NA_real_ 填充。

  • 第一个维度 (模拟次数) 取全部 (: 或留空)。

  • 第二个维度 (时间) 使用 test_idx 进行索引，这意味着只选择与当前折测试集对应的那些“时间片”。

  • 第三个维度 (站点数) 取全部。

通过赋值操作，当前折的模拟数据就被准确地填充到了 sim_full 数组的相应时间段内。循环结束后，sim_full 将包含整个研究时段所有站点、所有500次模拟的完整降水数据。

评估指标计算

# ==================================================================================
#  评估指标计算 - 定义和实施模型性能评估
# ==================================================================================
# 加载ROC曲线分析包
library(pROC)

# 定义湿日阈值（1mm，用于区分有无降水）
wet_threshold <- 1

message("\n--- 开始计算评估指标 ---")
# ------------------------------------------------------------------------------
#  数据验证 - 确保所有必要数据格式正确
# ------------------------------------------------------------------------------
# 验证站点坐标格式
if (!is.matrix(final_station_coords) || ncol(final_station_coords) != 2) {
  stop("无法加载'final_station_coords'或格式不正确")
}
# 验证结果矩阵格式
if (!is.matrix(obs) || !is.matrix(pred_mean) || !is.matrix(prob_for_auc) || !is.matrix(expected_value)) {
  stop("无法加载结果矩阵'obs', 'pred_mean', 'prob_for_auc'或'expected_value'")
}
# 验证维度一致性
if (nrow(final_station_coords) != ncol(obs) || ncol(obs) != ncol(pred_mean) ||
    ncol(obs) != ncol(prob_for_auc) || ncol(obs) != ncol(expected_value)) {
  stop("站点坐标数量与观测/预测/概率/期望值的站点数量不匹配")
}
message("评估数据加载完成")

  • 检查 final_station_coords（最终使用的站点坐标）是否是一个两列的矩阵

  • 检查 obs（观测值）、pred_mean（确定性预测）、prob_for_auc（概率预测）和 expected_value（模拟期望值）是否都是矩阵。

  • 检查这些矩阵的列数（代表站点数）是否与 final_station_coords 的行数（也代表站点数）一致，并且它们之间的列数是否也相互一致。

如果任一检查失败，就直接报错。

打印一些关于合并后的预测结果的基本描述性统计，初步了解数据分布

# ------------------------------------------------------------------------------
#  数据分布概览 - 了解预测结果的整体特征
# ------------------------------------------------------------------------------
message("\n数据分布概览:")
message("概率预测范围: ", paste(round(range(prob_for_auc, na.rm = TRUE), 3), collapse = " - "))
message("期望值范围: ", paste(round(range(expected_value, na.rm = TRUE), 2), collapse = " - "))
# 统计湿日数量
wet_idx <- obs >= wet_threshold & !is.na(obs) & !is.na(expected_value)
message("湿日数量: ", sum(wet_idx), " / 总有效日数: ", sum(!is.na(obs)))
message("湿日期望值概要:")
print(summary(expected_value[wet_idx]))

对应结果

--- 开始计算评估指标 ---
成功加载评估数据

数据分布概览：
概率预测范围：0 - 1
观测值范围：0 - 304.06
湿天数：277054 / 有效总天数：907542
湿天观测值汇总：
   最小值    第一四分位数    中位数     平均值     第三四分位数   最大值 
  0.000      1.716       3.654     5.255      6.744     304.062

这里我有疑问，为什么最小值是0？前面定义了湿天降雨量>= 1mm，那么 obs 的最小值确实应该是1，奇了怪了。

  • 第一四分位数 (1st Qu.: 1.716)：25%的湿天降水量小于等于1.716。

  • 中位数 (Median: 3.654)：50%的湿天降水量小于或等于3.654。这是湿天降水量的典型值。

  • 平均值 (Mean: 5.255)：湿天的平均降水量为5.255。平均值大于中位数，表明湿天降水量的分布是右偏的（即存在一些极大的降水值拉高了平均数）。

  • 第三四分位数 (3rd Qu.: 6.744)：75%的湿天降水量小于或等于6.744。

  • 最大值 (Max: 304.062): 极端降雨

各种降水统计指标函数定义

# ------------------------------------------------------------------------------
#  R01 - 湿日频率（降水日占总日数的比例）
# ------------------------------------------------------------------------------
calculate_R01 <- function(precip_series, threshold) {
  valid_days <- sum(!is.na(precip_series))
  if (valid_days == 0) return(NA)
  wet_days <- sum(precip_series >= threshold, na.rm = TRUE)
  return(wet_days / valid_days)
}

计算给定降水时间序列 precip_series 中，降水量大于或等于 threshold（通常是 wet_threshold）的天数占总有效天数（非NA天数）的比例。

# ------------------------------------------------------------------------------
#  SDII - 简单日降水强度指数
# ------------------------------------------------------------------------------
# 新版本：只计算模型判定为湿日（>0）的平均降水量，不使用1mm阈值
calculate_SDII_no_threshold <- function(precip_series) {
  wet_days_precip <- precip_series[precip_series > 0 & !is.na(precip_series)]
  if (length(wet_days_precip) == 0) return(NA)
  return(mean(wet_days_precip))
}

# 原始版本：用于观测数据，保持1mm阈值的一致性
calculate_SDII <- function(precip_series, threshold) {
  wet_days_precip <- precip_series[precip_series >= threshold & !is.na(precip_series)]
  if (length(wet_days_precip) == 0) return(NA)
  return(mean(wet_days_precip))
}

这里真的有说法，我不好讲，怎么说呢，要么是我个人没玩明白，要么就是原论文有问题，这里细节很多，深入能牵扯到gamma建模部分，我后面详细汇总好好写一下，这里是我代码的定时炸弹，或者glm可能就是这样的，这里后面再补充吧

calculate_SDII_no_threshold - 它计算 precip_series 中所有大于0（即模型预测有降水，即使量很小，未达到1mm的 wet_threshold）且非NA的那些天的平均降水量。

calculate_SDII - 计算 precip_series 中所有降水量大于或等于 threshold（通常是 wet_threshold，如1mm）且非NA的那些湿润日的平均降水量。

# ------------------------------------------------------------------------------
#  P98 - 极端降水指标（第98百分位数）
# ------------------------------------------------------------------------------
# 标准版本：使用阈值
calculate_P98 <- function(precip_series, threshold) {
  wet_days_precip <- precip_series[precip_series >= threshold & !is.na(precip_series)]
  if (length(wet_days_precip) < 2) return(NA)
  return(quantile(wet_days_precip, probs = 0.98, na.rm = TRUE, type = 8))
}

# 无阈值版本：用于模型输出
calculate_P98_no_threshold <- function(precip_series) {
  wet_days_precip <- precip_series[precip_series > 0 & !is.na(precip_series)]
  if (length(wet_days_precip) < 2) return(NA)
  return(quantile(wet_days_precip, probs = 0.98, na.rm = TRUE, type = 8))
}

calculate_P98 -计算给定降水序列 precip_series 中，湿润日（降水量 >= threshold）降水量的第98百分位数。

calculate_P98_no_threshold - 所有模型预测有降水（>0）的日子来计算第98百分位数

干湿状态转换概率计算函数:

# ------------------------------------------------------------------------------
#  转换概率计算 - 干湿状态转换
# ------------------------------------------------------------------------------
# 创建前后天干湿状态的数据框
create_transition_data <- function(precip_series, threshold) {
  is_wet <- precip_series >= threshold & !is.na(precip_series)
  if (length(is_wet) < 2) return(NULL)
  data.frame(current = is_wet[-1], previous = is_wet[-length(is_wet)])
}

# DW - 干转湿概率
calculate_DW <- function(precip_series, threshold) {
  trans_data <- create_transition_data(precip_series, threshold)
  if (is.null(trans_data)) return(NA)
  n_dry_days <- sum(!trans_data$previous)
  if (n_dry_days == 0) return(NA)
  n_dry_to_wet <- sum(!trans_data$previous & trans_data$current)
  return(n_dry_to_wet / n_dry_days)
}

# WW - 湿转湿概率（降水持续性）
calculate_WW <- function(precip_series, threshold) {
  trans_data <- create_transition_data(precip_series, threshold)
  if (is.null(trans_data)) return(NA)
  n_wet_days <- sum(trans_data$previous)
  if (n_wet_days == 0) return(NA)
  n_wet_to_wet <- sum(trans_data$previous & trans_data$current)
  return(n_wet_to_wet / n_wet_days)
}

create_transition_data: 辅助函数，用于创建分析干湿状态转换所需的数据结构。具体思路为:

此函数首先将输入的降水序列 precip_series 根据 threshold 转换为一个布尔向量 is_wet (TRUE代表湿，FALSE代表干)。然后，它创建一个数据框，其中 current 列是 is_wet 从第二个元素开始到末尾的序列，previous 列是 is_wet 从第一个元素开始到倒数第二个元素的序列。这样，数据框的每一行就代表了相邻两天的干湿状态（前一天和当天）。

calculate_DW - 干转湿 (Dry-to-Wet) 概率

calculate_WW - 湿转湿 (Wet-to-Wet) 概率 (降水持续性)

五大指标计算

calculate_AUC - ROC曲线下面积 (AUC)

# ------------------------------------------------------------------------------
#  AUC - ROC曲线下面积（评估概率预测的判别能力）
# ------------------------------------------------------------------------------
calculate_AUC <- function(probabilities, observed_binary) {
  valid_idx <- !is.na(probabilities) & !is.na(observed_binary)
  # 需要足够的样本和两类都存在
  if (sum(valid_idx) < 10 || length(unique(observed_binary[valid_idx])) < 2) {
    return(NA)
  }
  roc_obj <- tryCatch(
    roc(response = observed_binary[valid_idx], predictor = probabilities[valid_idx], quiet = TRUE),
    error = function(e) NULL
  )
  if (is.null(roc_obj)) return(NA)
  return(as.numeric(auc(roc_obj)))
}

输入是模型预测的降水发生概率 probabilities 和对应的真实观测二值序列 observed_binary (0代表干，1代表湿)。它首先筛选出有效的成对数据点，并检查样本量是否足够且两类事件（0和1）是否都存在。然后调用 pROC::roc() 生成ROC对象，再用 pROC::auc() 计算其面积。

calculate_COR_paper - Spearman相关系数

# ------------------------------------------------------------------------------
#  Spearman相关系数 - 评估预测与观测的一致性
# ------------------------------------------------------------------------------
calculate_COR_paper <- function(observed_precip, predicted_expected_value) {
  cor_val <- tryCatch(
    cor(observed_precip, predicted_expected_value, method = "spearman", use = "pairwise.complete.obs"),
    error = function(e) NA
  )
  return(cor_val)
}

计算观测降水量 observed_precip 和模型预测的期望降水量 predicted_expected_value 之间的Spearman等级相关系数。method = "spearman" 指定了相关系数的类型，use = "pairwise.complete.obs" 表示在计算时只使用成对的有效数据点（即对应位置上两者都非NA的数据）。

基于500次模拟的集合计算分布相关的评估指标

在定义了各种评估函数之后，此部分专门处理从500次蒙特卡洛模拟得到的降水集合数据 (sim_full)，并计算指标的集合平均预测值。

# ==================================================================================
#  基于500次模拟计算分布指标
# ==================================================================================
message("计算所有站点的评估指标...")
n_stations_eval <- ncol(obs)
metrics_results_all <- list()

# 二值化观测数据（用于AUC计算）
obs_binary <- ifelse(obs >= wet_threshold & !is.na(obs), 1L, 0L)
obs_binary <- apply(obs_binary, 2, as.integer)

message("从500次模拟中计算分布指标...")

obs_binary 是将整个时间序列的观测数据 obs 根据 wet_threshold 转换为0/1矩阵，主要用于后续计算AUC等需要二元输入的指标。

定义批量计算函数 (用于处理模拟集合 sim_full)

# ------------------------------------------------------------------------------
#  定义批量计算函数
# ------------------------------------------------------------------------------
# 通用计算函数：对每次模拟的每个站点应用指标函数
calc_by_member <- function(fun, threshold = wet_threshold) {
  t(apply(sim_full, 1, function(mat)          # 遍历每次模拟
    apply(mat, 2, fun, threshold = threshold) # 对每个站点计算指标
  ))                                           # 返回 500 × n_stations 矩阵
}

  • apply(sim_full, 1, ...): 对 sim_full 数组的第一个维度（即500次模拟）进行迭代。对于每一次模拟，它得到一个二维矩阵 mat，其维度是 [总天数, 站点数]，代表了该次模拟中所有站点在所有时间点的降水值。

  • apply(mat, 2, fun, threshold = threshold): 对这个二维矩阵 mat 的第二个维度（即站点）进行迭代。对每个站点的时间序列，调用传入的指标函数 fun (并传递 threshold) 来计算指标值。这会为当前模拟的每个站点生成一个指标值。

整个 apply (外层) 会返回一个矩阵，行是站点，列是模拟次数。通过 t() 转置，最终 calc_by_member 返回一个维度为 [模拟次数, 站点数] 的矩阵，其中每个元素是特定模拟次数下特定站点的指标值。

calc_SDII_by_member - SDII专用按成员计算函数 (无阈值版)

# SDII专用计算函数（不使用阈值）
calc_SDII_by_member <- function() {
  t(apply(sim_full, 1, function(mat)
    apply(mat, 2, calculate_SDII_no_threshold) # 使用无阈值版本
  ))
}

只有SDII的计算中，阈值为0，对于R01 P98阈值仍然还是按照1mm来计算

计算500次模拟的各项指标的集合分布

# ------------------------------------------------------------------------------
#  计算500次模拟的各项指标
# ------------------------------------------------------------------------------
# R01：湿日频率（仍需要1mm阈值定义湿日）
R01_sim <- calc_by_member(calculate_R01)

# SDII：简单日降水强度（不设阈值）
SDII_sim <- calc_SDII_by_member()

# P98：极端降水（使用1mm阈值保持与观测的一致性）
P98_sim <- t(apply(sim_full, 1, function(mat)
  apply(mat, 2, calculate_P98, threshold = wet_threshold)
))

计算基于模拟集合平均的指标预测值

# 计算500次模拟的集合平均作为最终预测值
SDII_pred_by_stn <- colMeans(SDII_sim, na.rm = TRUE)  # 各站点SDII预测值
P98_pred_by_stn <- colMeans(P98_sim, na.rm = TRUE)    # 各站点P98预测值

对上述得到的 _sim 矩阵（维度 [模拟次数, 站点数]）按列求平均值，得到每个站点关于这500次模拟的平均指标值。这可以被看作是这些指标的一个集合平均预测。

逐站点计算所有评估指标

遍历每个评估站点 (n_stations_eval 个)，并为每个站点计算之前定义的一系列评估指标

# ==================================================================================
#  逐站点计算所有评估指标
# ==================================================================================
for (i in seq_len(n_stations_eval)) {
  if (i %% 10 == 0) message("  处理站点 ", i, "/", n_stations_eval)

  # 提取当前站点的所有数据
  obs_s   <- obs[, i]                    # 观测值
  prob_s  <- prob_for_auc[, i]          # 降水概率
  cond_s  <- cond_int_for_cor[, i]      # 条件强度
  exp_s   <- expected_value[, i]        # 期望值
  obs_bin <- obs_binary[, i]            # 二值化观测
  pred_s <- pred_mean[, i]              # 确定性预测

从之前合并好的全局结果矩阵中，提取出当前站点 i 对应的列向量。

然后检查当前站点是否有足够的有效数据

  # ------------------------------------------------------------------------------
  #  数据充足性验证
  # ------------------------------------------------------------------------------
  valid_obs <- sum(!is.na(obs_s)) >= 30
  valid_prob <- sum(!is.na(prob_s)) >= 30
  valid_cond <- sum(!is.na(cond_s)) >= 30

  if (!(valid_obs && valid_prob && valid_cond)) {
    # 数据不足，返回NA值
    metrics_results_all[[i]] <- list(
      station_id = final_station_ids[i],
      lon = final_station_coords[i, 1],
      lat = final_station_coords[i, 2],
      R01_obs = NA, R01_pred = NA,
      SDII_obs = NA, SDII_pred = NA,
      P98_obs = NA, P98_pred = NA,
      DW_obs = NA, DW_pred = NA,
      WW_obs = NA, WW_pred = NA,
      AUC = NA, COR_paper = NA
    )
    next
  }

计算当前站点 i 的各项观测指标

使用之前定义的 calculate_... 系列函数，基于当前站点的观测序列 obs_s 和湿日阈值 wet_threshold (1mm) 来计算。

R01_o  <- calculate_R01(obs_s, wet_threshold)  # 观测的湿日频率
SDII_o <- calculate_SDII(obs_s, wet_threshold) # 观测的湿日平均强度 (使用1mm阈值)
P98_o  <- calculate_P98(obs_s, wet_threshold)  # 观测的湿日P98极端降水 (使用1mm阈值)
DW_o   <- calculate_DW(obs_s, wet_threshold)   # 观测的干转湿概率
WW_o   <- calculate_WW(obs_s, wet_threshold)   # 观测的湿转湿概率

计算当前站点 i 的各项预测指标

  # ------------------------------------------------------------------------------
  #  计算预测指标
  # ------------------------------------------------------------------------------
  R01_p  <- calculate_R01(pred_s, wet_threshold)  # R01仍使用1mm阈值
  SDII_p <- SDII_pred_by_stn[i]  # 使用500次模拟的均值
  P98_p  <- P98_pred_by_stn[i]   # 使用500次模拟的均值

  # 转换概率使用概率阈值0.5转为二值序列
  pred_binary <- ifelse(prob_s >= 0.5, 1L, 0L)
  DW_p <- calculate_DW(pred_binary, 1)
  WW_p <- calculate_WW(pred_binary, 1)

  # 计算判别能力指标
  AUC_v    <- calculate_AUC(prob_s, obs_bin)
  COR_v    <- calculate_COR_paper(obs_s, exp_s)

  • R01_p: 使用 glimpr 输出的确定性降水预测序列 pred_s 和1mm阈值来计算预测的湿日频率。 这里要改成500次蒙特卡洛采样得到

  • SDII_p 和 P98_p: 这两个指标的预测值直接采用了之前基于500次蒙特卡洛模拟的集合平均结果(SDII_pred_by_stn[i] 和 P98_pred_by_stn[i])。即使用模拟集合的平均表现作为其“预测值”，而不是仅仅依赖于确定性GLM的输出

  • DW_p 和 WW_p: 预测的干湿转换概率，这里计算思路有问题，首先将 glimpr 输出的逐日降水发生概率 prob_s 与0.5比较，得到一个二值化的预测发生序列 pred_binary (≥0.5则为1，否则为0)，这里不对，应该改成和历史数据进行比较，然后基于这个 pred_binary 序列（将其视为降水序列，其中1代表湿，0代表干，阈值设为1即可正确处理）来计算DW和WW。

  • AUC_v: 使用 glimpr 输出的降水发生概率 prob_s 和二值化的真实观测序列 obs_bin 来计算ROC曲线下面积。

  • COR_v: 使用原始的观测降水序列 obs_s 和基于500次模拟的集合平均期望降水序列 exp_s 来计算它们之间的Spearman等级相关系数。

保存当前站点 i 的所有计算得到的指标结果

将当前站点 i 的站点信息（ID，经纬度）以及所有计算得到的观测指标 (_o) 和预测指标 (_p)、AUC值、相关系数值，作为一个列表元素存入 metrics_results_all[[i]]。

  # ------------------------------------------------------------------------------
  #  保存站点结果
  # ------------------------------------------------------------------------------
  metrics_results_all[[i]] <- list(
    station_id = final_station_ids[i],
    lon = final_station_coords[i, 1],
    lat = final_station_coords[i, 2],
    R01_obs = R01_o, R01_pred = R01_p,
    SDII_obs = SDII_o, SDII_pred = SDII_p,
    P98_obs = P98_o, P98_pred = P98_p,
    DW_obs = DW_o, DW_pred = DW_p,
    WW_obs = WW_o, WW_pred = WW_p,
    AUC = AUC_v, COR_paper = COR_v
  )
}

计算相对偏差

# ==================================================================================
#  计算相对偏差 - 评估模型系统性偏差
# ==================================================================================
# 相对偏差计算函数
calculate_relative_bias <- function(pred_val, obs_val) {
  bias <- ifelse(is.na(obs_val) | is.na(pred_val) | obs_val == 0,
                 NA,
                 100 * (pred_val - obs_val) / obs_val)
  return(bias)
}

# 计算所有指标的相对偏差并整理结果
metrics_df_eval <- metrics_df_all_raw %>%
  mutate(
    R01_Bias = calculate_relative_bias(R01_pred, R01_obs),
    SDII_Bias = calculate_relative_bias(SDII_pred, SDII_obs),
    P98_Bias = calculate_relative_bias(P98_pred, P98_obs),
    DW_Bias = calculate_relative_bias(DW_pred, DW_obs),
    WW_Bias = calculate_relative_bias(WW_pred, WW_obs)
  ) %>%
  select(station_id, lon, lat, AUC, COR_paper,
         R01_Bias, SDII_Bias, P98_Bias, DW_Bias, WW_Bias)

相对偏差计算函数 calculate_relative_bias: 计算单个预测值 pred_val 相对于观测值 obs_val 的相对偏差百分比。公式为 100 * (预测值 - 观测值) / 观测值。为了避免除以零或对NA值进行无效计算，函数首先检查 obs_val 或 pred_val 是否为NA，或者 obs_val 是否为0。

计算所有指标的相对偏差并整理结果到 metrics_df_eval

绘制评估指标地图可视化

plotMetricMap_fixed 画图而已，

#降尺度 4
决策树