【Python】重点サンプリング（Importance Sampling)の理論と実装

 

本記事では、マルコフ連鎖を利用しない静的なモンテカルロと呼ばれる手法の一つである『重点サンプリング』を簡単に説明して、そのアルゴリズムをPythonで実装します。

 

重点サンプリング（Importance Sampling)の基本

 

まずは、重点サンプリング（Importance Sampling）の解説を行います。

 

サンプリングを用いた期待値評価

 

まず、対象となる離散確率分布を次のように定義します。

$$p(x) = \frac{1}{Z_{p}} \tilde{p}(x)$$

ここで、$Z_{p}$は規格化定数で、次の関係を満たします。

$$Z_{p} = \int dx \tilde{p}(x)$$

そして、この確率分布に従う確率変数$x$に関する関数$f(x)$の期待値を$\mathbb{E}_{p(x)}[f(x)]$とおくと、期待値は次のように定義されます。

$$\mathbb{E}_{p(x)}[f(x)] = \int dx p(x) f(x)$$

一般的に、この積分を実行できるケースはごくわずかです（正規分布等の特殊な確率分布の場合は積分が実行できる可能性があります）

また、グリットに分けて数値計算を行う手法も次元が大きくなると現実的な時間で実行不可能となります。仮に$n$次元変数に関する積分を10個のグリットに分けて行う場合、$10^{n}$の和を計算する必要があります泣

これらの問題を解決するために、モンテカルロ法のアイデアが有効です！

具体的に、モンテカルロ法のアイデアは以下の通りです。

具体的には以下のように近似を行います。

 

確率論の定理（中心極限定理）から、モンテカルロ法の期待値評価は、サンプル数$T \to \infty$の極限で真の期待値に収束します。

また、モンテカルロ法の期待値評価による誤差は、変数の次元に依存せずに減少していきます。

期待値の問題を『確率分布からのサンプリング問題』に変換することができました。

具体的な例を通して、サンプリングによる期待値近似を理解したい方は、モンテカルロ法のイメージのサイコロの例を一読してみてください。

＊本記事で『サンプリング』と言った場合は、数式で与えられた確率分布から具体的なサンプル（列）を生成することを意味しています。

 

重点サンプリングの目的

 

ここまで、確率分布$p(x)$に従う期待値が確率分布$p(x)$からのサンプリング問題により近似できることを紹介しました。

しかし、一般的に確率分布$p(x)$の規格化定数$Z_{p}$が簡単に計算できる場合は稀です。そのため、$p(x)$からのサンプリングもそこまで容易ではありません。

重点サンプリングは、規格化定数が未知でサンプリングが困難な$p(x)$に従う確率変数$x$の期待値やその規格化定数を別のサンプリングが容易な確率分布$q(x)$からサンプリングにより解決する手法です。

詳細は、以降説明していきます。その前に簡単にこれまでの流れをまとめておきます。

 

重点サンプリングの仕組み

 

重点サンプリングでは、まずサンプリングが容易かつ規格化定数が既知の確率分布$q(x)$を用意します。

$$q(x) = \frac{1}{Z_{q}} \tilde{q}(x)$$

ここで、$Z_{q}$は規格化定数で容易に計算可能であることを仮定します。以降、この分布$q(x)$を提案分布と呼ぶことにします。

また、この確率分布$q(x)$は以下の条件を満たす必要があります。

以降の解説では、確率分布$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 \\&= \mathbb{E}_{q(x)} \left[A(x) \frac{p(x)}{q(x)}  \right] \end{align}$$

この単純な式変形により、$p(x)$に従う確率変数$x$の$A(x)$の期待値を計算する問題が、$q(x)$に従う確率変数$x$の$A(x) p(x)/q(x)$の期待値を計算する問題に帰着されます。

さらに式変形を進めると次のように表せます。

