R双预测变量高斯过程回归出现Xi - Xj非数值参数报错如何修复

你之前定义的平方指数核与协方差函数仅支持单变量输入，当输入为二维预测变量数据框时，outer函数传入核函数的参数为数据框对象，无法直接执行减法运算，因此触发报错。

1 重写适配多维输入的核函数与协方差计算逻辑 #

原来的单变量核函数直接对两个标量做差平方，多维场景下需要计算两个输入点之间的欧氏距离平方：

# 适配多维输入的平方指数核
SE <- function(Xi, Xj, l) {
  # Xi、Xj为长度等于特征维度的向量
  exp(-0.5 * sum((Xi - Xj)^2) / l^2)
}
# 适配多维输入的协方差矩阵生成函数
cov <- function(X, Y, l=1) {
  # X、Y为n行d列、m行d列的输入矩阵/数据框，d为特征维度
  n <- nrow(X)
  m <- nrow(Y)
  cov_mat <- matrix(0, nrow = n, ncol = m)
  for (i in 1:n) {
    for (j in 1:m) {
      cov_mat[i,j] <- SE(as.numeric(X[i,]), as.numeric(Y[j,]), l)
    }
  }
  return(cov_mat)
}

2 调整二维预测变量的生成逻辑 #

如果要覆盖二维平面所有网格点做预测，需要生成网格组合而不是直接拼接两个等长序列：

# 生成二维网格点作为预测输入，适当降低长度减少计算量，50*50=2500个点计算协方差会比较慢
x_predict_1 <- seq(-5,5,len=20) 
x_predict_2 <- seq(-6,6,len=20)
x_predict <- expand.grid(x1 = x_predict_1, x2 = x_predict_2)
l <- 1

# 生成先验协方差矩阵，不会再报错
COV <- cov(x_predict, x_predict, l)

3 先验采样与可视化示例 #

# 采样20个先验函数
values <- MASS::mvrnorm(20, rep(0, length=nrow(x_predict)), COV)
# 整理数据格式
dat <- x_predict %>% 
  mutate(rep = list(1:nrow(values))) %>% 
  unnest(rep) %>% 
  mutate(value = as.vector(t(values)))

# 可视化示例：绘制其中1个采样的先验函数热力图
ggplot(dat %>% filter(rep == 1), aes(x=x1, y=x2, fill=value)) +
  geom_tile() +
  scale_fill_viridis_c() +
  labs(title = "二维输入下的高斯过程先验采样示例", x = "预测变量1", y = "预测变量2")

4 加入观测值的拟合逻辑调整 #

# 构造二维观测数据
obs <- data.frame(
  x1 = c(-4, -3, -1,  0,  2),
  x2 = c(-2, 0, 1, 2, -1),
  y = c(-2,  0,  1,  2, -1)
)

# 计算后验均值与协方差，加微小扰动避免矩阵奇异
cov_xx_inv <- solve(cov(obs[,c("x1","x2")], obs[,c("x1","x2")], l) + diag(1e-6, nrow(obs))) 
Ef <- cov(x_predict, obs[,c("x1","x2")], l) %*% cov_xx_inv %*% obs$y
Cf <- cov(x_predict, x_predict, l) - cov(x_predict, obs[,c("x1","x2")], l) %*% cov_xx_inv %*% cov(obs[,c("x1","x2")], x_predict, l)

# 后验采样
values_post <- MASS::mvrnorm(20, Ef, Cf)

# 整理后验均值与置信区间数据
gp <- x_predict %>% 
  mutate(
    Ef = as.numeric(Ef),
    sigma = 2*sqrt(diag(Cf))
  )

# 可视化后验均值热力图+观测点
ggplot(gp, aes(x=x1, y=x2, fill=Ef)) +
  geom_tile() +
  geom_point(data = obs, aes(x=x1, y=x2, size = y), color = "red") +
  scale_fill_viridis_c(name = "后验均值f(x)") +
  labs(x = "预测变量1", y = "预测变量2")

• 二维预测场景下预测点数量是两个维度长度的乘积，50*50=2500个点的协方差矩阵计算耗时会远高于单变量的50个点，测试阶段可以适当降低网格长度提升运行速度
• 计算观测值协方差矩阵时加入了1e-6的单位矩阵扰动，避免观测点重复或接近导致矩阵奇异无法求逆

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