ベイズマックスへようこそ。
このブログでは,心理学と統計が大好きなmutopsyが,
統計に関するあれやこれやを自由気ままに投稿します。

2019年01月01日

既知の確率質量関数からMAP推定値と最高密度区間 (HDI) を計算するためのRの自作関数

 あけましておめでとうございます。mutopsyです。年末にノロウイルスらしきものに感染してしまいましたがぶじ元気に2019年を迎えることができました。今年も何卒よろしくお願いいたします。新年早々ではありますが,必要に駆られてRの自作関数を作ったので共有します。

 既に分かっている確率質量関数を使ってMAP (最大事後確率) 推定値とHDI (highest density interval; 最高密度区間) を求める関数です。実現値として離散値kを返す既知の確率質量関数 (e.g., 二項分布) で事後分布が表現できる場合に,そのような事後分布のMAPとHDIを計算するのに使えます。ちなみに連続値を返す確率密度関数の場合でも,離散化すればこの関数を使って近似的にMAP推定値とHDIを求めることができるはずです。

 2019年1月2日追記:MAPとなる値やHDIの範囲が2つ以上ある場合 (例えば2峰分布など) にも対応できるように関数を修正しました。結果の出力の仕方も微修正。

Rコード01: 既知の確率質量関数からMAP推定値とHDIを計算する関数
hdimass <- function(prob.k,start.k=0,credMass=0.95){
  #警告メッセージ
  if(sum(prob.k)>1) print("警告:確率の和が1を越えています。")
  if(credMass<0|credMass>1) print("警告:credMassの値が正しくありません。")
  
  #HDIの計算の下準備
  credMass1m <- 1 - credMass
  accum<-1:length(prob.k)
  for(k in 1:length(prob.k)){
    if(k ==1){
      accum[k] <- sort(prob.k)[k]
    }else{
      accum[k] <- accum[k-1] + sort(prob.k)[k]
    }
  }
  Interval <- grep(TRUE,prob.k>max(sort(prob.k)[accum<=(min(abs(accum-credMass1m))+credMass1m)]))
  
  #HDIが複数ある場合やHDIが存在しない場合 (全範囲がHDI) を考慮
  if (length(Interval) <= 1){
    Check <- 0
    Interval <- 1:length(prob.k)
  } else{
    Check <- Interval[2:length(Interval)]-Interval[1:(length(Interval)-1)]-1
  }
  if(sum(Check)==0){
    Lower <- min(Interval) + (start.k-1)
    Upper <- max(Interval) + (start.k-1)
  } else{
    Lower <- c(min(Interval),Interval[grep(TRUE,Check!=0)+1]) + (start.k-1)
    Upper <- c(Interval[grep(TRUE,Check!=0)],max(Interval)) + (start.k-1)
  }
  
  #実際にその範囲の値となる確率
  ActualProb <- sum(prob.k[Interval])
  
  #MAPの計算
  MAP <- grep(TRUE,prob.k == max(prob.k)) + (start.k-1)
  
  #結果をリストにまとめて返す
  Result <- list(MAP=MAP,Lower=Lower,Upper=Upper,ActualProb=ActualProb)
  return(Result)
}

 第1引数 (prob.k) としてそれぞれのkが実現する確率のベクトルを,第2引数 (start.k) としてkの最小値 (デフォルトはstart.k=0) を与えることで,MAPと95% HDIの上限・下限を格納したリストが返り値として得られます。また,実際にHDIの範囲内の値をとる確率も返してくれます (離散分布なので95% HDIであっても厳密には95%にはなりませんが,95%に最も近い値になるはずです)。第1引数 (prob.k) には現実的な長さのベクトルを与えれば十分です。例えば N = 20,θ = 0.5 の二項分布であれば,長さが20もあれば総和がほぼ1になるので計算上は問題ありません (以下の例では一貫して長さ500のベクトルを与えていますがもっと少なくても結果は変わりません)。ちなみに第3引数であるcreadMassに(0,1)の範囲の値を与えれば (例えばcredMass = 0.99),異なる確度のHDIを出力することもできます (99% HDIなど)。

 計算の中身を簡単に説明します。まず最初にfor文で確率質量 (prob.k)を昇順にソートし,順番に累積確率を計算しています。次に,累積確率が閾値 (95% HDIなら5%)に最も近い点を探し,各々のkがHDIに含まれるか否かを論理値ベクトルで表現します。この論理値ベクトルを使って,(credMass)% HDIに含まれるkの値を抽出し,Intervalという内部の変数に格納しています。HDIが1つしかない場合は,Intervalの最小値と最大値ををHDIの下限・上限とみなしています。HDIが2つ以上ある場合は,Interval内で数値が飛び飛びになっている部分を見つけてその変化点をHDIの境界とみなしています。MAP推定値は最大事後確率を実現するkの値を求めることで計算しています。

 実際にこの関数を使って,さまざまなパラメータを持つ二項分布のMAPとHDIを計算してみました。

