Science

【Python】棄却サンプリングの実装と理論

【Python】棄却サンプリングの実装と理論

 

本記事では、マルコフ連鎖を利用しない例として『棄却サンプリング』と呼ばれるサンプリング手法を紹介します。

理解度を高めるために、Pythonによる実装を紹介します。

 

 

棄却サンプリング(Rejection Sampling)の解説

まずは、棄却サンプリング(Rejection Sampling)の解説を行います。

 

棄却サンプリングの設定

 

棄却サンプリングの設定は、サンプリングしたい確率分布の規格化定数は未知ですが、規格化されていない\(\tilde{p}(x)\)は既知であるような状況を考えます。

すなわち以下のような確率分布で、\(\tilde{p}(x)\)は既知ですが、、規格化定数\(Z\)は知ることができないような状況を考えます。

$$p(x) = \frac{1}{Z} \tilde{p}(x),~~~Z = \int p(x) dx$$

規格化定数が未知という状況は、複雑な分布になればなるほど生じる問題です。

 

棄却サンプリングのアルゴリズム

 

まず、棄却サンプリングでは、任意の\(x\)に対して以下の条件を満たし、サンプリングが容易な提案分布\(q(x)\)を決定します。

$$kq(x) \ge \tilde{p}(x)$$

すなわち、\(kq(x)\)が\(p(x)\)の全てを覆うようなケースです。

このような提案分布を使用し、棄却サンプリングのアルゴリズムは以下のようになります。

  1. \(q(x)\)から、\(x_{i}\)を生成
  2. 区間\([0, kq(x_{i})]\)の一様分布から乱数\(u_{i}\)を生成
  3. \(u_{i} \le \tilde{p}(x_{i})\)ならばサンプルを受理、\(u_{i} > \tilde{p}(x_{i})\)ならサンプルを棄却

③で受理したサンプルは\(p(x)\)に従います。

 

棄却サンプリングの実装(Python)

ここからは、棄却サンプリングを実装していきましょう。

まずは、使用するライブラリをインポートしてください。

import numpy as np
import matplotlib.pyplot as plt
from scipy import stats
plt.style.use('seaborn-dark')
np.random.seed(0)

 

 

サンプリングを目指す確率分布について

 

今回は以下の混合ガウス分布からのサンプリングを目指します!

$$p(x \mid \mu_{1}, \mu_{2}, \sigma_{1}^{2}, \sigma_{2}^{2}, \pi_{1}, \pi_{2}) = \pi_{1} \mathcal{N}(x \mid \mu_{1}, \sigma_{1}^{2}) + \pi_{2} \mathcal{N}(x \mid \mu_{1}, \sigma_{1}^{2})$$

ここで、\(\pi_{1}, \pi_{2}\)は、\(\pi_{1} + \pi_{2}\)を満たします。

棄却サンプリングでは、規格化定数が未知の状態でのサンプリングを行う手法なので、規格化されていない混合ガウス分布を以下のように定義します。

\begin{align}&\tilde{p}(x \mid \mu_{1}, \mu_{2}, \sigma_{1}^{2}, \sigma_{2}^{2}, \pi_{1}, \pi_{2}) \\ = &\pi_{1} \exp \left(-\frac{(x – \mu_{1})^{2}}{2 \sigma_{1}^{2}} \right) + \pi_{2} \exp \left(-\frac{(x – \mu_{1})^{2}}{2 \sigma_{1}^{2}} \right) \end{align}

規格化された分布と規格化されていない分布を以下にプロットします。

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

# unnormalized Mixed Gaussian
def mgauss(x, mu1=2.0, mu2=-2.0, scale1=1.0, scale2=1.0, pi=0.5):
    pi1 =pi
    pi2 =1-pi1
    a1 = - ((x - mu1)**(2)) / (2*scale1**(2))
    a2 = - ((x - mu2)**(2)) / (2*scale2**(2)) 
    return pi1*np.exp(a1)+pi2*np.exp(a2)


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')
plt.legend(fontsize=17)
plt.tick_params(axis='y', labelsize=20)
plt.tick_params(axis='x', labelsize=20)
plt.show()

 

<output>

規格化された分布と規格化されていない分布青色が規格化されていない分布で、オレンジ色が規格化された分布

 

 

提案分布の設定

 

今回棄却サンプリングに使用する提案分布は、以下のガウス分布です。

$$q(x) = \frac{1}{\sqrt{2 \pi}} \exp \left(-\frac{(x-\mu)^{2}}{2 \sigma^{2}}  \right)$$

まずは、\(q(x)\)が\(\tilde{p}(x)\)にぎりぎり覆い被さるような\(k\)を計算して、棄却サンプリングを実装します。

 

# 提案分布
def q(x, loc=0, scale=1):
    return stats.norm.pdf(x, loc=loc, scale=scale)

x = np.linspace(-10.0, 10.0, 500)

# 最適なkを計算
k_opt = np.max(mgauss(x)/q(x, scale=2))

fig, ax = plt.subplots(figsize=(10, 6))
ax.plot(x, mgauss(x), linewidth=3, label='unnormalized gaussian mixture')
ax.plot(x, k_opt*q(x, scale=2), linewidth=3, label='normalized gaussian')
plt.legend(fontsize=13)
plt.tick_params(axis='y', labelsize=20)
plt.tick_params(axis='x', labelsize=20)
plt.show()

 

<output>

 

