やりたいこと
・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()
