aMCMC for Horseshoe: algorithms

• • Exact Algorithm
  • • Step 1: sampling $\eta$
    • Step 2: sampling $\xi$
    • Step 3: sampling ${\sigma }^2$
    • Step 4: sampling $\beta$
  • Approximate Algorithm
  • Correlation and Traceplot

A limitation of Bayesian shrinkage: lack of computational efficiency under very-high dimension.

Framework
$$\begin{gathered} L(z|W\beta ,{\sigma }^2)=(2\pi {\sigma }^2{)}^{?\frac{N}{2}}{e}^{?\frac{1}{2{\sigma }^2}(z?W\beta {)}^{\prime }(z?W\beta )} \\ {\beta }_j|{\sigma }^2,{\eta }_j,\xi {～}_{iid}N(0,\frac{{\sigma }^2}{\xi {\eta }_j}),j=1,?,p \\ {\eta }^{?1/2}{～}_{iid}{C}^{+}(0,1),\xi ～C^{+}(0,1),{\sigma }^2～IG(\frac{w}{2},\frac{w}{2}) \end{gathered}$$

Exact Algorithm

A blocked Metropolis-within-Gibbs algorithm:

sample ${\eta }_j～p({\eta }_j|\xi ,\beta ,{\sigma }^2)$

MCMC $\xi$ with kernel $\log ({\xi }^{?})～N(\log (\xi ),s)$

sample ${\sigma }^2|\eta ,\xi ～IG(\frac{w+N}{2},\frac{w+z^{\prime }{M}_{\xi }^{?1}z}{2})$

sample $\beta |\eta ,\xi ,{\sigma }^2～N((W^{\prime }W+({\xi }^{?1}D{)}^{?1}{)}^{?1}{W}^{\prime }z,{\sigma }^2(W^{\prime }W+({\xi }^{?1}D{)}^{?1}{)}^{?1})$

where
$$\begin{gathered} D=diag({\eta }_j^{?1}),M_{\xi }=I_N+{\xi }^{?1}WD{W}^{\prime } \\ p(\xi |\eta )=|M_{\xi }{|}^{?1/2}(\frac{w}{2}+\frac{w+z^{\prime }{M}_{\xi }^{?1}z}{2}{)}^{?\frac{w+N}{2}}\frac{1}{\sqrt{\xi }(1+\xi )} \end{gathered}$$

Step 1: sampling $\eta$

$$\begin{gathered} p({\eta }_j|\xi ,\beta ,{\sigma }^2)\propto p(\beta |{\eta }_j,\xi ,{\sigma }^2)p({\eta }_j) \\ \propto \sqrt{{\eta }_j}{e}^{?\frac{{\beta }^2}{2{\sigma }^2{\eta }_j\xi }}p({\eta }_j) \end{gathered}$$

${\eta }_j^{?1/2}～C^{+}(0,1){=}_{d}X$,
$$\begin{gathered} p({\eta }_j)=\frac{2}{\pi }\frac{1}{1+({\eta }_j^{?1/2}{)}^2}|?\frac{1}{2}{\eta }_j^{?3/2}| \\ =\frac{1}{\pi }\frac{1}{{\eta }_j^{1/2}+{\eta }_j^{3/2}} \end{gathered}$$

$$p({\eta }_j|\xi ,\beta ,{\sigma }^2)\propto \sqrt{{\eta }_j}{e}^{?\frac{{\beta }^2}{2{\sigma }^2{\eta }_j\xi }}\frac{1}{{\eta }_j^{1/2}+{\eta }_j^{3/2}}=\frac{1}{1+{\eta }_j}{e}^{?\frac{{\beta }^2{\eta }_j\xi }{2{\sigma }^2}}$$

Rejection sampler (Appendix S1)
Given $?\in (0,1)$, this sampler is used to sample from
$$h_{?}(t)=C_{?}\frac{e^{??t}}{1+t},t>0$$

where
$$\begin{gathered} {\int }_{t>0}{h}_{?}(t)dt={\int }_{t>0}{C}_{?}\frac{e^{??t}}{1+t}dt=1?C_{?}=\frac{e^{??}}{E_1(?)} \\ {E}_1(?)={\int }_{?}^{+\infty }\frac{e^{?t}}{t}dt \end{gathered}$$

Algorithm
Step 1: Draw z and u independently
$z～h_L,u～U(0,1)$. To draw $z～h_L$，

Step 2: accept z as sample from h with acceptance rate
Here acceptance rate is $e^{?(f?f_L)(z)}$.

Step 2: sampling $\xi$

Transition: $\log ({\xi }^{?})～N(\log (\xi ),s)$
Acceptance rate is
$$\alpha =\min (\frac{g(\xi |{\xi }^{?})p({\xi }^{?}|\eta ,z)}{g({\xi }^{?}|\xi )p(\xi |\eta ,z)},1)$$

