はじめに

Machine Learning Advent Calenderの25日目の記事です。2015年も終わりですね

ネタに悩んでいたのですが，最近社内勉強会でトピックモデル本の輪読をやっている中で，
改めてサンプリングって難しいけど大事だよなぁと感じたので
いろんなサンプリングについて調べてまとめようと思います．

筆者について

Gunosyという会社のデータ分析部というところで，世の中に最適にニュースを届けるためのロジックを構築したり，ユーザ行動の分析をしたりしています．

あと東京大学の博士後期課程で業務の傍ら研究活動もやっています．

専門分野はコンテンツ評価・ユーザ行動分析で，そのために機械学習を用います．ユーザ側です．論文読んで実装したりはします．
統計については機械学習やデータ分析の文脈で学習・利用をしていますが，専門と名乗れるようなものではありません．

本エントリはそのような筆者によるものなので，理論の奥深いところや最新動向などについては触れられませんが，
サンプリングをあまり知らない人や，難しいと感じている方にとって有益なエントリになれば幸いです．

サンプリングとはなにか?

サンプリングとは分布$p(z)$からその分布に従うサンプル$Z^{(l)}=(z_1,.., z_l)$を得ることです
サンプルを大量に得ることができれば，分布から直接計算を行うことが難しい場合にもサンプルから答えをえることができます．
機械学習とサンプリングの関わりとしては，統計的学習があげられます．
統計的学習とは観測データの背後に潜むルールを統計的に記述し，自動的に獲得することです．
例えばトピックモデルは語が生成されるルールに対して背景にその文書のトピックが存在するというルールを記載し，トピックから語がどのように生成されるかというルールを獲得することで，文書のトピックを分析したり，分類を行ったりする統計的学習手法です．
このような場合においてベイズ推定を用いて背後にあるパラメータの事後確率から，予測分布を獲得しますが，その計算を解析的に行うことは多くの場合不可能です．
その場合においてもサンプルを大量に得ることで近似的な予測分布を得ることができるようになります．

シンプルなサンプリングアルゴリズム

サンプリングアルゴリズムとしてはMCMCが有名ですが，まずはサンプリングとはどのようなものかを考えるために棄却サンプリングを紹介します．

棄却サンプリングとは，サンプルを生成したい分布を覆うようなよく知られた分布から生成したサンプルを用いて，サンプリングを行う方法です．
予測分布$p(z)$のサンプルを提案分布$q(z)$から求めることを考えます．
このとき$p(z) < kq(z)$が満たされているものとします

PRML11.2の例を参考にコーシー分布を用いてガンマ分布を推定することを試みてみます．
まず目的のガウス分布を覆うようなコーシー分布を作ってみます．
今回は厳密に行うわけではないので，目視でざっくりです．

可視化にはJuliaを用いています．
$a=3.0, b=1.0$のガンマ分布を覆うコーシー分布は$\mu=2.0, \beta=1.8$に2.0を乗じました．

棄却サンプリングのコードを以下に示します．
（Juliaは最近始めたばかりなのであんまりいいコードではないです)

function rejection_sampling(proposal_distribution::Distribution, target_distribution::Distribution, k::Float64, M::Int)
    y = rand(proposal_distribution, M)
    u = rand(M)
    d = map( (x)->pdf(target_distribution, x) / (k*pdf(proposal_distribution, x)), y)
    s = Float64[]
    for i in 1:M
        if d[i] >= u[i]
            push!(s, k*y[i])
        end
    end
    s
end

棄却サンプリングは非常にシンプルです．

提案分布$q$に従う乱数$y$を生成する
[0,1]の一様分布に従う乱数$u$を生成する
$p(y)/(k*q(y)) > u$だった場合それを$u$を$p(z)$に従うサンプルとする．

棄却サンプリングでやりたいことは
$[0, kq(y)]$に従う一様乱数を生成し，それが$p(y)$を下回る場合にサンプルとして採用することです．
$[0, kq(y)]$に従う一様乱数は確率密度関数$kq(z)$の曲線の下側になります．
その乱数のうち，p(y)より大きい物を破棄することは，予測分布の曲線より上側のサンプルを破棄することです．
$[0, kq(y)]$に従う一様乱数は[0,1]の一様乱数に$kq(z)$を乗じるのと同義のため，
3のような条件で得ることができます．

上記の関数を実行すると以下の様な結果になりました．

10000回試行して3010個の乱数が得られたので，7割程度のサンプルは破棄されたといえます．

ヒストグラムで可視化してみます．

先ほど示したガンマ分布の確率密度関数と同じような曲線を描いているのがわかります．
このように異なる分布から求めたい分布のサンプルを得ることができます．

Markov chain Monte Carlo methods(MCMC)

