ギブスサンプリング：やってることのイメージ
統計の計算とかをやろうとすると、サンプリングという方法をとって計算をさせるという場面がよく起こります。
どういうときに使うのかというと、例えば確率密度関数に従う確率変数zを引数に取るある関数
の期待値が計算したい場合などです。この場合計算するべきは
の計算になるわけですが、実際に計算するとなると、この関数がきれいな形をしていないと解析的に解くことが難しくなります。要するに不定積分をして数値を代入するだけでは解けなくなります。
もしどうしてもこれを計算したかったら、近似的にこの積分に一致するように計算を行っていく必要が出てきます。
一番シンプルに考えた時に思いつくのは「積分の区間を均一に区切ってやって、その時の関数の値で近似を行う」という方法です。
しかしこの方法にも問題があります。それは確率変数がすごく大きな次元を取っていたときに、計算に必要な区切りの数が爆発してしまうという問題です。（次元の呪いと言います）
ですから全部の領域を区切っていく方法には限界があります。
そこで新たに考えられるのは、確率によってzをたくさん選んできて（これを
(iは何番目に選ばれたか）とあらわすことにします）、その
を使って各々で
を計算してやって、それの平均値をとるという方法です。数式にすると
という風になります。こうするともはや計算するときにすべての領域を見る必要がなくなっています。また感覚的にもわかるように、選ぶ数をどんどん増やしていくとだんだんと求めたい値に収束していく、という性質も持ち合わせています。しかしこの方法も問題がないわけではなく、今度は全空間をなめる必要がなくなった代わりに、確率からzを選ばなくてはならなくなります。
もし確率がすべての空間で同じである（一様分布）ならば、適当な乱数をその範囲内の値から0の間をとるように発生させて、それをzとして採用すればよいですがふつう世の中にあるような問題を解こうとすると、もっと複雑な確率分布に従った現象は多く存在していますから、どうにかしてそのような難しい確率分布にそってzを取り出してくることが必要になります。
このzを取り出してくる操作のことをサンプリングと言い、いろいろなサンプリングの方法が提案されています。
今回はそのうちの一つである、ギブスサンプリングというものを作ってみました。
ギブスサンプリングはサンプリング法の中でもマルコフ連鎖モンテカルロ法（MCMC)という部分に分類されます。
マルコフ連鎖モンテカルロは、サンプル点を直前に取り出したサンプルに依存して新しいサンプルを取得します。これがマルコフ連鎖という名前のもとになっています。その中でもギブスサンプリングでは、その新しい点を取り出す操作をある次元にだけにそって行います。
別の言い方をすると、例えば3次元であればx方向だけを変える、と決めれば直前のy,zはそのままにして、xだけを新しく選びます。そしてその選ばれたx'を使ってを新しいサンプルとして採用します。この操作をいろいろな次元にそって行っていくことで、たくさんのサンプル点を次々に生成していきます。なので一回のサンプリングでは、変数は一つしか変わりません。
ではそのｘはどう選ぶのかというと、次のxの条件付き確率によって選びます。
ギブスサンプリングが有効になるのは、この条件付き確率が簡単な場合です。なぜなら、結局は新しい点は、その条件付き確率（提案分布）から選んでこないといけませんから、この条件付き確率が難しい形をしていると意味がないからです。具体的な例で言うと、多次元のガウス分布などがあります。多次元ガウス分布はその条件付き確率も必ずガウス分布になることが知られていて、ガウス分布に従うサンプリングは、一様分布を用いて簡単に作成することができる（Marsaglia,Box-Muller法など）ので、一回の繰り返しは一様分布に従うサンプリングができれば簡単に行えます。ですからギブスサンプリングを用いれば、多次元ガウス分布によるサンプリングを行うことができるのです。*1*2
これを使って、2次元の場合にギブスサンプリングを行うプログラムをpythonで作成しました。
# coding: utf-8 __author__ = "nyk510" """ ギブスサンプリングを用いた, 二次元ガウス分布からのサンプリング """ import numpy as np import matplotlib.pyplot as plt import scipy.stats as stats def gibbs_sampling(nu, cov, sample_size): """ ギブスサンプリングを用いて与えられた共分散, 平均値を持つ 多次元ガウス分布からのサンプリングを行う関数 :param np.ndarray nu: 平均値 :param np.ndarray cov: 共分散 :param int sample_size: サンプリングする数 :return: :rtype: np.ndarray """ samples = [] n_dim = nu.shape[0] # start point of sampling start = [0, 0] samples.append(start) search_dim = 0 for i in range(sample_size): if search_dim == n_dim - 1: """ search dimension selection is cyclic. it can be replaced random choice. """ search_dim = 0 else: search_dim = search_dim + 1 prev_sample = samples[-1][:] A = cov[search_dim][search_dim - 1] / float(cov[search_dim - 1][search_dim - 1]) # A*Σ_yy = Σ_xy _y = prev_sample[search_dim - 1] # previous values of other dimension # p(x|y) ~ N(x|nu[x]+A(_y-nu[y]),Σ_zz) # Σ_zz = Σ_xx - A0*Σ_yx mean = nu[search_dim] + A * (_y - nu[search_dim - 1]) sigma_zz = cov[search_dim][search_dim] - A * cov[search_dim - 1][search_dim] sample_x = np.random.normal(loc=mean, scale=np.power(sigma_zz, .5), size=1) prev_sample[search_dim] = sample_x[0] samples.append(prev_sample) return np.array(samples) if __name__ == '__main__': # 2 dimension normal distribution nu = np.ones(2) covariance = np.array([[0.5, 0.5], [0.5, 3]]) # eig_values: 固有値 # eig_vectors: 固有ベクトル eig_values, eig_vectors = np.linalg.eig(covariance) average_eigs = np.average(eig_values) sample = gibbs_sampling(nu, covariance, 1000) fig, ax1 = plt.subplots(figsize=(6, 6)) ax1.scatter(sample[:, 0], sample[:, 1], marker="o", facecolor="none", alpha=1., s=30., edgecolor="C0", label="Samples" ) # 答え合わせ # scipy.stats を用いて多次元ガウス分布の確率密度関数を計算 multi_norm = stats.multivariate_normal(mean=nu, cov=covariance) X, Y = np.meshgrid(np.linspace(nu[0] - average_eigs * 2, nu[0] + average_eigs * 2, 1000), np.linspace(nu[1] - average_eigs * 2, nu[1] + average_eigs * 2, 1000)) Pos = np.empty(X.shape + (2,)) Pos[:, :, 0] = X Pos[:, :, 1] = Y Z = multi_norm.pdf(Pos) ax1.contour(X, Y, Z, colors="C0", label="True Probability Density Function") ax1.legend() ax1.set_title("Gibbs Sampling") fig.tight_layout() fig.savefig("figures/gibbs_sampling.png", dpi=150)
ちゃんと目標の分布に従ってサンプルが選ばれていることがわかります。
