サンプリング法メモ

確率 プログラミング

はじめに

ある分布に従った乱数を生成したいことがよくあったりする(一様分布に従う一様乱数や正規分布に従う正規乱数など)。
ベイズ統計学なんかだと、自然共役事前分布が使えないような複雑な分布の場合にMCMCで分布のサンプリングをして計算したりするので、結構使う場面は多い(はず)。
ちょっと調べてみたので、メモ。(基本、PRML第11章のはなし)

一様分布からの(疑似)乱数生成

一様乱数の生成の仕方について。少し。

線形合同法
  • X_{n+1} = (A*X_{n}+B) mod Mという形式で乱数Xを次々に生成していく方法
  • C言語のライブラリで実装されているrand関数
  • 多次元にすると均等に分布しなかったり、下位ビットがランダムじゃなかったり、といろいろ問題がある
xorshift

http://ja.wikipedia.org/wiki/Xorshift

  • xor演算とビットシフトだけで生成していく方法
  • 線形合同法よりも長周期でCのrand()を使うぐらいならこっちを使った方がよさそう(コードも短い)
メルセンヌ・ツイスター

http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/mt.html

  • これまでの疑似乱数生成アルゴリズムの欠点を克服した現状最強の疑似乱数生成の方法
    • 使えるならば使った方が良いと思うけど、欠点がないわけではないので、使うときは注意

サンプリング法

一様分布以外の(連続型の)確率分布からサンプリング(標本抽出)を行う方法について。
以降、[0,1]の一様乱数の生成アルゴリズムは与えられていると仮定。

逆関数法
  • 目的の確率分布の確率分布関数(累積分布関数)F(x)がわかっていて、その逆関数F^{-1}(z)が(容易に)求められる場合の方法
  • 一様分布する確率変数zがx=F^{-1}(z)によって変換されると、xは目的の確率分布に従う
Marsaglia法(Box-Mullar法,Ziggurat法)
  • 極座標を使った正規乱数の生成方法(？)
    • 指数分布に従う乱数でも出てくるのでMarsagliaさんがいろいろがんばったんだと思う
階段関数の近似によるサンプリング

http://www.sci.kagoshima-u.ac.jp/~ebsa/wakimoto01/pdf/ch04-04.pdf

  • 分布の定義域を小区分に分割し、一様乱数を生成する。小区分の高さをf_{k}と区間の幅の積の累積和を考え、その乱数値が含まれる区間のxの中央値を返す
    • 逆関数の離散バージョンみたいな感じ？
棄却サンプリング
  • 単純で標準的ではない比較的複雑な分布p(z)のサンプリングをしたいとき、zのすべての値についてk*q(z)>=p~(z)となる簡単にサンプリングできる分布q(z)を用意し、乱数z_0をq(z)から生成し、乱数u_0を[0,k*q(z_0)]の一様分布から生成する。もし、u_0>p~(z_0)なら破棄し、そうでないならu_0を乱数として採用する
    • kは定数
    • p~(z)は正規化された分布。p(z) = 1/Z p~(z), (Zは正規化定数)
  • 要は、分布を覆うような分布(提案分布)を使って、提案分布から生成した乱数の採用・不採用を決める
  • 提案分布の選択の仕方によっては、棄却するものが多くなったりして無駄が多い
適応的棄却サンプリング
  • 棄却サンプリングにおいて元の分布を覆うような提案分布(包絡分布)を決めるのは結構難しい
  • 分布p(z)の観測値に基づいてその場で包絡関数を構築していく方法
    • 導関数の評価を行う方法、Metropolis-Hastingsステップをおく方法(適応的棄却Metropolisサンプリング)、などがある
  • 上記同様、棄却するものが多くなったり、高次元には適さないなどの問題がある
重点サンプリング
  • この方法は、直接サンプリングが難しい分布場合にある関数の期待値を求める方法
    • なので、乱数生成という枠組みではない
  • 棄却サンプリングのように提案分布を利用して、重要度重みを考慮したすべてのサンプリング値から求める
Sampling-Importance-Resampling(SIR)
  • 提案分布q(z)や定数kの選び方は難しいので、重点サンプリングの重みを活用する
  • 2段階により構成され、1段階目では提案分布q(z)からL個のサンプルを抽出し、2段階目で各サンプルの重みを求め、最後に、重みを確率とするL個の離散分布からサンプリングする方法
メトロポリス法
  • Markov chain(一つ前の状態のみに依存する確率過程) Monte Carlo(乱数を使って数学的問題を解く)法の一つ
  • 対称な提案分布q(x)を定め、現在x_tにいるときサンプル候補x'を提案分布q(x)から生成する、もし一様乱数uがmin(1,p~(x')/p~(x_t) )>uならばx_{t+1}=x'とする、そうでないならばx'を破棄し、x_{t+1}=x_tとし提案分布q(x|x_{t+1})からサンプル候補を生成し、繰り返す
メトロポリス・ヘイスティングス法
  • メトロポリス法の改良
    • 提案分布が対称でない場合にも適用できる
  • 受理確率をmin(1, {p~(x')q_k(x_t|x')} / {p~(x_t)q_k(x'|x_t) )とする
    • kは考えられる可能な遷移集合の要素のラベル
    • 提案分布が対称の時( q(a|b)=q(b|a) )はメトロポリス法に一致する
ギブスサンプリング
  • メトロポリス・ヘイスティングス法の特別な場合(とみなせる)
  • 条件付き分布からサンプリングを繰り返す方法
  • 適当な初期値からはじめ、目的の確率分布に含まれるパラメータの一つ以外を固定した状態で残りの一つをサンプリングすることを繰り返していく
    • 例えばp(x,y)という分布ならば、初期値(x0,y0)->p(x,y0)からx1をサンプリング->p(x1,y)からy1をサンプリング->p(x,y1)から...
    • 初期値から安定したサンプリングになるまでの期間(burn-in,焼き入れ)があるため、一般的に、このburn-inの部分は破棄する
  • 条件付き確率分布が計算しやすい必要がある
スライスサンプリング
  • メトロポリス法はステップサイズに対して敏感
  • 分布の特徴に応じて自動的に調整する方法
ハイブリッドモンテカルロアルゴリズム

参考文献