稀疏数据分析：马蹄估计量及其理论性质

• 基本框架
• • half-Cauchy分布
  • 为什么它叫马蹄估计量
  • 后验均值、shrinkage与$\kappa$
  • $\theta$的边缘先验分布的阶
• 马蹄估计的一致性
• 收敛速率
• • Clarke-Barron(1990)引理
  • 马蹄估计的最优收敛速率

这是对The horseshoe estimator for sparse signal这篇论文的回顾，这篇论文在Modern Bayesian statistics与Bayesian Machine Learning领域比较重要，它提出了一种新的先验结构——horseshoe prior，基于这种先验得到的均值的后验估计在shrinkage上具有类似LASSO的性质，也就是保留数据中稀疏的信号、让噪声缩减为0。

基本框架

假设观察值是$y_i$，它服从正态分布
$$y_i{\sim }_{iid}N({\theta }_i,{\sigma }^2)$$

我们引入下面的分为三层的先验：
$$\begin{gathered} {\theta }_i|{\lambda }_i\sim N(0,{\lambda }_i^2) \\ {\lambda }_i|\tau \sim C^{+}(0,\tau ) \\ \tau |\sigma \sim C^{+}(0,\sigma ) \end{gathered}$$

其中$C^{+}(0,a)$是尺度参数为$a$的half-Cauchy分布。假设$\sigma$的先验是Jeffrey先验，即密度函数与$1/\sigma$成正比。

half-Cauchy分布

如果$X\sim C^{+}(0,a),a>0$，则称$X$服从尺度参数为$a$的half-Cauchy分布，它的密度函数是
$$f(x)=\frac{2a}{\pi (x^2+a^2)}$$

我们先验证一下归一性：
$${\int }_0^{\infty }f(x)dx={\int }_0^{\infty }\frac{2a}{\pi (x^2+a^2)}dx=\frac{2}{\pi }\arctan (x/a){|}_0^{\infty }=1$$

当然这个分布的期望也是不存在的
$${\int }_0^{+\infty }xf(x)dx={\int }_0^{\infty }\frac{2ax}{\pi (x^2+a^2)}dx=\frac{a}{\pi }\ln (x^2+a^2){|}_0^{+\infty }=+\infty$$

为什么它叫马蹄估计量

考虑${\lambda }_i$的边缘先验分布，
$$\begin{gathered} p({\lambda }_i,\tau ,\sigma )\propto \frac{\tau }{\pi ({\lambda }_i^2+{\tau }^2)}\frac{\sigma }{\pi ({\tau }^2+{\sigma }^2)}\frac{1}{\sigma }\propto \frac{\tau }{({\lambda }_i^2+{\tau }^2)({\tau }^2+{\sigma }^2)} \\ p({\lambda }_i)\propto {\int }_0^{\infty }{\int }_0^{\infty }\frac{\tau }{({\lambda }_i^2+{\tau }^2)({\tau }^2+{\sigma }^2)}d\sigma d\tau \end{gathered}$$

定义${\kappa }_i=1/(1+{\lambda }_i^2)$，这个量在Bayesian shrinkage中非常重要，我们在下一个小标题介绍它的意义，但我们可以先分析它的先验分布。现在我们只想做一点定性分析，了解一下${\kappa }_i$的先验的形状，所以简单起见假设$\sigma =\tau =1$，于是
$$\begin{gathered} p({\lambda }_i)=\frac{2}{\pi ({\lambda }_i^2+1)},{\lambda }_i>0 \\ p(k_i)=p({\lambda }_i({\kappa }_i))|{\lambda }_i^{\prime }|=\frac{1}{\pi }{x}^{-1/2}(1-x{)}^{-1/2},x\in (0,1) \end{gathered}$$

因此$k_i\sim Beta(1/2,1/2)$，懒得自己画图我就扒了百度百科的图，看$\alpha =\beta =0.5$（粉红色）那条，那就是我们${\kappa }_i$的先验分布，是不是非常像一个马蹄铁的形状，所以这种先验结构被称为马蹄先验，基于这种先验的贝叶斯方法被称为马蹄估计。

