本記事では、マルコフ連鎖を利用しない静的なモンテカルロと呼ばれる手法の一つである『重点サンプリング』の理論とPythonによる実装を紹介します。
重点サンプリング（Importance Sampling)の理論
まずは、重点サンプリング（Importance Sampling）の解説を行います。
重点サンプリングの設定
まずは、以下のような確率分布について考えます。
$$p(x) = \frac{1}{Z_{p}} \tilde{p}(x)$$
ここで、\(Z\)は規格化定数で、以下の関係を満たします。
$$Z_{p} = \int dx \tilde{p}(x)$$
一般に高次元になると規格化定数\(Z_{p}\)の計算は困難になることが多いです。
そのため、今回は\(\tilde{p}(x)\)は簡単に計算できるが、\(Z_{p}\)の計算は困難となる状況を想定します。
そして、重点サンプリングの目的は、評価可能な\(\tilde{p}(x)\)をもとに、規格化定数\(Z_{p}\)や\(p(x)\)の統計量\(A(x)\)の期待値を評価することです。
重点サンプリングの理論
まず重点サンプリングでは、サンプリングが容易かつ規格化定数が既知の以下のような確率分布\(q(x)\)を用意します。
$$q(x) = \frac{1}{Z_{q}} \tilde{q}(x)$$
以降、このような分布\(q(x)\)を提案分布と呼ぶことにします。
ここで、確率分布は以下の条件を満たす必要があります。
\(q(x)\)の条件
\(p(x) \neq 0\)となる任意の\(x\)に対して\(q(x) \neq 0\)となる。
以降の解説では、確率分布\(p(x)\)の期待値を\(\mathbb{E}_{p(x)}[\cdot]\)と表して、確率分布\(q(x)\)の期待値を\(\mathbb{E}_{q(x)}[\cdot]\)と表すことにします。
規格化定数と期待値の評価方法
重点サンプリングの基本的な考え方は、ある統計量\(A(x)\)の期待値が、以下のような式変形で、サンプリングが容易な確率分布\(q(x)\)のサンプル平均の形に変形できることです。
\begin{align} \mathbb{E}_{p(x)}[A(x)] &= \int A(x) p(x) dx \\ &= \int A(x) \frac{p(x)}{q(x)} q(x) dx \\&= \frac{Z_{q}}{Z_{p}} \int q(x) A(x) \frac{\tilde{p}(x)}{\tilde{q}(x)} dx \\ &=\frac{Z_{q}}{Z_{p}} \frac{1}{M} \sum_{m=1}^{M} A(x^{(m)})\frac{f(x_{q}^{(m)})}{g(x_{q}^{(m)})} \\&= \frac{Z_{q}}{Z_{p}} \frac{1}{M} \sum_{m=1}^{M} A(x^{(m)}) w^{(m)} \end{align}
ここで、\(w^{(m)}\)は以下のように定義しました。
$$w^{(m)} := \frac{f(x_{q}^{(m)})}{g(x_{q}^{(m)})}$$
一般的に\(w^{(m)}\)は重要度重み（importance weight）と呼ばれます。
あとは、規格化定数\(Z_{p}\)が計算できれば、\(A(x)\)の期待値が評価できます。
この後すぐに説明しますが、規格化定数\(Z_{p}\)も重要度重みを用いて評価することができます。
\(q(x)\)の『\(p(x) \neq 0\)となる任意の\(x\)に対して\(q(x) \neq 0\)となる』という条件は、重要度重みが適切に定義されるように要請しています。
規格化定数の計算
規格化定数は、重要度重みのサンプル平均により計算することができます。
実際、規格化定数の比は、以下のように変形できます。
\begin{align} \frac{Z_{p}}{Z_{q}} &= \frac{1}{Z_{q}} \int \tilde{p}(x) dx \\ &= \int \frac{\tilde{p}(x)}{\tilde{q}(x)} q(x) dx \\ &\approx \frac{1}{M} \sum_{m=1}^{M} w^{(m)} \end{align}
いま、\(Z_{q}\)は、既知なので規格化は以下のように計算することができます。
$$Z_{p} \approx \frac{Z_{q}}{M} \sum_{m=1}^{M} w^{(m)}$$
この結果を統計量\(A(x)\)の期待値評価の式に代入すると、期待値は以下のように評価することができます。
\begin{align} \mathbb{E}_{p}[A(x)] &\approx \frac{Z_{q}}{Z_{p}} \frac{1}{M} \sum_{m=1}^{M} A(x^{(m)}) w^{(m)} \\ &= \frac{\sum_{m=1}^{M} A(x^{(m)}) w^{(m)}}{\sum_{m=1}^{M} w^{(m)}}\end{align}
アルゴリズム的には、提案分布\(q(x)\)からサンプリングを行い、重要度重みと統計量の重要度重み付き和を計算すれば良いです。
ここまでの結果を含め、アルゴリズムを整理します。
重点サンプリングのアルゴリズム
重点サンプリングのアルゴリズムを紹介します。
重点サンプリングのアルゴリズム
- \(q(x)\)からサンプル\((x^{(m)})_{m=1, \ldots M}\)を生成
- 各サンプルに対して、重み\(w_{(m)}\)と統計量\(A(x^{(m)})\)を計算
得られた重み\((w^{(m)})_{m=1, \ldots, M}\)と統計量\((A(x^{(m)}))_{m=1, \ldots, M}\)を用いて以下を計算する
- 期待値
$$\mathbb{E}_{p}[A(x)] \approx \frac{\sum_{m=1}^{M} A(x^{(m)}) w^{(m)}}{\sum_{m=1}^{M} w^{m}}$$
- 規格化定数
$$Z_{p} \approx \frac{Z_{q}}{M} \sum_{m=1}^{M} w^{(m)}$$
この後は、このアルゴリズムをPythonで実装して挙動を確かめます。
重点サンプリングの具体例（Python）
ここでは、実際にPythonを用いて重点サンプリングを実装していきます。
今回使用する例は以下の一次元混合ガウス分布です。
$$\pi \mathcal{N}(0, 1) + (1- \pi) \mathcal{N}(4, 1)$$
重点サンプリングの設定では、規格化されていない確率分布が既知という状況なので、以下の関数のみ既知であることを仮定します。
このような状況下で平均と規格化定数\(Z\)を推定する問題を考えます。
今回は、正解が計算できて、平均は\(4\)となり、規格化定数\(Z\)は、\(\sqrt{2 \pi}\)となります。
提案分布をガウス分布として、以下の二つのガウス分布を用いて挙動を確かめます。
- 提案分布1 : \(\mathcal{N}(0, 1)\)
- 提案分布2 : \(\mathcal{N}(2, 4)\)
使用するライブラリをインポート
まずは、以下のように実装に使用するライブラリをインポートしてください。
import numpy as np
from matplotlib import pyplot as plt
# SEEDを固定
np.random.seed(0)
混合ガウス分布と提案分布を定義
まずは、規格化されていない混合ガウス分布とガウス分布を定義しておきます。
# target
def f(x, mu1=0.0, mu2=4.0, pi=0.5):
pi1 =pi
pi2 =1-pi1
a1 = - ((x - mu1)**(2)) / 2
a2 = - ((x - mu2)**(2)) / 2
return pi1*np.exp(a1)+pi2*np.exp(a2)
# proposal
def g(x, mu=0, scale=1):
return np.exp(-((x - mu)**2)/(2*(scale**2)))
定義した関数を用いて、ターゲットとなる混合ガウス分布と提案分布をプロットしてみましょう。
x=np.linspace(-10.0, 10.0, 1000)
fig, ax = plt.subplots(figsize=(16, 9))
ax.plot(x, f(x), linewidth=3, label='target dist')
ax.plot(x, g(x, mu=0, scale=1), linestyle='--', linewidth=3, label=r'proposal dist1 : $\mu=0, \sigma=1$')
ax.plot(x, g(x, mu=2, scale=2), linestyle='--', linewidth=3, label=r'proposal dist2 : $\mu=2, \sigma=2$')
plt.legend(fontsize=17)
plt.tick_params(axis='y', labelsize=20)
plt.tick_params(axis='x', labelsize=20)
plt.show()
<output>
提案分布2（\(\mathcal{N}(2, 4)\)）の方が目的となる確率分布\(p\)を全体的に覆っているのがわかります。
重点サンプリングのアルゴリズム
重点サンプリングのアルゴリズムを以下に示します。
def importance_sampling(n_samples=100, prop_mean=0, prop_scale=1):
'''重点サンプリング(Target : 多変量ガウス分布, 提案分布: ガウス分布)
Parameters
--------------
num : 提案分布から得るサンプル数
prop_mean : 提案分布(ガウス分布)の平均
prop_scale : 提案分布(ガウス分布)の標準偏差
Returns
---------
norm_term : 規格化定数の推定値
mean : 平均の推定値
'''
# 提案分布からのサンプル
samples = np.random.normal(loc=prop_mean, scale=prop_scale, size=n_samples)
# weightの計算
weights = f(samples)/g(samples, mu=prop_mean, scale=prop_scale)
# 規格化定数
norm_term = (np.sum(weights)/n_samples)*np.sqrt(2*np.pi*(prop_scale**2))
# 平均
mean = np.sum(samples*weights)/weights.sum()
return norm_term, mean
先ほど説明した平均\(0\), 標準偏差\(1\)のガウス分布と平均\(2\), 標準偏差\(2\)のガウス分布のケースでアルゴリズムの結果を示していきます。
重点サンプリングの結果
今回は、重点サンプリングによって得られた規格化定数と平均から真の値を引いた値のサンプル数依存性を確かめます。
各サンプルごとの重点サンプリングは以下のコードで実行できます。
loss_norms1, loss_norms2 = [], []
loss_means1, loss_means2 = [], []
for i in range(1, 5000, 10):
norm_term1, mean1 = importance_sampling(n_samples=i, prop_mean=0, prop_scale=1)
norm_term2, mean2 = importance_sampling(n_samples=i, prop_mean=2, prop_scale=2)
loss_norm1 = norm_term1- np.sqrt(2*np.pi)
loss_norm2 = norm_term2- np.sqrt(2*np.pi)
loss_mean1 = mean1 - 2
loss_mean2 = mean2 - 2
loss_norms1.append(loss_norm1)
loss_norms2.append(loss_norm2)
loss_means1.append(loss_mean1)
loss_means2.append(loss_mean2)
ここで、リストは以下を格納しています。
- 『loss_norm1』、『loss_mean1』 : \(q(x)\)として平均\(0\)、標準偏差\(1\)のガウス分布を利用した結果
- 『loss_norm2』、『loss_mean2』 : \(q(x)\)として平均\(2\)、標準偏差\(2\)のガウス分布を利用した結果
得られた規格化定数の結果をプロットします。
規格化定数の推定結果。横軸は重点サンプリングで使用したサンプル数を表す。
次に平均の推定結果を示します。
平均の推定結果。横軸は重点サンプリングのサンプル数を表す。
分布\(q(x)\)が似ている分布でないと正しい推定結果が得られていないことがわかります。
実は、分布\(q(x)\)が似ていない場合は、重要度重みの分散が大きくなり、推定精度が悪くなってしまいます。
実際、各提案分布に対して重みの分散を評価すると以下のようになります。
def eval_weight(n_samples=100, prop_mean=0, prop_scale=1):
samples = np.random.normal(loc=prop_mean, scale=prop_scale, size=n_samples)
# weightの計算
weights = f(samples)/g(samples, mu=prop_mean, scale=prop_scale)
return weights
weights1 = eval_weight(n_samples=5000, prop_mean=1, prop_scale=1)
weights2 = eval_weight(n_samples=5000, prop_mean=2, prop_scale=2)
print(f'【weight variance】 proposal dist1 : {weights1.var():.3f}, proposal dist2 : {weights2.var():.3f} ')
<output>
【weight variance】 proposal dist1 : 1440.412, proposal dist2 : 0.089
このように形状が類似していない提案分布1の場合、重みの分散が大きくなっています。
まとめ
本記事では、重点サンプリングの理論とPythonによる実装を紹介しました。
皆さんも本記事で紹介したコードのパラメータを変えて遊んでみてください。
Pythonを学習するのに効率的なサービスを紹介していきます。
まず最初におすすめするのは、Udemyです。
Udemyは、Pythonに特化した授業がたくさんあり、どの授業も良質です。
また、セール中は1500円定義で利用することができ、コスパも最強です。
下記の記事では、実際に私が15個以上の講義を受講して特におすすめだった講義を紹介しています。
他のPythonに特化したオンライン・オフラインスクールも下記の記事でまとめています。
自分の学習スタイルに合わせて最適なものを選びましょう。
また、私がPythonを学ぶ際に使用した本を全て暴露しているので参考にしてください。