日本語ではマルコフ連鎖モンテカルロ法と呼ばれます．
確率分布のマルコフ連鎖を元にサンプリングを行うアルゴリズムの総称です．
先ほど紹介した棄却サンプリングは生成するために多くの乱数を棄却しなくてはなりません．
また覆うような分布を作ることは変数が1つであれば作ることは容易ですが，変数が増えるほど生成することが困難担っていきます．
多次元の場合にも効率よく生成できるのがMCMCです．
本エントリでは代表的なMetropolis-Hastingsについて述べます．
よく知られているギブスサンプリングやHMCなども基本的には同じ考え方で作られているため，
こちらを理解することでぐっと理解しやすくなるはずです．

Metropolis Hastings

MCMCがなぜ優れているかを知るために最もシンプルなMetopolis-Hastingsを見てみましょう．
よく話に出てくるギブスサンプリングはこの特殊系であるといえます．
Metropolis Hastingにおいても提案分布を用いてサンプルが良い値であればその値を採用し，そうでなければ採用しません．
この良い値かどうかというのが棄却サンプリングでは確率密度関数の下側にいるかどうかでしたが，
Metropolis-Hastingsでは前に生成されたサンプルの値がどれだけ生成されやすい値かを元に，
次に生成された値がより生成されやすい値であればあるほど受理しやすく，
生成されにくい値であれば受理されにくくなります．（この受理されるか否かを決める確率を受理確率といいます）
このように前に生成された値を用いて決めるので，マルコフ連鎖モンテカルロといいます．
ちなみにモンテカルロとはモナコの首都でギャンブルが盛んな地域であり，
乱数を用いた手法にはよく出てくる名前です．

あるステップ$\tau$において，$z^{\tau}$という値が生成されたとき，次に分布$q_k(z|z^{\tau})$から生成される値$z^*$の受理確率を以下のように定義します.

A_k(z^*, z^{\tau}) = min (1, \frac{p(z^*)q_k(z^{\tau}|z^*)}{p(z^{\tau})q_k(z^*|z^{\tau})})

この式の意味を考えてみましょう．
$p(z)$は求めたい分布における生成確率です．
なので
$ \frac{p(z^\star)}{p(z^{\tau})} $は新たに生成された値が前に生成された値より，どれだけ生成されやすいかを表しています．
$q_k(z^\star|z^{\tau})$は今回生成された値が提案分布において生成される確率を表しています．
$q_k(z^{\tau}|z^*)$は次にまた前の値が生成される確率を表しています．
この値が大きければ逆にも遷移する可能性が高いということです．
この逆にも遷移できるということは、詳しい説明は省きますが（というかできませんが）マルコフ連鎖によって分布が推定できる条件になっているらしいです

ちなみに$q_k$がガウス分布のような対象な分布の場合はこの項を省くことができます．
この形はMetropolis法と呼ばれています．
Metropolis-Hasting法はMetropolis法を対象ではない分布を提案分布にしても用いることができるようにしたものですが，このようにしてみるとMetropolis-Hasting法の特別な形がMetropolis法に見えますね．

ではこちらも簡単に実装をしてみます．

function metropolis_sampling(target_distribution::Distribution, M::Int, burn_in::Int)
    m = M + burn_in
    x0 = randn()
    samples = Float64[]
    for i in 1:m
        x1 = x0 + randn()
        rate = pdf(target_distribution, x1)/pdf(target_distribution, x0)
        if rate < 1 && rate < rand()
            continue
        end
        x0 = x1
        if i < burn_in
            continue
        end
        push!(samples, x1)
    end
    samples
end

提案分布を以前に生成された値が平均値になるガウス分布にしました．
この場合は左右対称になるので受理関数の左側の部分は入りません．
burn_inというのはMCMCで必要になる要素です．
MCMCでは前の要素によってサンプルが生成されるため、初期の頃は初期値に依存して正確なサンプルが得られない状況にあります
そのためある程度捨てる区間が必要で，それをburn_inと呼びます

実行すると以下のようになりました

ガウス分布が生成できていることがわかると思います．

まとめ

本エントリではサンプリングの手法についてまとめました．
とっつきにくいところが多いサンプリングについてを自分なりにわかりやすくなるようまとめたつもりです．
お読みいただいた方の少しでもお役にたてれば幸いです．

宣伝

Gunosyではデータマイニング領域に関する勉強会を隔週でおこなっています．
元々大学のゼミのような活動を会社でも定期的にやることで勉強のペースを作りたいというモチベーションでやっており，まもなく100回を迎えます．
現在はトピックモデル本の輪読と論文紹介をやっておりますので，興味のある方は是非参加してみてください．

参考資料

パターン認識と機械学習下巻 11章
 データ解析のための統計モデリング入門
 トピックモデルによる統計的潜在意味解析
 Gibbs sampling(jags) VS HMC – Hamiltonian Monte Carlo(stan)
Juliaで学ぶ Hamiltonian Monte Carlo (NUTS 入り)
R による棄却サンプリング法（おれが）入門
[Stanで統計モデリングを学ぶ(2): そもそもMCMCって何だったっけ？]
(http://tjo.hatenablog.com/entry/2014/02/08/173324)