Rによる演習

演習問題のマグロの例を用いてRでも計算してみましょう．ただし，Ｒを使う前にまず手計算を行い，答の確認のために利用して下さい．

自分でRで式をつくって計算する場合

分散\(\sigma^2\)が既知のとき

例えば，信頼水準 95\(\%\)の時．

y <- c(42, 40, 50, 39, 40, 43, 42, 48, 43)
Ns <- length(y)   #データの数(n)
Mean <- mean(y)   #標本平均
sigma <- 6        #既知の標準偏差sigma

alpha <- 0.05
ZZ <- qnorm(1-alpha/2, mean=0, sd=1)
LB <- Mean - sigma/sqrt(Ns)*ZZ
UB <- Mean + sigma/sqrt(Ns)*ZZ
data.frame(LB,UB)
##         LB       UB
## 1 39.08007 46.91993

例えば，信頼水準 90\(\%\), 99\(\%\)の時も同様に計算してみて下さい．

分散\(\sigma^2\)が未知のとき

今度は\(\sigma\)を推定し，\(t\)分布を用いることになります．

y <- c(42, 40, 50, 39, 40, 43, 42, 48, 43)
Ns <- length(y)   #データの数(n)
Mean <- mean(y)   #標本平均
sigma.est <- sd(y)    #未知の標準偏差sigmaの推定値

alpha <- 0.05
TT <- qt(1-alpha/2, df=Ns-1)   #自由度 (n-1)のt分布のパーセント点
LB <- Mean - sigma.est/sqrt(Ns)*TT
UB <- Mean + sigma.est/sqrt(Ns)*TT
data.frame(LB,UB)
##        LB      UB
## 1 40.1497 45.8503

Rの関数”t.test”を用いて計算する方法

分散\(\sigma^2\)が未知のとき

この信頼区間はRの既存の関数 \(t.test\) を用いても容易に計算できる． この関数は仮説検定でも使える関数で，それについては次週扱います．

Res <- t.test(y, conf.level=0.95); Res
## 
##  One Sample t-test
## 
## data:  y
## t = 34.789, df = 8, p-value = 5.098e-10
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  40.1497 45.8503
## sample estimates:
## mean of x 
##        43
Res$conf.int
## [1] 40.1497 45.8503
## attr(,"conf.level")
## [1] 0.95

分散が既知の場合

分散が既知の場合の信頼区間は，実際にはほとんど使われません．分散が予め分かっているデータなどほとんどないからです．ただ，この授業でもそうですが，説明の順序と分かりやすいさから分散既知の場合を通して方法の概念を理解して，そして分散が未知の場合の実際的な解析方法を学びます．

Rで分散が既知の場合の関数がありませんでしたので，ここでは私が簡易的に作った関数を示しておきます．参考にしてみて下さい．

z.test <- function(y, mu, sigma, alpha){
   
Ns <- length(y)
Mean <- mean(y)
ZZ <- qnorm(1-alpha/2, mean=0, sd=1)

# For confidence interval
LB <- Mean - sigma/sqrt(Ns)*ZZ
UB <- Mean + sigma/sqrt(Ns)*ZZ
cat("# Mean = ", Mean, "\n")   #catは表示するコマンド，"\n"は改行コード
cat(paste("#", eval(1-alpha)*100, "percent confidence interval: "), c(LB,UB), "\n")

# For statistical test (これは今回は気にしないでください)
if(missing('mu')==FALSE){
   Zvalue <- sqrt(Ns)*(Mean-mu)/sigma
   cat(paste("#", "p-value = "), 2*(1-pnorm(abs(Zvalue))), "\n")
}

}
z.test(y, sigma=6, alpha=0.05)
## # Mean =  43 
## # 95 percent confidence interval:  39.08007 46.91993

図4のコード

# シミュレーションの設定
set.seed(2020) #乱数のシード(結果の再現性の為)
Nsim <- 10000  #シミュレーションの回数
Ns <- 50       #シミュレーションの各回におけるデータ数
mu <- 5        #シミュレーションで仮定した真の期待値パラメータmu
sigma <- 2     #シミュレーションで仮定した真の標準偏差パラメータsigma
CI <- array(NA, c(Nsim, 2)) #シミュレーションで計算された信頼区間を保存する箱（10000*2の配列形式)

# シミュレーションの実行（毎回データを50個生成し，それを基に信頼区間を計算する操作を10000回繰り返す)
for(i in 1:Nsim){
  y <- rnorm(Ns, mu, sigma)
  CI[i,] <- t.test(y, conf.level=0.95)$conf.int
}

# 信頼区間が真のmuを含んでいるかどうかをチェック
Judge <- (CI[,1] <= mu & mu <= CI[,2]) #信頼区間が真の値を含んでいるかどうかをTRUE or FALSEで判断
#含んでいる場合を1, 含んでいない場合を0とし，平均する事で含んでいる割合が計算される
CI.actual.level <- mean(Judge)         

# 図示(10000回分を表示するわけにいかないので，ランダムに100回分を抽出して表示)
Nuse <- 100
Use <- sample(1:Nsim, Nuse)　# 1～10000の中から100をランダムに取り出す

plot(0,0, type="n", xlab="mu", ylab="", xlim=c(mu-sigma, mu+sigma), ylim=c(0,Nuse))
#真のmuを含んでいない場合を赤で表示
for(i in 1:Nuse) lines(CI[Use[i],], c(i,i), col=2-Judge[Use[i]]) 
lines(c(mu,mu), c(0,Nuse), col="blue") #縦の線