后验均值、shrinkage与$\kappa$

现在来填一个小坑，$\kappa$为什么重要？我们可以做一点简单的推导来理解$\kappa$的含义，考虑非常简单的情况，固定$\tau =\sigma =1$，先验可以被简化为两层$$\begin{gathered} \theta |\lambda \sim N(0,{\lambda }^2) \\ \lambda \sim C^{+}(0,1) \end{gathered}$$

定义$\kappa =1/(1+{\lambda }^2)$，则${\lambda }^2=\frac{1-\kappa }{\kappa }$
$$\begin{gathered} p(k)=p(\lambda (\kappa ))|{\lambda }^{\prime }|=\frac{1}{\pi }{\kappa }^{-1/2}(1-\kappa {)}^{-1/2},\kappa \in (0,1) \\ p(\theta |\kappa )=\frac{1}{\sqrt{2\pi {\lambda }^2}}{e}^{-\frac{{\theta }^2}{2{\lambda }^2}}=\sqrt{\frac{\kappa }{2\pi (1-\kappa )}}{e}^{-\frac{\kappa {\theta }^2}{2(1-\kappa )}} \end{gathered}$$

我们尝试写出后验均值的表达式，如果$y\sim N(\theta ,1)$
$$\begin{gathered} E[\theta |y]={\int }_{-\infty }^{+\infty }\theta p(\theta |y)d\theta ={\int }_{-\infty }^{+\infty }\theta {\int }_0^{1}p(\theta ,\kappa |y)d\kappa d\theta \\ {=}_{Fubini}{\int }_0^{1}d\kappa {\int }_{-\infty }^{+\infty }\theta p(\theta ,\kappa |y)d\theta \end{gathered}$$

下面我们用Hierarchical技巧：
$$p(\theta ,\kappa |y)=p(\theta |\kappa ,y)p(\kappa |y)$$

于是
$$\begin{gathered} {\int }_0^{1}d\kappa {\int }_{-\infty }^{+\infty }\theta p(\theta ,\kappa |y)d\theta ={\int }_0^{1}p(\kappa |y)d\kappa {\int }_{-\infty }^{+\infty }\theta p(\theta |\kappa ,y)d\theta \\ ={\int }_0^{1}E[\theta |\kappa ,y]p(\kappa |y)d\kappa \end{gathered}$$

因为$\theta |\kappa \sim N(0,\frac{1-\kappa }{\kappa }),y|\theta \sim N(\theta ,1)$，这是一个经典的共轭分布族，直接用正态分布共轭分布族的结论：
$$E[\theta |\kappa ,y]=(1-\kappa )y$$

所以
$${\int }_0^{1}E[\theta |\kappa ,y]p(\kappa |y)d\kappa ={\int }_0^1(1-\kappa )yp(\kappa |y)d\kappa =(1-E[\kappa |y])y$$

综上
$$E[\theta |y]=(1-E[\kappa |y])y$$

这是$\theta$的后验均值，也是它贝叶斯估计，我们把$y$理解为信号的一个观测，$\theta$代表信号，$N(0,1)$是噪声，于是这个估计的作用是从观测中去掉噪声还原信号，不难发现
$$\frac{y-E[\theta |y]}{y}=E[\kappa |y]$$

也就是说我们基于马蹄估计还原的信号比观测更小，于是马蹄估计是一种shrinkage estimation，relative amount of shrinkage等于$E[\kappa |y]$，我们希望信号的$E[\kappa |y]$非常小，噪声的$E[\kappa |y]$非常大，这样我们就可以用这个方法把信号从噪声观测中还原出来，下文会分析马蹄估计理论上具有这个性质。

$\theta$的边缘先验分布的阶

$\theta$的边缘先验分布为
$$p(\theta )={\int }_0^{\infty }\frac{1}{\sqrt{2\pi }\lambda }{e}^{-\frac{{\theta }^2}{2{\lambda }^2}}\frac{2}{\pi (1+{\lambda }^2)}d\lambda$$

