基于Gaussian Process的加拿大人口插值绘图问题求助

数据获取代码 #

library(readr)
library(dplyr)

url <- "https://www150.statcan.gc.ca/t1/tbl1/en/dtl!downloadDbLoadingData-nonTraduit.action?pid=1710000901&latestN=0&startDate=19800101&endDate=20240701&csvLocale=en&selectedMembers=%5B%5B1%2C2%2C3%2C4%2C5%2C6%2C7%2C8%2C9%2C10%2C11%2C12%2C14%2C15%5D%5D&checkedLevels="
data <- data.frame(read_csv(url))
data <- data[data$GEO == "Canada",]
data$date_column <- as.Date(paste(data$REF_DATE, "-01", sep=""), format="%Y-%m-%d")
finaldf <- data |> 
  dplyr::select(date_column, VALUE) |>
  dplyr::arrange(date_column)

需求说明 #

基于上述季度人口数据，使用高斯过程（Gaussian Process）插值生成缺失的月度人口数据，并绘制包含不确定性区间的结果图。

原高斯过程代码存在的问题及修复建议 #

核心问题分析 #

你改编的代码存在以下关键问题，导致结果不符合预期：

• 核函数长度尺度与时间跨度不匹配
• 未对原始数据做标准化处理
• 预测点粒度不符合月度要求
• 可视化时未将数值日期转回日期格式
• 协方差矩阵计算存在数值稳定性风险

修复后的完整代码 #

library(MASS)
library(tidyverse)

set.seed(12345)

# 1. 处理日期与标准化数据
finaldf$date_numeric <- as.numeric(finaldf$date_column)
obs <- data.frame(x = finaldf$date_numeric, y = finaldf$VALUE)
# 标准化人口数据，解决数值尺度问题
obs$y_scaled <- scale(obs$y)

# 2. 定义平方指数核函数
SE <- function(Xi, Xj, l) {
  exp(-0.5 * (Xi - Xj)^2 / l^2)
}

cov_function <- function(X, Y, l) {
  outer(X, Y, SE, l)
}

# 3. 生成月度预测点（而非等长序列）
date_seq_monthly <- seq(min(finaldf$date_column), max(finaldf$date_column), by = "month")
x_predict <- as.numeric(date_seq_monthly)
# 调整长度尺度为5年（匹配人口增长的长期趋势）
l <- 365 * 5

# 4. 计算后验分布
jitter <- 1e-4  # 增大抖动项提升稳定性
cov_xx <- cov_function(obs$x, obs$x, l) + jitter * diag(length(obs$x))
cov_xx_inv <- ginv(cov_xx)  # 使用广义逆避免奇异问题

# 计算标准化后的后验均值与协方差
Ef_scaled <- cov_function(x_predict, obs$x, l) %*% cov_xx_inv %*% obs$y_scaled
Cf <- cov_function(x_predict, x_predict, l) - cov_function(x_predict, obs$x, l) %*% cov_xx_inv %*% cov_function(obs$x, x_predict, l)
Cf <- Cf + jitter * diag(nrow(Cf))

# 将后验均值还原为原始尺度
Ef <- Ef_scaled * attr(obs$y_scaled, "scaled:scale") + attr(obs$y_scaled, "scaled:center")

# 5. 生成后验样本
values_post <- mvrnorm(200, Ef, Cf)
dat_post <- data.frame(x = x_predict, t(values_post)) %>%
  tidyr::pivot_longer(-x, names_to = "rep", values_to = "value") %>%
  mutate(rep = as.numeric(as.factor(rep)))

# 6. 准备可视化数据（转回日期格式）
gp <- data.frame(
  date = as.Date(x_predict, origin = "1970-01-01"),
  Ef = Ef,
  sigma = 2 * sqrt(pmax(0, diag(Cf))) * attr(obs$y_scaled, "scaled:scale")  # 不确定性区间还原尺度
)
obs$date <- as.Date(obs$x, origin = "1970-01-01")

# 7. 绘制结果图
ggplot(dat_post) + 
  # 后验样本曲线
  geom_line(aes(x = as.Date(x, origin = "1970-01-01"), y = value, group = rep), 
            color = rgb(0.7, 0.1, 0.4), alpha = 0.2) +
  # 95%不确定性区间
  geom_ribbon(data = gp, aes(x = date, ymin = Ef - sigma, ymax = Ef + sigma), 
              fill = rgb(0.7, 0.1, 0.4), alpha = 0.3) +
  # 后验均值曲线
  geom_line(data = gp, aes(x = date, y = Ef), size = 1, color = "darkred") +
  # 原始季度数据点
  geom_point(data = obs, aes(x = date, y = y), size = 2, color = "black") +
  labs(title = "加拿大人口月度插值（高斯过程）", x = "日期", y = "人口数量") +
  theme_minimal()

备选方案：线性插值（na.approx） #

如果你只需要简单的月度插值，不需要不确定性估计，可以使用以下线性插值代码：

library(readr)
library(dplyr)
library(lubridate)
library(zoo)

temp_file <- tempfile(fileext = ".csv")

tryCatch({
    url <- "https://www150.statcan.gc.ca/t1/tbl1/en/dtl!downloadDbLoadingData-nonTraduit.action?pid=1710000901&latestN=0&startDate=19800101&endDate=20240701&csvLocale=en&selectedMembers=%5B%5B1%2C2%2C3%2C4%2C5%2C6%2C7%2C8%2C9%2C10%2C11%2C12%2C14%2C15%5D%5D&checkedLevels="
    download.file(url, temp_file, mode = "wb")
    
    data <- data.frame(read_csv(temp_file))
    
    data <- data[data$GEO == "Canada",]
    data$date_column <- as.Date(paste(data$REF_DATE, "-01", sep=""), format="%Y-%m-%d")
    final_df <- data[,c("date_column", "VALUE")]
    
    final_df$date_column <- as.Date(final_df$date_column)
    
    current_date <- floor_date(Sys.Date(), "month")
    date_seq <- seq(from = min(final_df$date_column),
                    to = current_date,
                    by = "month")
    
    monthly_df <- data.frame(date_column = date_seq)
    
    merged_df <- merge(monthly_df, final_df, by = "date_column", all.x = TRUE)
    merged_df$VALUE_interpolated <- na.approx(merged_df$VALUE, na.rm = FALSE)
    
}, error = function(e) {
    message("数据处理出错: ", e$message)
}, warning = function(w) {
    message("数据处理警告: ", w$message)
}, finally = {
    if(file.exists(temp_file)) {
        unlink(temp_file)
    }
})

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