如何在R中从第一性原理实现双变量密度计算与绘图？

原代码核心问题排查 #

• 遍历逻辑错误：原代码sapply(X, ...)是遍历矩阵的单个元素，而非逐行遍历二维坐标点，不符合密度计算要求
• 协方差矩阵运算错误：cov**0.5是对矩阵每个元素开平方，实际需要的是协方差矩阵行列式的平方根sqrt(det(cov))
• 运算符号误用：R 中矩阵乘法需使用%*%，*为元素级运算，会导致维度不匹配报错
• 变量命名冲突：原代码使用cov作为协方差矩阵变量名，覆盖了R内置的cov()协方差计算函数，存在潜在报错风险
• 向量运算维度问题：指数部分的矩阵乘法需要保证行向量×矩阵×列向量的维度匹配，最终输出标量密度值

修正后完整实现代码 #

library(matlib)
library(ggplot2)

# 生成2D网格坐标
number_cor <- 200 # 可根据精度需求调整为1000，数值越大计算耗时越长
x <- y <- seq(-5, 15, length.out = number_cor)
grid <- expand.grid(u = x, v = y)
X <- as.matrix(grid)

# 设置分布参数
d <- 2
p <- 0.75
# 若要得到偏态分布效果，可将mu设置为非零向量，例如 mu <- c(1,1)
mu <- c(0, 0)
cov_mat <- matrix(c(2, 0, 0, 2), nrow = 2, ncol = 2)
# 提前计算协方差矩阵的逆和行列式，避免循环中重复计算提升效率
det_cov <- det(cov_mat)
inv_cov <- inv(cov_mat)
# 求和上下界
lower <- 1
upper <- 10

# 单个坐标点的密度计算函数
mvgsn_pdf_single <- function(point) {
  k_seq <- seq(lower, upper, 1)
  sum_density <- 0
  for (k in k_seq) {
    weight <- p * ((1 - p) ^ (k - 1))
    norm_const <- 1 / (((2 * pi) ^ (d / 2)) * sqrt(det_cov) * (k ^ (d / 2)))
    center_diff <- point - k * mu
    exp_term <- exp((-1 / (2 * k)) * t(center_diff) %*% inv_cov %*% center_diff)
    sum_density <- sum_density + weight * norm_const * as.numeric(exp_term)
  }
  return(sum_density)
}

# 逐行计算所有网格点的密度值
grid$density <- apply(X, 1, mvgsn_pdf_single)

# 绘制双变量密度热图
ggplot(grid, aes(x = u, y = v, fill = density)) +
  geom_raster(interpolate = TRUE) +
  scale_fill_viridis_c() +
  labs(x = "u", y = "v", fill = "密度值") +
  theme_minimal()

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