[data analysis] 커널밀도추정(kernel density estimator)

커널밀도추정(kernel desity estimation)

내용

• 밀도추정(Density estimation)
• Histogram
• Kernel Density Estimation
• python을 적용한 kde의 계산

밀도추정(Density estimation)

확률, 통계에서 밀도(density)는 확률(probability)를 나타냅니다. 그러므로 밀도 추정은 관찰된 데이터에서 관찰치의 확률을 추정하는 것으로 밀도 함수의 추정치를 구성하는 것입니다. 즉, 다음 식과 같이 일정한 구간에 포함된 포인트의 개수(빈도수)나 면적등을 계산하는 것으로 관찰치의 밀도를 추정할 수 있습니다.

$$\begin{array}{}\text{(식 1)}& P(a<x<b)={\int }_a^{b}f(d)dx\mathrm{\forall }a<b\end{array}$$

식 1과 같이 추정치를 계산하기 위해서는 밀도함수인 f(x)를 알아야 합니다. 정규분포와 같이 기존의 분포함수를 가정하여 사용할 수 있습니다. 이러한 분포에 적용하기 위해서는 평균과 분산과 같이 모수를 가정할 수 있어야 합니다. 이러한 분석을 모수적 분석(parametic analysis)이라합니다. 반면에 모수를 가정할 수 없는 경우를 비모수 분석(nonparametic analysis)라고 하며 이 경우 커널 밀도 함수를 적용합니다.

특정한 밀도 함수의 가정없이 데이터의 분포에서 특정한 관찰치의 밀도를 계산하는 고전적이지만 자주 사용하는 방법으로 히스토그램이 있습니다.

Histogram

빈도 히스토그램 형태로 데이터를 그룹화하는 것은 다양한 추정 절차의 기초에 내재되어 있는 고전적 방법론입니다. 유용한 시각적 정보를 제공하여 데이터 표현 장치로 사용되었지만 밀도 추정 방법으로서 비모수 통계에서 근본적인 역할을 했습니다.

기본적으로 히스토그램은 빈 높이로 정의되는 계단 함수로, 빈 높이란 각 빈에 포함된 관측치의 비율을 빈 너비로 나눈 값과 같습니다.

분포 함수 $F_x$로부터 랜덤변수 $X_1,\cdots ,X_n$를 관찰하고 분포함수는 연속이라고 가정합니다. $x_1,\cdots ,x_n$을 랜덤변수의 실현에서 관측된 데이터 포인트라고 가정합니다. 시작 포인트 $x_0$에서 구간(bin) $I_j=[x_0+jh,x_0+(j+1)h]$로 정의합니다. 이 정의에서 $j=1,\cdots ,k$입니다.

$$\begin{array}{rl}\text{(식 2)}& P(X=I_j)& ={\int }_{I_j}f(x)dx=f(u)h& u\in I_j\end{array}$$

식 2는 연속적인 유한 함수에 대한 평균값 정리에서 따릅니다. 직관적으로, X가 구간 $I_j$에 속할 확률을 $I_j$의 관측치 비율로 근사할 수 있습니다. 즉,

