重点サンプリング

Pythonに慣れることを目的として、まずは比較的簡単なところからやっていきたい。ということで、重点サンプリングをまずはやる。

以下のコードにより、一様分布を提案分布として、重点サンプリングにてt分布が得られる。

PythonやRには任意の自由度のt分布から乱数を発生させる関数が用意されているため、こんなことをやる必要は全く無いが、そこは勉強のためなのでスルー

# -*- coding: utf-8 -*-
import matplotlib.pyplot as plt
import numpy as np
import seaborn as sns

#サンプリング数
N=100000
#自由度
nu=5

#ターゲット分布である自由度nuのt分布から乱数をN個発生させる
y=np.random.standard_t(nu, size=N)
plt.hist(y,bins=100,normed=True)
maxx=max(y)
minx=min(y)

#提案分布はターゲット分布の最大、最小区間を取る一様分布とする
x=np.random.rand(N)*(maxx-minx)+minx
plt.hist(x,bins=100,normed=True)

#ターゲット分布の確率密度÷提案分布の確率密度にてウェイトを計算
w=stats.t.pdf(x,nu)/(1/(maxx-minx))
#ウェイトを標準化する
w=w/sum(w)
#標準化したウェイトにてリサンプリング
k=np.random.choice(range(0,N),N,p=w)
x=x[k]
plt.hist(x,bins=100,normed=True)

#サンプリングされた分布のモーメントが理論値と合うか確認
np.mean(x)
np.std(x)
np.sqrt(nu)/(np.sqrt(nu-2))

これを実行すると、こんなグラフが得られる。

グラフの青は正解用にNumPyの乱数を使って生成したt分布、緑は提案分布の一様分布、赤が提案分布を重点サンプリングすることによって得られたt分布。
青と赤はほぼ重なっているため、想定通りサンプリングできたと考えていいだろう。
また、平均もほぼ0で、標準偏差も理論値と近いものとなっている。