カーネル密度推定
PRMLを読んでいます。とりあえず出てくる式はすべて自分で追って、変形とかもやってるのでなかなか進まなくてもどかしいですが、もう多分確率密度関数とかを復習することは無いだろうし…と思いながらやってます。
今日は二章のおわりまで行きました。
ノンパラメトリックの推定方法のところで、この間「外野フライをカーネル密度推定を用いて解析する」っていう論文を読んで気になってたカーネル密度推定のお話が出てきました。ざっくり言うと、その点の周りにいっぱい点があれば、そこの確率は高いよねっていうことです。発想としては至ってシンプルなんですね。直感的。
なんか可視化してみたら面白そうだったので、pythonの練習も兼ねて実装をしてみました。
やってるのは、一次元のガウス分布に従うサンプルから確率密度関数を推定するっていうものです。
import numpy as np import matplotlib.pyplot as plt from scipy.stats import multivariate_normal class KernelEstimation(): def __init__(self,trainingdata,h,D,kernel): self.trainingdata = trainingdata self.h = h self.D = D self.kernel = kernel def labeltext(self): return "kernel: "+self.kernel+" sigma: "+str(self.h) def _kernel(self,d): """ return kernel distance """ if self.kernel == "gausian": norm = np.linalg.norm(d) return np.exp(-norm*norm/(2.*self.h*self.h)) return 0 def _normalize(self): """ normalization parametor """ if self.kernel == "gausian": return 1/np.power((2*self.h*self.h*np.pi),(self.D/2.))/self.trainingdata.size else: return 0 def predict(self,x): """ predict by kernal estimate """ prob = np.zeros(x.size) for i in range(x.size): for k in range(self.trainingdata.size): prob[i] += self._kernel(x[i]-self.trainingdata[k]) return prob*self._normalize() rv = multivariate_normal(1,1) x = np.linspace(-5,5,100) plt.plot(x,rv.pdf(x),"k-") samples = rv.rvs(100)#sampleの生成 plt.hist(samples,normed=True, facecolor='green',alpha=0.3) kernel = KernelEstimation(samples,0.30,1,"gausian") plt.plot(x,kernel.predict(x),"b--",label=kernel.labeltext()) plt.legend() plt.show()
KernelEstimationというクラス作ってみました、完全に練習ですが少しは見通しが良くなってるかな?
sigmaの値を変えると推定の雰囲気が色々と変わるので面白いですね。