カーネル密度推定法の式(13)を使って、実際に確率分布を推定してみたいと思います。推定する対象は、棄却サンプリングを使って、適当な形状の分布を用いることにします。
早速結果です
青色の線がサンプリングした確率分布で、オレンジの線がカーネル密度推定法で推定した確率分布です。のパラメータは3水準で実験してみました。カーネル密度推定法で書いた通り、
の値を大きくとると推定される分布形状がなだからかになり、値を小さくとるとサンプル点に過敏に反応した形状となっています。せっかく棄却サンプリングを使ったので、もうちょっと面白い形状の分布からサンプリングすればよかった。
実際に使うときははどういう基準で決めればいいのかな。グラフで可視化できるようなものなら何とでもなりそうですが、多次元の場合だと難しい。
今回のコードです。
# カーネル密度推定 import numpy as np import matplotlib.pyplot as plt # サンプリングする分布（区間[-1,1]で積分してだいたい1になるように調整した） def p(x): return np.cos(0.5*np.pi*x)*np.abs(3*x-3)*np.abs(3*x-3)/13 # 棄却サンプリング def sampling(): k = 3 # 採択するまでループ while True: # 提案分布q(z)からサンプリング z = np.random.uniform(-1, 1) # [0,kq(z)]の一様分布からサンプリング u = k*np.random.uniform(0, 0.5) # 棄却するか判定 if p(z) > u: return z # ガウスカーネルを使ったカーネル密度推定 def kernel_density_estimation(x, sample): h = 0.1 return 1.0/(np.sqrt(2*np.pi*h**2))*np.exp(-(x-sample)**2/(2*h**2))/N # ランダムシードを固定 np.random.seed(0) # サンプル数 N = 50 # 棄却サンプリング samples = np.array([sampling() for i in range(N)]) # カーネル密度推定で確率分布を推定する dens = np.array([np.sum(np.vectorize(kernel_density_estimation)(i, samples)) for i in np.linspace(-1, 1, 100)]) x = np.linspace(-1, 1, 100) y = p(x) plt.title("h=0.1") plt.xlim(-1, 1) plt.plot(x, y) plt.plot(np.linspace(-1, 1, 100), dens) plt.show()
