Bayesian Fusion Estimation via t-Shrinkage

TL;DR

This paper introduces a Bayesian t-shrinkage prior for high-dimensional parameter modeling with unknown block structures, demonstrating superior performance over existing methods in simulations and real data analyses.

Contribution

It proposes a novel t-fusion prior that induces stronger shrinkage and posterior consistency, extending Bayesian fusion and clustering techniques.

Findings

t-fusion outperforms Bayesian Laplace-fusion in simulations and real data.

The method achieves better posterior convergence than Dirichlet process models.

t-fusion effectively induces posterior blocking in high-dimensional parameters.

Abstract

Shrinkage prior has gained great successes in many data analysis, however, its applications mostly focus on the Bayesian modeling of sparse parameters. In this work, we will apply Bayesian shrinkage to model high dimensional parameter that possesses an unknown blocking structure. We propose to impose heavy-tail shrinkage prior, e.g., $t$ prior, on the differences of successive parameter entries, and such a fusion prior will shrink successive differences towards zero and hence induce posterior blocking. Comparing to conventional Bayesian fused lasso which implements Laplace fusion prior, $t$ fusion prior induces stronger shrinkage effect and enjoys a nice posterior consistency property. Simulation studies and real data analyses show that $t$ fusion has superior performance to the frequentist fusion estimator and Bayesian Laplace-fusion prior. This $t$-fusion strategy is further developed to conduct a Bayesian clustering analysis, and simulation shows that the proposed algorithm obtains better posterior distributional convergence than the classical Dirichlet process modeling.

Figures

Tables

Table 1. Table 1: Comparison among Bayesian t-fusion, Bayesian Laplace fusion, and L 1 subscript 𝐿 1 L_{1} fusion. The report are based on average results from 100 replications. Refer to Section 4.1 for the definition of R 𝑅 R , W 𝑊 W and B 𝐵 B .

Methods	Bayesian t-fusion	Bayesian Laplace fusion	$L_1$ fusion
$L_2$ error of $\theta$	1.3243	2.2323	1.5916
Standard error	0.0602	0.0324	0.0421
$L_1$ error of $\theta$	9.6964	16.6296	11.6745
Standard error	0.3419	0.2725	0.3407
$W$ (the smaller the better)	0.05584	0.2182	0.2128
Standard error	0.0028	0.0037	0.0077
$B$ (the larger the better)	1.4243	0.6792	1.1302
Standard error	0.0684	0.0211	0.0398
$R$ (the smaller the better)	387.32	85.6	1360.9
Standard error	32.0628	6.2375	53.2685

Table 2. Table 2: Comparison among Bayesian t-shrinkage, Bayesian DP prior, and L 1 subscript 𝐿 1 L_{1} fusion. The report are based on average results from 100 replications. Refer to Section 4.2 for the definition of R 𝑅 R , W 𝑊 W and B ~ ~ 𝐵 \tilde{B} .

Methods	$t$-shrinkage (10)	$t$-shrinkage (9)	DP prior(13)	$L_1$ fusion (12)
posterior mean of $L_2^2$ error	20.559	20.5660	40.7029	—
Standard error	0.6885	0.4589	3.1209	—
$L_2$ error of $\theta$	4.1873	4.3925	3.9584	5.0039
Standard error	0.0872	0.0531	0.1621	0.0739
$L_1$ error of $\theta$	28.1099	33.4543	26.8122	43.8450
Standard error	0.7385	0.4845	1.7122	0.7909
$W$ (the smaller the better)	0.3517	0.4428	0.2183	0.4551
Standard error	0.0093	0.0063	0.0061	0.0062
$\tilde{B}$ (the larger the better)	1.4652	1.3392	1.3885	0.9855
Standard error	0.0218	0.0147	0.0419	0.0146
$R$ (the smaller the better)	1752.24	2294.76	—	1973.7
Standard error	31.5641	14.6102	—	11.9937

Equations

$$y_i={\theta }_i^{\ast }+{\epsilon }_i$$

$${\hat{\theta }}^{\mathrm{FL}}=\arg \min \left(\frac{\Vert y-\theta {\Vert }^2}{2}+{\lambda }_1\sum _{i=1}^n|{\theta }_i|+{\lambda }_2\sum _{i=2}^n|{\vartheta }_i|\right).$$

$${\hat{\theta }}^{\mathrm{F}}=\arg \min \left(\frac{\Vert y-\theta {\Vert }^2}{2}+\lambda \sum _{i=2}^n|{\vartheta }_i|\right)=\arg \min \left(\frac{\Vert y-\theta {\Vert }^2}{2}+\lambda \sum _{i=2}^n|{\theta }_i-{\theta }_{i-1}|\right),$$

$$\begin{array}{rl}{\sigma }^2& \sim \text{mbox}Inverse-Gamma(a_{\sigma },b_{\sigma }),{\theta }_1|{\sigma }^2\sim N(0,{\sigma }^2{\lambda }_1),\\ ({\theta }_i-{\theta }_{i-1})|{\sigma }^2& \sim \text{mbox}Laplace(\lambda /\sigma ),\text{mbox}foralli=2,\dots ,n\end{array}$$

$$({\theta }_i-{\theta }_{i-1})|{\sigma }^2,{\lambda }_i\sim N(0,{\lambda }_i{\sigma }^2),{\lambda }_i\sim \text{mbox}exp(-{\lambda }^2/2),$$

$$\begin{array}{rl}{\lambda }_i^{-1}\sim & \text{mbox}Inverse-Gaussian(\lambda \sigma /|{\theta }_i-{\theta }_{i-1}|,{\lambda }^2)\text{mbox}foralli=2,\dots ,n,\\ {\sigma }^2\sim & \text{mbox}Inverse-Gamma\left(a_{\sigma }+n,b_{\sigma }+\frac{\Vert y-\theta {\Vert }^2}{2}+\frac{{\theta }_1^2}{2{\lambda }_1}+\sum _{i=2}^n\frac{({\theta }_i-{\theta }_{i-1}{)}^2}{2{\lambda }_i}\right),\\ {\theta }_i\sim & N({\mu }_i,{\nu }_i),\end{array}$$

$$\begin{array}{rl}{\sigma }^2& \sim \text{mbox}Inverse-Gamma(a_{\sigma },b_{\sigma }),{\theta }_1|{\sigma }^2\sim N(0,{\sigma }^2{\lambda }_1),\\ ({\theta }_i-{\theta }_{i-1})|{\sigma }^2& \sim t_{df}(s\sigma ),\text{mbox}foralli=2,\dots ,n\end{array}$$

$$({\theta }_i-{\theta }_{i-1})|{\sigma }^2,{\lambda }_i\sim N(0,{\lambda }_i{\sigma }^2),{\lambda }_i\sim \text{mbox}Inverse-Gamma(a_t,b_t),$$

$${\lambda }_i\sim \text{mbox}Inverse-Gamma\left(a_t+1/2,b_t+\frac{({\theta }_i-{\theta }_{i-1}{)}^2}{2{\sigma }^2}\right)\text{mbox}foralli=2,\dots ,n.$$

$$\begin{split}&\int_{-|G^{*}|\log(n)/n^{2}}^{|G^{*}|\log(n)/n^{2}}f_{\vartheta}(x)dx\geq 1-n^{1+u},\mbox{ for some $u>0$, }\\ &-\log(\underline{\pi}_{\vartheta})=O(\log n),\quad\mbox{ where }\underline{\pi}_{\vartheta}=\min_{|x|\leq\max_{i}|\vartheta^{*}_{i}/\sigma^{*}|+1}f_{\vartheta}(x),\\ &-\log(\underline{\pi}_{\theta})=O(|G^{*}|\log n),\quad\mbox{ where }\underline{\pi}_{\theta}=\min_{|x|\leq|\theta^{*}_{1}/\sigma^{*}|+1}f_{\theta}(x),\\ &a_{\sigma}\log(1/b_{\sigma})+b_{\sigma}/\sigma^{*2}+(a_{\sigma}+2)\log(\sigma^{*2})=O(|G^{*}|\log n),\end{split}$$

$$\pi (\Vert \theta -{\theta }^{\ast }\Vert \ge M{\sigma }^{\ast }\sqrt{n}{ϵ}_n|y)\to 0,$$

$$\pi (\{|G(\theta ,\sigma )\backslash G^{\ast }|>\delta |G^{\ast }|\}|y)\to 0,$$

$$\begin{array}{rl}& \log (\underset{i}{\max }|{\vartheta }_i^{\ast }|/{\sigma }^{\ast })=O(\log (n)),\\ & {\theta }_1^{\ast 2}/({\lambda }_1{\sigma }^{\ast 2})+\log ({\lambda }_1)=O(\log n),\end{array}$$

$$P(|t_{df}(s)|\ge \sqrt{\log (n)/n})\approx 1/n.$$

