廿TT

譬如水怙牛過窓櫺 頭角四蹄都過了 因甚麼尾巴過不得

母比率の信頼区間確かめ算

A/Bテストのガイドライン:仮説検定はいらない(Request for Comments|ご意見求む) - 廿TT

以前この記事で母比率の95%信頼区間を出したんですが、ちょっとやり方が不味かったかなーと思いまして、やり直しました。

結論としては、とりあえず母比率の信頼区間を出すときは群馬大学の青木先生のコードを使ってください。
R -- 母比率の信頼区間

ぼくのやり方はこちら。
なにも考えずにqbinomが返してくる値をnで割ってる。

my_prop.conf_95<- function(r,n){
p <- r/n
upper <-qbinom(0.975,size=n,prob=r/n)/n
lower <-qbinom(0.025,size=n,prob=r/n)/n
list(lower.p=lower, upper.p=upper)
}

ぼくの書いた関数と青木先生の書いた関数で結果を比較します。

#青木先生のコード
source("http://aoki2.si.gunma-u.ac.jp/R/src/prop_conf.R", encoding="euc-jp")
CI1<- prop.conf(175, 500) # F 分布を使って正確な信頼区間を求める場合
CI2<- prop.conf(175, 500, approximation=TRUE) # 正規分布による近似を行う場合
#
#ぼくのやり方
myCI<-my_prop.conf_95(175, 500) #なにも考えずにqbinomが返してくる値をnで割ってる

出力結果をみたところ、ぼくのやり方(myCIのほう)は粗い近似になっちゃってるみたいです。

> CI1
$lower.p
[1] 0.3081869

$upper.p
[1] 0.3935989

> CI2
$lower.p
[1] 0.3094802

$upper.p
[1] 0.3928071

> myCI
$lower.p
[1] 0.308

$upper.p
[1] 0.392

ただ、元の目的がエラーバーで傾向を見るというものでしたので、グラフを書いて見たところ、

plot(c(p,p,p),xlab="",ylab="" ,pch=16, xaxt  = "n", bty = "n",
     ylim=range(c(unlist(CI1),unlist(CI2),unlist(myCI))))
arrows(1:3,c(CI1$upper.p,CI2$upper.p,myCI$upper.p),
       1:3,c(CI1$lower.p,CI2$lower.p,myCI$lower.p),
       length=.05,angle=90,code=3, col="blue3")
axis(side=1,at=1:3, labels=c("CI1", "CI2", "myCI"))

A/Bテストのガイドライン:仮説検定はいらない(Request for Comments|ご意見求む) - 廿TT
↑この記事自体の致命的な欠陥にはなってないんじゃないかなー、と思います。自分に甘いかな。
f:id:abrahamcow:20140210060047p:plain

「F 分布を使って正確な信頼区間を求める」という意味がわからない方は、とりあえずブラックボックスとして「F 分布を使って正確な信頼区間を求める」方を使ってください。

入門・演習 数理統計

入門・演習 数理統計

追記(2014年2月15日)

上のように書いたところ、すぐにコメントでご教示いただきました。
R デフォルトのbinom.test関数で、F 分布を使ったのと同じ正確な信頼区間が求められます。

> ans <- binom.test(175, 500, conf.level = 0.95)
> ans$conf.int
[1] 0.3081869 0.3935989
attr(,"conf.level")
[1] 0.95

ありがとうございました。