参考文献：
Density Estimation for Statistics and Data Analysis. B.W.Silverman
峰度定义：https://en.wikipedia.org/wiki/Kurtosis
偏度定义：https://en.wikipedia.org/wiki/Skewness
halvorsen 的博客：Numerical Integration in Python
文章目录
在博客核密度估计：带宽宽度选择——原理阐述中，我们讨论了最佳带宽 h h h 的选择原理，即使得 I M S E ( f ^ ) IMSE(\hat{f}) IMSE(f^) 最小， f ^ \hat{f} f^ 为经验概率密度函数，是从随机向量（样本） X = ( X 1 , X 2 , ⋯ , X n ) X=(X_1,X_2,\cdots,X_n) X=(X1,X2,⋯,Xn) 根据样本在各带宽下的频率估计而出的。
我们知道，在忽略
o
(
h
4
+
(
n
h
)
−
1
)
o(h^4+(nh)^{-1})
o(h4+(nh)−1) 的前提下，我们得到最佳宽度的表达式为：
h
o
p
t
=
κ
2
−
2
/
5
κ
1
/
5
{
∫
[
f
(
2
)
(
x
)
]
2
d
x
}
−
1
/
5
n
−
1
/
5
h_{opt} = \kappa_2^{-2/5} \kappa^{1/5} \{\int [f^{(2)}(x)]^2 dx \}^{-1/5} n^{-1/5}
hopt=κ2−2/5κ1/5{∫[f(2)(x)]2dx}−1/5n−1/5
其中：
κ
2
=
∫
k
(
v
)
v
2
d
v
κ
=
∫
k
2
(
v
)
d
v
\begin{aligned} \kappa_2 =& \int k(v) v^2 dv \\ \kappa = & \int k^2(v) dv \end{aligned}
κ2=κ=∫k(v)v2dv∫k2(v)dv
上述公式中，
k
(
v
)
k(v)
k(v) 为核函数，其常数表达式可为正态分布，在式中是已知的；
f
(
x
)
f(x)
f(x) 为待估计的概率密度函数，
f
(
2
)
(
x
)
f^{(2)}(x)
f(2)(x) 为其二阶导数；
n
n
n 为样本容量。
可以看出， h o p t h_{opt} hopt 与 f ( 2 ) ( x ) f^{(2)}(x) f(2)(x) 有关。
Rule of Thumb 方法
rule of thumb 方法也叫经验法，从上述分析中我们可以看出，由于
h
o
p
t
h_{opt}
hopt 与
f
(
2
)
(
x
)
f^{(2)}(x)
f(2)(x) 有关，而
f
(
x
)
f(x)
f(x) 是未知待求解的。但我们可以在求解窗口宽度
h
o
p
t
h_{opt}
hopt 时，将
f
(
x
)
f(x)
f(x) 假设为一个已知函数，通常为方差为
σ
2
\sigma^2
σ2 的正态分布，于是可以求解出：
∫
f
(
2
)
(
x
)
2
d
x
=
σ
−
5
∫
ϕ
(
2
)
(
x
)
2
d
x
=
3
8
π
−
1
/
2
σ
−
5
\int f^{(2)}(x)^2 dx = \sigma^{-5} \int \phi^{(2)}(x)^2 dx = \frac{3}{8}\pi^{-1/2}\sigma^{-5}
∫f(2)(x)2dx=σ−5∫ϕ(2)(x)2dx=83π−1/2σ−5
Φ
(
x
)
\Phi(x)
Φ(x) 为正态分布，
f
(
x
)
f(x)
f(x) 的均值可以通过变量代换消掉 ，最终可以估计出
∫
f
(
2
)
(
x
)
2
d
x
≈
0.212
σ
−
5
\int f^{(2)}(x)^2 dx \approx 0.212\sigma^{-5}
∫f(2)(x)2dx≈0.212σ−5。
假设我们采用了标准正态分布作为核函数，那么就可以求解出最佳的窗口宽度
h
o
p
t
h_{opt}
hopt 如下所示：
h
o
p
t
=
(
4
3
)
1
/
5
σ
n
−
1
/
5
=
1.06
σ
n
−
1
/
5
\begin{aligned} h_{opt}& = (\frac{4}{3})^{1/5} \sigma n^{-1/5} \\ &= 1.06 \sigma n^{-1/5} \end{aligned}
hopt=(34)1/5σn−1/5=1.06σn−1/5
如果随机变量确实服从正态分布，那么使用此种方法得出的最佳带宽是比较准确的。但如果数据并非服从正态分布，而是服从某种多峰分布，得使用这种方法，会造成核密度估计的过分平滑，实际分布的细节会被 cover 掉。
我们也可以从 h o p t h_{opt} hopt 的定义中看出这一点： f ( 2 ) ( x ) f^{(2)}(x) f(2)(x) 作为函数的二阶导数，它反映了函数的凹凸程度。在相同的方差下，正态分布的上凸或下凹的程度是明显大于多峰分布的1。因此，若假设 f ( x ) f(x) f(x) 是正态分布来求解 h o p t h_{opt} hopt，会由于 f ( 2 ) ( x ) f^{(2)}(x) f(2)(x) 较大，从而使得 h h h 过大，造成过平滑。
另一些变种的方法如下：
h
o
p
t
=
0.79
×
IQR
n
−
1
/
5
h_{opt} = 0.79 \times \text{IQR} n^{-1/5}
hopt=0.79×IQRn−1/5
和
h
o
p
t
=
0.9
min
(
σ
,
nIQR
)
n
−
1
/
5
h_{opt} = 0.9 \min(\sigma, \text{nIQR}) n^{-1/5}
hopt=0.9min(σ,nIQR)n−1/5
其中 IQR 为四分位距，nIQR 为标准四分位距，其取值大致为
n
I
Q
R
=
0.7413
×
I
Q
R
nIQR = 0.7413 \times IQR
nIQR=0.7413×IQR。
双峰分布情况下，ROT 方法的有效性讨论
我们用两个正态分布来构造一个双峰分布，其表达式如下所示：
y
=
0.5
×
(
1
2
π
exp
(
−
(
x
−
μ
1
)
2
2
)
+
1
2
π
exp
(
−
(
x
−
μ
2
)
2
2
)
)
y = 0.5\times \bigg( \frac{1}{\sqrt{2\pi}}\exp \big( \frac{-(x-\mu_1)^2}{2} \big) + \frac{1}{\sqrt{2\pi}}\exp \big( \frac{-(x-\mu_2)^2}{2} \big)\bigg)
y=0.5×(2π
1exp(2−(x−μ1)2)+2π
1exp(2−(x−μ2)2))
当
μ
1
=
−
2
,
μ
2
=
2
\mu_1=-2, \mu_2=2
μ1=−2,μ2=2 时，可以得出双峰分布的概率密度图像如下所示：
我们可以根据带宽宽度最优解公式
h
o
p
t
=
κ
2
−
2
/
5
κ
1
/
5
{
∫
[
f
(
2
)
(
x
)
]
2
d
x
}
−
1
/
5
n
−
1
/
5
h_{opt} = \kappa_2^{-2/5} \kappa^{1/5} \{\int [f^{(2)}(x)]^2 dx \}^{-1/5} n^{-1/5}
hopt=κ2−2/5κ1/5{∫[f(2)(x)]2dx}−1/5n−1/5，求出双峰分布的最佳带宽宽度
h
o
p
t
h_{opt}
hopt，并于使用经验公式
h
r
o
t
=
1.06
σ
n
−
1
/
5
h_{rot}= 1.06 \sigma n^{-1/5}
hrot=1.06σn−1/5 求出经验最佳带宽宽度，并求出两者的比（可以约掉公因数）：
ratio
=
h
o
p
t
h
r
o
t
=
(
∫
f
(
2
)
(
x
)
2
d
x
3
8
π
−
1
/
2
σ
−
5
)
−
1
/
5
\begin{aligned} \text{ratio} &= \frac{h_{opt}}{h_{rot}} \\ &=\bigg(\frac{\int f^{(2)}(x)^2 dx}{\frac{3}{8} \pi^{-1/2} \sigma^{-5}} \bigg)^{-1/5} \end{aligned}
ratio=hrothopt=(83π−1/2σ−5∫f(2)(x)2dx)−1/5
于是，我们遍历
∣
μ
1
−
μ
2
∣
∈
[
0
,
6
]
| \mu_1 - \mu_2 |\in[0, 6]
∣μ1−μ2∣∈[0,6] 的不同双峰分布，并求出
∣
μ
1
−
μ
2
∣
| \mu_1 - \mu_2 |
∣μ1−μ2∣ 的峰峰间距对应的
h
o
p
t
/
h
r
o
t
h_{opt}/h_{rot}
hopt/hrot：
可以看到在
∣
μ
1
−
μ
2
∣
∈
[
0
,
2
]
| \mu_1 - \mu_2 |\in[0,2]
∣μ1−μ2∣∈[0,2] 时，用经验公式求出的带宽宽度与最佳带宽宽度是比较相近的。这实际上是因为峰峰间距为
[
0
,
2
]
[0,2]
[0,2] 时，双峰分布的峰峰间距实际上并不明显，其概率密度函数的图像与正态分布是非常相近的，我们可以画出峰，峰间距为 2 时，双峰分布的概率密度图像如下所示：
可以看到与正态分布是非常相近的，此时双峰分布的双峰特点并不明显。因此当使用经验公式求取最佳宽度时，与实际的最佳宽度比较接近。当双峰分布的间距大于 2 时，很明显的，实际最佳宽度
h
o
p
t
h_{opt}
hopt 与使用经验公式求取的最佳宽度
h
r
o
t
h_{rot}
hrot 的差距变大，
h
r
o
t
−
h
o
p
t
h_{rot} - h_{opt}
hrot−hopt 的值，随着峰峰间距的增大而增大，此时使用核密度图，会造成过分平滑。
偏峰分布情况下，ROT 方法的有效性讨论
同样的我们考虑经验法求解的最佳宽度与实际最佳宽度的偏差，在偏锋分布的情况下又是怎么样的呢？我们用 ln-正态分布作为偏峰分布，其表达式如下所示：
y
=
1
σ
x
2
π
exp
(
−
ln
2
(
x
)
2
σ
2
)
y = \frac{1}{\sigma x \sqrt{2\pi}} \exp\Big( -\frac{\text{ln}^2(x)}{2\sigma^2} \Big)
y=σx2π
1exp(−2σ2ln2(x))
该函数在不同的标准差下，其图像如下所示：
根据偏度的求解公式，我们可以得出，标准差
σ
∈
(
0
,
0.7
]
\sigma\in(0,0.7]
σ∈(0,0.7] 时，偏度系数的取值范围为
[
0
,
3
)
[0, 3)
[0,3)。我们画出偏度系数的取值为
[
0
,
3
)
[0, 3)
[0,3) 的情况下，对应的
h
o
p
t
/
h
r
o
t
h_{opt}/h_{rot}
hopt/hrot 都取值，并绘制偏度系数与该比值构成的图像，如下所示：
可以看到当偏度系数增大时，用经验公式求取出来的最佳宽度与实际最佳宽度的偏差逐渐增大，导致使用核函数估计得出的图像过分平滑。
峰度系数不同情况下，ROT 方法的有效性讨论
t 分布与正态分布十分相似，但是 t 分布可以调整分布的峰值系数，而正态分布的峰度系数为 0。经验法求取最佳宽度时，是用正态分布去近似实际概率密度函数的，所以我们可以用峰度系数不同的 t 分布，来观察峰度对 ROT 方法有效性的影响。
当 t 分布的均值为 0，方差为 1 时，其概率密度函数的表达式如下所示：
y
=
Γ
(
(
n
+
1
)
/
2
)
π
n
Γ
(
n
/
2
)
(
1
+
x
2
/
n
)
−
(
n
+
1
)
/
2
y = \frac{\Gamma\Big((n+1)/2\Big)}{\sqrt{\pi n}\Gamma(n/2)}\Big( 1+x^2/n\Big)^{-(n+1)/2}
y=πn
Γ(n/2)Γ((n+1)/2)(1+x2/n)−(n+1)/2
不同峰度系数下的 t 分布与正态分布的概率密度图像如下所示：
可以看到峰度系数越大，则该函数的概率密度图的尾巴越“厚实”。当峰度为零时，此时*度非常的高，t 分布与正态分布十分接近。
我们画出偏度系数的取值为
[
0
,
120
]
[0, 120]
[0,120] 的情况下，对应的
h
o
p
t
/
h
r
o
t
h_{opt}/h_{rot}
hopt/hrot 都取值，并绘制偏度系数与该比值构成的图像，如下所示：
从上图可以看出，对于峰度系数不同的概率密度函数来说，采用经验法得出的最佳宽度与实际的最佳宽度，非常接近。
结束语
非常感谢各位能够静下心来读完本博客，也希望各位在阅读本博客时，能够边读边思考，有所收获。
若各位希望得到本博客中，分析和画图用到的源程序，请关注、点赞并私信我索取。
若各位有优秀的论文，并认为具有极高的复现价值，也可以推荐给博主。
再次感谢各位的耐心观看，希望本博客能对您的知识有所益助。
-
如下图所示：
图中的单峰分布为标准正态分布，图中的双峰分布为分别用正态分布 0.5 × N ( − 1 , 0.01 ) + 0.5 × N ( 1 , 0.01 ) 0.5\times N(-1,0.01)+0.5\times N(1,0.01) 0.5×N(−1,0.01)+0.5×N(1,0.01) 得出的概率密度函数。图中的概率密度是根据 [ − 10 , 10 ] [-10,10] [−10,10] ，间隔为 0.1 的采样数据，以及两个分布的概率密度函数画出来的。根据方差公式的计算分别得出图中正态函数和双峰分布的方差分别为：1.00 和 1.01，可以近似的认为两个分布的方差相同。从下图可以看出，很明显的正态分布的凸度，较双峰分布大。↩︎