Deep Karmaning

技術系の話から日常のことまで色々と書きます

Bayesian Bootstrapを試してみる

Bayes R

概要

twitterを見てみたらBayesian Bootstrapなるものが有用ということを知りました。

とはいえ初めて聞いてどんなものかよくわかっていないのですが、 色々とネットを検索したところRのbayesbootパッケージなるものを見つけたので今回は試してみます。

github.com

なお本記事は自分自身も理解が十分ではない部分がありますので、間違い等はあるかもしれませんがよろしくお願いいたします。

Bayesian Bootstrapとは

そもそもBayesian Bootstrapとは何なのでしょうか。

どうやら因果推論で有名なRubinによって1981年に紹介された手法のようで、論文は以下になります。

projecteuclid.org

正直軽く読んでも全くよくわからなかったので、bayesbootの作者の解説記事を見てみます。

www.sumsar.net

version oneの部分を主に参考にしたいのですが、Byesian Bootstrapのモデルとしては以下のようになるようです。

f:id:rf00:20180904150828p:plain

dは取りうるすべての値のベクトル（離散か連続）、 xは実際に得られたデータ、 πはdが得られる確率となるようです。

ここでd=(0,1)のときのBayesian Bootstrapは、x_i ~ Bernoulli(p)、pをp∼Beta(0,0) とするのと同じになります。

dが離散値や整数のカウントデータの場合は、各値が出る確率を考えることに違和感はないと思いますが。 dが連続値のときは違和感のある過程を置く必要があります。つまり3.11が出る確率はとか3.14が出る確率はみたいなことを仮定していく必要があるため奇妙な感じがします*1。

モデルとして以上のような形になるのですが、実際上bayesbootパッケージ内でBayesian Bootstrapはどのような手順で行っているのかを見ていきたいと思います。

www.sumsar.net

bayesbootの作者の実装上の手順としては以下のようです。

  • データ数と同じ次元のuniform Dirichlet分布から重みを生成する

uniform Dirichlet分布はGamma(1, 1)から生成することができ、Exponential(1) と同様になるとのことです。

例としてはこのように生成するようです。nがサンプルサイズでn1がリサンプリング数で、最後に合計が1になるように正規化していますね。

dirichlet_sample <- matrix( rexp(n * n1, 1) , ncol = n, byrow = TRUE)
dirichlet_sample <- dirichlet_sample / rowSums(dirichlet_sample)
  • 生成された重みに基づいてリサンプリングして統計量を計算する

先程のuniform Dirichlet分布で各サンプル毎に重みがついたので、この重みを考慮してリサンプリングを行い、そこで得られたサンプルで統計量を計算するということを行います。

例えばガンマ分布に基づいて生成されたダミーデータをもとにBayesian Bootstrapで中央値を求めるならばこんな感じです。

x <- rgamma(10, shape = 2, scale = 1)
n <- length(x)
n1 <- 3000
n2 <- 1000

dirichlet_sample <- matrix( rexp(n * n1, 1) , ncol = n, byrow = TRUE)
dirichlet_sample <- dirichlet_sample / rowSums(dirichlet_sample)

bb_median <- rep(NA, n1)
for(i in 1:n1) {
  data_sample <- sample(x, size = n2, replace = TRUE, prob = dirichlet_sample[i,])
  bb_median[i] <- median(data_sample)
}
quantile(bb_median, c(0.025, 0.975))

ここでやっているのは、元データに3000パターン(n1)の重みの付け方をして、1パターンごとに1000個サンプリング(n2)して中央値を求める、 それを全パターン同様に繰り返し中央値を求めその信頼区間を出すというようなイメージです。

実験

今までBayesian Bootstrapに関して見てきたのですが、実際、通常のBootstrapよりもどれくらいいいのかはよくわからないところです*2。

そこで今回は実験を実施してBayesian Bootstrapと通常のBootstrapの比較をしたいと思います。

実験の方針としては以下の記事を真似します。

hoxo-m.hatenablog.com

実験設定

それでは実験の設定を説明します。

元データのサンプルサイズを変化させながら通常のBootstrapとBayesian Bootstrapの期待値の信頼区間*3を求め、その信頼区間の間に真の期待値の値がどの程度入るかを検証していきます。

信頼区間は各サンプルサイズごとに100回繰り返します。つまり100回Bootstrap信頼区間を計算し、その100個の信頼区間の中で何回真の期待値が入ったかを比較することになります。

元データのサンプルサイズとしては、5, 10, 20, 50, 100, 100を用意しました。

元データの分布はガンマ分布で生成しました。パラメータは形状母数が2、尺度母数が1で期待値が2になります*4。

今回実験に使った全コードは以下です。

bayesian_bootstrap_experiments.r · GitHub

実験結果

実験の結果としては以下になりました。

f:id:rf00:20180904165432p:plain

横軸は元データのサンプルサイズで、縦軸は100回の試行回数のうちどの程度が信頼区間内に真の値が入っているかを示しています。赤がBayesian Bootstrapで、青が通常のBootstrapです。

これを見ると、サンプルサイズがかなり小さい時(n= 5, 10, 20)は通常のBootstrapが精度が高く、サンプルサイズが一定を超えたところでBayesian Bootstrapの精度が上回っています。

Bayesian Bootstrapが精度を上回ったときのサンプルサイズも50と小さいので、小さいときも有効であることを示唆していますが、小さすぎると駄目なようですね。

また十分に大きいサンプルサイズでもBayesian Bootstrapのほうが精度が良いので、小さすぎないとき以外*5はBayesian Bootstrapが有効なのかもしれません。

まとめ

今回はRのbayesbootパッケージを使ってでBayesian Bootstrapを試してみました。

実験結果からは通常のBootstrapに比べてBayesian Bootstrapは元データのサンプルサイズが小さ過ぎるときはよくないですが、 それでも十分に小さいサンプルサイズでも通常のBootstrapよりは精度が出ることが検証できました。

個人的には今後Bayesian Bootstrapは活用していきたいなと思いました。

それでは間違い等ありましたら指摘お願いいたします。

*1:とはいえ記事中には“all data as observed are discrete”ということが書かれていて、捉え方を変えると納得できる部分もあります

*2:JMPという統計ソフトウェアの解説ではサンプルサイズが小さくて不確実性が高いときが有効という解説もありますhttp://www.jmp.com/japan/support/help/13/ba-bootstrapping-2.shtml

*3:信頼区間としてはいろいろな種類がありますが、今回はパーセンタイル法を利用

*4:https://ja.wikipedia.org/wiki/%E3%82%AC%E3%83%B3%E3%83%9E%E5%88%86%E5%B8%83を参考のこと

*5:小さすぎないときがどの程度かは難しいですね。。。