$$\begin{array}{}\text{(식 3)}& P(X=I_j)\approx \frac{\mathrm{\#}\{x_i\in I_j\}}{n}\end{array}$$

식 2와 3을 적용하여 확률밀도 f(x)는 식 4로 추정될 수 있습니다.

$$\begin{array}{rl}\text{(식 4)}& \widehat{f_h}(x)& =\frac{\mathrm{\#}\{x_i\in I_j\}}{nh}=\frac{1}{nh}\sum _{i=1}^{n}I(x_i\in I_j)x\in I_{j}I(x_i\in I_j)& =\{\begin{array}{ll}1& \text{if}x\in I_j\\ 0& \text{oterwise}\end{array}\end{array}$$

히스토그램 추정치의 평활도는 평활화 매개변수 h(대역폭, bandwidth)에 의해 제어되며, 이는 모든 비모수 곡선 추정치가 공유하는 특성입니다. 작은 대역폭을 선택하면 들쭉날쭉한 추정치가 생성되는 반면, 더 큰 대역폭은 지나치게 평활화된 히스토그램 추정치를 생성하는 경향이 있습니다(Hardle, 1991 참조).

그림 1은 동일한 무작위 생성 데이터의 두 히스토그램의 예를 보여줍니다. 왼쪽의 히스토그램은 작은 대역폭으로 추정되어 결과적으로 많은 빈이 있는 반면, 오른쪽의 히스토그램은 큰 대역폭으로 계산되어 더 적은 수의 빈이 생성됩니다(h 선택에 대한 경험적 규칙은 Sturges의 규칙입니다: $h=1+{\log }_{2}n$).

그림 1. 히스토그램에서 bin의 영향

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
plt.rcParams['font.family'] ='NanumGothic'
plt.rcParams['axes.unicode_minus'] =False
import seaborn as sns
sns.set_style("darkgrid")

a=np.random.RandomState(3)
x=a.normal(size=500)
fig, axs=plt.subplots(1, 2, sharey=True)
sns.histplot(x=x, bins=50, stat="density", ax=axs[0])
sns.histplot(x=x, bins=10, stat="density", ax=axs[1])
plt.show()

히스토그램을 작성할 때는 대역폭뿐만 아니라 각 빈 에지의 시작점(위 그림의 각 bin의 모서리에 위치한 포인트)도 선택해야 합니다. 이러한 선택은 모양에 대해 생성할 수 있으며, 따라서 다른 추정치가 생성됩니다. 빈 에지 문제는 히스토그램의 단점으로, 커널 밀도 추정치와 같은 다른 추정치에서는 공유되지 않습니다. 또 다른 단점은 히스토그램 추정치가 일반적으로 매끄럽지 않아 노이즈로 인해 관찰되었을 수 있는 융기부를 표시한다는 것입니다.

Kernel Density Estimation

KDE의 기본 개념은 간단합니다. 샘플에서 위치 주변에 발생하는 데이터 포인터가 많을수록 해당위치에서 관찰이 발생할 가능성이 높아집니다. 이는 선형적으로 간격이 있는 포인트들에서 다양한 위치의 모든 관찰의 가중거리를 사용하는 것으로 구현됩니다. 이러한 개념을 식 5와 같이 나타낼 수 있습니다.

$$\begin{array}{rl}{\hat{p}}_n(x)& =\frac{1}{nh}\sum _{i=1}^{n}K\left(\frac{X-X_i}{h}\right)\\ X_i:& \text{고정된 위치}\\ \text{(식 5)}& K(x):& \text{Kernel function}h:& \text{Smoothing bandwidth}\end{array}$$

대표적으로 적용되는 kernel function은 다음과 같습니다.

kernel function의 예

kernel	함수
uniform	$k(x)=\frac{1}{2},|x|\le 1$
Epanechnikov	$k(x)=\frac{3}{4}(1-x^2),|x|\le 1$
Gaussian	$k(x)=\frac{1}{\sqrt{2\pi }}\exp (-\frac{1}{2}{x}^2)$

커널밀도 함수 중 gaussian에 대한 커널함수는 다음과 같습니다. 먼저 단변수에 대해 소개합니다.

$$\begin{array}{rl}\text{(식 6)}& x_i& =\{x_1,x_2,\cdots ,x_N\}\hat{f}:\text{KDE of PDF}\hat{f}(x,h)& =\frac{1}{N}\sum _{i=1}^N{K}_h(x-x_i)\end{array}$$

기본적으로 커널밀도추정 접근에서 각 데이터 포인트에 매끄러운 스케일링된 커널 함수를 중심화한 다음 평균을 계산합니다. 가징 일반적인 커널 중 하나는 가우시안 커널 입니다(식 7).

$$\begin{array}{}\text{(식 7)}& K(u)=\frac{1}{\sqrt{2\pi }}\exp (-\frac{u^2}{2})\end{array}$$

식 8의 K_h는 커널의 스케일된 버전입니다.

$$\begin{array}{}\text{(식 8)}& K_h(u)=\frac{1}{h}K\left(\frac{u}{h}\right)\end{array}$$

커널의 h는 bandwidth라고 하며 이 추정량의 매우 중요한 역할을 합니다.

표준정규분포를 기반으로 하는 가우시안 함수로부터 밀도함수를 직관적으로 생각해보면 0 점에서 정점을 이루는 함수는 식 9로 출발합니다.

$$\begin{array}{}\text{(식 9)}& \exp (-x^2)\end{array}$$

모든 구간 x에 대해 적분이 1이 되어야 하기 때문에 식 9는 식 10과 같이 수정됩니다.

$$\begin{array}{}\text{(식 10)}& K(x)=\frac{1}{\sqrt{2\pi }}\exp (-\frac{x^2}{2})\end{array}$$

이 두함수의 그래프는 그림 2와 같습니다. 결과적으로 식 10은 평균과 분산이 각각 0과 1인 가우시안 분포와 동일한 커널함수입니다.

그림 2. 가우시안 커널함수

def kernelFun(x):
    return np.exp(-x**2/2)/np.sqrt(2*np.pi)

x=np.linspace(-5, 5, 1000)
k=np.exp(-x**2)
k1=kernelFun(x)
sns.lineplot(x=x, y=k, label=r"$\exp(-x^2)$")
sns.lineplot(x=x, y=k1, label=r"$\frac{1}{\sqrt{2\pi}}\exp(-\frac{x^2}{2})$")
plt.legend()
plt.plot()

단일 데이터 포인트에서 x에 데이터 세트 $X=\{x_1,x_2,\cdots ,x_n\}$에서 i번째 포인트를 빼서 이동한 값을 커널에 적용할 수 있습니다. 또한 이 차이를 더 넓게 또는 더 좁게하기 위해 상수 h로 나누어 줍니다. 이 h값을 커널대역폭(kernel bandwidth)라고 합니다. 이 함수 역시 단위면적이 되어야 하므로 식 11과 같이 h로 나누어줍니다.

$$\begin{array}{}\text{(식 11)}& \frac{1}{h}K\left(\frac{x-x_i}{h}\right)\end{array}$$

다음 sympy를 사용하여 식 10과 첫번째 포인트 $x_1$에 대해 식 11을 적분한 결과를 나타낸 것입니다.

from sympy  import *

a=symbols("a")
f=exp(-a**2/2)/sqrt(2*pi)
F=f.integrate((a,-5, 5))
F.evalf(3)

1.0

h=0.3
x1=1
f1=exp(-((a-x1)/h)**2/2)/(h*sqrt(2*pi))
F1=f1.integrate((a,-5, 5))
F1.evalf(3)

1.0

첫번째 포인트 $x_1$에 대한 커널밀도함수는 식 12와 같습니다.

$$\begin{array}{}\text{(식 12)}& \frac{1}{h}K\left(\frac{x-x_1}{h}\right)\end{array}$$

식 12의 확률(pdf)는 그림 3과 같습니다.

그림 3. 지정 포인트 x₁=1인 가우시안 커널함수의 pdf.

h=0.3
x1=1
k2=kernelFun((x-x1)/h)/h
sns.lineplot(x=x, y=k2, label=r"$\frac{1}{h}k(\frac{x-x_1}{h})$")
plt.scatter(x1, 0, color="k", label=r"$x_1$")
plt.legend()
plt.show()

두번째 포인트에 대한 커널밀도함수는 식 5와 같습니다.

$$\begin{array}{}\text{(식 13)}& \frac{1}{h}K\left(\frac{x-x_2}{h}\right)\end{array}$$

식 12와 13의 두 점에 대한 단일 PDF는 각 PDF의 합이 됩니다(식 14).

$$\begin{array}{rl}f(x)& =\frac{1}{h}K\left(\frac{x-x_1}{h}\right)+\frac{1}{h}K\left(\frac{x-x_2}{h}\right)\\ \text{(식 14)}& & =\frac{1}{2h}[K\left(\frac{x-x_1}{h}\right)+K\left(\frac{x-x_2}{h}\right)]& =\frac{1}{2h}\sum _{i=1}^{2}K\left(\frac{x-x_i}{h}\right)\end{array}$$

그림 4. 지정 포인트 x₁와 x₂에 대한 가우시안 커널함수의 pdf.

식 14를 일반화하면 식 15와 같습니다.

$$\begin{array}{}\text{(식 15)}& f(x)=\frac{1}{nh}\sum _{i=1}^{n}K\left(\frac{x-x_i}{h}\right)\end{array}$$

python을 적용한 kde의 계산

scipy.stats.gaussian_kde(1차원객체)와 sklearn.neighbors.KernelDensity(2차원객체) 클래스를 사용하여 변수의 kde를 계산할 수 있습니다. 클래스 KernelDensity()의 경우

from scipy.stats import gaussian_kde
from sklearn.neighbors import KernelDensity

rng = np.random.RandomState(42)
X = rng.random_sample((100, 1))
X.shape

(100, 1)

X1=np.ravel(X)
kernel=gaussian_kde(X1)
y=kernel.pdf(X1)
y[:10]

array([0.94851699, 0.60328742, 0.93610674, 0.91989456, 1.05686655,
       1.05684708, 0.85001855, 0.81104559, 0.92045114, 0.93784738])

kde=KernelDensity(kernel='gaussian', bandwidth=0.1).fit(X)
logD=kde.score_samples(X)
logD[:10]

array([-0.06407906, -0.47322023, -0.04934331, -0.08301556,  0.10383486,
        0.10382813, -0.10005549, -0.1694874 , -0.08246769, -0.05395054])

sns.set_style("darkgrid")
fig, (ax1, ax2)=plt.subplots(1,2, figsize=(7, 3), gridspec_kw=dict(width_ratios=[3,3]), sharey=True)
sns.lineplot(x=X1, y=y, color="b", ax=ax1)
ax1.set_xlabel("X")
ax1.set_ylabel("kde")
ax1.set_title("gaussian_kde()")
sns.lineplot(x=X[:,0], y=np.exp(logD), color="r", ax=ax2)
ax2.set_xlabel("X")
ax2.set_title("KernelDensity()")
plt.show()

그림 5. gaussian_kde와 KernelDensity 클래스에 의한 kde.