如何在R中绘制椭球概率密度图？复现论文3D散点图遇阻

scatterplot3D包功能较为基础，无法实现目标图中的密度着色散点、透明椭球轮廓以及自定义视角等核心元素，推荐使用rgl（适合静态/交互式3D绘图）或plotly（侧重交互式）来复现，以下是具体实现步骤：

步骤1：生成椭球分布模拟数据 #

目标图对应多元正态分布样本，先生成符合椭球分布的数据集：

library(MASS)
library(rgl)

set.seed(123)
mean_vec <- c(0, 0, 0)
# 自定义协方差矩阵控制椭球形状
cov_mat <- matrix(c(2, 1, 0.5, 1, 3, 1, 0.5, 1, 2), nrow=3)
data <- mvrnorm(n=5000, mu=mean_vec, Sigma=cov_mat)
colnames(data) <- c("X", "Y", "Z")

步骤2：计算点的概率密度（用于颜色映射） #

通过核密度估计给每个数据点匹配密度值，作为颜色依据：

library(ks)
# 3D核密度估计
kd <- kde(x=data, h=Hpi(data))
# 给数据添加密度列
data$density <- predict(kd, x=data)

步骤3：绘制带密度着色的3D散点 #

用rgl实现散点的颜色渐变，调整点的大小和透明度：

# 用viridis配色匹配目标图的蓝-黄-红渐变
col_map <- viridis::viridis(n=100)
color_indices <- cut(data$density, breaks=100, labels=FALSE)
point_colors <- col_map[color_indices]

# 绘制散点，隐藏坐标轴和标签匹配目标图
plot3d(data$X, data$Y, data$Z, 
       col=point_colors, 
       size=1,
       alpha=0.6,
       xlab="", ylab="", zlab="",
       axes=FALSE)

步骤4：添加透明椭球轮廓 #

根据协方差矩阵生成椭球网格，添加线条式透明轮廓：

# 生成椭球网格函数
ellipse3d <- function(cov, center=c(0,0,0), level=0.95, segments=50) {
  eig <- eigen(cov)
  radii <- sqrt(eig$values * qchisq(level, 3))
  theta <- seq(0, 2*pi, length.out=segments)
  phi <- seq(0, pi, length.out=segments)
  x <- radii[1] * outer(cos(theta), sin(phi))
  y <- radii[2] * outer(sin(theta), sin(phi))
  z <- radii[3] * outer(rep(1, segments), cos(phi))
  rot <- eig$vectors
  coords <- cbind(c(x), c(y), c(z)) %*% rot
  coords <- coords + matrix(center, nrow=nrow(coords), ncol=3, byrow=TRUE)
  list(x=matrix(coords[,1], nrow=segments), 
       y=matrix(coords[,2], nrow=segments), 
       z=matrix(coords[,3], nrow=segments))
}

# 添加95%置信椭球（白色透明线条）
ell <- ellipse3d(cov_mat, center=mean_vec, level=0.95)
surface3d(ell$x, ell$y, ell$z, 
          col="white", alpha=0.2,
          front="lines", back="lines")

步骤5：调整视角匹配目标图 #

手动拖动调整或用参数固定视角：

# 可手动拖动窗口调整视角，满意后用rgl.viewpoint()获取参数
view3d(theta=30, phi=30, zoom=0.8)

步骤6：导出静态图片 #

rgl.snapshot("ellipsoid_density_plot.png", fmt="png")

替代方案：plotly交互式版本 #

如果需要交互式图表，用plotly更便捷：

library(plotly)

plot_ly(data, x=~X, y=~Y, z=~Z, 
        type="scatter3d", mode="markers",
        marker=list(color=~density, colorscale="Viridis", opacity=0.6, size=2)) %>%
  add_trace(type="mesh3d", x=ell$x, y=ell$y, z=ell$z, 
            opacity=0.2, color="white", showlegend=FALSE) %>%
  layout(scene=list(xaxis=list(visible=FALSE), yaxis=list(visible=FALSE), zaxis=list(visible=FALSE)))

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