nykergoto’s blog

機械学習とpythonをメインに、雑多な内容をとりとめなく扱うブログです。

カーネル密度推定

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というクラス作ってみました、完全に練習ですが少しは見通しが良くなってるかな?
f:id:dette:20150723191716p:plain
sigmaの値を変えると推定の雰囲気が色々と変わるので面白いですね。