$$\begin{array}{rl}{\sigma }^2& \sim \text{mbox}Inverse-Gamma(a_{\sigma },b_{\sigma }),{\theta }_{\hat{r}(1)}|{\sigma }^2\sim N(0,{\sigma }^2{\lambda }_1),\\ ({\theta }_{\hat{r}(i)}-{\theta }_{\hat{r}(i-1)})|{\sigma }^2,{\lambda }_i& \sim N(0,{\sigma }^2{\lambda }_i),{\lambda }_i\sim \text{mbox}Inverse-Gamma(a_t,b_t)\text{mbox}foralli=2,\dots ,n.\end{array}$$

$$\begin{array}{rl}{\lambda }_i\sim & \text{mbox}Inverse-Gamma\left(a_t+1/2,b_t+\frac{({\theta }_{\hat{r}(i)}-{\theta }_{\hat{r}(i-1)}{)}^2}{2{\sigma }^2}\right)\text{mbox}foralli=2,\dots ,n.\\ {\sigma }^2\sim & \text{mbox}Inverse-Gamma\left(a_{\sigma }+n,b_{\sigma }+\frac{\Vert y-\theta {\Vert }^2}{2}+\frac{{\theta }_{\hat{r}(1)}^2}{2{\lambda }_1}+\sum _{i=2}^n\frac{({\theta }_{\hat{r}(i)}-{\theta }_{\hat{r}(i-1)}{)}^2}{2{\lambda }_i}\right),\\ {\theta }_{\hat{r}(i)}\sim & N({\mu }_i,{\nu }_i),\end{array}$$

$$\pi (\Vert \theta \Vert \le \sqrt{M\log n}|y)<1/2,$$

$$\begin{split}&\mbox{update $\theta_{r(i)}$'s, $\lambda_{i}$'s and $\sigma^{2}$ according to (\ref{gibbs2})},\\ &\mbox{update the $r(\cdot)$ to be the rank statistic of the current sample $\theta$}.\end{split}$$

$$\begin{array}{rl}r& \sim \text{mbox}Uniformdistributionoverallpossiblepermutations,\\ {\sigma }^2& \sim \text{mbox}Inverse-Gamma(a_{\sigma },b_{\sigma }),{\theta }_{r(1)}|{\sigma }^2\sim N(0,{\sigma }^2{\lambda }_1),\\ ({\theta }_{r(i)}-{\theta }_{r(i-1)})|{\sigma }^2,{\lambda }_i& \sim N(0,{\sigma }^2{\lambda }_i),{\lambda }_i\sim \text{mbox}Inverse-Gamma(a_t,b_t)\text{mbox}foralli=2,\dots ,n.\end{array}$$

$$\pi (\Vert \theta -{\theta }^{\ast }\Vert \ge M{\sigma }^{\ast }\sqrt{n}{ϵ}_n|y)\to 0,$$

$$\hat{\theta }=\arg \min \Vert y-\theta {\Vert }^2/2+\lambda \sum _{i=2}^n|{\theta }_{\hat{r}(i)}-{\theta }_{\hat{r}(i-1)}|,$$

$$\begin{array}{rl}& y_i|{\theta }_i,{\sigma }^2\sim N({\theta }_i,{\sigma }^2),{\theta }_i\sim G,\\ & G\sim \text{mbox}Dirichlet-Process(N(0,\lambda ),a)\text{mbox}and{\sigma }^2\sim \text{mbox}Inverse-Gamma(a_{\sigma },b_{\sigma })\end{array}$$

$$\tilde{B}=\text{mbox}thebottom10\%quantileof\{|{\hat{\theta }}_i-{\hat{\theta }}_j|{\}}_{\{w_{ij}^{\ast }=0\}}.$$

$$\begin{array}{rl}& Pr({\chi }_d^2(\kappa )>d+\kappa +2x+\sqrt{(4d+8\kappa )x})\le \exp (-x),\text{mbox}and\\ & Pr({\chi }_d^2(\kappa )<d+\kappa -\sqrt{(4d+8\kappa )x})\le \exp (-x).\end{array}$$

$$Pr(X\ge k+1)\le 1-\Phi (\text{mbox}sign(k-nv)\{2nH(v,k/n){\}}^{1/2}),$$

$$P^{\ast }\left\{\frac{m(D_n)}{f^{\ast }(D_n)}\ge a_n\right\}\ge 1-a_n^{\prime },$$

$$\begin{array}{rl}& E^{\ast }\left(\pi (C_n\cup B_n)|D_n)\right)\le \frac{b_n+c_n}{a_n}+a_n^{\prime }+b_n^{\prime }.\end{array}$$

$$\begin{array}{rl}& E^{\ast }\pi (C_n\cup B_n)|D_n)=E^{\ast }\pi (C_n\cup B_n)|D_n)(1-\varphi (D_n))1_{{\Omega }_n}\\ +& E^{\ast }\pi (C_n\cup B_n)|D_n)(1-\varphi (D_n))(1-1_{{\Omega }_n})+E^{\ast }\pi (C_n\cup B_n)|D_n)\varphi (D_n)\\ \le & E^{\ast }\pi (C_n\cup B_n)|D_n)(1-\varphi (D_n))1_{{\Omega }_n}+E^{\ast }(1-1_{{\Omega }_n})+E^{\ast }\varphi (D_n)\\ \le & E^{\ast }\pi (C_n\cup B_n)|D_n)(1-\varphi (D_n))1_{{\Omega }_n}+b_n^{\prime }+a_n^{\prime }\\ \le & E^{\ast }\{m(D_n,C_n\cup B_n)/a_n{f}^{\ast }(D_n)\}(1-\varphi (D_n))+b_n^{\prime }+a_n^{\prime }.\end{array}$$

$$\begin{array}{rl}& E^{\ast }(1-\varphi (D_n))m(D_n,C_n\cup B_n)/f^{\ast }(D_n)={\int }_{C_n\cup B_n}{\int }_{\mathcal{X}}[1-\varphi (D_n)]f_{\theta }(D_n)d{D}_n\pi (\theta )d\theta \\ \le & {\int }_{C_n}{E}_{\theta }(1-\varphi (D_n))\pi (\theta )d\theta +{\int }_{B_n}{\int }_{\mathcal{X}}{f}_{\theta }(D_n)d{D}_n\pi (\theta )d\theta \le b_n+c_n.\end{array}$$

$$\begin{array}{rl}\varphi (y)=1\{& \Vert {\hat{\theta }}_G-{\theta }^{\ast }\Vert \ge \sqrt{n}{\sigma }^{\ast }{ϵ}_n\text{mbox}and|{\hat{\sigma }}_G^2-{\sigma }^{\ast 2}|>{\sigma }^{\ast 2}{ϵ}_n\\ & \text{mbox}forallG\supset G^{\ast },|G|\le (1+\delta )|G^{\ast }|\}\end{array}$$

Full text

Keywords: $t$-shrinkage prior; Bayesian fusion; Bayesian clustering; posterior consistency

1 Introduction

High dimensionality plays an important role in modern statistical applications such as genomics, image processing, finance and etc. An overview for the development of high dimensional analysis can be found in [47] and references therein. To overcome ill-posed problems that involve high dimensional parameters, one usually assumes that the true parameter value lies in a low dimensional subspace. To obtain such low dimensional estimation, the idea of regularization is commonly used, via penalized likelihood approaches or using informative prior specifications. Various penalty functions have been proposed for consistent frequentist estimation, including Lasso [44], SCAD [11], adaptive Lasso [53] and MCP [52]. For high dimensional Bayesian inferences, sparsity induced prior, such as spike-and-slab prior [17, 28, 32, 41, 51, 7, 38, 18], is widely used for model selection. Fused Lasso [45] considers another type of low dimensional embedding of a high dimensional parameter $\theta=(\theta_{i})_{i=1}^{p}$ where the successive differences $\vartheta_{i}=\theta_{i}-\theta_{i-1}$ are assumed to be sparse as well, in other words, there exists a consecutive block partition of $\theta_{i}$’s, such that $\theta_{i}$’s are constant within each block. Fused lasso method proposes a penalty function $\lambda_{1}\sum|\theta_{i}|+\lambda_{2}\sum|\vartheta_{i}|$ which consists of two terms that encourage sparsity among $\theta_{i}$’s and $\vartheta_{i}$’s respectively.

In this work, we consider the following Gaussian mean problem:

[TABLE]

where $\varepsilon_{i}$’s are iid normal error with unknown variance $\sigma^{2}$. Similarly to Fused Lasso applications, we also assume true parameter $\theta^{*}$ is blocky in the sense that there exists a partition $\{{\cal B}_{1}^{*},\dots,{\cal B}_{s}^{*}\}$ of $\{1,...,n\}$ such that $\theta^{*}_{i}$’s are constant for all $i\in{\cal B}_{k}^{*}$. Correspondingly, we define set $G^{*}=\{2\leq i\leq n:\vartheta_{i}^{*}:=\theta_{i}^{*}-\theta_{i-1}^{*}\neq 0\}$ whose number of elements is supposed to be much smaller than $n$. We are interested in conducting Bayesian structure recovery of $\theta^{*}$. Motivated by the $L_{1}$-fusion penalty used by Fused lasso estimator [45], as well as the development of the Bayesian lasso [34], [24] introduced Bayesian fused lasso by imposing independent Laplace priors on all successive differences. The implementation of Laplace shrinkage prior can significantly reduce the posterior sampling costs comparing to spike-and-slab modeling, and conceptually, the Laplace prior can be nicely interpreted as a Bayesian counterpart of $L_{1}$ penalty. However, many recent Bayesian theoretical developments show that in the context of sparse linear regression models, Laplace prior fails to achieve satisfactory posterior contraction [7, 5, 42]. It is believed that the posterior inconsistency of Laplace prior is due to its exponentially light tail, and [42] suggests to use heavy tail prior distribution for sparse linear regression models, which can induce sufficient Bayesian shrinkage effect and thereafter guarantee to recover the sparsity structure.

