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

2019年01月09日

二項分布モデルと幾何分布モデルで推定した確率パラメータの事後分布の比較

 こんにちは。mutopsyです。この記事では,同一のカウントデータを二項分布でモデリングした場合と幾何分布でモデリングしたときに,推定される確率パラメータの事後分布が一致することを説明します。間違いがあったらご指摘ください。

 次の2つの例を考えてみましょう。

  1. N箱のチョコボールを購入したら銀のエンゼルが1枚当たった。銀のエンゼルが当たる確率θ1を推定せよ。
  2. 銀のエンゼルが当たるまでチョコボールを購入したらN箱目で当たった。銀のエンゼルが当たる確率θ2を推定せよ。

θ1θ2はどちらも試行にかかわらず一定の値をとり,かつ各試行は独立であると仮定します。ここではチョコボールを例に挙げていますが,2つの事象 (e.g., 当たるか外れるか) のカウントデータであれば何でも構いません。どちらの文も,「銀のエンゼルがN箱中1箱当たった」という同じ結果を述べていますが,データの収集方法に違いがあります。1番目のケースではNがあらかじめ決まっていた状況で銀のエンゼルが1回当たったという結果が得られていますが,2番目のケースでは銀のエンゼルが1回当たることがあらかじめ決まっていた状況でN箱を要したという結果が得られているという違いです。この違いは確率パラメータの推定結果に影響を与えるのでしょうか。

 この例のようにNの決め方 (決まり方) が異なるケースでは,異なる統計モデルを想定するのが自然です。1番目の例のように,Nが所与のときに銀のエンゼルが当たった枚数 (K = 1) のデータが得られた場合には,二項分布を用いて

1 ~ Binomial(N, θ1)

と表現することができます。一方,2番目の例のように,銀のエンゼルが当たった枚数 (K = 1) が所与のときにNというデータが得られた場合には,幾何分布を用いて

N ~ Geometric(θ2)

と表現できます。幾何分布は,同一の確率パラメータを持つベルヌーイ分布に従う試行を独立に繰り返し行ったときに初めて成功する (今回の例では銀のエンゼルが当たる) のに要する回数の分布です。二項分布では銀のエンゼルが当たった回数 (K = 1) が説明されるのに対し,幾何分布では銀のエンゼルが当たるまでに要した試行回数 (N) が説明されるのがポイントです。

 では,2つのモデルでθ1θ2の推定結果がどうなるのかを確認してみましょう。今回は,N = {10, 20, 50, 100} の場合のθ1θ2をStanでベイズ推定してみます。4通りのN と2つのモデル (Is_Censored = 0なら二項分布モデル,Is_Censored = 1なら幾何分布モデル) の組み合わせを用いて8個のθをまとめて推定します。Stanコードは次の通り。rstan2.18.2,Rtools3.5,R3.5.2,Windows10で実行しました。

Stanコード01: 二項分布モデルと幾何分布モデルの比較
data {
  int N_ID; //標本の数 = 8
  int N[N_ID]; //試行回数 = {10,20,50,100,10,20,50,100}
  int K[N_ID]; //成功回数 (今回は1で固定;2以上の場合は負の二項分布になる)
  int Is_Censored[N_ID]; //打ち切られたデータか否か = {0,0,0,0,1,1,1,1}
}

parameters{
  real<lower=0,upper=1> theta[N_ID]; //当たる確率
}

model{
  for (i in 1:N_ID){
    if (Is_Censored[i] == 0){
      K[i] ~ binomial(N[i],theta[i]); //二項分布
    } else if (Is_Censored[i] == 1){
      N[i]-K[i] ~ neg_binomial(K[i],theta[i]/(1-theta[i])); //負の二項分布 (K = 1のとき幾何分布と一致)
    }
  }
}

幾何分布は負の二項分布の特殊なケースなのでモデル上は負の二項分布を使って表現しています (現時点ではStanには幾何分布用の関数はありません)。neg_binomialの第2引数にはθそのものではなくθのオッズを与える必要があります。また,左辺にはNではなくN-1を与えます (Stanで使われている負の二項分布の定義に合わせるため)。このモデルを走らせるためのRコードは以下の通り。

Rコード01: MCMCの実行
library(rstan)
N_ID <- 8
N <- c(10,20,50,100,10,20,50,100)
K <- rep(1,N_ID)
Is_Censored <- c(0,0,0,0,1,1,1,1)

#コンパイル
stanmodel <- stan_model(file="binomgeom.stan")

#データの準備
data <- list(N_ID = N_ID, N = N,
             K = K, Is_Censored = Is_Censored)

#サンプリング
fit <- sampling(stanmodel, data=data, seed=123,
                 chains=4, iter=250500, warmup=500, thin=1)

 各々の場合のθの事後分布とMAP推定値は以下の通り。

posterior_binomgeom.png
図1. 幾何分布モデルと二項分布モデルで推定したθの事後分布とそのMAP推定値。