它关于${\theta }^{-1}$的阶是我们比较关注的量，因为现代贝叶斯统计认为$\theta$的先验在$0$处的值越大，并且在尾部关于${\theta }^{-1}$的阶越小，用这种先验来做稀疏数据分析的效果就越好。

定理1 关于这个边缘先验分布$p(\theta )$，我们有下面两个结论：

${\lim }_{\theta \to 0}p(\theta )=\infty$

$\frac{K}{2}\log (1+\frac{4}{{\theta }^2})<p(\theta )<K\log (1+\frac{2}{{\theta }^2}),K=(2{\pi }^3{)}^{-1/2}$

证明
第一条，如果$\theta \to 0$，讨论积分
$$2K{\int }_0^{\infty }\frac{1}{\lambda }\frac{1}{(1+{\lambda }^2)}d\lambda =2K[\ln |x|-0.5\ln (1+x^2)]{|}_0^{+\infty }=\infty$$

所以${\lim }_{\theta \to 0}p(\theta )=\infty$。

第二条，做换元$u=1/{\lambda }^2$，
$$p(\theta )=K{\int }_0^{+\infty }\frac{1}{1+u}{e}^{-\frac{{\theta }^{2}u}{2}}du$$

再做换元$z=1+u$，
$$p(\theta )=K{e}^{\frac{{\theta }^2}{2}}{\int }_1^{\infty }\frac{1}{z}{e}^{-\frac{z{\theta }^2}{2}}dz=K{e}^{\frac{{\theta }^2}{2}}{E}_1(\frac{{\theta }^2}{2})$$

这里的$E_1$是一个特殊函数，它被称为exponential integral function，关于它有一个不等式
$$0.5e^{-t}\log (1+\frac{2}{t})<E_1(t)<e^{-t}\log (1+\frac{1}{t}),\forall t>0$$

所以
$$\frac{K}{2}\log (1+\frac{4}{{\theta }^2})<p(\theta )<K\log (1+\frac{2}{{\theta }^2}),K=(2{\pi }^3{)}^{-1/2}$$

当$\theta$足够大时，上界与下界关于${\theta }^{-1}$的阶都是2，因为
$$\begin{gathered} \log (1+\frac{2}{{\theta }^2})=2{\theta }^{-2}+o({\theta }^{-2}) \\ \log (1+\frac{4}{{\theta }^2})=4{\theta }^{-2}+o({\theta }^{-2}) \end{gathered}$$

马蹄估计的一致性

定理2 假设$p(|y-\theta |)$是给定样本数据为$y$的似然函数，这个式子表示分布的参数是$\theta$，并且它是一个位置参数；假设$\theta$的边缘先验分布为$p(\theta )$，它是一个混合高斯分布，即$\theta |\lambda \sim N(0,{\lambda }^2)$，其中$\lambda$的先验为$p(\lambda ),\lambda >0$，假设似然函数与边缘先验分布$p(\theta )$使$y$的边缘密度$m(y)$有界，定义下面三个量
$$\begin{gathered} m^{\ast }(y)=\int p(|y-\theta |)p^{\ast }(\theta )d\theta \\ {p}^{\ast }(\theta )={\int }_{\lambda >0}p(\theta |\lambda )p^{\ast }(\lambda )d\lambda \\ {p}^{\ast }(\lambda )={\lambda }^{2}p(\lambda ) \end{gathered}$$

则
$$E[\theta |y]=\frac{1}{m(y)}\frac{d}{dy}{m}^{\ast }(y)$$

如果$p(|y-\theta |)$是正态的，那么
$$E[\theta |y]=y+\frac{d}{dy}\log m(y)$$

评注
这个定理不但看上去比较长，而且也是一个非常重要的定理。在原文中，作者们还证明了如果$y\sim N(\theta ,1),\theta |\lambda \sim N(0,{\lambda }^2),\lambda \sim C^{+}(0,\tau )$，$\tau$为已知的常数，则$|y-E[\theta |y]|$有界，并且
$$\underset{|y|\to \infty }{\lim }\frac{d}{dy}\log m(y)\to 0$$

