Science

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

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

 

本記事では、マルコフ連鎖を利用しない静的なモンテカルロと呼ばれる手法の一つである『重点サンプリング』の理論とPythonによる実装を紹介します。

 

重点サンプリング(Importance Sampling)の理論

 

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

 

重点サンプリングの設定

 

まずは、以下のような確率分布について考えます。

 

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

 

ここで、\(Z\)は規格化定数で、以下の関係を満たします。

 

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

 

一般に高次元になると規格化定数\(Z_{p}\)の計算は困難になることが多いです。

そのため、重点サンプリングの目的は、評価可能な\(f(x)\)をもとに規格化定数や\(p(x)\)の統計量\(A(x)\)の期待値を評価することです。

 

 

重点サンプリングの理論

 

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

 

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

 

ここで、確率分布は以下の条件を満たす必要があります。

\(q(x)\)の条件

\(p(x) \neq 0\)となる任意の\(x\)に対して\(q(x) \neq 0\)となる。

 

まずは、\(A(x)\)の期待値を求める方法を説明します。

以降の解説では、確率分布\(p(x)\)の期待値を\(\mathbb{E}_{p}[\cdot]\)と表して、確率分布\(q(x)\)の期待値を\(\mathbb{E}_{q}[\cdot]\)と表すことにします。

 

規格化定数と期待値の評価方法

 

 

重点サンプリングの基本的な考え方は、ある統計量\(A(x)\)の期待値が、以下のような式変形で、サンプリングが容易な確率分布\(q(x)\)からのサンプル平均で変形できることです。

 