We find that the above phenomenon holds for Bayesian fusion estimation as well: imposing Laplace prior on $(\theta_{i}-\theta_{i-1})$ leads to a smoothly varying $\theta_{i}$ estimation rather than a blocky $\theta$, thus it fails to identify the blocking structure. Therefore, in this paper, we propose to use independent student-$t$ priors on successive differences $\theta_{i}-\theta_{i-1}$ for a Bayesian fusion problem. Our results show that such a simple $t$ fusion Bayesian modeling leads to very accurate posterior estimation. More importantly, comparing with Laplace prior or frequentist $L_{1}$ penalization, its performance on detecting the blocking structure is much better. The asymptotic posterior convergence induced by $t$ fusion prior is investigated as well. A related Bayesian work is [40] who proposed to use a Normal-Exponential-Gamma (NEG) prior for the successive differences. However, their Bayesian inference is only based on the maximum a posterior (MAP) estimator, while our application tries to fully utilize the whole posterior distributional information.

Furthermore, we consider a practically useful extension to Bayesian fusion estimation. Instead of assuming that $\theta^{*}$ has a consecutive blocking structure, it is more realistic to assume that $\theta^{*}$ possesses an unknown clustering structure. In other words, some unnecessarily consecutive $\theta_{i}^{*}$’s share the same value within an ${\cal B}_{s}^{*}$. In a broader scope, such a clustering problem can viewed as a simplest example of subgroup analysis where we assume that a subject-related parameter $\theta$ follows an unknown grouping structure. For example, in clinical trial studies, the treatment effects may vary across different subpopulations, but remain the same for the patients belonging to the same subpopulation. If one can correctly identify the subpopulation structure, then specific medical therapies can be prescribed for each subpopulation to maximize the treatment effectiveness. The existing Bayesian clustering analysis [49, 15, 31, 3] usually impose discrete priors on the clustering structure, along with a conditional prior on $\theta$ given specific clustering structure. In contrast, we propose to directly model the parameter $\theta$ via $t$ fusion prior. Such a prior specification allows a computationally efficient Gibbs posterior sampling algorithm. Our simulation shows that our procedure yields reasonable cluster structure recovery, and moreover it beats the usual Dirichlet process prior in terms of posterior contraction.

This paper is organized as follows. In Section 2, we study the Bayesian fusion problem with $t$ prior specification. We will present the posterior asymptotic result, and discuss its difference from the Laplace prior. In Section 3, we will use the $t$-fusion prior to solve the clustering problem. Several simulation studies and one real data application are presented in Section 4. At last, Section 5 provides more discussions and remarks. All technical proofs are postponed to the appendix section.

Throughout this work, the following notation is used. Given two positive sequences $\{a_{n}\}$ and $\{b_{n}\}$, $a_{n}\succ b_{n}$ means $\lim(a_{n}/b_{n})=\infty$ and $a_{n}\asymp b_{n}$ means $-\infty<\lim\inf(a_{n}/b_{n})\leq\lim\sup(a_{n}/b_{n})<\infty$. $\|x\|$ and $\|x\|_{1}$ denote $L_{2}$ and $L_{1}$ norms of vector $x$.

2 Bayesian fusion via $t$-shrinkage

2.1 Bayesian Modeling

Suppose we observe independent data $\{y_{i}:i=1,\dots,n\}$ following model (1). The indexing of the data has certain practical or scientific meaning, under which we can assume that the parameter vector $\theta^{*}$ is “stepwise”, in the sense that most of the successive differences $\vartheta_{i}=\theta_{i}-\theta_{i-1}$ are exactly 0. To induce the sparsity for both $\theta_{i}$’s and $\vartheta_{i}$’s, [45] proposed the following fused lasso estimator

[TABLE]

If one is not interested in pursuing the sparsity of $\theta$’s, then a fusion estimator [35] can be used

[TABLE]

for some tuning parameter $\lambda$. The above objective functions are both convex and fast computation algorithms are developed, e.g. [29, 45]. The penalty term $\lambda\sum_{i=2}^{n}|\theta_{i}-\theta_{i-1}|$ is can be interpreted as the negative logarithm of prior density used for Bayesian inferences, therefore, a nature Bayesian expansion to (2) is Laplace (double exponential) prior modeling [24, 40]. To account for the unknown variance parameter $\sigma^{2}$ and $\theta_{1}$, a convenient prior specification could be

[TABLE]

where Laplace($a$) denotes the distribution with cdf $f(x)\propto\exp(-a|x|)$.

According to [1], the above Laplace prior $\mbox{Laplace}(\lambda/\sigma)$ can be rewritten as a scale mixture of normal distributions:

[TABLE]

where exp$(a)$ denotes the exponential distribution $f(x)\propto\exp(-ax)$. This hierarchical representation for Laplace prior leads to a Gibbs sampling update that is similar to the Bayesian lasso [34]:

[TABLE]

where $\nu_{i}^{-1}=1/\sigma^{2}+1/\lambda_{i+1}\sigma^{2}+1/\lambda_{i}\sigma^{2}$ and $\mu_{i}=\nu_{i}(y_{i}/\sigma^{2}+\theta_{i+1}/\lambda_{i+1}\sigma^{2}+\theta_{i-1}/\lambda_{i}\sigma^{2})$, $\lambda_{n+1}$ and $\lambda_{0}$ are considered to be infinite, and $\mbox{Inverse-Gaussian}(a,b)$ denotes inverse Gaussian distribution with cdf $f(x)\propto x^{-3/2}\exp[-b(x-a)^{2}/(2a^{2}x)]$.

Despite the popularity of Laplace prior in many applications, recent Bayesian works [7, 5, 42] point out that, if we impose independent $\theta_{i}|\sigma^{2}\sim\mbox{Laplace}(\lambda/\sigma)$ in high dimensional sparse regression under Laplace priors, then the induced posterior has only a sub-optimal contraction rate, or even diverges. In other words, the posterior distribution of $\theta$ doesn’t contract into a small neighborhood around true value $\theta^{*}$ appropriately. For a blocky parameter $\theta$, we observe similar empirical results, as showed in the toy example in Section 2.3: the Laplace fusion prior fails to shrink the observations, which belongs to the same block, towards a same value. Hence, the resultant Bayesian estimate of $\theta$ doesn’t have a step-wise pattern at all.

Following the theoretical discovery of [42], we consider using a class of heavy tailed priors for the successive differences $\vartheta_{i}$’s. Specifically, this work will assign an $t$-shrinkage prior:

[TABLE]

where $t_{a}(b)$ denotes $t$-distribution with degree of freedom $a$ and scale parameter $b$. Note that the above $t$ distribution can be rewritten as an inverse-gamma scaled Gaussian mixture as

[TABLE]

where $a_{t}$, $b_{t}$ satisfy $df=2a_{t}$ and $s=\sqrt{b_{t}/a_{t}}$. Under this $t$-prior, the posterior distribution still allows a full conditional Gibbs sampler, where the update for $\theta_{i}$’s and $\sigma^{2}$ are exactly the same as in (4) and the update of $\lambda_{i}$’s follows

[TABLE]

To further understand the difference between the Laplace fusion prior and $t$ fusion prior, we compare their conditional prior $\pi(\theta_{i}|\theta_{i-1},\theta_{i+1},\sigma)$. Figure 1 plots the function $-\log[\pi(\theta_{i}|\theta_{i-1}=-1,\theta_{i+1}=1,\sigma=1)]$, up to a constant, for both prior specifications. It is clear to see that the conditional $t$-fusion prior allocates most of its prior mass at the two small neighborhoods centered at $\theta_{i-1}$ and $\theta_{i+1}$, given a sufficiently small scale parameter $s$. In other word, the prior introduces a strong shrinkage effect on $\theta_{i}$, towards either $\theta_{i-1}$ or $\theta_{i+1}$. Therefore, for all $i=2,\dots,n$, $\theta_{i}$ will merge with either $\theta_{i-1}$ or $\theta_{i+1}$ in the posterior distribution, which thereafter induces a posterior blocking structure. On the other hand, the conditional Laplace-fusion prior has a uniform prior density within the interval $[\theta_{i-1},\theta_{i+1}]$. Hence, it doesn’t encourage the posterior of $\theta_{i}$ to be grouped with either $\theta_{i-1}$ or $\theta_{i+1}$. It is worth to mention that the NEG-fusion prior [40] also has a similar plot pattern for its conditional prior density function. But a critical difference between $t$ prior and NEG prior, is that the density for NEG prior is non-differentiable at 0. Therefore, the MAP of NEG-fusion prior possesses an exact blocking structure, and [40] only use this MAP for Bayesian inferences rather than the whole posterior distributional information. But the $t$ prior is continuous differentiable everywhere (the functions displayed in the upper plot of Figure 1 is actually smooth at -1 and 1). Hence, its MAP doesn’t have blocky structure, and in this work we will utilize all the posterior samples for the Bayesian analysis.