Rコード02: パラメータと引数を変えながら二項分布で試してみた結果
> #(1) 成功確率0.5,試行回数20の二項分布 (95% HDI)
> hdimass(dbinom(0:500,20,0.5),start.k = 0)
$MAP
[1] 10

$Lower
[1] 6

$Upper
[1] 14

$ActualProb
[1] 0.9586105

> #(2) 成功確率0.2,試行回数20の二項分布 (95% HDI)
> hdimass(dbinom(0:500,20,0.2),start.k = 0)
$MAP
[1] 4

$Lower
[1] 1

$Upper
[1] 7

$ActualProb
[1] 0.9563281

> #(3) 成功確率0.6,試行回数100の二項分布 (95% HDI)
> hdimass(dbinom(0:500,100,0.6),start.k = 0)
$MAP
[1] 60

$Lower
[1] 51

$Upper
[1] 69

$ActualProb
[1] 0.948118

> #(4) 成功確率0.6,試行回数100の二項分布 (99% HDI)
> hdimass(dbinom(0:500,100,0.6),start.k = 0,credMass = 0.99)
$MAP
[1] 60

$Lower
[1] 48

$Upper
[1] 72

$ActualProb
[1] 0.9896389

> #(5) 成功確率0.5,試行回数20,k>=7の切断二項分布 (95% HDI)
> hdimass(dbinom(7:500,20,0.5)/sum(dbinom(7:500,20,0.5)),start.k = 7,credMass = 0.95)
$MAP
[1] 10

$Lower
[1] 7

$Upper
[1] 13

$ActualProb
[1] 0.9388129

> #(6) MAPとHDIが2つずつある混合二項分布 (95% HDI)
> hdimass(dbinom(0:500,100,0.3) * 0.5 + dbinom(0:500,100,0.7) * 0.5)
$MAP
[1] 30 70

$Lower
[1] 22 61

$Upper
[1] 39 78

$ActualProb
[1] 0.9501802
#↑95% HDIは[22,39]または[61,78]の範囲

 いい感じに計算できていると思います。

 最後に,この関数を作った経緯を少しだけ説明します。僕が今行っている分析で離散パラメータの推定を行っているのですが,MCMCの結果からMAP推定値とHDIを十分高い精度で計算するためには発生させるMCMCサンプルの数 (JAGSのn.iterの値) をかなり増やさなければならないことが分かりました。これは現実的ではなかったので,数値積分を利用して事後分布を近似的に計算し,その事後分布からMAP推定値とHDIを計算することにしました。その副産物として今回の関数ができたという訳です。Rを使って事後分布を近似的に求める方法についてもいつか記事で解説するかもしれません。

 という訳で,需要があるかどうかも怪しい誰得関数の紹介から2019年が始まりました。Enjoy!

【Rに関するTipsの最新記事】
posted by mutopsy at 16:52 | Rに関するTips