R在散点图上叠加双变量正态密度(椭圆)

R superimposing bivariate normal density (ellipses) on scatter plot

网站上也有类似的问题,但是我找不到这个看似很简单的问题的答案。我在 Old Faithful 数据集上混合了两个高斯:

if(!require("mixtools")) { install.packages("mixtools");  require("mixtools") }
data_f <- faithful
plot(data_f$waiting, data_f$eruptions)
data_f.k2 = mvnormalmixEM(as.matrix(data_f), k=2, maxit=100, epsilon=0.01) 
data_f.k2$mu # estimated mean coordinates for the 2 multivariate Gaussians
data_f.k2$sigma # estimated covariance matrix 

我只是想为均值向量 data_f.k2$mu 和协方差矩阵 data_f.k2$sigma 描述的模型的两个高斯分量叠加两个椭圆。得到类似的东西:

对于那些感兴趣的人,here 是创建上面图的 MatLab 解决方案。

您可以使用包 mixtools 中的 ellipse-函数。最初的问题是此函数从您的图中交换 x 和 y。我会尝试解决这个问题并更新答案。 (我会把颜色留给别人...)

plot( data_f$eruptions,data_f$waiting)
for (i in 1: length(data_f.k2$mu))  ellipse(data_f.k2$mu[[i]],data_f.k2$sigma[[i]])

如果你对颜色也感兴趣,你可以使用后验来得到合适的组。我用 ggplot2 做到了,但首先我使用@Julian 的代码展示了彩色解决方案。

# group data for coloring
data_f$group <- factor(apply(data_f.k2$posterior, 1, which.max))
# plotting
plot(data_f$eruptions, data_f$waiting, col = data_f$group)
for (i in 1: length(data_f.k2$mu))  ellipse(data_f.k2$mu[[i]],data_f.k2$sigma[[i]], col=i)

我的版本使用 ggplot2

# needs ggplot2 package
require("ggplot2")
# ellipsis data 
ell <- cbind(data.frame(group=factor(rep(1:length(data_f.k2$mu), each=250))), 
             do.call(rbind, mapply(ellipse, data_f.k2$mu, data_f.k2$sigma, 
                                   npoints=250, SIMPLIFY=FALSE)))

# plotting command
p <- ggplot(data_f, aes(color=group)) + 
  geom_point(aes(waiting, eruptions)) +
  geom_path(data=ell, aes(x=`2`, y=`1`)) +
  theme_bw(base_size=16)
print(p)

使用 mixtools 内部绘图功能:

plot.mixEM(data_f.k2, whichplots=2)