二項分布モデルで推定したθ1を半透明の赤色,幾何分布モデルで推定したθ2を半透明の青色で描画していますが,ほとんど重なっていて紫色にしか見えません。MAPも,多少の誤差はありますがほぼ一致しています。どうやら,どちらのモデルを使っても推定結果は一致ししそうです。

 両者が一致することは,事後分布を解析的に求めれば分かります。まず,二項分布の確率質量関数は,Nを所与として

posterior_binomgeom.png

と表されます。K = 1のとき,

posterior_binomgeom.png

となります。データが1つしかないときのθ1の尤度関数は

posterior_binomgeom.png

と表せます。最初のNθ1とは無関係な定数なので消去しています。最後に,θ1について無情報事前分布を仮定して尤度関数を正規化することで,次の事後分布が得られます。

posterior_binomgeom.png

 次に幾何分布ですが,確率質量関数は

posterior_binomgeom.png

と表されます。二項分布モデルの場合と同じように,これをθ2の関数 (尤度関数) とみなして正規化すると,

posterior_binomgeom.png

という事後分布が得られます。先ほどの二項分布モデルの事後分布と一致していることが確認できると思います。実はこの事後分布はどちらも (α, β) = (2, N) のベータ分布と一致します。二項分布と幾何分布の確率質量関数の違いは組み合わせ (C) を掛けるか否かの違いのみで,この部分は確率パラメータに依存しないため,事後分布を計算する過程で消去されるという訳です。ちなみに,データが複数ある場合やKが2以上の場合 (つまり幾何分布ではなく負の二項分布を使う場合) にも事後分布は同じパラメータを持つベータ分布 (α = 1 + 成功回数,β = 1 + 失敗回数) となり一致します (計算は省略)。

 という訳で,二項分布モデルを使っても幾何分布モデルを使っても確率パラメータの推定結果は一致することが分かりました。とはいえ,データの持つ情報 (データの収集方法の違い) を正確に表現にするためにもモデル上は両者を区別したほうが良さそうです。

posted by mutopsy at 15:22 | 統計

2019年01月04日

StanやJAGSの計算結果からパラメータのMAP推定値とHDIをまとめて得るためのRの自作関数

 こんにちは。mutopsyです。この記事では,StanやJAGSのfitオブジェクトを渡すとパラメータのMAP推定値とHDI (最高密度区間) をまとめて返してくれる自作関数を紹介します。

 StanやJAGSでパラメータの推定結果の要約を出力すると,通常は事後分布のEAP (i.e., 平均) や分位点の値が出力されますが,事後分布のMAP推定値やHDIが知りたい場合もありますよね。そんなわけで,そういう関数を作ってみました (もう既にあるかもしれないけど知ったこっちゃない)。自作関数といっても,HDIの計算にはHDIntervalというパッケージを利用しますが。関数は以下の通り。

Rコード01: MCMCサンプルからMAP推定値とHDIを計算する関数
library(HDInterval)
maphdi <- function(fit,credMass=0.99,include=TRUE,pars=NA){
  if(credMass<0|credMass>1) print("警告:credMassの値が正しくありません。")
  Class <- class(fit)[1]
  if(Class=="stanfit"){
    e <- as.data.frame(rstan::extract(fit))
  }else if(Class=="rjags"){
    e <- as.data.frame(fit$BUGSoutput$sims.list)
  }else if(Class=="data.frame"){
    e <- fit
  }else{
    print("エラー:未対応なオブジェクトです。")
    stop()
  }

  if(is.na(pars[1])) pars <- colnames(e)
  if (include==TRUE){
    e <- e[,pars]
  }else{
    e <- e[,!colnames(e) %in% pars]
  }
  if(length(pars)==1){
    MAP <- density(e)$x[which.max(density(e)$y)]
    HDI95 <- hdi(e)
    HDIAny <- hdi(e,credMass = credMass)
    Result <- data.frame(MAP = MAP, Lower95 = HDI95[1], Upper95 = HDI95[2],
                         LowerAny = HDIAny[1],UpperAny = HDIAny[2])
    rownames(Result) <- pars
  }else{
    MAP <- apply(e,2,function(z) density(z)$x[which.max(density(z)$y)])
    HDI95 <- apply(e,2,HDInterval::hdi)
    HDIAny <- apply(e,2,function(x) HDInterval::hdi(x,credMass = credMass))
    Result <- data.frame(MAP = MAP, Lower95 = HDI95[1,], Upper95 = HDI95[2,],
                         LowerAny = HDIAny[1,],UpperAny = HDIAny[2,])
  }

  colnames(Result)[4:5] <- c(
    paste("Lower",round(credMass*100,2),sep=""),
    paste("Upper",round(credMass*100,2),sep="")
  )
  return(Result)
  }

 試しにRからStanのfitオブジェクトを渡してみます。データとモデルは何でもよいので適当なものを渡します。

