からっぽのしょこ

読んだら書く!書いたら読む!同じ事は二度調べ(たく)ない

5.5:重点サンプリング【ゼロつく4のノート】

はじめに

 『ゼロから作るDeep Learning 4 ――強化学習編』の独学時のまとめノートです。初学者の補助となるようにゼロつくシリーズの4巻の内容に解説を加えていきます。本と一緒に読んでください。

 この記事は、5.5節の内容です。重点サンプリングの計算を確認します。

【前節の内容】

www.anarchive-beta.com

【他の記事一覧】

www.anarchive-beta.com

【この記事の内容】

5.5 方策オフ型と重点サンプリング

 前節までは、方策オン型で推定を行いました。次節では、方策オフ型で推定します。オフ型のモンテカルロ法では、重点サンプリングにより期待値を求めます。この節では、通常のサンプリング(オン型)と重点サンプリング(オフ型)による期待値計算を比較します。

5.5.2 重点サンプリング

 重点サンプリングでは、対象の確率分布とは別の確率分布のサンプルを用いて対象の確率分布に関する期待値を求めます。

・数式の確認

 まずは、重点サンプリングの計算式を導出します。

・期待値の近似式

 ある離散型の確率分布$\pi(x)$の期待値$\mathbb{E}_{\pi}[x]$を求めたいとします。確率変数$x$の期待値は、次の式で定義されます。

$$ \mathbb{E}_{\pi}[x] = \sum_x x \pi(x) $$

 $\pi(x)$による期待値であることを$\mathbb{E}_{\pi}[\cdot]$で表します。

 モンテカルロ法では、$\pi(x)$からサンプリングした$x^{(i)}$の標本平均で$\mathbb{E}_{\pi}[x]$を近似するのでした(5.2.1項)。

$$ \begin{aligned} x^{(i)} &\sim \pi(x) \qquad (i = 1, 2, \dots, n) \\ \mathbb{E}_{\pi}[x] &\simeq \frac{x^{(1)} + x^{(2)} + \cdots + x^{(n)}}{n} \end{aligned} $$

 ここで、$\sim$は左辺の変数が右辺の分布に従うことを、$\simeq$は近似を表します。

・重点サンプリングによる期待値の近似式

 重点サンプリングでは、別の確率分布$b(x)$からサンプリングして$\mathbb{E}_{\pi}[x]$を求めることを考えます。
 $\mathbb{E}_{\pi}[x]$の計算式に$\frac{b(x)}{b(x)} = 1$を掛けて変形します。

$$ \begin{align} \mathbb{E}_{\pi}[x] &= \sum_x x \pi(x) \\ &= \sum_x x \frac{b(x)}{b(x)} \pi(x) \\ &= \sum_x x \frac{\pi(x)}{b(x)} b(x) = \mathbb{E}_b \left[ x \frac{\pi(x)}{b(x)} \right] \tag{5.7} \end{align} $$

 $\pi(x)$による期待値$\mathbb{E}_{\pi}[\cdot]$を、$b(x)$による期待値$\mathbb{E}_b[\cdot]$で表現できました。
 $x$の係数の項について

$$ \rho(x) = \frac{\pi(x)}{b(x)} $$

とおきます。

