「統計学関連なんでもあり」の過去ログ--- 045

No.16897　二項分布による検出力　　【赤岳】　2012/05/10(Thu) 21:08

いつも大変お世話になっております。
Rを使い始めて1ヵ月ほどの初心者です。以下のことをご教示いただけませんでしょうか。
不具合の発現率が1％です。
200個を調べたとき，この不具合が少なくとも1件発生する確率は，
Power<- 1-dbinom(0, 200, 0.01)
で求められるかと思います。

また，200個見たとき，不具合が3件超発生する確率を求めると，
a<- dbinom(0, 200, 0.01)
b<- dbinom(1, 200, 0.01)
c<- dbinom(2, 200, 0.01)
d<- dbinom(3, 200, 0.01)
PowerA<- 1-(a+b+c+d)
と考えています。この考えでよろしいでしょうか。
また，これを普遍的にするため，関数を作成しました。

PowerB<-function(n,pt,rate)
{
aa<-dbinom(0,pt,rate)
cc<-0
for(i in 1:n){
bb<-dbinom(n,pt,rate)
cc<-cc+bb
}
return(1-(cc+aa))
}
＃nは不具合の件数，ptは不具合を調べる母数，rateは不具合の発生率
PowerB（3, 200, 0.01) #これが解。ただし，n<0

ところが，上記PowerAとBを実行すると，
> PowerA
[1] 0.141966
> PowerB(3,200,0.01)
[1] 0.3219543
となり，どうしてもAとBが一致しません。
一体，どこが間違っているのでしょうか。
ご教示お願いします。

No.16898　Re: 二項分布による検出力　　【青木繁伸】　2012/05/10(Thu) 21:51

プログラム中の
bb<-dbinom(n,pt,rate)
は
bb<-dbinom(i,pt,rate)
でしょう。
<- の前後や , の後や { の前などに適切に空白を置いて，インデントもちゃんとやれば，プログラムが見やすくなり，バグにも気づきやすくなるのではないかと思います。
a <- dbinom(0, 200, 0.01)
b <- dbinom(1, 200, 0.01)
c <- dbinom(2, 200, 0.01)
d <- dbinom(3, 200, 0.01)
PowerA <- 1 - (a + b + c + d)
PowerA
PowerB <- function(n, pt, rate) {
    aa <- dbinom(0, pt, rate)
    cc <- 0
    for (i in 1:n) {
        bb <- dbinom(i, pt, rate)
        cc <- cc + bb
    }
    return(1 - (cc+aa))
}
PowerB(3, 200, 0.01)  #これが解。ただし，n<0 
また，dbinom(0, pt, rate)を特別扱いする必要はないのだから，
PowerB <- function(n, pt, rate) {
    cc <- 1
    for (i in 0:n) {
        cc <- cc - dbinom(i, pt, rate)
    }
    return(cc)
}
とすればよいし，
> 1-sum(dbinom(0:3, 200, 0.01))
[1] 0.141966
> sum(dbinom(4:200, 200, 0.01))
[1] 0.141966
でよいし，何よりかんより，pbinom があるのだから，
> pbinom(3, 200, 0.01, lower.tail=FALSE)
[1] 0.141966
> 1-pbinom(3, 200, 0.01)
[1] 0.141966
で十分でしょう。

dbinom(0:200, 200, 0.01) をやって，どこからどこまでを加えたらどうなるかを見てみると悩まなくて済むでしょう。

No.16899　Re: 二項分布による検出力　　【赤岳】　2012/05/10(Thu) 22:47

青木先生，
早速のご教示ありがとうございます。
iは0から指定できるのも，pbinom関数があるのも，sumを用いるのも，初めて知りました。
これだけのコマンドでしたが，今日半日かかっていました。
ほんとにありがとうございました。