2.2 Posterior contraction of Bayesian $t$-fusion

In this section, we study the theoretical performance induced by the $t$-fusion prior specification (5). Our theoretical investigation follows the framework of [42], which studies the posterior convergence rate of coefficient $\beta$ in high dimensional sparse regression models $y=X\beta+\varepsilon$. Note the model (1) can also be represented as a sparse linear regression, where the design matrix $X$ is a $n$ by $n$ matrix whose lower triangle entries are all 1 and $\beta=(\theta_{1},\vartheta_{2},\dots,\vartheta_{n})$ is a unknown sparse vector. The following theorem studies the general posterior convergence properties given an independent prior over all $\vartheta_{i}$’s.

Theorem 2.1** (Posterior consistency).**

Assume that $|G^{*}|\prec n/\log(n)$, and prior specification follows that $\sigma^{2}\sim\mbox{Inverse-Gamma}(a_{\sigma},b_{\sigma})$, $\theta_{1}$ and $\vartheta_{i}$’s are conditionally independent given $\sigma^{2}$ with prior density $\pi(\theta_{1},\vartheta_{i}^{\prime}s|\sigma)\propto(1/\sigma)^{n}f_{\theta}(\theta_{1}/\sigma)\prod_{i=2}^{n}f_{\vartheta}(\vartheta_{i}/\sigma)$. Furthermore, if

[TABLE]

then there exist a constant $M$, and $\epsilon_{n}\asymp\sqrt{|G^{*}|\log n/n}$, such that the posterior distribution satisfies

[TABLE]

where the convergence holds in probability or in $L_{1}$ w.r.t. the probability measure of $y$.

The proof closely follows the Theorem A.1 in [42], and for the sake of readability, the proof is provided in the Appendix. The first inequality of sufficient condition set (6) requires that the prior imposed on $\vartheta_{i}$’s is highly concentrated around zero, such that it induces sufficient Bayesian shrinkage effect for those $\vartheta_{i}$’s whose true value is 0. The rest inequalities of (6) essentially require that the prior density at true parameter is at least of order $e^{-cn\epsilon_{n}^{2}}$ for some $c$, and this helps to prevent over-shrinkage for those $\vartheta_{i}$’s whose true value is not 0. Similar conditions, which need the prior to be “thick” at true parameter values, are regularly used in Bayesian literature [17, 23, 12, 13]. Given the concrete forms for prior density $f_{\theta}$ and $f_{\vartheta}$, the second and third inequalities of (6) are equivalent to some upper bound constraints on the magnitude of $\theta_{1}^{*}$ and $\max_{i}|\vartheta_{i}^{*}|$’s (see e.g., Corollary 2.1). The fourth inequality of (6) trivially holds for any fixed $a_{\sigma}$ and $b_{\sigma}$ if $\sigma^{*2}$ is assumed to be a constant. If the unknown error variance is supposed to be varying with respect to $n$, e.g., the studies of Gaussian sequence models [19] commonly assume that $\sigma^{*2}\propto n^{-1}$, then one can choose a fixed $a_{\sigma}$ and $b_{\sigma}\asymp n^{-\kappa}$ for some $\kappa>0$. Under such a choice, condition (6) holds as long as $-k\log n\leq\log(\sigma^{*2})<K\log n$ for some positive constant $K$.

The above result states that almost all the posterior mass contract into a neighborhood of $\theta^{*}$ with a radius $M\sigma^{*}\sqrt{n}\epsilon_{n}$, that is, the posterior convergence rate is of an order $\sigma^{*}\sqrt{|G^{*}|\log n}$. Note that if the partition index set $G^{*}$ were known, the oracle rate of contraction turns out to be $O(\sigma^{*}\sqrt{|G^{*}|})$. Hence, the Bayesian shrinkage achieves the ideal risk up to a logarithmic term in $n$. In frequentist literature, Theorem 2.7 of [35] showed that the convergence rate of Fused Lasso is no larger than $O(\sigma^{*}\sqrt{|G^{*}|\log|G^{*}|})$. However, this rate is not directly comparable since an additional minimal signal strength condition, which ensures that $G^{*}$ can be fully recovered in probability, is imposed in this paper.

The posterior distribution of $\vartheta$ are always continuous, and doesn’t directly provide Bayesian inferences for the block structure, or equivalently, the unknown $G^{*}$. The following result characterizes the some asymptotic performance of posterior block partition via discretization.

Theorem 2.2** (Posterior selection).**

Assume the conditions of Theorem 2.1 hold, denote $G(\theta,\sigma)=\{i:|\vartheta_{i}/\sigma|<\epsilon_{n}/n\}$, then the posterior of $G(\theta,\sigma)$ satisfies

[TABLE]

for some fixed constant $\delta$, where the convergence holds in probability and in $L_{1}$.

Therefore, the posterior distribution of $G(\theta,\sigma)$, which is induced by the posterior of ($\theta,\sigma^{2})$ and mapping $G(\cdot,\cdot)$, can be treated as a discrete posterior distribution for the unknown $G^{*}$ and used for Bayesian block partition selection. Theorem 2.2 essentially claims that the number false positive selections of via discretization $G(\cdot,\cdot)$ is bounded in posterior probability. Such a result is comparable to the model selection behavior of Bayesian Dirichlet-Laplace shrinkage [5]. It is worth to note that if certain minimal signal strength condition holds as well, i.e., $\min_{\vartheta_{i}^{*}>0}|\vartheta_{i}^{*}|$ is bounded away from zero, with addition assumptions, one can derive an even stronger posterior selection consistency result that $\pi(G(\theta,\sigma)=G^{*}|y)\rightarrow^{p}1$, following the proof of Theorem 2.3 in [42]. Readers of interests can easily derive such posterior selection consistency by themselves.

It is not difficult to verify the imposed conditions for the proposed $t$ fusion shrinkage prior (5).

Corollary 2.1**.**

When $f_{\theta}$ is the cdf of normal distribution $N(0,\lambda_{1})$ and $f_{\vartheta}$ is the cdf of a $t$-distribution with scale parameter $s$, then the first three inequalities of (6) hold when

[TABLE]

and $s=n^{-c}$ for some sufficiently large $c$.

The above results hold not only for the scaled $t$ prior, but also for any other choice of $f_{\vartheta}$, as long as $f_{\vartheta}$ has a polynomially decaying tail. Note that the above results can not be generalized to light tailed distributions such as Laplace-fusion prior, since it will lead to an unrealistic sufficient condition that $\max_{i}|\vartheta_{i}^{*}|/\sigma^{*}=o(1)$.

2.3 Bayesian posterior inference

In this section, we will illustrate the posterior behavior of Bayesian $t$-fusion by a toy example, and compare it to Bayesian Laplace-fusion. We will also discuss other issues related to hyperparameter choice.

A simulation data was generated with $n=100$ and $\sigma^{*}=0.5$. The data and underlying true parameter value are plotted in Figure 2. Three estimation procedures are considered: 1) $L_{1}$ fusion estimation (2) where the tuning parameter are selected by cross validation; 2) Bayesian Laplace prior (3) with $\lambda=\sqrt{2\log(n)}$ and $\lambda_{1}$=5; 3) Bayesian t prior (5) with $a_{t}=2$, $b_{t}=0.005$, $a_{\sigma}=0.5$, $b_{\sigma}=0.5$ and $\lambda_{1}=5$. The posterior samples are obtained by Gibbs sampler for 2000 iterations. The frequentist estimator and boxplots of Bayesian marginal posterior distributions are displayed in Figure 2. Note that for the sake of readability of the figure, we don’t draw the outliers in these boxplots.

The comparison shows that the $L_{1}$ fusion penalty leads to a strictly sparse $\widehat{\vartheta}$ estimation and the estimated $\widehat{\theta}$ does have a blocking structure, but there is a mild over-partition issue. Due to the nature of Bayesian shrinkage prior, the two Bayesian approaches only produce continuous posterior samples of $\vartheta_{i}$. Comparing with the Laplace prior, the advantages of $t$ prior are quite obvious. Its posterior concentration is better, i.e., shorter boxes and whiskers for the boxplots, and posterior mean is much more close to the true red curve. Although its posterior estimation is not exactly blocky, one can visually identify the three blocks. In contrast, the Laplace prior generates larger posterior variance, and its posterior center smoothly fluctuates around the true step function. As discussed at the end of Section 2.1, Laplace fusion prior itself merely has shrinkage effect for the successive differences, and hence the data fluctuation carries over to the posterior distribution of $\theta$. In this case, it is not clear to recover the underlying block structure.

