Deep Karmaning

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

Bayesian Bootstrapをjuliaでやってみる

概要

先日Rのパッケージを使ってBayesian Bootstrapをやってみました。

rf00.hatenablog.com

今回はより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))

f:id:rf00:20180924120410p:plain

このような結果です。

ついでにmedianで可視化してみましょう。最後の引数だけmedianに変わって、

histogram(bayesian_bootstrap_statistics(iris, index_sample, :SepalLength, median))

f:id:rf00:20180924121112p:plain

以上の結果が出ました。

事前分布を変えた場合どのようになるかを最後に試してみましょう*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))

以下のような分布になり、

f:id:rf00:20180924121647p:plain

mean(bayesian_bootstrap_statistics(iris, index_sample, :SepalLength, mean))

分布の平均は6.134599999999999でした。あまりGamma(1, 1)と変わらない結果がでました。

まとめ

本記事では改めてBayesian Bootstrapをjuliaで試してみました。

自分の理解ではBayesian Bootstrapはサンプリングする際に事前の分布を設定するのが通常のBootstrapとの違いかなという理解です。分布自体は主観性が入るのを考慮すれば何でもいいのかなと思う次第です。ただし、一様ディリクレ分布((Gamma(1, 1)))以外を使う理由はあまりなさそうだとは思います。

今後はサンプリングしたデータを使って重回帰などを行ってパラメータの分布を取るということをやってみたいです。今回はver 1.0.0のjuliaでGLMパッケージが利用できずそこまでできませんでした。

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

*1:Easy Bayesian Bootstrap in R - Publishable Stuff

*2:現実的に事前分布を変更したいと思うことはなかなかないと思いますが

*3:Gamma分布以外ももちろん使えます