分布的估计

有了服从某一分布的若干样本之后,如果我们能推出这个分布,那这对我们的接下来的所有操作都是十分有益的。因此从样本估计分布是一个十分重要的话题。

如果读者对统计有一些基本的认识,应该不难想到,如果我们的样本足够多,那么在绘制直方图的时候每个柱形都可以足够细(在保持每个柱形都有一定的样本量的前提下),此时我们的直方图就是其概率密度分布的一个良好的近似。一个更加数学化的版本就是,我们使用经验累计概率密度函数$F_n(x)\overset{\Delta}{=}\frac{1}{n}\sum_{i=1}^{n}\boldsymbol{1}(X_i\leq x)$来近似原分布的累计概率密度函数$F(x)\overset{\Delta}{=}\mathbb{P}(X\leq x)$。而强大数定律保证了,当$n$充分大时,$F_n(x)\overset{\mathrm{a.s.}}{\rightarrow} F(x),\,\forall x\in\mathbb{R}$。因此,我们可以使用

$$\widehat{p}_n(x)\overset{\Delta}{=}\frac{F_n(x+h)-F_n(x-h)}{2h}$$

来近似我们的目标

$p(x)\approx\frac{F(x+h)-F(x-h)}{2h}$

对于$\widehat{p}_n(x)$,我们对其稍作整理

$$\widehat{p}_n(x)=\frac{1}{2nh}\sum_{i=1}^{n}\boldsymbol{1}(x-h<X_i\leq x+h)\overset{\Delta}{=}\frac{1}{n}\sum_{i=1}^{n}\frac{1}{h}K\left(\frac{X_i-x}{h}\right)$$

其中$K(u)\overset{\Delta}{=}\frac{1}{2}\boldsymbol{1}(-1<u\leq 1)$。事实上这就是一个核函数。

估计的误差

让我们仔细地考察我们$\widehat{p}_n(x)$的计算中$h$的选取。显然,要让估计的准度足够高,我们需要选取更小的$h$。但过小的$h$同样会带来问题——过小的$h$可能意味着只有极少数的点被统计到,这也就导致我们估计的方差会增大。这只是我们定性的分析,让我们定量的写下这一过程。

对于某一点$x_0$处的密度估计,我们定义

$$\mathrm{MSE}(x_0)\overset{\Delta}{=}\mathbb{E}_p[(\widehat{p}_n(x_0)-p(x_0))^2]$$

通过经典的代数变形,我们可以将其分解为

$$\mathrm{MSE}(x_0)=\underbrace{(\mathbb{E}_p[\widehat{p}_{n}(x_0)]-p(x_0))^2}_{\mathrm{bias}(x_0)^2}+\underbrace{\mathbb{E}_p[(\widehat{p}_{n}(x_0)-\mathbb{E}_p[\widehat{p}_{n}(x_0)])^2]}_{\mathrm{Var}(\widehat{p}_n(x_0))}$$

方差项

设概率密度函数$p(x)\leq p_\max<\infty$,且$K:\mathbb{R}\rightarrow\mathbb{R}$满足$\int K^2(u)\,\mathrm{d}u<\infty$。那么我们有

$$\mathrm{Var}(\widehat{p}_n(x_0))\leq\frac{p_\max\int K^2(u)\,\mathrm{d}u}{nh}=C_1\frac{1}{nh}$$

为了证明这一式子,首先考察

$$\widehat{p}_n(x_0)-\mathbb{E}_p[\widehat{p}_n(x_0)]=\frac{1}{nh}\sum_{i=1}^{n}\{K(\frac{X_i-x}{h})-\mathbb{E}_p[K(\frac{X_i-x}{h})]\}\overset{\Delta}{=}\frac{1}{nh}\sum_{i=1}^{n}\eta_i(x_0)$$

不难发现,$\eta_i(x_0),i=1,2,\ldots,n$独立同分布,且均值为$0$,方差为

$$\begin{align}\mathrm{Var}(\eta_i(x_0))&=\mathbb{E}_p[\eta_i^2(x_0)]\\&\leq\mathbb{E}_p\left[K^2\left(\frac{X_i-x_0}{h}\right)\right]\\&=\int K^2\left(\frac{z-x_0}{h}\right)p(z)\,\mathrm{d}z\\&=\int K^2(u)p(x_0+uh)h\,\mathrm{d} u\\&\leq p_\max h\int K^2(u)\,\mathrm{d}u\end{align}$$

因此

$$\mathrm{Var}(\widehat{p}_n(x_0))=\mathrm{Var}(\widehat{p}_n(x_0)-\mathbb{E}_p[\widehat{p}_n(x_0)])=\frac{1}{n^2h^2}\sum_{i=1}^{n}\mathrm{Var}(\eta_i(x_0))\leq\frac{p_\max\int K^2(u)\,\mathrm{d}u}{nh}$$

偏差项

在讨论这一项前我们先定义区间$T$上的Holder class $\Sigma(\beta,L),\,\beta,L\in\mathbb{R}^{+}$。其是一系列函数的集合,这些函数满足两个性质

  • 该函数$\lfloor\beta\rfloor$阶可导
  • $\lvert f^{\lfloor\beta\rfloor}(x)-f^{\lfloor\beta\rfloor}(x^\prime) \rvert\leq L\lvert x-x^\prime \rvert^{\beta-\lfloor\beta\rfloor},\,\forall x,x^\prime\in T $