Note that for both (3) and (5) prior specifications, their posterior distributions highly depend on the choice of hyperparameter, especially the scale parameters of the Laplace (parameter $1/\lambda$) and $t$ prior (parameter $s$) distribution. To understand the influence of the scale parameter, we increase and decrease the scale parameter by a factor of $\sqrt{10}$, for both Laplace and $t$ priors. The corresponding posteriors are displayed in Figure 3.

From Figure 3 and its comparison with Figure 2, we see that, for $t$ prior specification, smaller scale hyperparameter leads to smaller posterior variation, and vise versa. Choosing an overly large scale parameter weakens the shrinkage effect, thus it fails to shrink the $\theta_{i}$’s that belong to the same block towards same value, and the corresponding posterior estimation is not flat within blocks. On the other hand, choosing an overly small scale hyperparameter, although yields very strong shrinking and grouping effect, may potentially over-partition the data. In conclusion, the scale parameter of $t$ fusion prior controls the aggressiveness of posterior block partition. For Laplace prior, both increasing or decreasing scale parameter can not remedy its poor behavior demonstrated in Figure 2. The choice of scale parameter only affects the smoothness of posterior of $\theta$. Smaller scale parameter leads to smoother non-parametric estimation which is similar to smooth spline regression; larger scale parameter leads to a rather rugged estimation. This phenomenon can be partially explained by Figure 1. As the scale of the Laplace prior decreases (i.e., $\lambda$ increases), the conditional prior of $\theta_{i}$ given $(\theta_{i-1},\theta_{i+1})$ acts more and more like a uniform distribution over $[\theta_{i-1},\theta_{i+1}]$. Hence, conditional on $(\theta_{i-1},\theta_{i+1})$, the posterior of “local” total variation, $|\theta_{i+1}-\theta_{i}|+|\theta_{i}-\theta_{i-1}|$, always decreases to the minimum value $|\theta_{i+1}-\theta_{i-1}|$ when the scale parameter decreases to zero, regardless of the observation value $y_{i}$.

As showed in the previous toy example, the hyperparameter value plays an important role for the performance of Bayesian posterior inferences. The theoretical suggestion, i.e., the first inequality of (6), gives a very small scale parameter, which in practice leads to severe over-partition issue under moderate $n$. Unlike frequentist high dimensional penalization estimation whose tuning parameter is usually determined by cross validation or selection criterion such as EBIC [9, 10], choosing a proper prior hyperparameter posts many difficulty for Bayesian statisticians. Conventional choices include imposing a hyper-prior on the hyperparameter [48, 8] or empirical Bayes [36]. For high dimensional sparse GLMs, [28] suggested choosing a hyperparameter such that the posterior mean and posterior mode are close. Several Bayesian works [40, 4] considered only using the sparse MAP as the Bayesian estimator, and thus abandon the whole posterior distributional information. This strategy somehow reduces the Bayesian computation to a frequentist optimization problem, and then the hyperparameter can be determined by EBIC criterion. It is beyond the scope of this work to theoretically study how to select an appropriate hypereparameter, i.e., the scale parameter of the $t$ prior. An empirical suggestion based on authors’ experience is to choose the scale parameter of (5) such that

[TABLE]

Although it doesn’t quite satisfy the conditions in Corollary 2.1 which suggests $s=n^{-c}$ for some large $c$, but in practice, it indeed yields a reasonable and stable Bayesian performance. Note that our prior choice $a_{t}=2$, $b_{t}=0.005$ approximately satisfies (7).

3 Bayesian $t$-shrinkage for Bayesian clustering

In this section, we would like to extend the applications of Bayesian $t$ fusion shrinkage to Bayesian clustering problem, where the parameter $\theta$ in model (1) is assumed to follow an unknown clustering structure. In other words, the observations $y_{i}$’s are not organized in a sensible order, and we no longer assume that the true cluster only consists of consecutive indexes.

A full Bayesian clustering analysis, or subgroup identification, usually imposes a prior on the clustering structure, including the number of clusters and how to partition observations into these clusters. For example, one could consider that the cluster structure arises from a tree splitting process, and impose certain prior on the tree nodes [3, 15]. The resultant posterior sampling hence usually requires an inefficiently reversible-jump Metropolis-Hastings move to travel across different clustering models in the sampling space. Another common approach is to model a mixture distribution via nonparametric priors such as Dirichlet process or Chinese restaurant process [33]. These mentioned approaches enforce explicit clustery posterior samples by directly applying discrete prior on subgroup structure. In this section, we will show how to use Bayesian $t$-shrinkage to induce implicit posterior clustering structure.

To implement shrinkage prior, we need to formulate the problem in a proper statistical modeling such that its parameter possesses certain sparsity feature. Conceptually, imposing shrinkage priors on $\theta_{i}-\theta_{i-1}$ shall not work for this clustering problem since $\{\theta_{i}-\theta_{i-1}\}_{i=2}^{n}$ is no longer a sparse vector. However, if some meaningful permutation $r$ of $\theta_{i}$ is known such that $\{\theta_{r(i)}-\theta_{r(i-1)}\}_{i=2}^{n}$ is indeed a sparse vector, for example, $\theta_{r(i)}$ is ascending or descending, then the problem will reduce to Bayesian shrinkage fusion studied in the previous section.

However, in practice, such $r(\cdot)$ is generally not available. Thus one simple solution would be to substitute $r(\cdot)$ with some estimator $\widehat{r}()$. For model (1), a trivial estimator is the rank statistic of the responses $y_{i}$, i.e., $y_{\widehat{r}(i)}=y_{(i)}$ where $y_{(i)}$ denotes the order statistics of $y_{i}$’s. Thus this leads to hybrid Bayesian modeling, depending on a frequentist pilot estimator $\widehat{r}$:

[TABLE]

The corresponding posterior Gibbs sampler follows:

[TABLE]

where $\nu_{i}^{-1}=1/\sigma^{2}+1/\lambda_{i+1}\sigma^{2}+1/\lambda_{i}\sigma^{2}$ and $\mu_{i}=\nu_{i}(y_{\widehat{r}(i)}/\sigma^{2}+\theta_{\widehat{r}(i+1)}/\lambda_{i+1}\sigma^{2}+\theta_{\widehat{r}(i-1)}/\lambda_{i}\sigma^{2})$. The idea of taking advantage of a pilot estimator $\widehat{\theta}$ and its rank statistic ${\widehat{r}(i)}$ has also been implemented in the two-stage frequentist approach [20].

It is worth to mention that posterior consistency Theorems 2.1 and 2.2 for Bayesian shrinkage fusion don’t apply to the above Bayesian modeling (8) even when $\theta_{\widehat{r}(i)}$ is in a ascending order, because $\varepsilon_{\widehat{r}(i)}$ are no longer iid error observations. Such a data dependent prior (8) can also be interpreted as a misspecified Bayesian modeling [22], where the data is the order statistics $y_{(i)}$ from a mixture distribution, and we model $y_{(i)}$’s as independent Gaussian variables whose mean function $E(y_{(i)})$ is a step function.

Although the prior modeling (8) is quite natural and intuitive, there are several issues. First, comparing with the iid observations, it is more difficult to to identify the underlying clustering structure of sorted observations. For example, it is much more difficult to visually identify the 3-cluster structure in the toy data showed in Figure 4 than the toy data showed in Figure 2. Secondly, (8) fails to take account of uncertainty of the estimator $\widehat{r}$, and the strict monotonicity of $y_{\widehat{r}(i)}$ will always carry over to the posterior of $\theta_{\widehat{r}(i)}$. Thirdly, over-clustering will occur. To understand this, let us consider that all data $y_{i}$’s have the same mean, i.e., there is only one cluster. But the mean function of the sorted data, $E(y_{(i)})$, is actually a strictly increasing function. From the perspective of Bayesian misspecification, our posterior of $\theta_{\widehat{r}(i)}$’s should contract towards the best step function, in the sense of minimizing the Kullback–Leibler divergence [22], rather than towards the constant function $E(y_{i})$. This causes the over-clustering, and furthermore, the posterior consistency of $t$ shrinkage prior established in Corollary 2.1 doesn’t hold anymore.

Theorem 3.1**.**

Assume that all $y_{i}\sim N(0,1)$ with known variance $\sigma^{*}=1$. If the prior (8) is used (except that there is no necessity to impose a prior on $\sigma^{2}$) and the scale parameter satisfies $-\log(b_{t})\asymp\log n$, then in high probability

[TABLE]

i.e., the rate-$\sqrt{\log n}$ posterior consistency fails.

The negative result in Theorem 3.1 motivates us to propose an adaptive pseudo Bayesian shrinkage approach. Instead of using a fixed estimator $\widehat{r}$, we update the “working” order $r$ of $\theta$ over iterations. To be specific, we modify the above Gibbs sampling iteration as:

[TABLE]

In other words, we update $r(\cdot)$ to the rank statistic of the newly obtained $\theta$.

There are a couple of rationales behind the algorithm (10). First, the update of $r$ potentially allows algorithm (10) to incorporate the uncertainty of rank statistic into the posterior sampling. Heuristically, if the current $r(\cdot)$ is close to the true ranking, (10) will shuffle the ordering of $\theta$ within each cluster, hence $y_{r(i)}$ acts like independent samples yielded by random permutation within cluster, rather than order statistics. Secondly, the sampling algorithm (10) can be somehow connected with a full Bayesian modeling:

[TABLE]

Under (11), the conditional posterior of $r(\cdot)$ follows $\pi(r|\theta,\lambda_{i}^{\prime}s,y)\propto\exp\{-\sum_{i=2}^{n}(\theta_{r(i)}-\theta_{r(i-1)})^{2}/\lambda_{i}^{2}\sigma^{2}\}$, from which sampling is expensive. However, when most of $\lambda_{i}^{2}$s are tiny, the distribution of $\pi(r|\theta,\lambda_{i}^{\prime}s,y)$ will be highly concentrated around its MAP which is approximately the rank statistic of $\theta$. Therefore, updating $r$ to the rank statistic of current sample $\theta$, as in (10), can be viewed as a convenient sampling of $\pi(r|\theta,\lambda_{i}^{\prime}s,y)$ under the uniform hyperprior of $r$. And such a hyperprior does benefit the posterior asymptotic performance, at least it remedies the posterior inconsistency described in Theorem 3.1 when there is only one underlying cluster.

Theorem 3.2**.**

Assume that all $y_{i}\sim N(\theta_{0},\sigma^{2})$, i.e. $\theta^{*}$ is a vector of $\theta_{0}$’s. If prior (11) is used with $b_{t}=n^{-c}$ for some sufficiently large $c$, and $\theta_{0}^{2}/(\lambda_{1}\sigma^{*2})+\log(\lambda_{1})=O(\log n)$ then there exist constant $M$ and $\epsilon_{n}\asymp\sqrt{\log n/n}$, such that

[TABLE]

where the convergence holds in probability and $L_{1}$.

Further investigations will be pursued to assess the posterior contraction induced by (11) when the true parameter $\theta^{*}$ has multiple-cluster structure. The minimal cluster difference $\min_{\{(i,j):\theta_{i}\neq\theta_{j}\}}|\theta_{i}-\theta_{j}|$ shall play an crucial role for the Bayesian asymptotics. If the minimal cluster difference is bounded by constant, then even the best Bayes classifier makes as many as $O(n)$ misclassifications, which leads to $\sqrt{n}$-rate $L_{2}$ error. And we conjecture that when the minimal cluster difference is large than $\sqrt{M\log n}$ for some constant $M$, the posterior distribution induced by (11) leads to satisfactory Bayesian consistency.

It is worth to mention that the $r$-update step in algorithm (10) may critically change the stochastic stability of the algorithm and that the ergodicity of Markov chains generated by (10) is still unclear. Therefore, in our application we use it with caution. In all our simulations and toy examples, we initialize with $r=\widehat{r}$, and the $r$-update step is only implemented every other 20 iterations after certain burn-in period.

In frequentist literature, another popular approach to cluster $\theta_{i}$’s is to impose a pairwise difference penalty as $\sum_{i\neq j}p(|\theta_{i}-\theta_{j}|)$ for some penalty $p(\cdot)$[30, 39, 21]. Although it is tempting to develop a Bayesian counterpart, i.e., using a prior of form $\pi(\theta)\propto\prod_{i\neq j}\pi(\theta_{i}-\theta_{j})$ with sparsity induced $\pi$, there are several problems. First, the pairwise differences $\{\theta_{i}-\theta_{j}\}_{i\neq j}$ must satisfy triangle inequality, hence the prior specification $\pi(\theta)\propto\prod_{i\neq j}\pi(\theta_{i}-\theta_{j})$ is actually not an independent prior over all pairwise difference $\{\theta_{i}-\theta_{j}\}_{i\neq j}$, and it will be difficult to characterize the effect of such dependency in the prior distribution. Second, even under clustering structure, $\{\theta_{i}-\theta_{j}\}_{i\neq j}$ are actually not sparse. For example, if there are two balanced groups among $\theta_{i}$’s, then half of the pairwise differences are nonzero. Imposing a sparse prior on a non-sparse system will lead to severe overshrinkage problem, and our simulation shows that such prior shrinks all $\theta_{i}$ towards $\bar{y}$. Third, such pairwise difference prior will also heavily increase the computational burden of posterior sampling.

To illustrate and compare the performance of (9) and (10), a simple toy example is conducted. We simulate a data set with $n=100$ data points that belong to 3 underlying subgroups. For both (9) and (10), we choose $a_{t}=2$, $b_{t}=0.005$, $a_{\sigma}=b_{\sigma}=0.5$ and $\lambda_{1}=5$. In addition, we compare them with frequentist $L_{1}$ fusion

[TABLE]

and Bayesian Dirichlet process modeling

[TABLE]

with $\lambda=5$, $a=0.1$ and $a_{\sigma}=b_{\sigma}=0.5$. The posterior sampling algorithms of Dirichlet process modeling are discussed in [33]. The frequentist estimator and posterior boxplots of $\theta_{i}$s under different priors are plotted in Figure 4. This figure clearly shows that frequentist $L_{1}$ fusion estimator fails to yield clustering structure for monotone sorted $y_{\widehat{r}(i)}$’s, especially for the middle portion of the data, one always have $\widehat{\theta}_{\widehat{r}(i)}=y_{(i)}$. Dirichlet process induces a quite reasonable posterior clustering, however its posterior concentration is not good, in the sense that the posterior variance is quite large. In the opposite, the $t$ modeling (8) has a strong prior concentration, which contributes to a better posterior variation, and encourages posterior clustering. However, consistent to our previous arguments, a severe over-clustering occurs, and the posterior identifies more than 6 clusters. At last, the posterior obtained by algorithm (10) combines both advantages of DP prior and (8). Comparing with DP prior, it has much smaller posterior variance; comparing with prior (8), the over-clustering issue is greatly alleviated.

4 Simulation and data anlaysis

4.1 Bayesian fusion simulations

In our first simulation studies, we consider model (1) with $\sigma^{*}=0.5$ and true parameter $\theta^{*}\in\mathbb{R}^{100}$ having three consecutive blocks. These three blocks are ${\cal B}_{1}^{*}=\{1,\dots,b_{1}\}$, ${\cal B}_{2}^{*}=\{b_{1}+1,\dots,b_{1}+b_{2}\}$, ${\cal B}_{3}^{*}=\{b_{1}+b+2+1,\dots,b_{1}+b_{2}+b_{3}\}$ where $(b_{1},b_{2},b_{3})\sim\mbox{multinomial}(100,(1/3,1/3,1/3))$, and $\theta_{i}^{*}=0$, 2 or 4 within each block respectively. We compare the performance among Bayesian $t$ fusion (5), Bayesian Laplace fusion (3) and frequentist $L_{1}$ fusion (2), based on 100 replications. The choices of the hyperparameters are same to the toy example discussed in Section 2.3.

To compare the accuracy of parameter estimations, we calculate the $L_{2}$ and $L_{1}$ estimation errors, $\|\widehat{\theta}-\theta^{*}\|_{2}$ and $\|\widehat{\theta}-\theta^{*}\|_{1}$, for frequentist and Bayes estimator (posterior mean). To assess the performance of block structure recovery, we consider several comparison criteria. We define the “adjacency” matrix as $\Sigma=(\omega_{ij})=(1\{\theta_{i}=\theta_{j}\})$, and corresponding estimation error of $\Sigma$: $R=\|\Sigma^{*}-\widehat{\Sigma}\|_{1}=\sum_{i,j}|\widehat{\omega}_{ij}-\omega_{ij}^{*}|$. For $L_{1}$ fusion estimator, $\widehat{\Sigma}$ is trivially estimated by $\widehat{\omega}_{ij}=1\{\widehat{\theta}^{\rm F}_{i}=\widehat{\theta}^{\rm F}_{j}\}$. For Bayesian shrinkage approaches, since their posterior samples are continuous, it is necessary to “sparsify” the continuous posterior in order to retrieve a sparse structure estimation for $\Omega$. In the literature, this is usually done by 1) hard thresholding [27, 43, 6, 16], or 2) decoupling shrinkage and selection methods [14, 50]. All these mentioned approaches depend solely on the posterior scaling, and donot take into account of the degree of prior concentration. Therefore, following the suggestion made in [42], we estimate $\widehat{\omega}_{ij}=1\{|\widehat{\theta}_{i}-\widehat{\theta}_{j}|/\widehat{\sigma}\leq\pi_{1/2n}\}$, where $\widehat{\theta}_{i}$ and $\widehat{\sigma}$ are the Bayes estimator, $\pi_{1/2n}$ denotes the $(1-1/2n)$ quantile of the prior distribution imposed on successive difference $\vartheta_{i}/\sigma$. This choice tries to increases the robustness of Bayesian inference: if an overly small or large scale parameter $s$ is used, the estimation for $\omega_{i,j}$ will adapt accordingly. For $t$ shrinkage (5), $\pi_{1/2n}=st_{1/2n}$ where $s$ is the scale parameter, $t_{1/2n}$ is the $1-1/2n$ quantile of $t$ distribution with df$=2a_{t}$; for Laplace shrinkage (3), $\pi_{1/2n}=\log(n)/\lambda$. Since different priors have different corresponding $\pi_{1/2n}$ values, and comparison solely based on the value of $R$ may not be completely fair. Hence, we also consider the following two statistics: $W=\mbox{Average}_{\{w_{ij}^{*}\neq 0\}}|\widehat{\theta}_{i}-\widehat{\theta}_{j}|$ denotes within-block average variation, and $B=\min_{\{w_{ij}^{*}=0\}}|\widehat{\theta}_{i}-\widehat{\theta}_{j}|$ denotes the between-block separation. A larger $B$ and smaller $W$ indicate a better block identification performance. The results are summarized in Table 1.

