kernel density estimation是在概率论中用来估计未知的密度函数,属于非参数检验方法之一,由Rosenblatt (1955)和Emanuel Parzen(1962)提出,又名Parzen窗(Parzen window)
本文翻译自英国雷丁大学(Reading University)Xia Hong老师的讲义材料
1 概率密度函数
连续概率函数 p ( x ) p(x) p(x)的数学定义满足以下特性:
- x x x介于 a a a、 b b b两点之间的概率为 P ( a < x < b ) = ∫ a b p ( x ) d x P(a<x<b)=\int_a^bp(x)dx P(a<x<b)=∫abp(x)dx
- x是非负实数
- 概率函数的积分: ∫ − ∞ ∞ p ( x ) d x = 1 \int_{-\infty}^\infty p(x)dx=1 ∫−∞∞p(x)dx=1
最常见的概率函数是高斯函数(Gaussian Function,又称为正态分布): p ( x ) = 1 2 π σ exp ( − ( x − c ) 2 2 σ 2 ) p(x)=\frac{1}{\sqrt{2\pi\sigma}}\exp(-\frac{(x-c)^2}{2\sigma^2}) p(x)=2πσ1exp(−2σ2(x−c)2),这里 c c c是均值, σ \sigma σ是标准差
拓展到对于向量 x \mathbf{x} x,非负函数 p ( x ) p(\mathbf{x}) p(x)有以下特性:
- x \mathbf{x} x在区域 R R R里的概率为 P = ∫ R p ( x ) d x P=\int_R p(\mathbf{x})d\mathbf{x} P=∫Rp(x)dx
- 概率函数的积分为 ∫ p ( x ) d x = 1 \int p(\mathbf{x})d\mathbf{x}=1 ∫p(x)dx=1
2 密度估计
密度估计:给定的一系列数量为
n
n
n的样本
x
1
,
⋅
⋅
⋅
,
x
n
\mathbf{x}_1, \cdot\cdot\cdot,\mathbf{x}_n
x1,⋅⋅⋅,xn,可以估计密度函数
p
(
x
)
p(\mathbf{x})
p(x),从而根据任意新样本
x
\mathbf{x}
x可以得到输出
p
(
x
)
p(\mathbf{x})
p(x)
大部分未知密度函数估计方法的基本思想都很简单,主要是依赖于样本落在区域
R
R
R的概率
P
P
P,即有
P
=
∫
R
p
(
x
)
d
x
P=\int_R p(\mathbf{x})d\mathbf{x}
P=∫Rp(x)dx
假设区域
R
R
R很小,
P
(
x
)
P(\mathbf{x})
P(x)在区域内波动很小,上式可以写做
P
=
∫
R
p
(
x
)
d
x
≈
p
(
x
)
∫
R
d
x
=
p
(
x
)
V
P=\int_R p(\mathbf{x})d\mathbf{x}\approx p(\mathbf{x})\int_R d\mathbf{x}=p(\mathbf{x})V
P=∫Rp(x)dx≈p(x)∫Rdx=p(x)V,这里
V
V
V是区域
R
R
R的“量”(二维即为面积)
从另一方面看,假设
n
n
n个样本
x
1
,
⋅
⋅
⋅
,
x
n
\mathbf{x}_1, \cdot\cdot\cdot,\mathbf{x}_n
x1,⋅⋅⋅,xn都是独立且服从概率密度函数
p
(
x
)
p(\mathbf{x})
p(x),且
n
n
n个样本中有
k
k
k个落在区域
R
R
R里面,则有
P
=
k
/
n
P=k/n
P=k/n,因此
p
(
x
)
p(\mathbf{x})
p(x)的估计式为
p
(
x
)
=
k
/
n
V
p(\mathbf{x})=\frac{k/n}{V}
p(x)=Vk/n
3 Parzen窗密度估计
考虑 R R R是中心在 x \mathbf{x} x的超立方体(例如二维平面),令 h h h为超立方体的边缘长度,所以对于二维平面有有 V = h 2 V=h^2 V=h2,对于三维立体有 V = h 3 V=h^3 V=h3
引入
ϕ
(
x
i
−
x
h
)
=
{
1
∣
x
i
k
−
x
k
∣
h
<
=
1
/
2
,
k
=
1
,
2
0
o
t
h
e
r
w
i
s
e
\phi(\frac{\mathbf{x}_i-\mathbf{x}}{h})=\left\{ \begin{aligned} 1\quad& \frac{|x_{ik}-x_{k}|}{h}<=1/2, k=1,2 \\ 0\quad& otherwise \end{aligned} \right.
ϕ(hxi−x)=⎩⎨⎧10h∣xik−xk∣<=1/2,k=1,2otherwise
Parzen概率密度公式(二维)为
p
(
x
)
=
k
/
n
V
=
1
n
∑
i
=
1
n
1
h
2
ϕ
(
x
i
−
x
h
)
p(\mathbf{x})=\frac{k/n}{V}=\frac{1}{n}\sum_{i=1}^n {\frac{1}{h^2}\phi(\frac{\mathbf{x}_i-\mathbf{x}}{h})}
p(x)=Vk/n=n1∑i=1nh21ϕ(hxi−x),
ϕ
(
x
i
−
x
h
)
\phi(\frac{\mathbf{x}_i-\mathbf{x}}{h})
ϕ(hxi−x)即为窗函数
我们归纳这个思想并拓展到其他Parzen窗密度估计法中
例如,如果使用高斯函数,对于一维有:
p
(
x
)
=
1
n
∑
i
=
1
n
1
2
π
σ
exp
(
−
(
x
i
−
x
)
2
2
σ
2
)
p(x)=\frac{1}{n}\sum_{i=1}^{n}{\frac{1}{\sqrt{2\pi\sigma}}\exp(-\frac{(x_i-x)^2}{2\sigma^2})}
p(x)=n1∑i=1n2πσ1exp(−2σ2(xi−x)2),这是对
n
n
n个将数据点作为中心的高斯函数的简单求平均,公式中的
σ
\sigma
σ需要再做确定
例子:
给定一个系列的5个数据点
x
1
=
2
x_1=2
x1=2,
x
2
=
2.5
x_2=2.5
x2=2.5,
x
3
=
3
x_3=3
x3=3,
x
4
=
1
x_4=1
x4=1,
x
5
=
6
x_5=6
x5=6,参数
σ
=
1
\sigma=1
σ=1,中心
x
=
3
x=3
x=3的高斯函数作为窗函数,求出Parzen概率密度估计(pdf)
解答:
1
2
π
exp
(
−
(
x
1
−
x
)
2
2
)
=
1
2
π
exp
(
−
(
2
−
3
)
2
2
)
=
0.2420
\frac{1}{\sqrt{2\pi}}\exp(-\frac{(x_1-x)^2}{2})=\frac{1}{\sqrt{2\pi}}\exp(-\frac{(2-3)^2}{2})=0.2420
2π1exp(−2(x1−x)2)=2π1exp(−2(2−3)2)=0.2420
1
2
π
exp
(
−
(
x
2
−
x
)
2
2
)
=
1
2
π
exp
(
−
(
2.5
−
3
)
2
2
)
=
0.3521
\frac{1}{\sqrt{2\pi}}\exp(-\frac{(x_2-x)^2}{2})=\frac{1}{\sqrt{2\pi}}\exp(-\frac{(2.5-3)^2}{2})=0.3521
2π1exp(−2(x2−x)2)=2π1exp(−2(2.5−3)2)=0.3521
1
2
π
exp
(
−
(
x
3
−
x
)
2
2
)
=
1
2
π
exp
(
−
(
3
−
3
)
2
2
)
=
0.3989
\frac{1}{\sqrt{2\pi}}\exp(-\frac{(x_3-x)^2}{2})=\frac{1}{\sqrt{2\pi}}\exp(-\frac{(3-3)^2}{2})=0.3989
2π1exp(−2(x3−x)2)=2π1exp(−2(3−3)2)=0.3989
1
2
π
exp
(
−
(
x
4
−
x
)
2
2
)
=
1
2
π
exp
(
−
(
1
−
3
)
2
2
)
=
0.0540
\frac{1}{\sqrt{2\pi}}\exp(-\frac{(x_4-x)^2}{2})=\frac{1}{\sqrt{2\pi}}\exp(-\frac{(1-3)^2}{2})=0.0540
2π1exp(−2(x4−x)2)=2π1exp(−2(1−3)2)=0.0540
1
2
π
exp
(
−
(
x
5
−
x
)
2
2
)
=
1
2
π
exp
(
−
(
6
−
3
)
2
2
)
=
0.0044
\frac{1}{\sqrt{2\pi}}\exp(-\frac{(x_5-x)^2}{2})=\frac{1}{\sqrt{2\pi}}\exp(-\frac{(6-3)^2}{2})=0.0044
2π1exp(−2(x5−x)2)=2π1exp(−2(6−3)2)=0.0044
因此,
p
(
x
=
3
)
=
(
0.2420
+
0.3521
+
0.3989
+
0.0540
+
0.0044
)
/
5
=
0.2103
p(x=3)=(0.2420 + 0.3521 + 0.3989+0.0540 + 0.0044)/5 = 0.2103
p(x=3)=(0.2420+0.3521+0.3989+0.0540+0.0044)/5=0.2103
下面用图形化语言表示Parzen窗,每个数据点密度函数(虚线)对于最终的概率密度函数(实线)有相同的贡献度