[R program]正規分布を綺麗な感じで描く

ggplot2を使って正規分布を綺麗に描いてみた。

正規分布の累積確率、たまに忘れそうになるんだよね。。。

ggplot2のgeomリストはこちらです。

http://had.co.nz/ggplot2/

参考になる本はこちら（前にも紹介したけど）。

Rコードはこちら↓

矢印をlinesで描いちゃったんだけど、arrowsを使った方が楽だったね。。。

library(ggplot2)

#---数式で記述
gauss.dencity <- function(x) 1/sqrt(2*pi)*exp(-x^2/2)
plot(gauss.dencity, -5, 5, xaxs="i",yaxs="i", bty="l", ylim=c(0, 0.5), las=1)

#---ダミーデータで描く
x4 <- seq(-5, -4, by=0.01)
x3 <- seq(-3.99, -3, by=0.01)
x2 <- seq(-2.99, -2, by=0.01)
x1 <- seq(-1.99, -1, by=0.01)
x0 <- seq(-0.99, 0, by=0.01)
x99 <- seq(0.01, 5, by=0.01)

xy4 <- data.frame(x=x4, y=gauss.dencity(x4), range=rep("~-4", length(x4)))
xy3 <- data.frame(x=x3, y=gauss.dencity(x3), range=rep("~-3", length(x3)))
xy2 <- data.frame(x=x2, y=gauss.dencity(x2), range=rep("~-2", length(x2)))
xy1 <- data.frame(x=x1, y=gauss.dencity(x1), range=rep("~-1", length(x1)))
xy0 <- data.frame(x=x0, y=gauss.dencity(x0), range=rep("~0", length(x0)))
xy99 <- data.frame(x=x99, y=gauss.dencity(x99), range=rep("all", length(x99)))

xy <- rbind(xy4, xy3, xy2, xy1, xy0, xy99)

#ggplot(xy , aes(x, y, col=range)) + geom_line() + ylim(0, 0.5)
ggplot(xy , aes(x, y, fill=range)) + geom_area(alpha=0.3) + ylim(0, 0.5)

lines(c(-0.6, -0.6), c(0.425, 0.47))
lines(c(-0.6, -0.7), c(0.425, 0.435))
lines(c(-0.6, -0.5), c(0.425, 0.435))
text(-0.6, 0.49, "P< 0.5")

lines(c(-1.5, -1.5), c(0.25, 0.4))
lines(c(-1.5, -1.6), c(0.25, 0.26))
lines(c(-1.5, -1.4), c(0.25, 0.26))
text(-1.5, 0.42, paste("P< ", round(pnorm(-1), digits=2), sep=""))

lines(c(-2.3, -2.3), c(0.1, 0.2))
lines(c(-2.3, -2.4), c(0.1, 0.11))
lines(c(-2.3, -2.2), c(0.1, 0.11))
text(-2.3, 0.22, paste("P< ",  round(pnorm(-2), digits=3), sep=""))

lines(c(-3.3, -3.3), c(0, 0.1))
lines(c(-3.3, -3.4), c(0, 0.01))
lines(c(-3.3, -3.2), c(0, 0.01))
text(-3.3, 0.12, paste("P< ",  round(pnorm(-3), digits=4), sep=""))

text(1.5, 0.5, adj=0, paste("Φ(0.4) = ", round(qnorm(0.4), digits=2), sep=""))
text(1.5, 0.475, adj=0, paste("Φ(0.3) = ", round(qnorm(0.3), digits=2), sep=""))
text(1.5, 0.45, adj=0, paste("Φ(0.2) = ", round(qnorm(0.2), digits=2), sep=""), col=2)
text(1.5, 0.425, adj=0, paste("Φ(0.1) = ", round(qnorm(0.1), digits=2), sep=""))
text(1.5, 0.4, adj=0, paste("Φ(0.05) = ", round(qnorm(0.05), digits=2), sep=""))
text(1.5, 0.375, adj=0, paste("Φ(0.025)= ", round(qnorm(0.025), digits=2), sep=""), col=2)

#ついでにサンプルサイズ
text(0.5, 0.25, "sample size \n (a=0.05, b=0.2) \n = 2(1.96+0.84)^2*(s/d)^2 \n = 2*15.68*(s/d)^2 \n (= 32)", adj=0)



#---密度で記述
x <- rnorm(1000)
x.Lab <- data.frame(x=x, range=ifelse(x < -4, "~-4", ifelse(x < -3, "~-3", ifelse(x < -2, "~-2", ifelse(x < -1, "~-1", ifelse(x < 0, "~0", "all"))))))
ggplot(x.Lab, aes(x, fill=range)) + geom_density(alpha=0.3)




#おまけ
library(MASS)
data(quine)
q <- ggplot(quine, aes(Days, fill=Age))
q <- q + geom_density(alpha=0.3)
print(q)