The simulation results show that the Bayesian $t$ fusion yields the most accurate estimation in terms of both $L_{2}$ and $L_{1}$ errors. This Bayes estimator also induces the best clustering result, as it has the largest between-group separation and smallest within-group variation. In comparison, the frequentist fusion estimator has a worse accuracy. It is well known that $L_{1}$ penalty introduces estimation bias [11], but unlike the LASSO penalty which introduces a bias of $\lambda$, fuse-$L_{1}$ penalty introduces only a bias as small as $\lambda/n_{i}$ ([35]) where $n_{i}$ is the number of elements in each block. Hence, its estimation performance is still acceptable. The $L_{1}$ fusion also has a much larger $R$ statistics. This is consistent to observation from our toy example, that $L_{1}$ fusion has a mild over-partition issue. Bayesian Laplace fusion, on the other hand, has a much worse estimation behavior due to its insufficient shrinkage effect. As discussed in Section 2.3, the Laplace fusion tends to yield smooth changing $\widehat{\theta}_{i}$’s, thus it has a much smaller between-block separation. Note that Laplace shrinkage does obtain the smallest $R$ statistic, but this doesn’t imply that it has a good performance on structure recovery. It has a small $R$ statistics because the Laplace-fusion prior doesn’t have strong prior concentration and its corresponding $\pi_{1/2n}$ is much larger than $t$ fusion.

4.2 Bayesian clustering simulations

In our second simulation studies, we consider model (1) with $\sigma^{*}=0.5$ and true parameter $\theta^{*}\in\mathbb{R}^{100}$ having three unknown clusters. With equal probability, $\theta_{i}^{*}=0$, 2, or 4. We aim to compare different approaches including Bayesian $t$ modeling (9) and (10), Dirichlet process prior (13) and the frequentist bCARDS estimation (12) using $L_{1}$ fusion [20]. Besides the comparison of $L_{2}$ and $L_{1}$ errors of the Bayesian/frequentist estimator, we also compare the posterior mean of squared $L_{2}$ error, $E_{\pi(\theta|y)}\|\theta-\theta^{*}\|^{2}$, among the Bayesian approaches. Note that the $L_{2}$ error of Bayes estimator only tells how good is the posterior center, while the posterior mean of squared $L_{2}$ error tells how good is distributional posterior contraction.

Note that the simulation setting seems similar to our previous study in Section 4.1, but as mentioned in Section 3, clustering problem is much more difficult than fusion problem. To see that, assume there are three balanced clusters for $\theta$, let $z_{1}$ be the largest observation in the 0’s cluster and $z_{2}$ be the smallest observation in the 2’s cluster. By the well known result [37], $E(z_{1})\approx 1.04>0.96\approx E(z_{2})$. Therefore, misclassification will always happen for those extreme observations, and $B=\min_{\{w_{ij}^{*}=0\}}|\widehat{\theta}_{i}-\widehat{\theta}_{j}|$ is always around 0 for all methods, therefore, the comparison of $B$ statistic values are meaningless. Because of that, to assess how well these methods identify and separate unknown clusters, we re-define the $B$ statistic as

[TABLE]

The simulation results are summarized in Table 2.

The comparison shows that DP prior yields the best point estimator in terms of estimation error, and the adaptive $t$-shrinkage method (10) gives a slightly worse point estimator. But in terms of the posterior contraction, both $t$-shrinkage approaches have much smaller posterior mean of squared $L_{2}$ error. This result is consistent to the insight obtained from the toy example in Section 3 that LP prior leads to larger posterior variance. The three Bayesian approaches have approximately the same performance for between-cluster separation ($\tilde{B}$ statistic), while the DP prior has a smaller within-cluster variation ($W$ statistic). The frequentist estimator (12) has the worst performance in almost every aspect. In summary, DP prior does yield the best Bayes point estimation, but the adaptive $t$-shrinkage (10) has the best Bayesian contraction, and its posterior mean, as a point estimator, has reasonable performance for estimation and clustering structure recovery. Furthermore, by comparing the results between $t$-shrinkage method (10) and (9), we conclude that the $r$-update step implemented in algorithm (10) does improve the performance of $t$-shrinkage in every aspect for this clustering problem.

4.3 Real data set analysis

In this section, we consider a real comparative genomic hybridization (CGH) dataset [46]. The dataset is available from the R package cghFLasso and it contains n=990 observations. We are interested in a fusion estimation for the mean trend in order to detect the “hot-spot” region, i.e., the corresponding genes have extra DNA copies. We apply both Bayesian $t$-shrinkage fusion (5) and frequentist $L_{1}$ fusion (2). For the Bayesian application, the hyperparameter is chosen as $a_{t}=2$ and $b_{t}=0.0007$, following the suggestion of (7). Both point estimations are plotted in Figure 5. Frankly speaking, the two estimations share similar pattern, and it is difficult to judge which one is definitely better. To see the slight advantage of $t$-shrinkage, a closer comparison is displayed in the same figure as well, where we compare the two estimators within the segment [80,120]. It is obvious that both approaches induce exactly the same blocking structure. But the difference is that, for $t$-shrinkage prior, the $\widehat{\theta}_{i}$ within each block is very close to the sample mean of the block; but for $L_{1}$ fusion, there is a large bias between $\widehat{\theta}_{i}$ and the sample mean of the block.

5 Final remarks

In this work, we study the Bayesian inference for a vector parameter $\theta$ which has a unknown blocking structure, using Bayesian shrinkage prior. For the ease of representation, this work focuses on the application of $t$-fusion prior, but our theorems actually holds for any polynomially decaying prior distribution. We demonstrate that a simple $t$-fusion prior leads to satisfactory posterior contraction, and is powerful to recover the blocking structure. We recommend not to use Laplace fusion prior since it can only induce a smoothly varying posterior estimation. Although this work mainly focuses on the normal sequence models, but the presented $t$-fusion modeling and the posterior asymptotic results can be easily extended to more complicated regression models such as $y=X\beta+Z\theta+\epsilon$ for some blocky parameter $\theta$.

We also extend the use of shrinkage prior to a more general clustering problem. To the best of authors’ knowledge, this is the first attempt in literature to use Bayesian shrinkage to recover unknown clustering structure. The basic idea is to find a pilot order $\widehat{r}$, and then fuse all pairs of neighbors (that are determined by $\widehat{r}$) via shrinkage priors. An adaptive $r$-update step is further incorporated to improve the clustering performance. Simulations show that the proposed algorithms (9) and (10) have reasonable performance, and their posterior contraction is much better than conventional Dirichlet process modeling. Further theoretical investigation is necessary to understand their asymptotic behaviors. In practice, this idea can be generalized to clustering problem of multi-dimensional data as well, i.e., $y_{i}=\theta_{i}\in\mathbb{R}^{d}$. Since we can not rank multivariate vectors, the possible alternative is to construct a minimum spanning tree (MST)[26], and impose shrinkage fusion prior on the pair of $\theta_{i}$ that is connected by a edge in the MST.

Acknowledgement

This work is in memory of Prof. Jayanta Ghosh who has jointly supervised the first PhD student of the second author. Qifan Song’s research is sponsored by NSF DMS-1811812. Guang Cheng’s research is sponsored by NSF DMS-1712907, DMS-1811812, DMS-1821183, and Office of Naval Research, (ONR N00014-18-2759).

Appendix A Appendix

First, let us state some useful lemmas.

Lemma A.1** (Lemma 1 of [25]).**

Let $\chi^{2}_{d}(\kappa)$ be a chi-square distribution with degree of freedom $d$, and noncentral parameter $\kappa$, then we have the following concentration inequality

[TABLE]

Lemma A.2** (Theorem 1 of [54]).**

Let $X$ be a Binomial random variable $X\sim\mbox{B}(n,v)$. For any $1<k<n-1$

[TABLE]