Calculate $p(\xi |\eta ,z)$ (only normalization parameter left)
$$\begin{gathered} p(\xi |\eta ,z)\propto ?p(z|\beta ,{\sigma }^2)p(\beta |{\sigma }^2,\eta ,\xi )p({\sigma }^2)p(\xi )d\beta d{\sigma }^2 \\ \propto |M_{\xi }{|}^{?1/2}(\frac{w}{2}+\frac{1}{2}{z}^{\prime }{M}_{\xi }z{)}^{?\frac{N+w}{2}}\frac{1}{(\xi {)}^{1/2}+(\xi {)}^{3/2}}=p(\xi |\eta ) \end{gathered}$$

Calculate
$$\frac{g(\xi |{\xi }^{?})p({\xi }^{?}|\eta ,z)}{g({\xi }^{?}|\xi )p(\xi |\eta ,z)}=\frac{\frac{1}{\sqrt{2\pi s}\xi }{e}^{?\frac{(\log \xi ?\log {\xi }^{?}{)}^2}{2s}}p({\xi }^{?}|\eta )}{\frac{1}{\sqrt{2\pi s}{\xi }^{?}}{e}^{?\frac{(\log {\xi }^{?}?\log \xi {)}^2}{2s}}p(\xi |\eta )}=\frac{p({\xi }^{?}|\eta ){\xi }^{?}}{p(\xi |\eta )\xi }$$

Step 3: sampling ${\sigma }^2$

$$\begin{gathered} p({\sigma }^2|\eta ,\xi )\propto \int p(z|\beta ,{\sigma }^2)p(\beta |{\sigma }^2,\eta ,\xi )p({\sigma }^2)d\beta \\ =IG(\frac{w+N}{2},\frac{w+z^{\prime }{M}_{\xi }^{?1}z}{2}) \end{gathered}$$

Step 4: sampling $\beta$

$$\begin{gathered} p(\beta |\eta ,\xi ,{\sigma }^2,z)\propto p(z|\beta ,{\sigma }^2)p(\beta |{\sigma }^2,\eta ,\xi ) \\ =N(z|W^{\prime }\beta ,{\sigma }^2)N(\beta |0,{\sigma }^2{\xi }^{?1}D) \\ =N(\beta |(W^{\prime }W+({\xi }^{?1}D{)}^{?1}{)}^{?1}{W}^{\prime }z,{\sigma }^2(W^{\prime }W+({\xi }^{?1}D{)}^{?1}{)}^{?1}) \end{gathered}$$

Bhattacharya (2016)'s algorithm

Main computational bottleneck: $M_{\xi }=I_N+{\xi }^{?1}WD{W}^{\prime }$, cost $N^{2}p$ (if $N>p$, see Hahn et al (2018) otherwise use aMCMC).

Approximate Algorithm

key idea: $\beta$ is sparse, so not every column in $W$ worth computing. This may help reduce computational cost.

Observation: ${\beta }_j{～}_{iid}N(0,{\sigma }^2{\xi }^{?1}{\eta }_j^{?1})$, if ${\beta }_j$ shrunk to 0, $\xi {\eta }_j$ must be large;
$${\xi }^{?1}WD{W}^{\prime }=\sum _{j=1}^p{\xi }^{?1}{\eta }_j^{?1}{w}_j{w}_j^{\prime }$$

where $w_i$ is the i-th column of $W$, and in this case, it has no contribution. To avoid computing with these columns, import hard-threshold

$$D_{\delta }=diag({\eta }_j^{?1}{1}_{{\xi }_{max}^{?1}{\eta }_j^{?1}>\delta }),{\xi }_{max}=\max (\xi ,{\xi }^{?})$$

Approximate MCMC:

sample ${\eta }_j～p({\eta }_j|\xi ,\beta ,{\sigma }^2)$

MCMC $\xi$ with kernel $\log ({\xi }^{?})～N(\log (\xi ),s)$, ${\xi }_{max}=\max (\xi ,{\xi }^{?})$, replace $M_{\xi }$ by $M_{\xi ,\delta }$, $$M_{\xi ,\delta }=I_N+{\xi }^{?1}{W}_S{D}_S{W}_S^{\prime },S=\{j:{\xi }_{max}^{?1}{\eta }_j^{?1}>\delta \}$$and compute its inverse by Woodbury’s identity $$M_{\xi ,\delta }^{?1}=I_N?W_S(\xi D_S^{?1}+W_S^{\prime }{W}_S{)}^{?1}{W}_S^{\prime }$$

sample ${\sigma }^2|\eta ,\xi ～IG(\frac{w+N}{2},\frac{w+z^{\prime }{M}_{\xi ,\delta }^{?1}z}{2})$

sample $\beta |\eta ,\xi ,{\sigma }^2～N((W^{\prime }W+({\xi }^{?1}{D}_{\delta }{)}^{?1}{)}^{?1}{W}^{\prime }z,{\sigma }^2(W^{\prime }W+({\xi }^{?1}{D}_{\delta }{)}^{?1}{)}^{?1})$

Computational cost:
Exact $O(N^3)+O(N^{2}p)$
Approximate $O(s_{\delta }^{2}N)+O(Np),s_{\delta }=|S|$

Correlation and Traceplot

總結(jié)

以上是生活随笔為你收集整理的aMCMC for Horseshoe: algorithms的全部內(nèi)容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯(cuò)，歡迎將生活随笔推薦給好友。