这个结果说明了马蹄估计的一致性。这个结果的证明并不复杂，它的要点在于如何计算$\frac{d}{dy}\log m(y)$这个值，其中
$$m(y)=K{\int }_0^{\infty }{e}^{-\frac{y^2}{2(1+{\tau }^2{\lambda }^2)}}\frac{1}{\sqrt{1+{\lambda }^2{\tau }^2}}\frac{1}{1+{\lambda }^2}d\lambda$$

原文的思路是做换元，$z=\frac{1}{1+{\tau }^2{\lambda }^2}$，然后把积分凑成${\Phi }_1$函数，这个函数被称为Humbert series，
$${\Phi }_1(a,b,c;x,y)=\frac{\Gamma (c)}{\Gamma (a)\Gamma (c-a)}{\int }_0^1{t}^{a-1}(1-t{)}^{c-a-1}(1-xt{)}^{-b}{e}^{yt}dt$$

根据Gordy (1998)等式，Humbert series可以用Kummer第一类函数$M(a,b;x)$表示，
$${\Phi }_1(a,b,c;x,y)=e^x\sum _{n=0}^{\infty }\frac{(a{)}_n(b{)}_n}{(c{)}_n}\frac{y^n}{n!}M(c-a,c+n,-x)$$

Kummer第一类函数$M(a,b;x)$满足下面的等式(Slater 1960 Chapter 4)，
$$M(a,b;x)=\{\begin{array}{l}\frac{\Gamma (a)}{\Gamma (b)}{e}^x{x}^{a-b}(1+O(x^{-1})),x>0\\ \frac{\Gamma (a)}{\Gamma (b-a)}(-x{)}^{-a}(1+O(x^{-1})),x<0\end{array}$$

基于这几个特殊函数进行计算即可验证一致性。

证明
有两个有用的等式
$$\begin{gathered} \frac{d}{dy}p(y-\theta )=-\frac{d}{d\theta }p(y-\theta ) \\ {\lambda }^2\frac{d}{d\theta }N(\theta |0,{\lambda }^2)=-\theta N(\theta |0,{\lambda }^2) \end{gathered}$$

第一个等式就是Schwartz定理，因为$p$是概率密度，已经是一阶导数了，而二阶导数与求导次序无关，所以第一个等式成立；第二个等式就是对正态分布求导。

所以对于
$$E[\theta |y]=\int \theta p(\theta |y)d\theta$$

我们用前文使用过的Hierarchical技巧，
$$\begin{gathered} p(\theta |y)=\frac{1}{m(y)}p(y|\theta )p(\theta |\lambda )p(\lambda ) \\ =\frac{1}{m(y)}p(y-\theta )N(\theta |0,{\lambda }^2)p(\lambda ) \end{gathered}$$

于是
$$\begin{gathered} \int \theta p(\theta |y)d\theta =\int \frac{1}{m(y)}\theta p(y-\theta )N(\theta |0,{\lambda }^2)p(\lambda )d\theta \\ =\frac{1}{m(y)}\int \frac{d}{dy}p(y-\theta )N(\theta |0,{\lambda }^2)p^{\ast }(\lambda )d\theta =\frac{1}{m(y)}\frac{d}{dy}{m}^{\ast }(y) \end{gathered}$$

收敛速率

作为理论研究者，我们除了关心统计模型的渐近性质（比如一致性）外，我们还想了解模型的非渐近性质(non-asymptotic properties)，比如收敛速率，下面给出的是马蹄估计的后验分布趋近于真实分布的速率。

Clarke-Barron(1990)引理

用${\theta }_0$表示参数的真实值，$p(y|\theta )$表示样本的概率密度，用$L(p_1,p_2)$表示从$p_2$到$p_1$的Kullback-Leibler divergence，
$$L(p_1,p_2)=E_{p_1}[\log (p_1/p_2)]$$

$\forall ϵ>0$，定义${\theta }_0$的Kullback-Leibler信息邻域为
$$A_{ϵ}=\{\theta :L(p_{{\theta }_0},p_{\theta })\le ϵ\}$$

用${\mu }_n$表示后验分布，$n$为样本数，$\mu$表示先验概率分布，密度函数的贝叶斯估计为
$${\hat{p}}_n=\int p_{\theta }{\mu }_n(d\theta )$$

定义贝叶斯估计的Cesaro average risk为
$$R_n=\frac{1}{n}\sum _{j=1}^{n}L(p_{{\theta }_0},{\hat{p}}_j)$$

Clarke-Barron(1990)引理
如果先验满足$\forall ϵ>0$，$\mu (A_{ϵ})>0$，则
$$R_n\le ϵ-\frac{1}{n}\mu (A_{ϵ})$$

马蹄估计的最优收敛速率

定理4 假设先验满足$\forall ϵ>0$，$\mu (A_{ϵ})>0$

如果$\theta$的边缘先验分布连续有上界且支撑集包含${\theta }_0$的至少一个邻域，则$R_n=O(n^{-1}\log n)$

马蹄估计量的Cesaro average risk具有下面的性质: 如果${\theta }_0=0$，$R_n=O(n^{-1}(\log n-b\log \log n))$，其中$b$是一个常数；如果${\theta }_0\ne 0$，$R_n=O(n^{-1}\log n)$

证明
第一条以及第二条${\theta }_0\ne 0$，取${ϵ}_n=1/n$，如果$\forall ϵ>0$，$\mu (A_{ϵ})>0$，则
$$\mu (A_{{ϵ}_n})={\int }_{A_{{ϵ}_n}}p(\theta )d\theta$$

这里的$A_{{ϵ}_n}$是Kullback-Leibler信息邻域
$$\begin{gathered} A_{{ϵ}_n}=\{\theta :L(p_{{\theta }_0},p_{\theta })\le \frac{1}{n}\}=\{\theta :E_{p_{{\theta }_0}}[\log (p_{{\theta }_0}/p_{\theta })]\le \frac{1}{n}\} \\ \subset \{\theta :|\theta -{\theta }_0|\le \frac{1}{\sqrt{n}}\} \end{gathered}$$

因为$p(\theta )$有上界，$\exists C>0$，$p(\theta )<C$
$$\mu (A_{{ϵ}_n})\le {\int }_{{\theta }_0-1/\sqrt{n}}^{{\theta }_0+1/\sqrt{n}}Cd\theta =2C/\sqrt{n}=O(1/\sqrt{n})$$

根据Clarke-Barron(1990)引理，
$$R_n\le \frac{1}{n}-\frac{\log (C{n}^{-1/2})}{n}=O(\log n/n)$$

第二条${\theta }_0=0$，根据定理一，
$$p(\theta )\ge K\log (1+4{\theta }^{-2})$$

所以
$$\mu (A_{ϵ})\ge K{\int }_0^{\sqrt{ϵ}}\log (1+4{\theta }^{-2})d\theta$$

做换元$u=1/{\theta }^2$，
$$\begin{gathered} \mu (A_{ϵ})\ge K{\int }_{4/ϵ}^{\infty }\frac{\log (1+u)}{u^{3/2}}du \\ =K{ϵ}^{1/2}\log (1+4{ϵ}^{-1})+2K{\int }_{4/ϵ}^{\infty }\frac{1}{u^{1/2}(1+u)}du \end{gathered}$$

取$ϵ=1/n$，根据Clarke-Barron(1990)引理，可以得到$R_n=O(n^{-1}(\log n-b\log \log n))$。

总结

以上是生活随笔为你收集整理的稀疏数据分析：马蹄估计量及其理论性质的全部内容，希望文章能够帮你解决所遇到的问题。

如果觉得生活随笔网站内容还不错，欢迎将生活随笔推荐给好友。