基于均值与逆方差绘制的mu1、mu2置信椭圆呈Y形问题求助

我用已有的均值和逆方差，写了代码绘制mu1和mu2的95%置信椭圆，但输出图形是Y形（如下）：

library(car)
library(SIBER)
ellipse<-function(mu1,mu2)
{
  30*1.155461e-05*(1860.500-mu1)^2+ 30*4.11493e-07*(8354.133-mu2)^2
  - 60*1.98504e-06*(1860.500-mu1)*(8354.133-mu2)
}
mu1<-seq(1650,2100,50)
mu2<-seq(7000,9500,50)
z <- t(sapply(mu1,ellipse,mu2))
contour(mu1,mu2,z, 
        levels=6.91, 
        drawlabels=T, 
        axes=T,
        frame.plot=T, 
        xlab="mu1 (stf)",
        ylab="mu2 (ben strnegth)",
        main = "A 95% confidence ellipse for mu = (mu1,mu2)")
points(1860.5,8354.133)
segments(0,8354.133,1860.500,8354.133)
segments(1860.500,0,1860.500,8354.133) 

问题根源 #

1. 二次型矩阵非正定：置信椭圆的核心是正定二次型(x-μ)^T Σ⁻¹(x-μ)=χ²，如果逆方差矩阵Σ⁻¹存在负特征值，就会画出双曲线（你看到的Y形就是双曲线的一支），而非椭圆。
2. 临界值错误：95%置信水平下，二维变量对应的卡方临界值是qchisq(0.95, df=2)=5.991，你用的6.91不符合统计要求，也会干扰图形形态。
3. 网格精度不足：步长设置为50过大，会导致图形边缘粗糙，放大异常形态。

修复方案 #

方案1：用现成函数直接绘制（推荐） #

借助car包内置的ellipse函数，无需手动构建方程，避免计算错误：

library(car)

# 定义均值向量
mu <- c(1860.5, 8354.133)
# 从代码中提取逆方差矩阵Σ⁻¹
inv_Sigma <- matrix(
  c(30*1.155461e-05, -30*1.98504e-06,
    -30*1.98504e-06, 30*4.11493e-07),
  nrow = 2, byrow = TRUE
)
# 计算协方差矩阵Σ
Sigma <- solve(inv_Sigma)

# 绘制基础散点图
plot(mu[1], mu[2], 
     xlim = c(1650, 2100), 
     ylim = c(7000, 9500),
     xlab = "mu1 (stf)",
     ylab = "mu2 (ben strength)",
     main = "95%置信椭圆",
     pch = 16)
# 添加95%置信椭圆
car::ellipse(mu, Sigma, level = 0.95, add = TRUE)
# 添加辅助参考线
abline(v = mu[1], h = mu[2], lty = 2)

方案2：修正手动方程 #

如果要坚持手动计算，需确保二次型正定，并修正临界值与网格精度：

ellipse <- function(mu1, mu2) {
  inv_Sigma <- matrix(
    c(30*1.155461e-05, -30*1.98504e-06,
      -30*1.98504e-06, 30*4.11493e-07),
    nrow = 2
  )
  x <- c(mu1, mu2)
  mu <- c(1860.5, 8354.133)
  # 正确计算二次型
  t(x - mu) %*% inv_Sigma %*% (x - mu)
}

# 缩小步长，提升图形平滑度
mu1 <- seq(1650, 2100, 10)
mu2 <- seq(7000, 9500, 10)
# 用outer生成网格数据，更高效
z <- outer(mu1, mu2, Vectorize(ellipse))

contour(mu1, mu2, z,
        levels = qchisq(0.95, df = 2), # 正确的95%临界值
        drawlabels = TRUE,
        axes = TRUE,
        frame.plot = TRUE,
        xlab = "mu1 (stf)",
        ylab = "mu2 (ben strength)",
        main = "95%置信椭圆"
)
points(1860.5, 8354.133, pch = 16)
abline(v = 1860.5, h = 8354.133, lty = 2)

额外检查 #

运行eigen(inv_Sigma)$values查看逆方差矩阵的特征值，如果有负数，说明你提供的逆方差数据本身有误，需要重新核对计算过程——正定矩阵的特征值必须全为正，才能生成椭圆。

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