Bayesian Bootstrapをjuliaでやってみる
概要
先日Rのパッケージを使ってBayesian Bootstrapをやってみました。
今回はよりBayesian Bootstrapの理解を深めるためjuliaで、パッケージを使わずBayesian Bootstrap的なことをやってみたいと思います。
実験
まずは必要なパッケージとデータを読み込みます。データはirisを使います。
using Distributions using StatsBase using RDatasets using DataFrames using Plots using Random iris = dataset("datasets", "iris") iris = iris[shuffle(1:end), :] #ランダムに並び替える
次に関数を定義していきます。
まずはサンプリング用の関数です。この関数にはデータとパラメータ、サンプルの事前分布を適当に与えそれを重みとしてブートストラップサンプリングを実施します。返り値としてはデータのインデックスになるので、この返り値を活用して統計量を計算に使う設計です。
n1は重み付けのパターンを何個生成するか、n2は統計量を計算する際に使うサンプルサイズ、dでサンプルの事前分布を指定します。
function bayesian_bootstrap_sampling(data, n1, n2, d; index_mode::Bool = false) if(index_mode) sample_pool_data = 1:size(data)[1] sampled_data = Int64.(zeros(n1, n2)) else sample_pool_data = copy(data) sampled_data = zeros(n1, n2) end for i = 1:n1 weight = rand(d, size(data)[1]) weight = pweights(weight / sum(weight)) sampled_data[i, :] = sample(sample_pool_data, weight, n2, replace = true, ordered = false) end return sampled_data end
次に先程のサンプリング関数の返り値を使って、指定した統計量を求める関数です。
function bayesian_bootstrap_statistics(data, sample, val, f) ret = zeros(n1) for i = 1:size(sample)[1] ret[i] = f(data[sample[i, :], val]) end return ret end
それではこれまでの関数を使って結果を確認してみます。まずは1000個の重み付きのパターンを生成し、各パターンの中から5個サンプリングします。事前分布はGamma(1, 1)でn次元の一様ディリクレ分布と一致する分布にしました*1。データは50サンプルの中からサンプリングします。
n1 = 1000 n2 = 5 d = Gamma(1, 1) index_sample = bayesian_bootstrap_sampling(iris[1:50, :], n1, n2, d, index_mode = true)
これでサンプリングができました。この関数の返り値は以下のように重み付けのパターン数(n1)×統計量に使うサンプル数(n2)だけのデータのインデックスになります。
1000×5 Array{Int64,2}:
9 17 17 6 40
8 42 39 15 20
9 24 31 41 22
21 2 2 39 1
5 46 47 35 28
8 33 26 7 26
49 32 23 34 33
27 31 18 15 13
43 31 4 4 50
41 17 26 40 22
9 23 30 42 49
28 10 45 8 28
42 46 27 33 8
⋮
10 39 48 49 12
21 17 12 26 39
11 21 18 44 41
44 40 9 47 4
27 45 4 34 45
5 27 40 36 42
36 45 11 19 33
28 17 10 35 8
31 34 47 34 34
49 6 36 21 4
38 6 3 4 33
5 43 29 49 42
統計量を求めてみましょう。既に定義したbayesian_bootstrap_statisticsを使ってmeanを求めます。
bayesian_bootstrap_statistics(iris, index_sample, :SepalLength, mean)
3つ目の引数で統計量を求める対称の変数を指定し、4つ目の引数で求めたい統計量を指定しています。
結果は以下のようになりました。
1000-element Array{Float64,1}:
5.3
6.26
6.319999999999999
5.839999999999999
6.180000000000001
6.4
6.06
5.88
5.959999999999999
6.34
5.659999999999999
6.58
6.08
⋮
6.44
6.199999999999999
6.2
5.92
6.2
6.14
6.300000000000001
6.640000000000001
6.779999999999999
6.5200000000000005
6.16
6.42
重み付けのパターン分の統計量が求まりました。更にこの平均を見てみましょう。
mean(bayesian_bootstrap_statistics(iris, index_sample, :SepalLength, mean))
6.150060000000001となって全データで平均を出した場合の5.843333333333333よりは大きい値になっています。可視化をすると、
histogram(bayesian_bootstrap_statistics(iris, index_sample, :SepalLength, mean))
このような結果です。
ついでにmedianで可視化してみましょう。最後の引数だけmedianに変わって、
histogram(bayesian_bootstrap_statistics(iris, index_sample, :SepalLength, median))
以上の結果が出ました。
事前分布を変えた場合どのようになるかを最後に試してみましょう*2。適当に分布をGamma(10, 100)にしました*3。
n1 = 1000 n2 = 5 d = Gamma(10, 100) index_sample = bayesian_bootstrap_sampling(iris[1:50, :], n1, n2, d, index_mode = true) histogram(bayesian_bootstrap_statistics(iris, index_sample, :SepalLength, mean))
以下のような分布になり、
mean(bayesian_bootstrap_statistics(iris, index_sample, :SepalLength, mean))
分布の平均は6.134599999999999でした。あまりGamma(1, 1)と変わらない結果がでました。
まとめ
本記事では改めてBayesian Bootstrapをjuliaで試してみました。
自分の理解ではBayesian Bootstrapはサンプリングする際に事前の分布を設定するのが通常のBootstrapとの違いかなという理解です。分布自体は主観性が入るのを許容できるのであれば何でもいいのかなと思います。ただし、一様ディリクレ分布*4以外を使う理由が発生することはあまりないと思います。
今後はサンプリングしたデータを使って重回帰などを行ってパラメータの分布を取るということをやってみたいです。今回はver 1.0.0のjuliaでGLMパッケージが利用できずそこまでできませんでした。
それでは間違い等ありましたらご指摘お願いいたします。
*1:Easy Bayesian Bootstrap in R - Publishable Stuff
*2:現実的に事前分布を変更したいと思うことはなかなかないと思いますが
*3:Gamma分布以外ももちろん使えます
*4: Gamma(1, 1)
