餡子付゛録゛

ソフトウェア開発ツールの便利な使い方を紹介。

信頼区間がマイナスに入らない二項分布の区間推定

感染症の罹患経験など、無作為抽出のyes/noの観測値は二項分布に従うことが知られていますが、その場合も正規分布を用いて区間推定を行うことが多いです*1中心極限定理正規分布に漸近することを利用しているのですが、稀にしか観測されない事象で、全体の生起数が10にも満たない場合は不適切な推定になりがちです。区間の片方の端がマイナスになります。

例えば厚生労働省新型コロナウイルスへの抗体保有率の調査の場合、東京都で観測数が1971名、抗体保有者が2名となり、正規分布から区間推定を行うと0.0010(95%信頼区間-0.00039~0.0024)と言う結果になります。2項分布の生起確率でマイナスはあり得ないので、これはちょっとまずいです。プラスの区間になるように、計算精度を上げましょう。

ベイズの定理を使って、事前分布と尤度関数の積を1に正規化して事後確率を求める方法で、信頼区間を計算するための確率分布をつくります*2。生起確率を求めるときの2項分布の尤度関数は、観測数と生起数を固定した2項分布の確率密度関数となることに注意すると、計算機にやらせる分にはそんなに難しくない作業になります。

n <- 1971 # 観測数
m <- 2 # 生起数

p <- m/n # 生起確率
v <- n*p*(1-p) # 観測数の分散

# 正規化するための尤度関数の積分
s <- integrate(function(p){
    # pを求めるときの2項分布の尤度関数は、観測数と生起数を固定した2項分布の確率密度関数となる
    dbinom(m, n, p)
  }, 0, 1, rel.tol=1e-15)$value

# 確率密度関数
pdf <- function(p){
  dbinom(m, n, p)/s
}

# 累積分布関数
cdf <- function(ps){
  sapply(ps, function(p){
    integrate(function(a){
      pdf(a)
    }, 0, p)$value
  })
}

# 閾値(2分探索で計算)
threshold <- function(a){
  l <- 0
  h <- 1
  m <- 0.5
  for(i in 1:100){
    m <- (l + h)/2
    v <- cdf(m) - a
    if(v == 0){
      break;
    } else if(v < 0){
      l <- m
    } else {
      h <- m
    }
  }
  m
}

# 95%信頼区間とする
a <- 0.025
sprintf("%.4f(%.0f%%信頼区間%.5f~%.5f)", p, 100*(1-2*a), threshold(a), threshold(1-a))

f:id:uncorrelated:20200617045800p:plain

0.0010(95%信頼区間0.00031~0.00366)と計算できました。なお、正規分布からの区間推定は以下で計算できるので、mが50ぐらいのときはほぼ変化なしな事を確認すると、安心して使いやすいと思います。

se <- sqrt(v)/n # 生起確率の標準誤差
rng <- se*qnorm(a, lower.tail=FALSE)
sprintf("%.4f(%.0f%%信頼区間%.5f~%.5f)", p, 100*(1-2*a), p - rng, p + rng)

*1: 計算機の利用が一般的でなかった頃の教科書で、付属している表から紙と鉛筆で計算できる方法が紹介されたためだと思います。

*2:統計哲学的には信用区間と言った方がよいかもです。