\begin{align} \mathbb{E}_{p}[A(x)] &= \int  A(x) p(x) dx \\ &= \int A(x) \frac{p(x)}{q(x)} q(x) dx \\&= \mathbb{E}_{q}\left[A(x) \frac{p(x)}{q(x)} \right] \\ &\approx \frac{1}{M} \sum_{m=1}^{M} A(x^{(m)}) \frac{p(x^{(m)})}{q(x^{(m)})} \\ &=\frac{Z_{q}}{Z_{p}} \frac{1}{M} \sum_{m=1}^{M} A(x^{(m)})\frac{f(x^{(m)})}{g(x^{(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^{(m)})}{g(x^{(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 f(x) dx \\ &= \int \frac{f(x)}{g(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)\)からサンプリングを行い、重みや統計量の標本平均をとれば良だけです。

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

 

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

 

重点サンプリングのアルゴリズムを紹介します。

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

以下を\(M\)回繰り返す。

  1. \(q(x)\)からサンプル\(x^{(m)}\)を生成する
  2. 重み\(\tilde{w}_{m}\)と統計量\(A(x^{(m)})\)を計算する

得られた重み\(\tilde{w}_{1}, \ldots, \tilde{w}_{M}\)と統計量\(A(x^{(1)}), \ldots, A(x^{(M)})\)を用いて以下を計算する

  • 期待値

$$\mathbb{E}_{p}[A(x)] \approx \frac{\sum_{m=1}^{M} A(x^{(m)} \tilde{w}_{m}}{\sum_{m=1}^{M} \tilde{w}_{m}}$$

  • 規格化定数

$$Z_{p} \approx \frac{Z_{q}}{M} \sum_{m=1}^{M} \tilde{w}_{m}$$

 

以降、このアルゴリズムを簡単に具体例に対して適用します。

 

重点サンプリングの具体例(Python)

 

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

今回使用する例は分散が\(1\)の一次元混合ガウス分布です。

$$\pi \mathcal{N}(\mu_{1}, 1) + (1- \pi) \mathcal{N}(\mu_{2}, 1)$$

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

$$\pi \exp \left(- \frac{(x-\mu_{1})}{2}\right) + (1- \pi)\exp \left(- \frac{(x-\mu_{2})}{2}\right)$$

 

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

今回の場合、平均は\(\mu_{1} + \mu_{2}\)となり、規格化定数\(Z\)は、\(1/\sqrt{2 \pi}\)となります。

提案分布をガウス分布として、この値を重点サンプリングにより近似的に評価していきます。

 

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

 

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

import numpy as np
from matplotlib import pyplot as plt
import scipy.stats as stats

# SEEDを固定
np.random.seed(0)

 

混合ガウス分布を定義

 

まずは、規格化されている混合ガウス分布、規格化されていない混合ガウス分布とガウス分布を定義しておきます。

# Mixed Gaussian 
def mgauss_target(x, mu1=2.0, mu2=-2.0, pi=0.5):
    pi1 =pi
    pi2 =1-pi1
    return pi1*stats.norm.pdf(x, loc=mu1, scale=1.0)+pi2*stats.norm.pdf(x, loc=mu2, scale=1.0)

# unnormalized Mixed Gaussian
def mgauss(x, mu1=2.0, mu2=-2.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)

# Gaussian
def q(x, loc=0, scale=1):
    return stats.norm.pdf(x, loc=loc, scale=scale)

 

下記を入力し、これらの関数をプロットしてみましょう。

x=np.linspace(-10.0, 10.0, 1000)
fig, ax = plt.subplots(figsize=(16, 9))
ax.plot(x, mgauss(x), linewidth=3, label='unnormalized gaussian mixture')
ax.plot(x, mgauss_target(x), linewidth=3, label='normalized gaussian mixture')
ax.plot(x, q(x, loc=0, scale=1), linestyle='--', linewidth=3, label=r'normalized gaussian : $\mu=0, \sigma=1$')
ax.plot(x, q(x, loc=2, scale=2), linestyle='--', linewidth=3, label=r'normalized gaussian : $\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>

 

数値実験に使用する関数

今回は、上記にプロットしたように提案分布として、平均\(0\), 標準偏差\(1\)と平均\(2\), 標準偏差\(2\)のガウス分布を使用し、比較検討していきます。

 

 

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

 

重点サンプリングのアルゴリズムを以下に示します。

def importance_sampling(num=100, prop_mean=0, prop_scale=1):
    '''重点サンプリング(Target : 多変量ガウス, 提案分布: ガウス分布)
    
    Parameters
    num : 提案分布から得るサンプル数
    prop_mean : 提案分布(ガウス分布)の平均 
    prop_scale : 提案分布(ガウス分布)の標準偏差

    Output
    norm_term : 規格化定数の推定値
    mean : 平均の推定値
    '''
    # 提案分布からのサンプル
    samples = np.random.normal(loc=prop_mean, scale=prop_scale, size=num)
    # weightの計算
    weights = mgauss(samples)/q(samples,loc=prop_mean,scale=prop_scale)
    samples_eval = np.dot(samples, weights)
    
    # 規格化定数
    norm_term = np.sum(weights)/num
    # 平均
    mean = samples_eval.sum()/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(num=i, prop_mean=0, prop_scale=1)
    norm_term2, mean2 = importance_sampling(num=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)\)が似ていない場合は、重要度重みの分散が大きくなり、推定精度が悪くなってしまいます。

 

まとめ

 

本記事では、重点サンプリングの理論とPythonによる実装を紹介しました。

皆さんも本記事で紹介したコードのパラメータを変えて遊んでみてください。

ABOUT ME
努力のガリレオ
【運営者】 : 東大で理論物理を研究中(経歴)東京大学, TOEIC950点, NASA留学, カナダ滞在経験有り, 最優秀塾講師賞, オンライン英会話講師試験合格, ブログと独自コンテンツで収益6桁達成 【編集者】: イングリッシュアドバイザーとして勤務中(経歴)中学校教諭一種免許取得[英語],カナダ留学経験あり, TOEIC650点
【大学生・大学院生限定】Amazon Prime Studentが圧倒的におすすめ

 

『Amazon Prime Student』は、大学生・大学院生限定のAmazon会員制度です。

Amazonを使用している方なら、必ず登録すべきサービスといっても過言ではありません…

主な理由は以下の通りです。

  1. Amazon Prime』のサービスを年会費半額で利用可能
  2. 本が最大10%割引
  3. 文房具が最大20%割引
  4. 日用品が最大15%割引
  5. お急ぎ便・お届け日時指定便が使い放題
  6. 6ヶ月間無料で使用可能

 

特に専門書や問題集をたくさん買う予定の方にとって、購入価格のポイント10%還元はめちゃめちゃでかいです!

少なくとも私は、Amazon Prime Studentを大学3年生のときに知って、めちゃめちゃ後悔しました。

専門書をすでに100冊以上買っていたので、その10%が還元できたことを考えると泣きそうでした…ww

より詳しい内容と登録方法については下記を参考にしてください。

 

Prime Studentの特典・登録方法【学生限定最強特典!】Prime Studentの特典・登録方法【学生限定最強特典!】 学生にとって超おすすめな『Prime Student』のメリット...

 

 

登録も退会もめちゃめちゃ簡単なので、6ヶ月の無料体験期間だけは経験してみても損はないと思います。