##KDEとは？
”カーネル密度推定（カーネルみつどすいてい、英: kernel density estimation）は、統計学において、確率変数の確率密度関数を推定するノンパラメトリック手法のひとつ（Wikipedia）”とされており、機械学習などで様々に応用されています。KDEはscipyやseaborn、pandasに実装されている為、これらを使えば簡単にKDEをプロットすることができますが、実際のところ何をしてるのかイマイチよく分かっていないブラックボックスとして使ってしまっていました。本記事では、KDEを一から実装しつつ、パラメータの意味などについても整理してみたいと思います。
##KDEの定義
${x_1,x_2,...,x_n}$を独立かつ同一な分布に従うサンプルとして、KDEではこれらのサンプルの従う分布の確率密度関数fの近似$\hat{f_h(x)}$を以下の式によって推定します。
\hat{f_h(x)}=\frac{1}{nh}\sum_{i=1}^{n}K(\frac{x-x_i}{h})
ここで$h$はバンド幅、$K(v)$はカーネル関数と呼ばれます。カーネル関数としてガウス分布が有名ですが、他の確率密度関数を使うこともできます。pythonでの実装は非常にシンプルです。
#sampleはKDEの対象データ
def kde(x,sample,band_width,kernel):
n=len(sample)
return np.sum([1/(n*band_width)*kernel((x-sample[i])/band_width) for i in range(n)])
KDEの結果を確率密度関数として解釈できる理由としては、各カーネル関数が確率密度関数であることから理解できます。カーネル関数は非負なので、それらの和も非負になります。また、$\hat{f_h(x)}$の積分についても、以下に示す様に$1$に等しくなります（途中で置換積分と、確率密度関数の定義（全範囲での積分が$1$）を使っています）。
\int^{\infty}_{-\infty}\hat{f_h(x)}dx
=\int^{\infty}_{-\infty}\frac{1}{nh}\sum_{i=1}^{n}K(\frac{x-x_i}{h})dx\\
=\frac{1}{nh}\sum_{i=1}^{n}\int^{\infty}_{-\infty}K(\frac{x-x_i}{h})dx
=\frac{1}{nh}\sum_{i=1}^{n}\int^{\infty}_{-\infty}K(t)hdt\\
=\frac{1}{n}\sum_{i=1}^{n}1
=\frac{1}{n}*n
=1
では、KDEは実際のところ何をやっているのでしょう…？KDEではカーネル関数の総和を確率密度関数の推定として採用します。その様子を見える為に、正規分布から発生させたサイズ$15$のサンプルについてKDEを行い、それぞれのカーネル関数の寄与を図示してみましょう。
import numpy as np
import matplotlib.pyplot as plt
import scipy.stats as stats
import random
random.seed(10)
np.random.seed(1)
sample_size=15
mean=3 #サンプルを取る正規分布の平均
sigma=2 #サンプルを取る正規分布の分散
#正規分布の定義
def normal(x,mean,sigma):
return 1/np.sqrt(2*np.pi*sigma**2)*np.exp(-((x-mean)/sigma)**2)
#サンプルの発生
y=list(np.random.normal(loc=mean,scale=sigma,size=int(sample_size)))
sample=sorted(y)
fig,ax=plt.subplots(nrows=2,figsize=(20,15))
#サンプルを取った正規分布の描画
ax1=ax[0]
x_kernel=np.linspace(np.min(y),np.max(y))
ax1.plot(x_kernel,[normal(x_kernel[i],mean,sigma) for i in range(len(x_kernel))],label='source',linewidth=5,color='cyan')
#ガウスカーネル関数の定義
def normal_kernel(x):
return 1/np.sqrt(2*np.pi)*np.exp(-x**2/2)
#Scott's ruleによるバンド幅の設定（後述）
band_width=np.sqrt(np.var(sample,ddof=1)*((sample_size)**(-1/5))**2)
cmap=["#"+''.join([random.choice('0123456789ABCDEF') for j in range(6)]) for i in range(len(sample))]
#各カーネル関数へ(x-x_i)/hを代入し、1/nhを掛けた関数の描画
output=[]
for i in range(len(sample)):
y_kernel=1/(sample_size*band_width)*normal_kernel((x_kernel-sample[i])/band_width)
output.append(y_kernel)
ax1.plot(x_kernel,y_kernel,color=cmap[i])
ax1.vlines(sample[i],0,0.1) #x_iの点の描画
ax1.legend(fontsize=30)
ax2=ax[1]
#それぞれのカーネル関数の寄与を描画
ax2.stackplot(x_kernel,output,colors=cmap)
#KDEの結果を描画
x=np.linspace(np.min(sample),np.max(sample))
y=[kde(x[i],sample,band_width,normal_kernel) for i in range(len(x))]
ax2.plot(x,y,color='k',linewidth=10,label='kde',linestyle='dashed')
#サンプルを取った正規分布の描画
y=[normal(x[i],mean,sigma) for i in range(len(x))]
ax2.plot(x,y,color='cyan',linewidth=10,linestyle='dashed',label='source')
ax2.legend(fontsize=30)
若干見にくいですが、上の図からそれぞれのサンプル点の周囲に正規分布の形状のグラフが発生しており、下の図からこれらのグラフが集中している場所（元の正規分布の平均付近）ではより多くのグラフが寄与することで、KDEの結果に元の正規分布のピークが再現されていることが分かります。このことから、KDEは”データが集中しているところには、元の分布のピークがあったんだろう”というアイデアに基づいていると考えられます。
##KDEのパラメータ（バンド幅・カーネル関数）
###バンド幅について
KDEではカーネル関数の変数へはある点$x$と、$i$番目のサンプル点$x_i$の差をバンド幅$h$で割ったものが代入されています。このことから、同じ$x-x_i$の値であっても、バンド幅が大きい程カーネル関数へと代入される値は小さくなることが分かります。ガウス分布に代表されるカーネル関数は原点（$v=0$、すなわち$x=x_i$の点）から離れるほど値が小さくなることから、バンド幅を大きくすることで、$x_i$から離れた値が効いてくる様になります。最終的にこれらのカーネル関数の総和がKDEによる確率密度関数の推定結果となることから、大きなバンド幅はより滑らかな推定になることが期待できます。バンド幅を変化させた時の$\frac{1}{nh}K(\frac{x-x_i}{h})$の様子を見てみましょう。
#3種類のバンド幅のオプション
h=[0.05,0.3,1]
#サンプル数５
sample_size=5
y=list(np.random.normal(loc=mean,scale=sigma,size=int(sample_size)))
sample=sorted(y)
x_kernel=np.linspace(np.min(y)-1,np.max(y)+1,500)
cmap=["#"+''.join([random.choice('0123456789ABCDEF') for j in range(6)]) for i in range(len(sample))]
fig,ax=plt.subplots(nrows=3,figsize=(20,15))
#それぞれのバンド幅についてグラフを描く
for k in range(3):
h_k=h[k]
ax_k=ax[k]
for i in range(len(sample)):
y_kernel=[1/(sample_size*h_k)*normal_kernel((x_kernel[j]-sample[i])/h_k) for j in range(len(x_kernel))]
ax_k.plot(x_kernel,y_kernel,color=cmap[i])
ax_k.vlines(sample[i],0,np.max(y_kernel)*1.1)
ax_k.set_title(f'bandwidth={h_k}',fontsize=20)
これらのグラフより、バンド幅が広いほどより個々のカーネル関数の重なりが大きくなり、これらの総和で表現されるKDEの結果も滑らかなものとなると考えられます。では、適切なバンド幅ははどの様に設定すればよいのでしょうか？scipyでも使われている"Rule-of-thumb"として、Scott's-ruleとSilverman's-ruleが挙げられます。こちらのサイトを参考とさせて頂いたところ、一次元データの場合はそれぞれのルールに基づいてバンド幅は以下の様に推定されます。
Scott's\quad rule\\
h=n^{\frac{-1}{5}}*\sqrt{\hat{\sigma}}\\
Silverman's\quad rule\\
h=(\frac{3}{4}n)^{-\frac{1}{5}}*\sqrt{\hat{\sigma}}\\
*\hat{\sigma}は対象データの不偏分散、nは対象データの数
これらのルールに基づいて決定したバンド幅を使用した場合と使用しない場合のそれぞれについて、KDEの結果を比較してみましょう。
#Scott's ruleによるバンド幅
h_scott=np.sqrt(np.var(sample,ddof=1)*((sample_size)**(-1/5))**2)
#Silverman's ruleによるバンド幅
h_silverman=np.sqrt(np.var(sample,ddof=1)*((3/4*sample_size)**(-1/5))**2)
h=[0.05,0.3,3]
sample_size=300
y=list(np.random.normal(loc=mean,scale=sigma,size=int(sample_size)))
sample=sorted(y)
x=np.linspace(np.min(y)-1,np.max(y)+1,100)
fig,ax=plt.subplots(nrows=5,figsize=(20,20))
for k in range(3):
h_k=h[k]
ax_k=ax[k]
for i in range(len(sample)):
y_kde=[kde(x[i],sample,h_k,normal_kernel) for i in range(len(x))]
ax_k.plot(x,y_kde,color='red')
y=[normal(x[i],mean,sigma) for i in range(len(x))]
ax_k.plot(x,y,color='k')
ax_k.set_title(f'bandwidth={h_k}',fontsize=20)
#Scott's ruleによるKDEの描画
ax[3].plot(x,[kde(x[i],sample,h_scott,normal_kernel) for i in range(len(x))])
ax[3].plot(x,[normal(x[i],mean,sigma) for i in range(len(x))],color='k')
ax[3].set_title(f'Scott: {h_scott:3f}',fontsize=20)
#Silverman's ruleによるKDEの描画
ax[4].plot(x,[kde(x[i],sample,h_silverman,normal_kernel) for i in range(len(x))])
ax[4].plot(x,[normal(x[i],mean,sigma) for i in range(len(x))],color='k')
ax[4].set_title(f'Silverman: {h_silverman:3f}',fontsize=20)
バンド幅が小さい場合にはギザギザとした近似であり、バンド幅が大きすぎると滑らかな近似になりすぎてしまうことが分かります。また、Scott's rule、Silverman's ruleのいずれについても、それなりに良い近似を与えていることが観察できます。
###カーネル関数について
ガウスカーネルは広く使われており、scipyのKDEではgaussian_kdeとして実装されています。しかしながら、他にも利用可能なカーネル関数はいくつかあることから、これらについて比較・検討してみたいと思います。サンプルデータとして、正規分布の混合分布を用います。
mean1=3
mean2=15
sigma1=2
sigma2=4
sample_size=3000
def normal(x,mean,sigma):
return 1/np.sqrt(2*np.pi*sigma**2)*np.exp(-((x-mean)/sigma)**2)
def normal_kernel(x):
return 1/np.sqrt(2*np.pi)*np.exp(-x**2/2)
sample=[]
for i in range(sample_size):
idx=np.random.randint(1,3)
if idx==1:
sample.append(np.random.normal(loc=mean1,scale=sigma1))
else:
sample.append(np.random.normal(loc=mean2,scale=sigma2))
sample=sorted(sample)
fig,ax=plt.subplots()
ax.hist(sample,density=True,bins=100)
ax.set_xlim(np.min(sample),np.max(sample))
x=np.linspace(0,np.max(sample))
ax.plot(x,0.5*normal(x,mean1,sigma1)+0.5*normal(x,mean2,sigma2))
- ガウスカーネル
非常によく知られた標準正規分布をカーネル関数に用いるものです。scipyのgaussian_kdeによるKDEの結果と、実装したKDEの結果を比較してみましょう。
#Scott's ruleによるバンド幅の設定
band_width=np.sqrt(np.var(sample,ddof=1)*((sample_size)**(-1/5))**2)
#実装したKDEの結果
x=sorted(np.linspace(np.min(sample),np.max(sample),200))
y=[kde(x[i],sample,band_width,normal_kernel) for i in range(len(x))]
#ScipyによるKDEの結果
kernel_scipy=stats.gaussian_kde(sample)
y_scipy=kernel_scipy(sample)
fig,ax=plt.subplots(figsize=(10,7))
ax.plot(x,y,label='kde')
ax.plot(sample,y_scipy,linestyle='dashed',label='scipy')
ax.legend(fontsize=15)
どうやらほぼ一致している様です。カーネル関数の形と、もともとの確率密度関数の形状を比較してみましょう。
fig,ax=plt.subplots(nrows=2,figsize=(15,10))
#ガウスカーネルの描画
ax1=ax[0]
x=np.linspace(-5,5)
y=normal_kernel(x)
ax1.plot(x,y)
ax2=ax[1]
x=np.linspace(0,np.max(sample),200)
band_width=np.sqrt(np.var(sample,ddof=1)*((sample_size)**(-1/5))**2)
y=[kde(x[i],sample,band_width,normal_kernel) for i in range(len(x))]
ax2.plot(x,y,label='kde')
ax2.plot(x,0.5*normal(x,mean1,sigma1)+0.5*normal(x,mean2,sigma2),label='source')
ax2.legend(fontsize=15)
ピークの位置は上手く拾えていますが、広がりについてはもう一声という雰囲気ですね。
2.Top-hatカーネル
一様分布を採用したものです。カーネル関数の形がシルクハット（top hat）に似ていることから、この名前がついていると思われます。
def tophat_kernel(x):
if -1<x<1:
return 0.5
else:
return 0
fig,ax=plt.subplots(nrows=2,figsize=(15,10))
ax1=ax[0]
x=np.linspace(-5,5,200)
y=[tophat_kernel(x[i]) for i in range(len(x))]
ax1.plot(x,y)
ax2=ax[1]
x=np.linspace(0,np.max(sample),200)
band_width=np.sqrt(np.var(sample,ddof=1)*((sample_size)**(-1/5))**2)
y=[kde(x[i],sample,band_width,tophat_kernel) for i in range(len(x))]
ax2.plot(x,y,label='kde')
ax2.plot(x,0.5*normal(x,mean1,sigma1)+0.5*normal(x,mean2,sigma2),label='source')
ax2.legend(fontsize=15)
ガウスカーネルと比べてギザギザとしていますが、ピークの高さと位置はよく合っていそうな雰囲気です。
3.Linearカーネル
線形な関数をカーネル関数へと採用したものです。
def linear_kernel(x):
if 0<abs(x)<1:
return 1-abs(x)
else:
return 0
fig,ax=plt.subplots(nrows=2,figsize=(15,10))
ax1=ax[0]
x=np.linspace(-5,5,200)
y=[linear_kernel(x[i]) for i in range(len(x))]
ax1.plot(x,y)
ax2=ax[1]
x=np.linspace(0,np.max(sample),200)
band_width=np.sqrt(np.var(sample,ddof=1)*((sample_size)**(-1/5))**2)
y=[kde(x[i],sample,band_width,linear_kernel) for i in range(len(x))]
ax2.plot(x,y,label='kde')
ax2.plot(x,0.5*normal(x,mean1,sigma1)+0.5*normal(x,mean2,sigma2),label='source')
ax2.legend(fontsize=15)
Top-hatより滑らかですし、ガウスカーネルと比べて高さの合致も良いので、個人的には試したカーネルの中で一番良さそうな雰囲気を受けました。
4.Exponentialカーネル
指数関数を採用したものです。
def exponential_kernel(x):
return 0.5*np.exp(-abs(x))
fig,ax=plt.subplots(nrows=2,figsize=(15,10))
ax1=ax[0]
x=np.linspace(-5,5,200)
y=[exponential_kernel(x[i]) for i in range(len(x))]
ax1.plot(x,y)
ax2=ax[1]
x=np.linspace(0,np.max(sample),200)
band_width=np.sqrt(np.var(sample,ddof=1)*((sample_size)**(-1/5))**2)
y=[kde(x[i],sample,band_width,exponential_kernel) for i in range(len(x))]
ax2.plot(x,y,label='kde')
ax2.plot(x,0.5*normal(x,mean1,sigma1)+0.5*normal(x,mean2,sigma2),label='source')
ax2.legend(fontsize=15)
ガウスカーネルにソックリですね。正直違いが分かりませんでした。
5.Epanechnikovカーネル
この論文でカーネル関数の例として扱われています。定義は以下の様になっています。
K(v)=\frac{3}{4\sqrt{5}}(1-\frac{x^2}{5}) \quad (|x| \leq\sqrt{5})
def epanechnikov_kernel(x):
if 0<abs(x)<np.sqrt(5):
return (3/(4*np.sqrt(5)))*(1-(x**2)/5)
else:
return 0
fig,ax=plt.subplots(nrows=2,figsize=(15,10))
ax1=ax[0]
x=np.linspace(-5,5,200)
y=[epanechnikov_kernel(x[i]) for i in range(len(x))]
ax1.plot(x,y)
ax2=ax[1]
x=np.linspace(0,np.max(sample),200)
band_width=np.sqrt(np.var(sample,ddof=1)*((sample_size)**(-1/5))**2)
y=[kde(x[i],sample,band_width,exponential_kernel) for i in range(len(x))]
ax2.plot(x,y,label='kde')
ax2.plot(x,0.5*normal(x,mean1,sigma1)+0.5*normal(x,mean2,sigma2),label='source')
ax2.legend(fontsize=15)
これもガウスカーネルと違いがよく分かりませんでした。何か数学上便利な性質はありそうですが。
##おわりに
本記事ではKDEについてnumypyを使って実装し、その中身やパラメータの意味について整理してみました。バンド幅については経験則的な方法で十分と思われますが、カーネル関数はガウスカーネルとほぼ変わらない挙動を示すもの（Exponentinal、Epanechnikov）、ガウスカーネルと異なった挙動を示すもの（Top-hat、Linear）があったことから、ガウスカーネルで望ましい結果が得られない場合には異なる種類のカーネル関数を試すのも良いかもしれません。ここまでお読みいただき、ありがとうございました。何か誤り等ございましたら、コメントにてお伝え頂けますと幸いです。よろしくお願いいたします。