平均値の信頼区間のシミュレーションプロット

Rで，信頼区間を書いていいプロットするコードのメモ。
layout関数で，上手くプット画面を分割して，さも1枚のグラフかのように見せかけるのがポイント。
けっこう可愛い感じのグラフを描画できる。

# 
# 縦向きの，信頼区間の図。
# ポイントは，レイアウトとparの使い方。
# 
layout(matrix(c(1,2),ncol=1), heights=c(1/5,4/5))
par(mar=c(0,3,1,1))
mu <- 0
sigma <- 1
n <- 5
size <- 100

curve(dnorm(x,mu,sigma/sqrt(n)),from=-3, to = 3,axes=F, xlab = "",ylab="",lty =1, col =  rgb(1,0,1,alpha = 0.7),lwd = 2)
arrows(mu, 0, mu, dnorm(mu,mu,sigma/sqrt(n)), lwd=2, angle = 0,lty =1, col =  rgb(1,0,1,alpha = 0.7))

x <- matrix(rnorm(size*n,mean = mu,sd = sigma), ncol=n)
x_mean <- rowMeans(x)
upp <- x_mean + 1.96*sigma/sqrt(n)
low <- x_mean - 1.96*sigma/sqrt(n)
is_sig <- ((upp < 0) + (low < 0)) !=1

par(mar=c(3,3,0,1))
plot(x_mean,1:size, xlim = c(-3,3),col =ifelse(is_sig, rgb(1,0,0,alpha=0.7),  rgb(0,0,1,alpha=0.7)), pch = 1- is_sig ,xlab="x平均")
arrows(x_mean, 1:size, upp, 1:size, lwd=2, angle = 90,length = 0.03,col =ifelse(is_sig, rgb(1,0,0,alpha=0.7),  rgb(0,0,1,alpha=0.7)) , lty=1 + is_sig)
arrows(x_mean, 1:size, low, 1:size, lwd=2, angle = 90,length = 0.03,col =ifelse(is_sig, rgb(1,0,0,alpha=0.7),  rgb(0,0,1,alpha=0.7)) , lty=1 + is_sig)
abline(v = mu , lwd = 2,lty=2,col =  rgb(1,0,1,alpha = 0.7))

出来上がる図は次のようなもの。

#
# 横向きの，信頼区間の図。
#

layout(matrix(c(2,1),nrow=1), widths =c(4/5,1/5))
par(mar=c(3,0,1,1)) # 下，左，上，右の順で，余白を設定する。上が0なので，いい感じでくっつけられる。
mu <- 0
sigma <- 1
n <- 10
size <- 100

# 右に90度，図を転置させる。
t_dnorm <- function(x, mu, sigma){
  cbind(x, dnorm(x,mu,sigma)) %*% matrix(c(0,-1,1,0),ncol=2,byrow=T)
}
plot(t_dnorm(seq(from=-3, to =3, by=0.1), mu, sigma/sqrt(n)),type="l",axes=F, xlab = "",ylab="",lty =1, col =  rgb(1,0,1,alpha = 0.7),lwd = 2)
arrows(0, mu, dnorm(mu,mu,sigma/sqrt(n)), mu, lwd=2, angle = 0,lty =3, col =  rgb(1,0,1,alpha = 0.7))


x <- matrix(rnorm(size*n,mean = mu,sd = sigma), ncol=n)
x_mean <- rowMeans(x)
upp <- x_mean + 1.96*sigma/sqrt(n)
low <- x_mean - 1.96*sigma/sqrt(n)
is_sig <- ((upp < 0) + (low < 0)) !=1

par(mar=c(3,3,1,0)) # この余白の設定で，絶妙に図をくっつけるのがポイント。
plot(x_mean, ylim = c(-3,3),col =ifelse(is_sig, rgb(1,0,0,alpha=0.7),  rgb(0,0,1,alpha=0.7)), pch = 1- is_sig ,xlab="x平均")
arrows(1:size,x_mean, 1:size, upp, lwd=2, angle = 90,length = 0.03,col =ifelse(is_sig, rgb(1,0,0,alpha=0.7),  rgb(0,0,1,alpha=0.7)) , lty=1 + is_sig)
arrows(1:size,x_mean, 1:size, low,lwd=2, angle = 90,length = 0.03,col =ifelse(is_sig, rgb(1,0,0,alpha=0.7),  rgb(0,0,1,alpha=0.7)) , lty=1 + is_sig)
abline(h = mu , lwd = 2,lty=3,col =  rgb(1,0,1,alpha = 0.7))

出来上がる図は次のようなもの。