信頼区間を実感する

データ分析基礎ノート

正規分布からのランダムサンプリング

\(N(165,10)\)から10個のランダムサンプリング

ex<-rnorm(10,165,10)
m<-mean(ex)
ci<-1.96*sd(ex)/sqrt(10)
c(m-ci,m,m+ci)

## [1] 165.3364 171.2051 177.0738

ランダムサンプリングの繰り返し

これを100回繰り返す

rdata<-replicate(100,{
  ex<-rnorm(10,165,10)
  m<-mean(ex)
  ci<-1.96*sd(ex)/sqrt(10)
  c(m-ci,m,m+ci)
})
rdata<-data.frame(1:100,t(rdata))
colnames(rdata)<-c("trials","minci","mean","maxci")
head(rdata)

##   trials    minci     mean    maxci
## 1      1 156.5805 162.1288 167.6770
## 2      2 156.8756 163.6460 170.4164
## 3      3 157.9112 164.0962 170.2812
## 4      4 163.8681 169.1025 174.3368
## 5      5 162.3142 167.1485 171.9829
## 6      6 156.9508 164.4389 171.9269

作図

matplot(rdata$trials, rdata[,c("minci","mean","maxci")], type="l")
segments(0,165,100,165)

正答率

sum(rdata$minci<165 & rdata$maxci > 165)/100

## [1] 0.93

ggplot2でかっこよく作図（インストールされてなければinstall.packages("ggplot2")）

library(ggplot2)
ggplot(data = rdata, aes(x = mean ,y = trials)) + 
  geom_point() + 
  geom_errorbarh(aes(xmin=minci, xmax=maxci)) +
  geom_vline(xintercept = 165, color="red")

関数化（要gglpot2）．m: 正規分布の平均，s: 正規分布の標準偏差，n: サンプルサイズ，r: 繰り返し回数

CItest<-function(m,s,n,r) {
  rdata<-replicate(r,{
    ex<-rnorm(n,m,s)
    m<-mean(ex)
    ci<-1.96*sd(ex)/sqrt(n)
    c(m-ci,m,m+ci)
  })
  rdata<-data.frame(1:r,t(rdata))
  colnames(rdata)<-c("trials","minci","mean","maxci")
  sr<-sum(rdata$minci<m & rdata$maxci > m)/r
  cat(paste("正答率: ",sr))
  ggplot2::ggplot(data = rdata, aes(x = mean ,y = trials)) + 
    geom_point() + 
    geom_errorbarh(aes(xmin=minci, xmax=maxci)) +
    geom_vline(xintercept = m, color="red")
}
CItest(160,20,10,500)

## 正答率:  0.928