やりたいこと

・2D散布図を描画
・プロットの分布を2D正規分布で近似
・プロットの存在確率を示す確率楕円を描く

0. サンプルデータ

回転楕円みたいな分布を作る。

set.seed(13)
dat0 <- matrix(rnorm(1000), , 2)
dat0[,1] <- 2 * dat0[,1]

sita <- pi/3
rot <- matrix(c(cos(sita), -sin(sita), sin(sita), cos(sita)),2)

dat <- data.frame(dat0 %*% rot)

# 1. 確率楕円の算出

ellipse::ellipseを使って、楕円の中心をcentre、描きたい確率範囲をlevelで与える。

ellipse(cov(dat), centre=apply(dat, 2, mean), level=0.95)

2. 描画

楕円データの準備

level <- seq(0.5, 0.95, by=0.05)

datE <- NULL
for(i in level){
  datE <-rbind(datE, cbind(i, ellipse(cov(dat), centre=apply(dat, 2, mean), level=i)))
}
datE <- data.frame(datE)
datE$i <- factor(datE$i)

描画

ggplot()+
  layer(
    data = dat,
    mapping = aes(x=X1, y=X2),
    geom = "point",
    stat = "identity",
    position = "identity"
  )+
  layer(
    data = datE,
    mapping = aes(x=X1, y=X2, color=i),
    geom = "path",
    stat = "identity",
    position = "identity"
  )+
  scale_color_manual(values=grey(seq(0, 200, length=10)/255))+
  theme_bw()

layerを久しぶりに使ったらposition表記が必須になったようで戸惑った。

参考

書いたあとこんなのを見つけて…。
何かパラメータいじったら楕円になりそうな気がしなくもなくもないですね〜。

ggplot(dat, aes(X1, X2))+
  geom_density2d(n=10)+
  geom_point()