分布的估计

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

如果读者对统计有一些基本的认识,应该不难想到,如果我们的样本足够多,那么在绘制直方图的时候每个柱形都可以足够细(在保持每个柱形都有一定的样本量的前提下),此时我们的直方图就是其概率密度分布的一个良好的近似。一个更加数学化的版本就是,我们使用经验累计概率密度函数\(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\)阶可导

如果此时再加上约束\(\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\}\)。可以证明这样修正后误差不会更差。(当然,这样修正后其实也不是一个合法的概率密度函数,但作为单点处的估计来说是足够好的)