$$ \mathbb{E}_{\pi}[x] = \mathbb{E}_b[x \rho(x)] \tag{5.7'} $$

 $\rho(x)$を、$\mathbb{E}_{\pi}[x]$と$\mathbb{E}_b[x]$を対応付ける重みとみなします。

 よって、「$b(x)$のサンプル$x^{(i)}$」と「$x^{(i)}$の重み$\rho(x^{(i)})$」の積の標本平均により$\mathbb{E}_{\pi}[x]$を近似できます。

$$ \begin{aligned} x^{(i)} &\sim b(x) \qquad (i = 1, 2, \dots, n) \\ \mathbb{E}_{\pi}[x] &= \mathbb{E}_b[x \rho(x)] \simeq \frac{ x^{(1)} \rho(x^{(1)}) + x^{(2)} \rho(x^{(2)}) + \cdots + x^{(n)} \rho(x^{(n)}) }{ n } \end{aligned} $$

 以上で、重点サンプリングの式が得られました。

・プログラムで確認

 次は、通常のサンプリングと重点サンプリングを行って、期待値の近似を求めます。

 利用するライブラリを読み込みます。

# ライブラリを読み込み
import numpy as np
import matplotlib.pyplot as plt


・サンプリング

 モンテカルロ法により期待値を推定します。

 期待値を求めたい確率分布を設定します。

# 確率分布の設定
x = np.array([1, 2, 3]) # 確率変数
pi_x = np.array([0.1, 0.1, 0.8]) # 確率

 確率変数としてとり得る値$x$をx、対応する確率$\pi(x)$をpi_xとして値を指定します。

 設定した確率分布をグラフで確認します。

# 確率分布の棒グラフを作成
plt.figure(figsize=(8, 6))
plt.bar(x=x, height=pi_x) # 確率分布
plt.xticks(ticks=x)
plt.ylim(0, 1)
plt.xlabel('x')
plt.ylabel('probability')
plt.title('$\pi(x)$', fontsize=20)
plt.grid()
plt.show()

対象の確率分布


 $\pi(x)$の真の期待値を計算します。

# 真の期待値を計算
E_pi = np.sum(x * pi_x)
print(E_pi)
2.7


 モンテカルロ法により期待値の推定値を求めます。

# サンプルサイズを指定
N = 100

# サンプルの受け皿を初期化
samples = []

# サンプリング
for _ in range(N):
    # 確率分布πに従いサンプルを生成
    s = np.random.choice(a=x, p=pi_x)

    # サンプルを格納
    samples.append(s)

# 標本平均と標本分散を計算
mean = np.mean(samples)
var = np.var(samples)
print(mean)
print(var)
2.65
0.4475000000000001

 np.random.chioce()で乱数を生成して、リストsamplesに格納していきます。確率変数の引数ax、確率の引数ppi_xを指定します。

 標本平均(期待値の推定値)の推移をグラフで確認します。

# サンプルサイズ(試行回数)を作成
n = np.arange(1, N+1)

# 標本平均の推移を作図
plt.figure(figsize=(8, 6)) # 図の設定
plt.plot(n, np.cumsum(samples)/n, 
         label='sample mean') # 標本平均の推移
plt.hlines(y=E_pi, xmin=0, xmax=N, 
           color='red', linestyle='--', label='true mean') # 真の平均
plt.xlabel('iteration')
plt.ylabel('mean')
plt.suptitle('Sampling', fontsize=20)
plt.title('$x^{(i)} \sim \pi(x)$', loc='left')
plt.grid()
plt.legend()
plt.show()

標本平均の推移

 サンプルの累積和をnp.cumsum()で計算し対応する標本数で割って、試行ごとに標本平均を求めます。

 サンプルの度数をグラフで確認します。

# 出現回数を集計
x_smps, x_cnts = np.unique(samples, return_counts=True)

# サンプルのヒストグラムを作成
plt.figure(figsize=(8, 6))
plt.bar(x=x_smps, height=x_cnts, 
        label='sample') # サンプルの度数
plt.vlines(x=E_pi, ymin=0.0, ymax=x_cnts.max()*1.1, 
           color='red', linestyle='--', label='true mean') # 真の平均
plt.xticks(ticks=x)
plt.xlabel('x')
plt.ylabel('count')
plt.suptitle('Sampling', fontsize=20)
plt.title('$x^{(i)} \sim \pi(x)$', loc='left')
plt.grid()
plt.legend()
plt.show()

サンプルのヒストグラム


 ランダムな処理を行うため、実行する度に結果が異なります。そこで、複数回行った結果を確認します。

# シミュレーション回数を指定
runs = 200

# サンプルサイズを指定
N = 100

# 全てのシミュレーションのサンプルの受け皿を初期化
all_samples = np.empty((runs, N))

# 繰り返しシミュレーション
for run in range(runs):
    # サンプルの受け皿を初期化
    samples = []

    # サンプリング
    for _ in range(N):
        # 確率分布πに従いサンプルを生成
        s = np.random.choice(a=x, p=pi_x)

        # サンプルを格納
        samples.append(s)
    
    # 結果を格納
    all_samples[run] = samples
    
    # 標本平均と標本分散を計算
    mean = np.mean(samples)
    var = np.var(samples)
    print('run '+str(run+1)+' : mean='+str(np.round(mean, 2)) + ', var='+str(np.round(var, 2)))
run 1 : mean=2.71, var=0.43
run 2 : mean=2.63, var=0.47
run 3 : mean=2.71, var=0.43
run 4 : mean=2.7, var=0.43
run 5 : mean=2.56, var=0.59
(省略)
run 196 : mean=2.6, var=0.54
run 197 : mean=2.69, var=0.41
run 198 : mean=2.61, var=0.58
run 199 : mean=2.69, var=0.43
run 200 : mean=2.73, var=0.4

 モンテカルロ法による期待値の近似計算をrunsに指定した回数繰り返します。

 それぞれの推移を重ねて描画します。

# サンプルサイズ(試行回数)を作成
n = np.arange(1, N+1)

# 標本平均の推移を作図
plt.figure(figsize=(8, 6)) # 図の設定
for run in range(runs):
    plt.plot(n, np.cumsum(all_samples[run])/n, 
             alpha=0.1) # 標本平均の推移
plt.hlines(y=E_pi, xmin=0, xmax=N, 
           color='red', linestyle='--', label='true mean') # 真の平均
plt.xlabel('iteration')
plt.ylabel('mean')
plt.suptitle('Sampling', fontsize=20)
plt.title('$\pi=('+', '.join([str(np.round(val, 2)) for val in pi_x])+')$', loc='left')
plt.grid()
plt.legend()
plt.ylim(0, 10)
plt.show()

標本平均の推移

 どの結果も、真の期待値に近付いているのが分かります。

 次の内容と比較しやすいように、縦軸の幅を広くとって描画しています。

・重点サンプリング

 続いて、重点サンプリングにより期待値を推定します。

 サンプリングを行う確率分布を設定して、グラフを確認します。

# 確率分布の設定
b_x = np.array([1/3, 1/3, 1/3]) # 確率

# 確率分布の棒グラフを作成
plt.figure(figsize=(8, 6))
plt.bar(x=x, height=b_x) # 確率分布
plt.xticks(ticks=x)
plt.ylim(0, 1)
plt.xlabel('x')
plt.ylabel('probability')
plt.title('$b(x)$', fontsize=20)
plt.grid()
plt.show()

サンプルリング用の確率分布


 重点サンプリングにより期待値の推定値を求めます。

# サンプルサイズを指定
N = 100

# インデックス用の値を作成
idx = np.arange(len(b_x))

# サンプルの受け皿を初期化
samples = []

# サンプリング
for _ in range(N):
    # 確率分布bに従いサンプルを生成
    i = np.random.choice(a=idx, p=b_x)
    s = x[i]

    # 重みを計算
    rho = pi_x[i] / b_x[i]

    # 重み付けしたサンプルを格納
    samples.append(s * rho)

# 標本平均と標本分散を計算
mean = np.mean(samples)
var = np.var(samples)
print(mean)
print(var)
2.7330000000000005
10.294011000000003

 np.random.chioce()の確率変数の引数axのインデックス、確率の引数pb_xを指定します。生成したインデックスを使って、xから値を取り出します。
 また、対応する重みを計算して、重み付けした値を保存します。

 期待値の推定値の推移を確認します。

# サンプルサイズ(試行回数)を作成
n = np.arange(1, N+1)

# 標本平均の推移を作図
plt.figure(figsize=(8, 6)) # 図の設定
plt.plot(n, np.cumsum(samples)/n, 
         label='sample mean') # 標本平均の推移
plt.hlines(y=E_pi, xmin=0, xmax=N, 
           color='red', linestyle='--', label='true mean') # 真の平均
plt.xlabel('iteration')
plt.ylabel('mean')
plt.suptitle('Importance Sampling', fontsize=20)
plt.title('$x^{(i)} \sim b(x)$', loc='left')
plt.grid()
plt.legend()
plt.show()

重み付き標本平均の推移


 重み付けしたサンプルの度数を確認します。

# サンプルを重み付け
w_x = pi_x / b_x * x

# 出現回数を集計
x_smps, x_cnts = np.unique(samples, return_counts=True)

# サンプルのヒストグラムを作成
plt.figure(figsize=(8, 6))
plt.bar(x=x_smps, height=x_cnts, 
        width=0.25, label='sample') # サンプルの度数
plt.vlines(x=E_pi, ymin=0.0, ymax=x_cnts.max()*1.1, 
           color='red', linestyle='--', label='true mean') # 真の平均
plt.xticks(ticks=w_x)
plt.xlabel('$x \\rho(x)$')
plt.ylabel('count')
plt.suptitle('Importance Sampling', fontsize=20)
plt.title('$x^{(i)} \sim b(x)$', loc='left')
plt.grid()
plt.legend()
plt.show()

重み付きサンプルのヒストグラム


 重点サンプリングを繰り返し行います。

# シミュレーション回数を指定
runs = 200

# サンプルサイズを指定
N = 100

# インデックス用の値を作成
idx = np.arange(len(b_x))

# 全てのシミュレーションのサンプルの受け皿を初期化
all_samples = np.empty((runs, N))

# 繰り返しシミュレーション
for run in range(runs):
    # サンプルの受け皿を初期化
    samples = []

    # サンプリング
    for _ in range(N):
        # 確率分布bに従いサンプルを生成
        i = np.random.choice(a=idx, p=b_x)
        s = x[i]

        # 重みを計算
        rho = pi_x[i] / b_x[i]

        # 重み付けしたサンプルを格納
        samples.append(s * rho)
    
    # 結果を格納
    all_samples[run] = samples
    
    # 標本平均と標本分散を計算
    mean = np.mean(samples)
    var = np.var(samples)
    print('run '+str(run+1)+' : mean='+str(np.round(mean, 2)) + ', var='+str(np.round(var, 2)))
run 1 : mean=3.08, var=10.89
run 2 : mean=2.81, var=10.39
run 3 : mean=3.08, var=10.88
run 4 : mean=2.32, var=9.27
run 5 : mean=2.48, var=9.58
(省略)
run 196 : mean=2.74, var=10.25
run 197 : mean=2.7, var=10.0
run 198 : mean=2.75, var=10.24
run 199 : mean=2.61, var=9.93
run 200 : mean=3.15, var=10.93


 結果を確認します。

# サンプルサイズ(試行回数)を作成
n = np.arange(1, N+1)

# 標本平均の推移を作図
plt.figure(figsize=(8, 6)) # 図の設定
for run in range(runs):
    plt.plot(n, np.cumsum(all_samples[run])/n, 
             alpha=0.1) # 標本平均の推移
plt.hlines(y=E_pi, xmin=0, xmax=N, 
           color='red', linestyle='--', label='true mean') # 真の平均
plt.xlabel('iteration')
plt.ylabel('mean')
plt.suptitle('Importance Sampling', fontsize=20)
plt.title('$b=('+', '.join([str(np.round(val, 2)) for val in b_x])+')$', loc='left')
plt.grid()
plt.legend()
plt.ylim(0, 10)
plt.show()

重み付き標本平均の推移

 真の期待値を近似できているのが分かります。ただし、先ほどよりも散らばり具合が大きいです。次は、分散を小さくすることを考えます。

5.5.3 分散を小さくするには

 重点サンプリングにおいて分散が大きくなる理由を考えます。

・数式の確認

 まずは、真の分散と重点サンプリングによる分散を数式で確認します。

・分散の計算式

 分散は、次の式で定義されます。

$$ \mathbb{V}[x] = \mathbb{E}\Bigl[ (x - \mathbb{E}[x])^2 \Bigr] $$

 この式を展開して整理します。

$$ \begin{aligned} \mathbb{V}[x] &= \mathbb{E}\Bigl[ x^2 - 2 x \mathbb{E}[x] + (\mathbb{E}[x])^2 \Bigr] \\ &= \mathbb{E}[x^2] + - 2 \mathbb{E}[x] \mathbb{E}[x] + \mathbb{E} \Bigl[ (\mathbb{E}[x])^2 \Bigr] \\ &= \mathbb{E}[x^2] + - 2 (\mathbb{E}[x])^2 + (\mathbb{E}[x])^2 \\ &= \mathbb{E}[x^2] + - (\mathbb{E}[x])^2 \end{aligned} $$

 期待値の性質$\mathbb{E}[x + y] = \mathbb{E}[x] + \mathbb{E}[y]$、$\mathbb{E}[a x] = a \mathbb{E}[x]$、$\mathbb{E}[a] = a$を用いて変形しています。期待値$\mathbb{E}[x]$は定数です。
 または、確率変数と確率分布の積の総和の形に戻して、次のようにも求められます。

$$ \begin{aligned} \mathbb{V}[x] &= \sum_x \Bigl( x - \mathbb{E}[x] \Bigr)^2 p(x) \\ &= \sum_x \Bigl\{ x^2 - 2 x \mathbb{E}[x] + (\mathbb{E}[x])^2 \Bigr\} p(x) \\ &= \sum_x x^2 p(x) - 2 \mathbb{E}[x] \sum_x x p(x) + (\mathbb{E}[x])^2 \sum_x p(x) \\ &= \mathbb{E}[x^2] - 2 \mathbb{E}[x] \mathbb{E}[x] + (\mathbb{E}[x])^2 \\ &= \mathbb{E}[x^2] - (\mathbb{E}[x])^2 \end{aligned} $$

 総和の性質$\sum_x a x = a \sum_x x$、離散型確率分布の性質$\sum_x p(x) = 1$を用いて変形しています。

 分散は、「2乗の期待値$\mathbb{E}[x^2]$」と「期待値の2乗$(\mathbb{E}[x])^2$」の差で計算できるのが分かりました。

・重点サンプリングの分散

 重点サンプリングでは、サンプリング用の確率分布$b(x)$を用いて重み付きの期待値を求めるのでした。

$$ \mathbb{E}_{\pi}[x] = \mathbb{E}_b[x \rho(x)] $$

 重み$\rho(x)$は、2つの分布の比でした。

$$ \rho(x) = \frac{\pi(x)}{b(x)} $$

 重み付き変数$x \rho(x)$の分散を考えます。
 分散の計算式に$x \rho(x)$を代入します。

$$ \mathbb{V}_b[x \rho(x)] = \mathbb{E}_b \Bigl[ (x \rho(x))^2 \Bigr] - \Bigl( \mathbb{E}_b[x \rho(x)] \Bigr)^2 $$

 前の項は

$$ \begin{aligned} \mathbb{E}_b \Bigl[ (x \rho(x))^2 \Bigr] &= \sum_x \left( x \frac{\pi(x)}{b(x)} \right)^2 b(x) \\ &= \sum_x x^2 \frac{\pi(x) \pi(x)}{b(x) b(x)} b(x) \\ &= \sum_x x^2 \frac{\pi(x)}{b(x)} \pi(x) \\ &= \sum_x x^2 \rho(x) \pi(x) = \mathbb{E}_{\pi}[x^2 \rho(x)] \end{aligned} $$

と変形できます。
 また、後の項は$\mathbb{E}_{\pi}[x]$に置き換えられるので、それぞれ$\pi(x)$による期待値の項に置き換えます。

$$ \mathbb{V}_b[x \rho(x)] = \mathbb{E}_{\pi}[x^2 \rho(x)] - (\mathbb{E}_{\pi}[x])^2 $$

 重点サンプリングの分散$\mathbb{V}_b[x \rho(x)]$を、$\pi(x)$の分散

$$ \mathbb{V}_{\pi}[x] = \mathbb{E}_{\pi}[x^2] - (\mathbb{E}_{\pi}[x])^2 $$

と比較すると、前の項に関して$\rho(x)$倍されています。
 つまり、$\rho(x)$が大きいほど分散が大きくなります。

・プログラムで確認

 次は、真の分散と重点サンプリングの分散を計算した後に、5.5.2項の重点サンプリングを行います。

・期待値と分散の計算

 $\pi(x), b(x)$の値と重みの関係をグラフで確認します。

# 真の確率を指定
pi = 0.5

# サンプリング確率の値を作成
b_vals = np.linspace(0.0, 1.0, num=101)[1:] # 0を除く

# 重みを計算
rho_vals = pi / b_vals

# サンプリング確率と重みの関係を作図
plt.figure(figsize=(8, 6))
plt.plot(b_vals, rho_vals, label='$\\rho(x)$') # 重み
plt.vlines(x=pi, ymin=0.0, ymax=rho_vals.max(), 
           color='red', linestyle='--', label='$\pi(x)$') # 真の確率
plt.xlabel('$b(x)$')
plt.ylabel('$\\rho(x)$')
plt.title('$\\rho(x) = \\frac{\pi(x)}{b(x)}$', fontsize=20)
plt.grid()
plt.legend()
plt.ylim(0, 5)
plt.show()

確率分布と重みの関係

 $\pi(x) > b(x)$のとき$\rho(x) > 1$、$\pi(x) = b(x)$のとき$\rho(x) = 1$、$\pi(x) < b(x)$のとき$\rho(x) < 1$になります。
 ただし、$x$は複数の値をとります。例えば、$\pi(1) = 0.5, \pi(2) = 0.5$のとき、$b(1) = 0.8$とすると$\rho(1) = 0.625$となり1より小さくなりますが、$b(2) = 0.2$なので$\rho(2) = 2.5$となり、重点サンプリングの分散は真の分散より大きくなります。
 よって、$\pi(x)$と$b(x)$の値が近いほど分散も近くなるのが分かります。

 $\pi(x)$を設定して、期待値と分散を計算します。

# 確率分布の設定
x = np.array([1, 2, 3]) # 確率変数
pi_x = np.array([0.1, 0.1, 0.8]) # 確率

# 期待値を計算
E_pi = np.sum(x * pi_x)
print(E_pi)

# 2乗の期待値を計算
E_pi2 = np.sum(x**2 * pi_x)
print(E_pi2)

# 分散を計算
V_pi = E_pi2 - E_pi**2
print(V_pi)
2.7
7.7
0.40999999999999925


 続いて、$b(x)$を設定して、重み付けした期待値と分散を計算します。

# 確率分布の設定
#b_x = pi_x.copy()
b_x = np.array([1.0/3.0, 1.0/3.0, 1.0/3.0]) # 確率
#b_x = np.array([0.2, 0.2, 0.6]) # 確率

# 重みを計算
rho_x = pi_x / b_x

# 期待値を計算
E_b = np.sum(x * rho_x * b_x)
print(E_b)

# 2乗の期待値を計算
E_b2 = np.sum(x**2 * rho_x * pi_x)
print(E_b2)

# 分散を計算
V_b = E_b2 - E_b**2
print(V_b)
2.7
17.43
10.139999999999999

 pi_xb_xを近付けると分散が小さくなる(真の分散に近付く)のを確認できます。また、期待値の推定時の標本分散(samplesの分散)と近い値なのも分かります。

・重点サンプリング

 対象の確率分布pi_xに近付けたb_xを設定します。

# 確率分布の設定
b_x = np.array([0.2, 0.2, 0.6])


 5.5.2項のコードで、繰り返し重点サンプリングを行います。

run 1 : mean=2.61, var=2.69
run 2 : mean=2.68, var=2.62
run 3 : mean=2.94, var=2.32
run 4 : mean=2.91, var=2.43
run 5 : mean=2.7, var=2.56
(省略)
run 196 : mean=2.68, var=2.53
run 197 : mean=2.97, var=2.27
run 198 : mean=2.74, var=2.49
run 199 : mean=2.8, var=2.5
run 200 : mean=2.86, var=2.46


 結果を確認します。

重み付き標本平均の推移

 分散が小さくなっているのが分かります。

 この節では、重点サンプリングを確認しました。次節では、重点サンプリングを用いて行動価値関数を推定します。

参考文献


おわりに

 5章の内容はこれで終了ですが、ここで終わるともったいないです。この内容が分かれば難しくないので、付録Aもやってしまいましょう。

【次節の内容】

www.anarchive-beta.com