スライスサンプリングでお手軽サンプリング
ベイズ法などで用いるサンプリング法のなかで，スライスサンプリングというのがあります．原論文は以下になります．
SLICE SAMPLING (Neal, Annals of Statistics 2003)
サンプリング法で代表的なのはMetropolis-HastingやGibbs Samplingですが，これらにはそれぞれ問題点があります．例えばMetropolis法の問題点として，実際の問題に適応する場合，どのような提案分布を選んだら良いのか分からない，ということが挙げられます．よく使われるのは，現在のサンプル点を中心としたGauss分布ですが，分散をどう選べば良いか，ということなど，自明ではありません．提案分布が適切でない場合，サンプルの棄却率が高くなり，アルゴリズムの効率が落ちます．一方，この分散を小さくしすぎると，一度のサンプルで動く距離が小さくなり，これも効率の低下につながります．またGibbs Samplingの場合，サンプルしたい一変数を，他の変数で条件付けた分布からのサンプルを考えるわけですが，この分布が簡単に求まらない場合もあります．
このような問題に対処するため，スライスサンプリングは，提案分布を自分で考えなくても，現在のサンプル点に応じて自動でサンプル幅（上の分散に対応）を選んでくれるサンプリング法です．つまり，難しいことを考えなくてもどんな問題にも適応可能なわけです．そういう意味で，お手軽なサンプリング法と言えるでしょう．
さて，このスライスサンプリングですが，今回簡単に試せるライブラリを見つけたので紹介します．C++のコードになりますが，Mark JohnsonのAdaptor Grammarのパッケージの中に入っています．
Software available from Mark Johnson, Dept of Computing, Macquarie University, Sydney, Australia
実験
ちょっと簡単に試してみました．
メインで呼び出すメソッドは，
template <typename F, typename LogF, typename Uniform01> F slice_sampler1d(const LogF& logF0, //!< log of function to sample F x, //!< starting point Uniform01& u01, //!< uniform [0,1) random number generator F min_x = -std::numeric_limits<F>::infinity(), //!< minimum value of support F max_x = std::numeric_limits<F>::infinity(), //!< maximum value of support F w = 0.0, //!< guess at initial width unsigned nsamples=1, //!< number of samples to draw unsigned max_nfeval=200) //!< max number of function evaluations
です．logF0というのは，求めたい分布の，正規化されていないxの関数値の，対数を返す関数オブジェクトです．xには，1つ前のステップで求めた（つまり現在の）サンプル点を，u01には，[0,1)を返す乱数生成オブジェクト，を入れます．求めたい分布の関数で，対数を指定しないといけないのは，まあそういう決まりということで・・・．
例えば，適当な比率になっている2つの混合Gauss分布からのサンプルを求めるコードは，以下のようになります．
#include <iostream> #include <iomanip> #include "slice-sampler.h" using namespace std; struct LogGaussianMixture { LogGaussianMixture(double mean1, double mean2, double variance, double pi) : mean1(mean1), mean2(mean2), variance(variance), pi(pi) {} double operator()(double x) const { double y1 = exp(-(x - mean1) * (x - mean1) / (2 * variance)) + 1e-100; double y2 = exp(-(x - mean2) * (x - mean2) / (2 * variance)) + 1e-100; return log(pi * y1 + (1 - pi) * y2); } double mean1; double mean2; double variance; double pi; }; struct RandGen { double operator()() const { return (double)rand() / RAND_MAX; } }; int main(int argc, char *argv[]) { LogGaussianMixture log_gm(0, 10, 3, 0.2); RandGen rand_gen; double x = 0; int sample_size = 1000000; for (int i = 0; i < sample_size; ++i) { x = slice_sampler1d(log_gm, x, rand_gen); cout << setprecision(5) << x << " "; } return 0; }
ちょっと注意なのが，関数値にlogをかける前に，微小な値を足しています．これは，関数値が0になった場合に，log(0) = -∞となり，ライブラリ中のassertでこけるためですｗ．
出力したファイルをRで描画してみたのが，以下になります．
具体的なアルゴリズム
これまでサンプリングの中身には全く触れませんでしたが，あまりにブラックボックスだと信用出来ないと思うので，簡単に解説します．
以下，一次元の場合で考えますが，多次元でも同じです．
今，ある分布p(x)がサンプルをとりたい分布だとして，p(x)が正規化されていない，関数f(x)∝p(x)が手元にあるとします．これはベイズで事後分布を求めるときに自然に出てくる話です．
このとき，スライスサンプリングは，以下の手順でp(x)からサンプリングを行います．
- x を適当に初期化
- 以下を繰り返す
- 現在のxの値x0での，関数の値f(x0)を計算する．
- [0, f(x0))の一様分布から，値yをサンプリングする．
- そのyの値に対して，y <= f(x) となる，xの範囲Sを求める．（これは，現在のyの値で関数f(x)を水平に"スライス"し，f(x)と交わる点を求めることに相当）
- 範囲Sの一様分布からxをサンプルし，これを新しいxとする．
これは結局何をしているかというと，y = f(x)をグラフに書いて，f(x)より下にある部分の点(x,y)を，移動させるわけです．この(x,y)のうち，xだけを取り出すと，確率f(x)に沿った値でサンプルが得られる，という話になります．
ここで問題になるのが，y <= f(x)となるxの範囲Sを求める部分で，これはf(x)の関数形が複雑な場合自明ではありません．これに対処するため，論文中には2つのアプローチが紹介されています．両方の共通点は，最初に現在のサンプル点であるx0の回りの小さい区間を考え，それを何らかの枠組みで，Sを近似するほど大きく（大きすぎた場合は小さく）していくということです．Mark Johnsonの実装では，このうち，区間を（確率的に左右どちらかの方向に）倍に倍に増やしていく，doubling procedureというのに基づいています．
上の議論は，もともとの論文中での主張で，主に連続分布を対象にしています．一方，離散分布の場合でもスライスサンプリングは別の意味で有効性が示されていますが，こちらは[機械学習] スライスサンプリング - tsubosakaの日記などを参照してください．
もちろん実際の問題に対しては，もっと良いサンプリング法が個別に存在することもあると思いますが，とりあえず手軽に試せるということで，サンプルしたいパラメータの分布が（正規化されてなくても）書き下せるときは，応急処置的にこれを使って対処してみる，というのは良いのではないかと思います．このパッケージは一次元のサンプリングしか対応していませんが，多次元の場合，自明な方法はGibbs Samplingに組み込む（GibbsもNLPなどでよく出てくる自然な分布だと簡単に定式化出来ますが，一変数のサンプルが自明でない場合もある）ことで対応できます．また，ハイパーパラメータのサンプルなどであれば，多くの場合一次元ですが，尤度関数が非常に複雑ということが多いため，このサンプラーが有効だと思います．