$$\begin{align} \mathbb{E}_{q(x)} \left[A(x) \frac{p(x)}{q(x)}  \right] &= \int q(x)　A(x) \frac{p(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{\tilde{p}(x^{(m)})}{\tilde{q}(x^{(m)})} \\&= \frac{Z_{q}}{Z_{p}} \frac{1}{M} \sum_{m=1}^{M} A(x^{(m)}) w(x^{(m)})\end{align}$$

ここで、$w^{(m)}$は以下のように定義しました。

$$w(x^{(m)}) := \frac{\tilde{p}(x^{(m)})}{\tilde{q}(x^{(m)})}$$

一般的に$w^{(m)}$は重要度重み（importance weight）と呼ばれます。

また, $x^{(m)}$は、$p(x)$ではなく$q(x)$からのサンプル点です。

あとは、規格化定数$Z_{p}$が計算できれば、$A(x)$の期待値が評価できます。

この後すぐに説明しますが、規格化定数$Z_{p}$も重要度重みを用いて評価することができます。

 

規格化定数の計算

 

規格化定数は、重要度重みのサンプル平均により計算することができます。

実際、規格化定数の比は、以下のように変形できます。

$$\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(x^{(m)}) \end{align}$$

いま、$Z_{q}$は、既知なので規格化は以下のように計算することができます。

$$Z_{p} \approx \frac{Z_{q}}{M} \sum_{m=1}^{M} w(x^{(m)})$$

この結果を統計量$A(x)$の期待値評価の式に代入すると、期待値は以下のように評価することができます。

$$\begin{align} \mathbb{E}_{p(x)}[A(x)] &\approx \frac{Z_{q}}{Z_{p}} \frac{1}{M} \sum_{m=1}^{M} A(x^{(m)}) w(x^{(m)}) \\ &= \frac{\sum_{m=1}^{M} A(x^{(m)}) w(x^{(m)})}{\sum_{m=1}^{M} w(x^{(m)})}\end{align}$$

アルゴリズム的には、提案分布$q(x)$からサンプリングを行い、重要度重みと統計量の重要度重み付き和を計算すれば良いです。

ここまでの結果を含め、アルゴリズムを整理します。

 

重点サンプリングのアルゴリズム

 

 

この後は、このアルゴリズムをPythonで実装して挙動を確かめます。

 

 

重点サンプリングの具体例（Python）

 

ここでは、実際にPythonを用いて重点サンプリングを実装していきます。

今回使用する例は以下の一次元混合ガウス分布です。

$$\pi \mathcal{N}(0, 1) + (1- \pi) \mathcal{N}(4, 1)$$

重点サンプリングの設定では、規格化されていない確率分布が既知という状況なので、以下の関数のみ既知であることを仮定します。

 

このような状況下で平均と規格化定数$Z$を推定する問題を考えます。

今回は、正解が計算できて、平均は$4$となり、規格化定数$Z$は、$\sqrt{2 \pi}$となります。

提案分布をガウス分布として、以下の二つのガウス分布を用いて挙動を確かめます。

 

 

使用するライブラリをインポート

 

まずは、以下のように実装に使用するライブラリをインポートしてください。

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)

 

ここで、リストは以下を格納しています。

 

得られた規格化定数の結果をプロットします。

 

次に平均の推定結果を示します。

 

分布$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学習を効率化させるサービス

 

Pythonを学習するのに効率的なサービスを紹介していきます。

まず最初におすすめするのは、Udemyです。

Udemyは、Pythonに特化した授業がたくさんあり、どの授業も良質です。

また、セール中は1500円定義で利用することができ、コスパも最強です。

下記の記事では、実際に私が15個以上の講義を受講して特におすすめだった講義を紹介しています。

 

【最新】UdemyでおすすめのPythonコース｜東大生が厳選！10万を超える講座があるUdemyの中で、Pythonに関係する講座を厳選しました。また、本記事では、Udemyを使用しながらPythonをどのような順番で勉強するべきかを紹介しました。ぜひ参考にしてください。...

 

他のPythonに特化したオンライン・オフラインスクールも下記の記事でまとめています。

 

【最新】Pythonに強いプログラミングスクール7選｜東大生が厳選Pythonの流行と共に、Pythonに強いプログラミングスクールが増えてきました。本記事では、特にPythonを効率的に学ぶことができるプログラミングスクールを経験をもとに厳選して、内容を詳しく解説しています。...

 

自分の学習スタイルに合わせて最適なものを選びましょう。

また、私がPythonを学ぶ際に使用した本を全て暴露しているので参考にしてください。

 

• 【レベル別】Pythonを学ぶための本を厳選しました｜入門〜上級者まで
• 【レベル別】Pythonで機械学習を学ぶための本｜最短実務応用を目指せ！