如果此时再加上约束$\int f(x)\,\mathrm{d}x=1, f\ge 0$, 我们将其记为$\mathcal{P}(\beta,L)$

同时定义$l$阶核函数,其满足性质

  • $\int K(u)\,\mathrm{d}u$=1
  • $\forall j=1,2,\ldots,l,\,\int u^jK(u)\,\mathrm{d}u=0$
  • $\int u^{l+1}K(u)\mathbb{d}u\neq 0$

我们在本文的结尾构造了一个满足定义的核函数,这说明了其存在性。

在此基础上,我们开始我们的推导。假设$p\in\mathcal{P}(\beta,L)$,并令$K$为$l=\lfloor\beta\rfloor$阶核且满足$\int\lvert u\rvert^\beta\lvert K(u)\rvert\,\mathrm{d}u<\infty$。此时则有

$$\lvert\mathbb{E}_p[\widehat{p}_n(x_0)]-p(x_0)\rvert\leq \frac{L}{l!}\int\lvert u\rvert^\beta\lvert K(u)\rvert\,\mathrm{d}u\cdot h^\beta=C_2h^\beta$$

证明的核心步骤如下,其中第四个等式使用了泰勒展开的拉格朗日余项

$$\begin{align}\mathbb{E}_p[\widehat{p}_n(x_0)]-p(x_0)&=\frac{1}{h}\int K\left(\frac{x-x_0}{h}\right)p(x)\,\mathrm{d}x-p(x_0)\\&=\int K(u)p(x_0+uh)\,\mathrm{d}u-p(x_0)\\&=\int K(u)[p(x_0+uh)-p(x_0)]\,\mathrm{d}u\\&=\int K(u)(uh)^l/l! [p^{(l)}(x_0+\tau u h)]\,\mathrm{d}u\\&=\int K(u)(uh)^l/l! [p^{(l)}(x_0+\tau u h)-p^{(l)}(x_0)]\,\mathrm{d}u\\&\leq\int K(u) (uh)^l/l! \cdot L\lvert u h\rvert^{\beta-\lfloor{\beta\rfloor}}\,\mathrm{d}u\\&=\frac{Lh^\beta}{l!}\int u^\beta K(u) \,\mathrm{d}u\leq C_2h^\beta \end{align}$$

综上,我们有

$$\mathrm{MSE}\leq C_2^2h^{2\beta}+\frac{C_1}{nh}$$

最小化该式子,有$h_n^\star=\left(\frac{C_1}{2\beta C_2^2}\right)^{1/(2\beta+1)}n^{-1/(2\beta+1)}$,此时

$$\mathrm{MSE}=O(n^{-2\beta/(2\beta+1)})$$

这给出了随着样本量$n$的增加,$\mathrm{MSE}$的渐进趋势。

核函数的构造

之前我们只给出了核函数的定义,那么这样的函数是否存在呢?这里我们给出了一个构造。

取$\{\varphi_m(\cdot)\}_{m=0}^{\infty}$是区间$[-1,1]$上的勒让德多项式,即

$$\varphi_0(x)=1/\sqrt{2},\,\varphi_m(x)=\sqrt{\frac{2m+1}{2}}\frac{1}{2^m m!}\frac{\mathrm{d}^m}{\mathrm{d}x^m}[(x^2-1)^m]$$

可以证明

$$\int_{-1}^{1}\varphi_m(u)\varphi_k(u)\,\mathrm{d}u=\delta_{mk}=\left\{\begin{align}1,\,\mathrm{if}\,m=k\\0,\,\mathrm{if}\,m\neq k\end{align}\right.$$

则$K(u)=\sum_{m=0}^{l}\varphi_m(0)\varphi_m(u)\boldsymbol{1}(\lvert u\rvert\leq 1)$是一个$l$阶核。这是因为对于$u^j,\,j=0,1,\ldots,l$,其总可以表示各正交基的线性组合,即$u^j=\sum_{q=0}^j b_{qj}\varphi_q(u)$。此时

$$\int u^jK(u)\,\mathrm{d}u=\sum_{q=0}^j\sum_{m=0}^lb_{qj}\varphi_q(u)\varphi_m(0)\varphi_m(u)\,\mathrm{d}u=\sum_{q=0}^{j}b_{qj}\varphi_q(0)=u^j\lvert_{(u=0)}=\left\{\begin{align}1,\,\mathrm{if}\,j&=0\\0,\,\mathrm{if}\,j&=1,\ldots,l\end{align}\right.$$

另外,值得注意的是,显然我们的核函数$K$不会是恒正的。这可能会导致一个问题:我们的$\widehat{p}_n$将可能出现负数,因此不再是一个合法的概率密度函数。为此,一个简单的修正就是取$\max\{\widehat{p}_n(x),0\}$。可以证明这样修正后误差不会更差。(当然,这样修正后其实也不是一个合法的概率密度函数,但作为单点处的估计来说是足够好的)