where $\Phi$ is the cumulative distribution function of standard Gaussian distribution and $H(v,k/n)=(k/n)\log(k/nv)+(1-k/n)\log[(1-k/n)/(1-v)]$.

The next lemma is a refined result of Lemma 6 in [2]:

Lemma A.3**.**

Let $f^{*}$ be the true probability density of data generation, $f_{\theta}$ be the likelihood function with parameter $\theta\in\Theta$, and $E^{*}$, $E_{\theta}$ denote the corresponding expectation respectively. Let $B_{n}$ and $C_{n}$ be two subsets in parameter space $\Theta$, and $\phi_{n}$ be some test function satisfying $\phi_{n}(D_{n})\in[0,1]$ for any data $D_{n}$. If $\pi(B_{n})\leq b_{n}$, $E^{*}\phi(D_{n})\leq b_{n}^{\prime}$, $\sup_{\theta\in C_{n}}E_{\theta}(1-\phi(D_{n}))\leq c_{n}$, and furthermore,

[TABLE]

where $m(D_{n})=\int_{\Theta}\pi(\theta)f_{\theta}(D_{n})d\theta$ is the margin probability of $D_{n}$. Then,

[TABLE]

Proof.

Define $\Omega_{n}$ to be the event of $(m(D_{n}))/(f^{*}(D_{n}))\geq a_{n}$, and $m(D_{n},C_{n}\cup B_{n})=\int_{C_{n}\cup B_{n}}\pi(\theta)f_{\theta}(D_{n})d\theta$. Then

[TABLE]

By Fubini theorem,

[TABLE]

Combining the above inequalities leads to the conclusion. ∎

Proof of Theorem 2.1 and 2.2 Let $G=\{g_{1},g_{2},\dots,g_{d}\}$ be a generic subset of $\{2,\dots,n\}$, and it also represents potential a $(d+1)$-group structure of $\theta$ as $\{\{1,\dots,g_{1}\},\{g_{1}+1,\dots,g_{2}\},\dots,\{g_{d}+1,\dots,n\}\}$. Given $G$ and its corresponding blocking structure, $\widehat{\theta}_{G}(y)$ denotes the estimation of $\theta$ based on block mean, i.e. $\widehat{\theta}_{G,j}(y)=\sum_{i=g_{j}+1}^{g_{j+1}}y_{i}/(g_{j+1}-g_{j})$ for all $g_{j}+1\leq j\leq g_{j+1}$, and $\widehat{\sigma}^{2}_{G}(y)=\|y-\widehat{\theta}_{G}\|^{2}/(n-|G|-1)$.

To prove the posterior contraction, we will apply Lemma A.3. We define the following testing function

[TABLE]

for some $\delta>0$, and define $C_{n}$ and $B_{n}$ as:

[TABLE]

Note that when $G\supset G^{*}$, $\|\widehat{\theta}_{G}(y)-\theta^{*}\|^{2}\sim\sigma^{*2}\chi^{2}_{|G|+1}$ and $\|y-\widehat{\theta}_{G}\|^{2}\sim\sigma^{*2}\chi^{2}_{n-|G|-1}$, thus by Lemma A.1, we have that

[TABLE]

for some constant $c_{1}$, since $|G|=O(|G^{*}|)\prec n\epsilon_{n}^{2}$ and $\epsilon_{n}\prec 1$. Therefore,

[TABLE]

as long as $n\epsilon^{2}/[|G^{*}|\log n]$ is sufficiently large.

For any $(\theta,\sigma^{2})\in C_{n}$ satisfying $\|\theta-\theta^{*}\|\leq M\sqrt{n}\sigma^{*}\epsilon_{n}$ and $\sigma^{2}/\sigma^{*2}\geq(1-\epsilon_{n})/(1+\epsilon_{n})$, we define $\widehat{G}=\{i:\theta_{i}-\theta_{i-1}\geq\sigma\epsilon_{n}/n^{2}\}\cup G^{*}$ ( hence $|G|\leq(1+\delta)|G^{*}|$), thus

[TABLE]

for some $c_{2}$ given a large $M$, where the second inequality is due to the fact that $\|\widehat{\theta}_{\widehat{G}}(\theta)-\theta\|\leq\sqrt{n}\sigma\epsilon_{n}$ when $\theta\in B_{n}$.

For any $(\theta,\sigma^{2})\in C_{n}$ satisfying $\sigma^{2}/\sigma^{*2}<(1-\epsilon_{n})/(1+\epsilon_{n})$ or $\sigma^{2}/\sigma^{*2}>(1+\epsilon_{n})/(1-\epsilon_{n})$,

[TABLE]

for some $c_{2}^{\prime}$, where the noncentral parameter $\lambda<n\epsilon_{n}^{2}\prec(n-|G|-1)\epsilon_{n}$.

Combining the results from the previous two paragraph, it is easy to obtain that

[TABLE]

Now we consider the marginal posterior density of data $y$. With probability $P(\|\varepsilon\|\leq 2\sqrt{n}\sigma^{*})$ (which converges to 1),

[TABLE]

for some constant $c_{3}^{\prime}$. Besides,

[TABLE]

for some constant $c_{3}^{\prime\prime}$. Thus

[TABLE]

for some $c_{3}$, where $c_{3}$ can be sufficiently small when $n\epsilon_{n}^{2}/[|G^{*}|\log n]$ is large enough.

At last, we study the prior probability of set $B_{n}$. Due to the prior independence of $\vartheta_{i}^{\prime}$s, $\pi(B_{n})=\pi[\mbox{Bin}(n-1-|G^{*}|,p)>(\delta)|G^{*}|]$, where $p\leq(1/n)^{1+u}$. By Lemma A.2,

[TABLE]

for some $c_{4}$. Combine results (15), (16), (17) and (18), and we apply Lemma A.3 to get the posterior consistency result that

[TABLE]

given a sufficient large constants $\delta$ and $n\epsilon_{n}/|G^{*}|\log n$.

Proof of Theorem 3.1

Consider $y=\theta+d$, where $\theta_{i}^{*}\equiv 0$ for all $i$ and error $d$ is order statistics of standard normal variables, i.e. the density of $d$ is $f(d)=n!\prod\phi(d_{i})1(d_{1}\leq d_{2},\dots,d_{n-1}\leq d_{n})$ and $\phi$ denotes the standard normal density. The prior of $\theta$ follows $\pi(\theta)=\pi_{1}(\theta_{1})\prod_{i=2}^{n}\pi_{t,s}(\theta_{i}-\theta_{i-1})$, where $\pi_{\lambda_{1}}$ is the density of $N(0,\lambda_{1})$, and $\pi_{t,s}$ is the density of $t$ distribution with tiny scale parameter satisfying $-\log s\asymp\log n$, i.e. conditions in Corollary 2.1 holds, and we consider the misspecified posterior of form $\pi(\theta|D_{n})=\exp\{-(y-\theta)^{2}/2\}\pi(\theta)$.

Define $\mu\in\mathbb{R}^{n}$ as $\mu_{i}=0$ for all $1\leq i\leq k=3n/4$, and $\mu_{i}=Z_{0.25}/2$ for $i>k$ where $Z_{0.25}$ is the right 25% quantile of standard normal distribution, thus $\|\mu-\theta^{*}\|^{2}\asymp n$.

Let $\Delta\theta$ be any vector such that $\|\Delta\theta\|^{2}\leq M\log n$. Then

[TABLE]

And

[TABLE]

for some positive constant $c$ given sufficiently large $n$, where the inequalities above hold since $y_{n}\geq y_{n-1}\dots\geq y_{k+1}$, and $y_{k+1}\approx Z_{0.25}$ with high probability, due to large sample empirical quantile theory.

Combining the above two results, we have that the posterior density satisfies $\pi(\mu+\Delta\theta|D_{n})\gg\pi(\theta^{*}+\Delta\theta|D_{n})$ for any $\|\Delta\theta\|^{2}\leq M\log n$ with high probability. Therefore, more posterior mass is distributed within the $\sqrt{M\log n}$-radius ball centered at $\mu$ than at the true parameter $\theta^{*}$.

Proof of Theorem 3.2

The proof of this theorem is quite similar to the proof of Theorem 2.1 and 2.2. We define the same testing function as in the proof of theorem 2.1 and 2.2, and define the following two sets:

[TABLE]

Using the same arguments, one can still establish exponential separation results (15) and (16).

To establish (17), we notice that

[TABLE]

for some constant $c_{3}^{\prime}$ and

[TABLE]

This ensures (17).

As for the prior probability of $B_{n}$, if the scale parameter for the $t$ distribution is sufficiently small, i.e. $s=n^{-w}$ for some large $w$ and $\int_{\pm\epsilon_{n}/n^{2}}\pi_{t,s}(x)dx\geq 1-1/n^{1+u}$ for some sufficiently large $u$ where $\pi_{t,s}$ denotes the $t$ density function with scale parameter $s$, then for any ranking $r$,

[TABLE]

This hence implies that $-\log(\pi(B_{n}))\geq u\log n$.