提案分布と目的の分布黄色が提案分布で青色が規格化されていない混合ガウス分布

 

 

棄却サンプリングのアルゴリズム

 

いよいよ、棄却サンプリングアルゴリズムを実装します。

実装は以下のようになります。

def rejection_sampling(times=100, k=k_opt):
    '''Rejection Sampling

    Parameters
    times : 使用するサンプル数
    k : Rejection Samplingの条件を満たすようなk

    Return
    samples : acceptされたサンプル
    '''
    samples = []
    accepted_num = 0

    for i in range(times):
        x = np.random.normal(0, 2.0)
        u = np.random.uniform(0, k*q(x, scale=2))

        if u <= mgauss(x):
            samples.append(x)
            accepted_num += 1
    print(f'ACCEPTANCE RATE : {accepted_num/times}')
    return samples, accepted_num/times

 

 

棄却サンプリングの結果

 

棄却サンプリングで得られた結果と規格化されている混合ガウス分布を比較してみました。

samples = rejection_sampling(times=10000)
fig, ax = plt.subplots(figsize=(16, 9))
ax.hist(samples, bins=50, density=True, label='accepted sample')
ax.plot(x, mgauss_target(x), linewidth=3, label='unnormalized gaussian mixture')
plt.legend(fontsize=17)
plt.tick_params(axis='y', labelsize=20)
plt.tick_params(axis='x', labelsize=20)
plt.show()

 

<output>

rejection-samplingの結果

 

分布の形状の類似度と受理率の関係

 

以下のように、分布の類似度を\(k\)の値を大きくして小さくした場合、どのように受理されるサンプルが変わるかを確認してみます。

x = np.linspace(-10.0, 10.0, 500)
fig, ax = plt.subplots(figsize=(10, 6))
ax.plot(x, mgauss(x), linewidth=3, label='target')
for k in ks:
    ax.plot(x, k*q(x, scale=2), linewidth=3, label=f'k = {k:.1f}')
plt.legend(fontsize=13)
plt.tick_params(axis='y', labelsize=20)
plt.tick_params(axis='x', labelsize=20)
plt.show()

 

<output>

異なる高さの提案分布

このときの受理率は以下の様になります。

類似度と受理率の関係

この様に提案分布の類似度が異なると受理されるサンプルが大幅に下がります。

低次元の場合、このような分布を見つけることは簡単ですが、高次元になるとその様な分布を探すのは困難です。

 

まとめ

 

本記事では棄却サンプリングを利用して、規格化定数が未知の確率分布からサンプリングする方法を紹介しました。

後半では、具体例を示すため、Pythonによる実装をまとめました。

 

 

ABOUT ME
努力のガリレオ
【運営者】 : 東大で理論物理を研究中(経歴)東京大学, TOEIC950点, NASA留学, カナダ滞在経験有り, 最優秀塾講師賞, オンライン英会話講師試験合格, ブログと独自コンテンツで収益6桁達成 【編集者】: イングリッシュアドバイザーとして勤務中(経歴)中学校教諭一種免許取得[英語],カナダ留学経験あり, TOEIC650点
【必見】オンライン英会話を利用しないのは時代遅れ?

 

 

英語学習でオンライン英会話を利用しないのは、時代遅れです…

 

『スピーキングができないから、オンライン英会話はまだまだ後回し…』と思っている方は大間違いです。

 

最近のオンライン英会話は、教材が豊富で、もはや参考書・問題集は不要になりつつあります…

また、多くのオンライン英会話が最新メソッドを導入して、スピーキングのみならず、リーディング・リスニング・ライティングの能力も上がる教材を導入しています。

さらに、効率的な記憶ツールや学習カウンセリングのサービスが提供されることもあります!!

 

特に初心者の方は、下記の記事で紹介しているオンライン英会話をチェックしてみてください。

 

【東大生が厳選】初心者にオススメなオンライン英会話10選実は、最適なオンライン英会話はユーザーによって異なります。本記事では、複数のオンライン英会話を経験した私が、他社と比較しつつ特徴を明確にして本当にオススメできるオンライン英会話を10個厳選しました。...

 

中には、2週間の無料体験期間を提供しているオンライン英会話もあります。

英語学習が滞っている方や英語学習を習慣化できていない方は無料体験期間で良いので挑戦してみてください。

何事も行動することが大切です。

 

また、多くのオンライン英会話がTOEICやTOEFL等の外部英語試験の対策レッスンを提供してきています。

中には、追加料金なしで、本番さながらの対策や面接対策を行うことができるオンライン英会話もあります。

 

特にTOEIC対策に強いオンライン英会話は下記を参考にしてください。

 

【徹底比較】TOEIC対策に強いオンライン英会話5選|東大生が厳選!オンライン英会話はTOEICのスコアを効率的に上げる最強のツールです。本記事では、TOEICにオンライン英会話が効果的な理由とTOEIC対策に強いオンライン英会話を厳選しました。詳しい内容は記事を参考にしてください。...

 

TOEFL対策に強いオンライン英会話に関しては下記を参考にしてください。

 

【徹底比較】TOEFL対策に強いオンライン英会話5選|東大生が厳選!TOEFL対策にオンライン英会話は有効です。事実、私もTOEFL対策のためにオンライン英会話を使用してTOEFL ibtの4技能、TOEFL itpのスコアを爆上げすることができました。この記事では、私が思うTOEFL対策に最適なオンライン英会話を紹介しました。...