총 면적이 1인 히스토그램은 확률변수의 밀도를 보여준다 할 수 있다. 이것을 매끄럽게 표현하기 위해 고안한 방법이 커널 밀도 추정(kernel density estimation)이라 생각할 수 있다.

📌커널함수

커널 밀도함수 추정을 위해서 먼저 커널함수에 대해 알아보자. 커널함수란 적분값이 1이고, 원점대칭인 양함수를 의미한다. 이를 수식으로 나타내면 아래와 같다.

$$\int\limits_{-\infty}^{\infty} K(u)=1$$

$$K(-u)=K(u) \ for \ all \ values \ of \ u.$$

📌커널함수 밀도 추정

커널밀도함수(Kernel density estimation, KDE)는 커널함수를 바탕으로 확률밀도를 추정하는 것을의미한다. $x_1, \cdots , x_n$을 미지의 밀도함수 $f(x)$로부터의 임의 표본이라 하면, $\hat{f}(x)$는 다음과 같다.

$$\hat{f}(x)= {1\over n}\sum\limits_{i=1}^n {1\over b}K({{x-x_i}\over b}) $$

$K(z)$는 커널함수를 의미하고, $b>0$를 띠너비(bandwidth)라 한다. 띠너비는 작아질수록 커널함수는 뾰족해진다는 특징이있다. 가우시안커널(방사형기저함수)의 수식은 아래와 같다.

$$ K(z)={1\over{\sqrt{2\pi}}} e^{- {{z^2}\over 2}}$$

수식을 간단히 정리하면 관측자료마다 커널함수를 생성한 후, 커널 함수를 모두 더하고 자료수로 나누어 산출된다.

🌈Python 구현

import numpy as np
from sklearn.datasets import load_iris 
import pandas as pd
import random
import matplotlib.pyplot as plt

iris = load_iris() # sample data load
df = pd.DataFrame(data=iris.data, columns=iris.feature_names)
df['Species'] = iris.target
df['Species'] = df['Species'].map({0:"setosa", 1:"versicolor", 2:"virginica"})

def kernel(z):
    return (1/np.sqrt(2*np.pi)*np.exp(-(z**2)/2))
    
def kernel_density(df,b=.3304):
    f_hat=0
    for x_i in df:
        f_hat+=kernel((x-x_i)/b)/(len(df)*b)
    return f_hat

plt.hist(df['sepal length (cm)'],8,density=True)
plt.plot(x,kernel_density(df['sepal length (cm)'],b=.2871746))
plt.xlim(4,8)

🌈R 구현

library(MASS)
x=sort(runif(1000,-(min(iris$Sepal.Length)+1),
                    max(iris$Sepal.Length)+1))
x

#gaussian kernel
kernel=function(z){
    return(1/sqrt(2*pi)*exp(-(z^2)/2))
}
kernel_density=function(df,b=.3304){
    f_hat=0
    for(x_i in df){
        f_hat=f_hat+kernel((x-x_i)/b)/(length(df)*b)
    }
    return(f_hat)
}

hist(iris$Sepal.Length,freq=F,xlim=c(4,8))
lines(x,
     kernel_density(
         df=iris$Sepal.Length,b=.2871746),
         xlim=c(4,8),lty=3,col=2)
lines(density(iris$Sepal.Length),lty=2,col=3)

🌈2차원 확장

$$\hat{f}(x,y)= {1\over n}\sum\limits_{i=1}^n {1\over b_1}K({{x-x_i}\over b_1}) {1\over b_2}K({{x-x_i}\over b_2}) $$

x=iris$Sepal.Length
y=iris$Sepal.Width
nx <- length(x)
lims=c(range(x),range(y))
n <- rep(n, length.out = 2L)
gx <- seq.int(lims[1L], lims[2L], length.out = n[1L])
gy <- seq.int(lims[3L], lims[4L], length.out = n[2L])

h=c(bandwidth.nrd(x), bandwidth.nrd(y))
h <- h/4
ax <- outer(gx, x, "-")/h[1L]
ay <- outer(gy, y, "-")/h[2L]

z <- tcrossprod(
    matrix(dnorm(ax), , nx), 
    matrix(dnorm(ay), , nx))/(nx * h[1L] * h[2L])
z