Rコード02: 関数を使ってみる
>fit.eg <- sampling(stanmodel, data=data, seed=123,
                 chains=4, iter=10500, warmup=500, thin=1)
> maphdi(fit.eg)
              MAP    Lower95    Upper95    Lower99    Upper99
beta.1  1.8502036  1.3540767  2.3442492  1.1987086  2.5187496
beta.2  0.6478270  0.5188389  0.7805164  0.4767612  0.8295815
beta.3  0.7103451  0.5996639  0.8201469  0.5627496  0.8524817
beta.4 -0.5622637 -0.7994808 -0.3030058 -0.8891769 -0.2392025
sigma   0.3135227  0.2814639  0.3548575  0.2728315  0.3686322
lp__   97.7211912 93.6796922 99.1662181 91.7998798 99.3115815

 こんな風に,Stanのfitオブジェクトをmaphdi()に引数として渡すだけで各パラメータのMAP推定値と95% HDI・99% HDIを計算して出力することができます。Stan2.18, rstan2.18.2, Rtools3.5, R3.5.2, Windows10で動作確認済み。JAGSにも対応しています (JAGS4.3.0, R2jags0.5-7で確認済み)。

 また,fitオブジェクトの代わりに,MCMCサンプルが格納されたデータフレームを渡してもMAP推定値とHDIを返してくれます。

Rコード03: MCMCサンプルを直接渡す
> e <- as.data.frame(rstan::extract(fit.eg))
> head(e)
    beta.1    beta.2    beta.3     beta.4     sigma     lp__
1 2.194799 0.5785987 0.6255656 -0.4086454 0.3094092 97.46736
2 1.726557 0.6609587 0.7787883 -0.7231942 0.3058624 97.31513
3 2.053161 0.6021407 0.6534262 -0.3992853 0.3425830 96.40984
4 1.941408 0.6503287 0.6399921 -0.4178822 0.2976153 97.93040
5 2.518750 0.4448252 0.7676119 -0.7428853 0.3209268 91.90027
6 1.758346 0.6890341 0.6522970 -0.4046868 0.2985565 97.78386
> maphdi(e)
              MAP    Lower95    Upper95    Lower99    Upper99
beta.1  1.8502036  1.3540767  2.3442492  1.1987086  2.5187496
beta.2  0.6478270  0.5188389  0.7805164  0.4767612  0.8295815
beta.3  0.7103451  0.5996639  0.8201469  0.5627496  0.8524817
beta.4 -0.5622637 -0.7994808 -0.3030058 -0.8891769 -0.2392025
sigma   0.3135227  0.2814639  0.3548575  0.2728315  0.3686322
lp__   97.7211912 93.6796922 99.1662181 91.7998798 99.3115815

 引数のcredMassに0〜1の範囲の値を与える事で,99% HDIの代わりに任意の確度のHDIを出力することもできます。95% HDIは必ず出力されるようになっています。

Rコード04: 任意のHDIを出力
> maphdi(e,credMass = 0.50)
              MAP    Lower95    Upper95    Lower50    Upper50
beta.1  1.8502036  1.3540767  2.3442492  1.7011924  2.0385575
beta.2  0.6478270  0.5188389  0.7805164  0.6060299  0.6964437
beta.3  0.7103451  0.5996639  0.8201469  0.6729364  0.7490473
beta.4 -0.5622637 -0.7994808 -0.3030058 -0.6429855 -0.4713128
sigma   0.3135227  0.2814639  0.3548575  0.3030061  0.3280369
lp__   97.7211912 93.6796922 99.1662181 96.7857008 98.5912001

 引数のparsにパラメータの名前を与えれば,特定のパラメータの結果のみを出力することもできます。また,include = FALSE とすると,指定したパラメータ以外の結果を出力することもできます。

Rコード05: 一部のパラメータのみ計算

> maphdi(fit.eg,pars=c("beta.1","beta.2"))
            MAP   Lower95   Upper95   Lower99   Upper99
beta.1 1.850204 1.3540767 2.3442492 1.1987086 2.5187496
beta.2 0.647827 0.5188389 0.7805164 0.4767612 0.8295815
> maphdi(fit.eg,pars=c("beta.3"))
             MAP   Lower95   Upper95   Lower99   Upper99
beta.3 0.7103451 0.5996639 0.8201469 0.5627496 0.8524817
> maphdi(fit.eg,include=FALSE,pars=c("beta.1","beta.2"))
              MAP    Lower95    Upper95    Lower99    Upper99
beta.3  0.7103451  0.5996639  0.8201469  0.5627496  0.8524817
beta.4 -0.5622637 -0.7994808 -0.3030058 -0.8891769 -0.2392025
sigma   0.3135227  0.2814639  0.3548575  0.2728315  0.3686322
lp__   97.7211912 93.6796922 99.1662181 91.7998798 99.3115815

 以上です。Enjoy!

posted by mutopsy at 19:10 | Rに関するTips

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!

posted by mutopsy at 16:52 | Rに関するTips