Bayesian Fusion Estimation via t-Shrinkage
Qifan Song Guang Cheng


empty





Bayesian Fusion Estimation via t-Shrinkage
    Qifan Song and Guang Cheng
Qifan Song is an assistant professor in Department of Statistics, Purdue University. Email: qfsong@purdue.edu;
Guang Cheng is a professor in Department of Statistics, Purdue University. Email: chengg@purdue.edu.

    
===================================================================================================================================================================================================================================================





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.



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




§ 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 <cit.> 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 <cit.>, SCAD <cit.>, adaptive Lasso <cit.> and MCP <cit.>. For high dimensional Bayesian inferences, sparsity induced prior, such as spike-and-slab prior <cit.>, is widely
used for model selection. 
Fused Lasso <cit.> considers another type of low dimensional embedding of a high dimensional parameter θ=(θ_i)_i=1^p where the successive differences ϑ_i = θ_i-θ_i-1 are assumed to be sparse as well, in other words, there exists a consecutive block partition of θ_i's, such that θ_i's are constant within each block.
Fused lasso method proposes a penalty function λ_1∑|θ_i|+λ_2∑ |ϑ_i|
which consists of two terms that encourage sparsity among θ_i's and ϑ_i's respectively.

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

    y_i=θ_i^*+ε_i

where ε_i's are iid normal error with unknown variance σ^2. Similarly to Fused Lasso applications, we also assume true parameter θ^* is blocky in the sense that 
there exists a partition {_1^*,…, _s^*}
of {1,...,n} such that θ^*_i's are constant
for all i∈_k^*. Correspondingly, we define set G^*={2≤ i≤ n: ϑ_i^*:=θ_i^*-θ_i-1^*≠ 0 } whose number of elements is supposed to be much smaller than n.
We are interested in conducting Bayesian structure recovery of θ^*.
Motivated by the L_1-fusion penalty used by Fused lasso estimator <cit.>, as well as the development of the Bayesian lasso <cit.>, <cit.> 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 <cit.>.
It is believed that the posterior inconsistency of Laplace prior is due to its exponentially light tail, and
<cit.> 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 (θ_i-θ_i-1) leads to a smoothly varying θ_i estimation rather than a blocky θ, thus it fails to identify the blocking structure.
Therefore, in this paper, we propose to use independent student-t priors on successive differences θ_i-θ_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 <cit.> 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 θ^* has a consecutive blocking structure, it is more realistic to assume that θ^* possesses an unknown clustering structure. In other words, some unnecessarily consecutive θ_i^*'s share the same value within an _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 θ 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 <cit.> usually impose discrete priors on the clustering structure, along with a conditional prior on θ given specific clustering structure. In contrast, we propose to directly model the parameter θ 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 <ref>, 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 <ref>, we will use the t-fusion prior to solve the clustering problem. 
Several simulation studies and one real data application are presented 
in Section <ref>.
At last, Section <ref> 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≻ b_n means lim(a_n/b_n)=∞ and a_n≍ b_n means -∞<liminf(a_n/b_n)≤limsup(a_n/b_n)<∞.
x and x_1 denote L_2 and L_1 norms of vector x.




§ BAYESIAN FUSION VIA T-SHRINKAGE



 §.§ Bayesian Modeling

Suppose we observe independent data {y_i: i=1,…,n} following model (<ref>). The indexing of the data has certain practical or scientific meaning, under which we can assume that the parameter vector
θ^* is “stepwise", in the sense that most of the successive differences
ϑ_i=θ_i-θ_i-1 are exactly 0.
To induce the sparsity for both θ_i's and ϑ_i's, <cit.> proposed the following fused lasso estimator

    θ^ FL=min(y-θ^2/2+λ_1∑_i=1^n|θ_i|+λ_2∑_i=2^n|ϑ_i|).

If one is not interested in pursuing the sparsity of θ's, then a fusion estimator <cit.> can be used

    θ^ F=min(y-θ^2/2+λ∑_i=2^n|ϑ_i|)=min(y-θ^2/2+λ∑_i=2^n|θ_i-θ_i-1|),

for some tuning parameter λ. The above objective functions are both convex and fast computation
algorithms are developed, e.g. <cit.>. The penalty term λ∑_i=2^n|θ_i-θ_i-1| is can be interpreted as 
the negative logarithm of prior density used for Bayesian inferences, therefore, a nature Bayesian expansion to (<ref>) is Laplace (double exponential) prior modeling <cit.>. To account for the unknown variance parameter σ^2 and θ_1, a convenient prior specification could be 

    σ^2   ∼(a_σ,b_σ),  θ_1|σ^2∼ N(0, σ^2λ_1),
    
    (θ_i-θ_i-1)|σ^2    ∼(λ/σ), i=2,…, n

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

According to <cit.>, the above Laplace prior (λ/σ) can be rewritten as a scale mixture of normal distributions:

    (θ_i-θ_i-1)|σ^2,λ_i ∼ N(0, λ_iσ^2),  λ_i∼(-λ^2/2),

where (a) denotes the exponential distribution f(x)∝exp(-ax).
This hierarchical representation for Laplace prior leads to a Gibbs sampling 
update that is similar to the Bayesian lasso <cit.>:

    λ_i^-1∼   (λσ/|θ_i-θ_i-1|, λ^2) i=2,…, n,
    σ^2∼   (a_σ+n,b_σ+y-θ^2/2+θ_1^2/2λ_1
    +∑_i=2^n(θ_i-θ_i-1)^2/2λ_i),
    θ_i∼   N( μ_i,ν_i),

where 
ν_i^-1=1/σ^2+1/λ_i+1σ^2+1/λ_iσ^2 and 
μ_i = ν_i(y_i/σ^2+θ_i+1/λ_i+1σ^2+θ_i-1/λ_iσ^2), λ_n+1 and λ_0 are considered to be infinite, and 
(a,b) denotes inverse Gaussian distribution with cdf
f(x)∝ x^-3/2exp[-b(x-a)^2/(2a^2x)].

Despite the popularity of Laplace prior in many applications, recent Bayesian works <cit.> point 
out that, if we impose independent
θ_i|σ^2∼(λ/σ) 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 
θ doesn't contract into a small neighborhood around true value θ^* appropriately. For a blocky parameter θ, we observe similar empirical results,
as showed in the toy example in Section <ref>: 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 θ doesn't have a step-wise pattern at all.

Following the theoretical discovery of <cit.>, we consider using a class of heavy tailed 
priors for the successive differences ϑ_i's. Specifically, this work will
assign an t-shrinkage prior:

    σ^2   ∼(a_σ,b_σ),  θ_1|σ^2∼ N(0, σ^2λ_1),
    
    (θ_i-θ_i-1)|σ^2    ∼ t_df(sσ), i=2,…, n

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

    (θ_i-θ_i-1)|σ^2,λ_i ∼ N(0, λ_iσ^2),  λ_i∼(a_t, b_t),

where a_t, b_t satisfy df=2a_t and s=√(b_t/a_t).
Under this t-prior, the posterior distribution still allows 
a full conditional Gibbs sampler, where the update for θ_i's
and σ^2 are exactly the same as in (<ref>) and the
update of λ_i's follows

    λ_i∼(a_t+1/2,b_t+(θ_i-θ_i-1)^2/2σ^2)
    i=2,…, n.


To further understand the difference between the Laplace fusion prior and 
t fusion prior, we compare their conditional prior 
π(θ_i|θ_i-1,θ_i+1, σ). Figure <ref>
plots the function -log[π(θ_i|θ_i-1=-1,θ_i+1=1, σ=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 
θ_i-1 and θ_i+1, given a sufficiently small scale parameter
s. In other word, the prior introduces a strong shrinkage effect on θ_i, towards either θ_i-1 or θ_i+1. 
Therefore, for all i=2,…,n, θ_i will merge with either θ_i-1 or θ_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 [θ_i-1,θ_i+1].
Hence, it doesn't encourage the posterior of θ_i to be grouped with either θ_i-1 or θ_i+1.
It is worth to mention that the NEG-fusion prior <cit.>
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 <cit.> 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 <ref> 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.







 §.§ Posterior contraction of Bayesian t-fusion

In this section, we study the theoretical performance induced by the t-fusion prior specification (<ref>). Our theoretical investigation follows the framework of <cit.>, which studies the posterior convergence rate of coefficient β in high dimensional sparse regression models y=Xβ+ε.
Note the model (<ref>) 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 β=(θ_1,ϑ_2,…,ϑ_n) is a unknown sparse vector.
The following theorem studies the general posterior convergence properties given an independent prior over all ϑ_i's.

Assume that |G^*|≺ n/log(n), and prior specification follows that 
σ^2∼(a_σ,b_σ), θ_1 and ϑ_i's are conditionally independent given σ^2 with prior density π(θ_1,ϑ_i's|σ)∝ (1/σ)^nf_θ(θ_1/σ)∏_i=2^n f_ϑ(ϑ_i/σ).
Furthermore, if

    ∫_-|G^*|log(n)/n^2^|G^*|log(n)/n^2f_ϑ(x)dx≥ 1- n^1+u,
       -log(π_ϑ)=O(log n),  π_ϑ = min_|x|≤max_i|ϑ^*_i/σ^*|+1 f_ϑ(x),
       -log(π_θ)=O(|G^*|log n),  π_θ = min_|x|≤ |θ^*_1/σ^*|+1 f_θ(x),
       a_σlog(1/b_σ)+b_σ/σ^*2+(a_σ+2)log(σ^*2)=O(|G^*|log n),

then there exist a constant M, and ϵ_n≍√(|G^*|log n/n), such that the posterior distribution satisfies

    π(θ-θ^*≥ Mσ^*√(n)ϵ_n|y ) → 0,

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 <cit.>, and for the sake of readability, the proof is provided in the Appendix.
The first inequality of sufficient condition set (<ref>) requires that
the prior imposed on ϑ_i's is highly concentrated around zero, such that it induces sufficient Bayesian shrinkage effect for those ϑ_i's whose true value is 0. The rest inequalities of (<ref>) essentially require that
the prior density at true parameter is at least of order e^-cnϵ_n^2 for some c, and this helps to prevent over-shrinkage for those ϑ_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 <cit.>. Given the concrete forms for prior density f_θ and f_ϑ, the second and third inequalities of (<ref>) are equivalent to some upper bound constraints on the magnitude of
θ_1^* and max_i|ϑ_i^*|'s (see e.g., Corollary <ref>). The fourth inequality of (<ref>) trivially holds for any fixed a_σ and b_σ if σ^*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 <cit.> commonly assume that σ^*2∝ n^-1, then one can choose a fixed a_σ and 
b_σ≍ n^-κ for some κ>0. Under such a choice, condition (<ref>) holds as long as -klog n≤log(σ^*2)<Klog n for some positive constant K.

The above result states that almost all the posterior mass contract into
a neighborhood of θ^* with a radius Mσ^*√(n)ϵ_n, that 
is, the posterior convergence rate is of an order σ^*√(|G^*|log n). Note that if
the partition index set G^* were known, the oracle rate of contraction turns out to be O(σ^*√(|G^*|)). Hence, the Bayesian shrinkage achieves the ideal risk up to a logarithmic term in n. In frequentist literature, Theorem 2.7 of <cit.> showed that the convergence rate of Fused Lasso is no larger than O(σ^*√(|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 ϑ 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.

Assume the conditions of Theorem <ref> hold, denote
G(θ,σ)={i:|ϑ_i/σ|<ϵ_n/n}, then the posterior of G(θ,σ)
satisfies

    π({|G(θ,σ)\ G^*|>δ|G^*|  } |y)→ 0,

for some fixed constant δ, where the convergence holds in probability and in L_1.


Therefore, the posterior distribution of G(θ,σ), which is induced by the posterior of (θ, σ^2) and mapping G(·,·), can be treated as a discrete posterior distribution for the unknown G^* and used for Bayesian block partition selection. Theorem <ref> essentially claims that the number false positive selections of via discretization G(·,·) is bounded in posterior probability. Such a result is comparable to the model selection behavior of Bayesian Dirichlet-Laplace
shrinkage <cit.>. It is worth to note that if certain minimal signal strength condition holds as well, i.e., min_ϑ_i^*>0|ϑ_i^*| is bounded away from zero, with addition assumptions, one can derive an even stronger posterior selection consistency result that π(G(θ,σ)=G^*   |y)→^p 1, following the proof of Theorem 2.3 in <cit.>. 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 (<ref>).

When f_θ is the cdf of normal distribution N(0,λ_1) and f_ϑ is the cdf of a t-distribution with scale parameter s, then the first three inequalities of
(<ref>) hold when

    log(max_i|ϑ_i^*|/σ^*)=O(log(n)),
       θ_1^*2/(λ_1σ^*2)+log(λ_1)=O(log n),

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_ϑ, as long as f_ϑ 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|ϑ_i^*|/σ^*=o(1). 



 §.§ 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 σ^*=0.5. The data and underlying true parameter value are plotted in Figure <ref>.
Three estimation procedures are considered: 1) L_1 fusion estimation (<ref>) where 
the tuning parameter are selected by cross validation; 2) Bayesian Laplace prior (<ref>) with λ = √(2log (n)) and λ_1=5; 3) Bayesian t prior (<ref>)
with a_t=2, b_t=0.005, a_σ=0.5, b_σ=0.5 and λ_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 <ref>.
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 ϑ estimation and the estimated θ 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 ϑ_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 <ref>, Laplace fusion prior itself merely has shrinkage effect for the successive differences, and hence the data fluctuation carries over to the posterior distribution of θ. In this case, it is not clear to recover the underlying block structure.

Note that for both (<ref>) and (<ref>) prior specifications, their posterior distributions highly depend on the choice of hyperparameter, especially the scale parameters of the Laplace (parameter 1/λ) 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 √(10), for both Laplace and t priors. The corresponding posteriors are displayed
in Figure <ref>.



From Figure <ref> and its comparison with Figure <ref>, 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 θ_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 <ref>. The choice of scale parameter only affects the smoothness of posterior of θ. 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 <ref>.
As the scale of the Laplace prior decreases (i.e., λ increases), the 
conditional prior of θ_i given (θ_i-1,θ_i+1) acts more and more like a uniform distribution over [θ_i-1,θ_i+1].
Hence, conditional on (θ_i-1,θ_i+1), the posterior of “local" total variation, |θ_i+1-θ_i|+|θ_i-θ_i-1|, always
decreases to the minimum value |θ_i+1-θ_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 (<ref>), 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 <cit.>, 
choosing a proper prior hyperparameter posts many difficulty for Bayesian statisticians.
Conventional choices include imposing a hyper-prior on the hyperparameter 
<cit.> or empirical Bayes <cit.>.
For high dimensional sparse GLMs, <cit.> suggested choosing a hyperparameter such that the posterior mean and posterior mode are close. Several Bayesian works <cit.>
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 (<ref>) such that 

    P(|t_df(s)|≥√(log (n)/n)) ≈ 1/n.

Although it doesn't quite satisfy the conditions in Corollary <ref>
 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 (<ref>).

 


§ 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 θ in model
(<ref>) 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
<cit.>. 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 <cit.>.
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 θ_i-θ_i-1 shall not work for this clustering problem since {θ_i-θ_i-1}_i=2^n is no longer a sparse vector. However, if some meaningful permutation r of θ_i is known
such that {θ_r(i)-θ_r(i-1)}_i=2^n is indeed a sparse vector, for example, θ_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(·) is generally not available. Thus one simple solution would be to substitute r(·) with some estimator r(). 
For model (<ref>), a trivial estimator is the rank statistic of the responses y_i, i.e., y_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 r:

    σ^2   ∼(a_σ,b_σ),  θ_r(1)|σ^2∼ N(0, σ^2λ_1),
    
    (θ_r(i)-θ_r(i-1))|σ^2,λ_i    ∼ N(0, σ^2
    λ_i),   λ_i∼(a_t,b_t)  i=2,…, n.

The corresponding posterior Gibbs sampler follows:

    λ_i∼   (a_t+1/2,b_t+(θ_r(i)-θ_r(i-1))^2/2σ^2)
    i=2,…, n.
    σ^2∼   (a_σ+n,b_σ+y-θ^2/2+θ_r(1)^2/2λ_1
    +∑_i=2^n(θ_r(i)-θ_r(i-1))^2/2λ_i),
    θ_r(i)∼   N( μ_i,ν_i),

where ν_i^-1=1/σ^2+1/λ_i+1σ^2+1/λ_iσ^2 and 
μ_i = ν_i(y_r(i)/σ^2+θ_r(i+1)/λ_i+1σ^2+θ_r(i-1)/λ_iσ^2).






The idea of taking advantage of a pilot estimator θ and its rank statistic r(i) has also been implemented in the two-stage frequentist approach <cit.>.

It is worth to mention that posterior consistency Theorems <ref> and <ref> for Bayesian shrinkage fusion don't apply to the above Bayesian modeling (<ref>)
even when θ_r(i) is in a ascending order, because ε_r(i) are no longer iid error observations. 
Such a data dependent prior (<ref>) can also be interpreted 
as a misspecified Bayesian modeling <cit.>, 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 (<ref>) 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 <ref> than the toy data showed in Figure <ref>. 
Secondly, (<ref>) fails to take account of uncertainty of the estimator r, and the strict monotonicity of y_r(i) will always
carry over to the posterior of θ_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 θ_r(i)'s should contract towards the best step function, in the sense of minimizing the Kullback–Leibler divergence <cit.>,
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
 <ref> doesn't hold anymore.

Assume that all y_i∼ N(0,1) with known variance σ^*=1. If the prior (<ref>) is used (except that there is no necessity to impose a prior on σ^2) and the scale parameter satisfies -log(b_t)≍log n, then in high probability

    π(θ≤√(Mlog n)|y)<1/2,

i.e., the rate-√(log n) posterior consistency fails.


The negative result in Theorem <ref> motivates us to propose an adaptive pseudo Bayesian shrinkage approach. Instead of using a fixed estimator r, we update the “working" order r of θ over iterations. To be specific, we modify the above Gibbs sampling iteration as:

    ,
       .
 
In other words, we update r(·) to the rank statistic of the newly obtained θ. 

There are a couple of rationales behind the algorithm (<ref>).
First, the update of r potentially allows algorithm (<ref>) to incorporate the uncertainty of rank statistic into the posterior sampling. Heuristically, if the current r(·) is close to the true ranking, (<ref>) will shuffle the ordering of θ 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 (<ref>) can be somehow connected with a full Bayesian modeling:

    r   ∼,
    σ^2   ∼(a_σ,b_σ),  θ_ r(1)|σ^2∼ N(0, σ^2λ_1),
    
    (θ_ r(i)-θ_ r(i-1))|σ^2,λ_i    ∼ N(0, σ^2
    λ_i),   λ_i∼(a_t,b_t)  i=2,…, n.

Under (<ref>), the conditional posterior of r(·) follows 
π(r|θ,λ_i's, y)∝exp{-∑_i=2^n(θ_ r(i)-θ_ r(i-1))^2/λ_i^2σ^2}, from which sampling is expensive. However, when most of λ_i^2s are tiny, the distribution of π(r|θ,λ_i's, y) will be highly concentrated around its MAP
which is approximately the rank statistic of θ. Therefore, updating r to the rank statistic of current sample θ, as in (<ref>), can be viewed as 
a convenient sampling of π(r|θ,λ_i'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 <ref> when there is only one underlying cluster.


Assume that all y_i∼ N(θ_0,σ^2), i.e. θ^* is a vector
of θ_0's. If prior (<ref>) is used with b_t = n^-c for some sufficiently large c, and
θ_0^2/(λ_1σ^*2)+log(λ_1)=O(log n)
then there exist constant M and ϵ_n≍√(log n/n), such that

    π(θ-θ^*≥ Mσ^*√(n)ϵ_n|y)→ 0,

where the convergence holds in probability and L_1.

Further investigations will be pursued to assess the posterior contraction induced by (<ref>)
when the true parameter θ^* has multiple-cluster structure.
The minimal cluster difference min_{(i,j):θ_i≠θ_j}|θ_i-θ_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 √(n)-rate L_2 error. And we conjecture that when
the minimal cluster difference is large than √(Mlog n) for some constant M, the posterior distribution induced by (<ref>) leads to satisfactory Bayesian consistency.

It is worth to mention that the r-update step in algorithm (<ref>)
may critically change the stochastic stability of the algorithm and that the ergodicity of Markov chains generated by (<ref>) is still unclear.
Therefore, in our application we use it with caution. In all our simulations and toy examples, we initialize with r=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 θ_i's is 
to impose a pairwise difference penalty as ∑_i≠ jp(|θ_i-θ_j|) for some penalty p(·)<cit.>.
Although it is tempting to develop a Bayesian counterpart, i.e., using a prior of form
π(θ) ∝∏_i≠ jπ(θ_i-θ_j) with sparsity induced π, there are several problems.
First, the pairwise differences {θ_i-θ_j}_i≠ j
must satisfy triangle inequality, hence the prior specification π(θ) ∝∏_i≠ jπ(θ_i-θ_j)
is actually not an independent prior over all pairwise difference {θ_i-θ_j}_i≠ j, and it will be difficult to characterize the effect of such dependency in the prior distribution.
Second, even under clustering structure, {θ_i-θ_j}_i≠ j are actually not sparse. For example, if there are two balanced groups among θ_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 θ_i towards y̅. 
Third, such pairwise difference prior will also heavily increase the computational burden of posterior sampling.

To illustrate and compare the performance of (<ref>) and (<ref>), a simple toy example is conducted. 
We simulate a data set with n=100 data points that belong to 3 underlying subgroups. 
For both (<ref>) and (<ref>), we choose a_t=2, b_t=0.005, a_σ=b_σ=0.5 and λ_1=5.
In addition, we compare them with
frequentist L_1 fusion

    θ= miny-θ^2/2+λ∑_i=2^n|θ_r(i)-θ_r(i-1)|,

and Bayesian Dirichlet process modeling

    y_i|θ_i,σ^2∼ N(θ_i, σ^2),  θ_i∼ G,
       G∼(N(0, λ), a)σ^2∼(a_σ,b_σ)

with λ=5,  a=0.1 and a_σ=b_σ=0.5.
The posterior sampling algorithms of Dirichlet process modeling are discussed in
<cit.>.
The frequentist estimator and posterior boxplots of θ_is under different priors are plotted in Figure <ref>. 
This figure clearly shows that frequentist L_1 fusion estimator fails to yield clustering structure for monotone sorted y_r(i)'s, especially for the middle portion of the data, one always have θ_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 (<ref>) 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 (<ref>) combines both advantages of DP prior and (<ref>). Comparing with DP prior, it
has much smaller posterior variance; comparing with prior (<ref>),
the over-clustering issue is greatly alleviated.







§ SIMULATION AND DATA ANLAYSIS




 §.§  Bayesian fusion simulations



In our first simulation studies, we consider model 
(<ref>) with σ^*=0.5 and true parameter θ^*∈^100 having three consecutive blocks. These three blocks are
_1^* ={1,…, b_1}, _2^*={b_1+1,…, b_1+b_2}, 
_3^*={b_1+b+2+1,…, b_1+b_2+b_3} where
(b_1,b_2,b_3)∼(100, (1/3,1/3,1/3)),
and θ_i^* =0, 2 or 4 within each block respectively.
We compare the performance among Bayesian t fusion (<ref>), Bayesian Laplace fusion (<ref>)
and frequentist L_1 fusion (<ref>), based on 100 replications. The choices of the hyperparameters are same to the toy example discussed in Section <ref>.

To compare the accuracy of parameter estimations, we calculate the L_2 and L_1 estimation errors, θ-θ^*_2 and
θ-θ^*_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 Σ=(ω_ij)=(1{θ_i=θ_j}), and corresponding 
estimation error of Σ: R=Σ^*-Σ_1=∑_i,j
|ω_ij-ω_ij^*|.
For L_1 fusion estimator, Σ is trivially estimated by
ω_ij=1{θ^ F_i=θ^ 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 Ω.
In the literature, this is usually done by 1) hard thresholding
<cit.>, 
or 2) decoupling shrinkage and selection methods <cit.>.
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 <cit.>, we estimate
ω_ij=1{|θ_i-θ_j|/σ≤π_1/2n}, where
θ_i and σ are the Bayes estimator, 
π_1/2n denotes the (1-1/2n) quantile of the
prior distribution imposed on successive difference ϑ_i/σ.
This choice tries to increases the robustness of Bayesian inference: if an
overly small or large scale parameter s is used, the estimation for ω_i,j will adapt accordingly.
For t shrinkage (<ref>), π_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 (<ref>), π_1/2n=log(n)/λ.
Since different priors have different corresponding π_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=_{w_ij^*≠0}|θ_i-θ_j| denotes within-block average variation, and
B=min_{w_ij^*=0}|θ_i-θ_j| denotes 
the between-block separation.
A larger B and smaller W indicate a better block identification performance. 
The results are summarized in Table <ref>.




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
<cit.>, but unlike the LASSO penalty which introduces a bias of λ, fuse-L_1 penalty introduces only a bias as small as λ/n_i (<cit.>) 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 
<ref>, the Laplace fusion tends to yield smooth changing θ_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 π_1/2n is much larger than t fusion. 



 §.§  Bayesian clustering simulations

In our second simulation studies, we consider model 
(<ref>) with σ^*=0.5 and true parameter θ^*∈^100 having three unknown clusters.
With equal probability, θ_i^*=0, 2, or 4. 
We aim to compare different approaches including Bayesian t modeling (<ref>) and (<ref>),
Dirichlet process prior (<ref>) and the frequentist bCARDS estimation (<ref>) using L_1 fusion <cit.>.
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_π(θ|y)θ-θ^*^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 <ref>, but as mentioned in Section <ref>, clustering problem is much more difficult than fusion problem. 
To see that, assume there are three balanced clusters for θ, 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 <cit.>, E(z_1)≈ 1.04 > 0.96 ≈ E(z_2). Therefore, misclassification will always happen for those extreme observations, and 
B=min_{w_ij^*=0}|θ_i-θ_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 

    B̃={|θ_i-θ_j|}_{w_ij^*=0}.

The simulation results are summarized in Table <ref>.




The comparison shows that DP prior yields the best point estimator in terms of estimation error, and the adaptive t-shrinkage method (<ref>) 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 <ref> that LP prior leads to larger posterior variance.
The three Bayesian approaches have approximately the same performance for between-cluster
separation (B̃ statistic), while the DP prior has a smaller
within-cluster variation (W statistic).
The frequentist estimator (<ref>) has the worst performance in almost every aspect.
In summary, DP prior does yield the best Bayes point estimation, but the adaptive t-shrinkage (<ref>) 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 (<ref>) and (<ref>), we conclude that the r-update step implemented in algorithm (<ref>)
does improve the performance of t-shrinkage in every aspect for this clustering problem.



 §.§  Real data set analysis

In this section, we consider a real comparative genomic hybridization
(CGH) dataset <cit.>. The dataset is available from the R package  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 (<ref>) and frequentist
L_1 fusion (<ref>).
For the Bayesian application, the hyperparameter is chosen as 
a_t=2 and b_t=0.0007, following the suggestion of (<ref>).
Both point estimations are plotted in Figure <ref>.
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 θ_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 
θ_i and the sample mean of the block.







§ FINAL REMARKS

In this work, we study the Bayesian inference for a vector parameter θ 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β+Zθ+ϵ for some blocky parameter θ. 

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 r, and then fuse all pairs of neighbors (that are determined by r) via shrinkage priors. An adaptive r-update step is further incorporated to improve the clustering performance.
Simulations show that the proposed algorithms (<ref>) and (<ref>) 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=θ_i∈^d.
Since we can not rank multivariate vectors, the possible alternative is to construct a minimum spanning tree (MST)<cit.>, and impose shrinkage fusion prior on the pair of θ_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

First, let us state some useful lemmas.

	Let χ^2_d(κ) be a chi-square distribution with degree of freedom d, and noncentral parameter κ, then we have the following concentration inequality
	
    Pr( χ^2_d(κ)>d+κ+2x+√((4d+8κ)x))≤exp(-x), 
       Pr( χ^2_d(κ)<d+κ-√((4d+8κ)x))≤exp(-x).




        Let X be a Binomial random variable X∼(n,v). For any 1<k<n-1
	
    Pr(X≥ k+1)≤ 1- Φ((k-nv){2nH(v, k/n)}^1/2),

	where Φ 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 <cit.>:

	Let f^* be the true probability density of data generation, f_θ be the likelihood function
	with parameter θ∈Θ, and E^*, E_θ denote the corresponding expectation
	respectively. Let B_n and C_n be two subsets in parameter space Θ, and ϕ_n be some
	test function satisfying 
	ϕ_n(D_n)∈[0,1] for any data D_n. If π(B_n)≤ b_n, E^*ϕ(D_n)≤ b_n',
	sup_θ∈ C_nE_θ(1-ϕ(D_n))≤ c_n, and  furthermore,  
	
    P^*{m(D_n)/f^*(D_n)≥ a_n }≥ 1-a_n',

	where m(D_n)=∫_Θπ(θ)f_θ(D_n)dθ is the margin probability of D_n. Then,
	
    E^*(π(C_n∪ B_n)|D_n))≤b_n+c_n/a_n+a_n'+b_n'.



	Define Ω_n to be the event of (m(D_n))/(f^*(D_n))≥ a_n,
	and m(D_n, C_n∪ B_n) = ∫_C_n∪ B_nπ(θ)f_θ(D_n)dθ. Then
	
    E^*π(C_n∪ B_n)|D_n) = E^*π(C_n∪ B_n)|D_n)(1-ϕ(D_n))1_Ω_n
    
    	+   E^*π(C_n∪ B_n)|D_n)(1-ϕ(D_n))(1-1_Ω_n)+ E^*π(C_n∪ B_n)|D_n)ϕ(D_n)
    ≤    E^*π(C_n∪ B_n)|D_n)(1-ϕ(D_n))1_Ω_n+E^*(1-1_Ω_n)+ E^*ϕ(D_n)
    ≤    E^*π(C_n∪ B_n)|D_n)(1-ϕ(D_n))1_Ω_n+b_n'+a_n'
    ≤    E^*{m(D_n, C_n∪ B_n)/a_nf^*(D_n)}(1-ϕ(D_n))+b_n'+a_n'.

	
	By Fubini theorem,
	
    E^*(1-ϕ(D_n))m(D_n, C_n∪ B_n)/f^*(D_n) = ∫_C_n∪ B_n∫_[1-ϕ(D_n)]f_θ(D_n) dD_nπ(θ)dθ
    ≤   ∫_C_n E_θ(1-ϕ(D_n))π(θ)dθ+∫_ B_n∫_f_θ(D_n) dD_nπ(θ)dθ≤ b_n+c_n.

	
	Combining the above inequalities leads to the conclusion.

Proof of Theorem <ref> and <ref>
0.1in
Let G={g_1,g_2,…,g_d} be a generic subset of {2,…, n}, and it also represents potential a (d+1)-group structure of θ
as {{1,…, g_1},{ g_1+1,…, g_2}, …, {g_d+1,…, n }}.
Given G and its corresponding blocking structure, θ_G(y) denotes the estimation of θ based on block mean, i.e. θ_G,j(y) = ∑_i=g_j+1^g_j+1y_i/(g_j+1-g_j) for all g_j+1≤ j≤ g_j+1,
and σ^2_G(y)=y-θ_G^2/(n-|G|-1).

To prove the posterior contraction, we will apply Lemma <ref>. We define the following testing function

    ϕ(y) = 1{   θ_G-θ^*≥√(n)σ^*ϵ_n
        |σ^2_G-σ^*2|>σ^*2ϵ_n
       G⊃ G^*,
        |G|≤ (1+δ)|G^*|}

for some δ>0, and define C_n and B_n as:

    C_n={θ: θ-θ^*≤ M√(n)σ^*ϵ_n, (1-ϵ_n)/(1+ϵ_n)<σ^2/σ^*2<(1+ϵ_n)/(1-ϵ_n)}^c\ B_n,
       B_n={θ: {ϑ_i's}_i∉ G^* , δ|G^*|
     σϵ_n/n}.


Note that when G⊃ G^*, θ_G(y)-θ^*^2∼σ^*2χ^2_|G|+1 and y-θ_G^2∼σ^*2χ^2_n-|G|-1, thus by Lemma <ref>,
we have that 

    P(θ_G(y)-θ^*≥√(n)σ^*ϵ_n
    |σ^2_G-σ^*2|>σ^*2ϵ_n)
    ≤exp{-c_1nϵ_n^2},

for some constant c_1, since |G|=O(|G^*|)≺ nϵ_n^2 and ϵ_n≺ 1.
Therefore,

    E_(θ^*,σ^*2)ϕ(y)≤n-1  (1+δ)|G^*|exp{-c_1nϵ_n^2}=exp{-c_1'nϵ_n^2},

as long as nϵ^2/[|G^*|log n] is sufficiently large.

For any (θ,σ^2)∈ C_n satisfying θ-θ^*≤ M√(n)σ^*ϵ_n and 
σ^2/σ^*2≥(1-ϵ_n)/(1+ϵ_n), we define G={i:θ_i-θ_i-1≥σϵ_n/n^2}∪ G^* ( hence |G|≤ (1+δ)|G^*|), thus

    P_(θ,σ^2)(  θ_G(y)-θ^*≤√(n)σ^*ϵ_n)
      = P_(θ,σ^2)(  θ_G(y)-θ_G(θ) + θ_G(θ)-θ^*≤√(n)σ^*ϵ_n)
    ≤   P_(θ,σ^2)(θ_G(y)-θ_G(θ)≥θ_G(θ)-θ^*-√(n)σ^*ϵ_n) 
    ≤    P_(θ,σ^2)(θ_G(y)-θ_G(θ)≥ M√(n)σ^*ϵ_n-√(n)σϵ_n-√(n)σ^*ϵ_n) 
    ≤   P(χ_|G+1|^2≥[√(1-ϵ_n/1+ϵ_n) (M-1)-1] nϵ_n^2)≤exp{-c_2nϵ_n^2}

for some c_2 given a large M, where the second inequality is due to the fact that θ_G(θ)-θ≤√(n)σϵ_n when θ∈ B_n. 

For any (θ,σ^2)∈ C_n satisfying σ^2/σ^*2<(1-ϵ_n)/(1+ϵ_n)
or σ^2/σ^*2>(1+ϵ_n)/(1-ϵ_n),

    P_(θ,σ^2)(   |σ^2_G-σ^*2|<σ^*2ϵ_n)
      = P_(θ,σ^2)(   |y-θ_G^2/σ^*2(n-|G|-1)-1|<ϵ_n)
    ≤    P_(θ,σ^2)( 1-ϵ_n<  y-θ_G^2/σ^*2(n-|G|-1)<1+ϵ_n)
    ≤   P_(θ,σ^2)( |y-θ_G(y)^2/σ^2-(n-|G|-1)|>(n-|G|-1)ϵ_n)
    
      =   P_(θ,σ^2)( |χ_n-|G|-1^2(λ)-(n-|G|-1)|>(n-|G|-1)ϵ_n)≤exp{-c_2'nϵ_n^2}
 for some c_2', 
where the noncentral parameter λ<nϵ_n^2≺(n-|G|-1)ϵ_n.

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

    sup_(θ,σ^2)∈ C_nE_(θ,σ^2)[1-ϕ(y)]≤max{exp(-c_2nϵ_n^2),exp(-c_2'nϵ_n^2)}.


Now we consider the marginal posterior density of data y. With probability P(ε≤2√(n)σ^*) (which converges to 1), 

    m(y)/f^*(y)=∫_σ^2∫_θσ^*nexp{-θ^*-θ+ε^2/σ^2}
    dθ dσ^2
    /σ^nexp{-ε^2/σ^*2}
    ≥   ∫_σ^2∫_θexp{-θ^*-θ^2/σ^2
    -2(θ^*-θ)^Tε/σ^2+ε^2/σ^*2-ε^2/σ^2-nlog(σ/σ^*)}π(θ,σ^2)dθ dσ^2
    
    
    
    
    ≥   π(max{|θ_1-θ_1^*|, |ϑ_i-ϑ_i^*|}/σ≤  |G^*|log n/n^2, 0≤σ^2-σ^*2≤σ^*2|G^*|log n/n)exp{-c_3'|G^*|log n }

for some constant c_3'. Besides,

    π(max{|θ_1-θ_1^*|, |ϑ_i-ϑ_i^*|}/σ≤  |G^*|log n/n^2, 0≤σ^2-σ^*2≤σ^*2 |G^*|log n/n)
    ≥   π_σ(σ^*2)*O(σ^*2|G^*|log n/n)* π_θ*O( |G^*|log n/n^2)  
       *π_ϑ^|G^*| [ |G^*|log n/n^2)]^|G^*| * [π_ϑ({-|G^*|log n/n^2, |G^*|log n/n^2})]^n
    
    =   exp{-c_3”|G^*|log n}.
for some constant c_3”. 
Thus 

    m(y)/f^*(y)≥exp{-(c_3'+c_3”)|G^*|log n}
        =exp{-c_3nϵ_n^2},

for some c_3, where c_3 can be sufficiently small when nϵ_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 ϑ_i's,
π(B_n) = π[(n-1-|G^*|,p)>(δ)|G^*| ], where p≤ (1/n)^1+u. By Lemma <ref>,

    π(B_n) ≤exp{ -c_4δ |G^*|log n}

 for some c_4.
Combine results (<ref>), (<ref>), (<ref>) and (<ref>), and we apply Lemma <ref>
to get the posterior consistency result that 

    π(B_n∪ C_n|y) →^p 0,

given a sufficient large constants δ and nϵ_n/|G^*|log n.
 
 0.2in
Proof of Theorem <ref>

Consider y=θ+d, where θ_i^*≡ 0 for all i and error d is order statistics of standard normal variables, i.e. the density of d is f(d)=n!∏ϕ(d_i)1(d_1≤ d_2,…,d_n-1≤ d_n) and ϕ denotes the standard normal density. The prior of θ follows
π(θ)=π_1(θ_1)∏_i=2^nπ_t,s(θ_i-θ_i-1),
where π_λ_1 is the density of N(0,λ_1), and π_t,s is the density of t distribution with tiny scale parameter satisfying -log s≍log n, i.e. conditions in Corollary <ref> holds, and we consider the misspecified posterior of form
π(θ|D_n)=exp{-(y-θ)^2/2}π(θ).

Define μ∈^n as μ_i=0 for all 1≤ i≤ k= 3n/4, and μ_i=Z_0.25/2 for i> k where Z_0.25 is the right 25% quantile of standard normal distribution, thus μ-θ^*^2≍ n.

Let Δθ be any vector such that Δθ^2≤ Mlog n.
Then 

    -log(π(μ+Δθ)/π(θ^*+Δθ))
    =-log(π_t,s(Z_0.25/2+Δθ_k)/π_t,s(Δθ_k)) = O( -log s) = O(log n).

And

    log(exp{-(y-μ-Δθ)^2/2}/exp{-(y-θ^*-Δθ)^2/2})
    =[(y-θ^*-Δθ)^2-(y-μ-Δθ)^2]/2
    
    =   1/2∑_i=k+1^n[(y_i-Δθ_i)^2-(y_i-Z_0.25/2-Δθ_i)^2]=1/2∑_i=k+1^n[(y_i-Δθ_i)Z_0.25-Z_0.25^2/4]
    ≥   1/2[(y_k+1:n_1-√(nMlog n/4))Z_0.25-nZ_0.25^2/16]≥ cn,

for some positive constant c given sufficiently large n, where the inequalities above hold
since y_n≥ y_n-1…≥ y_k+1, and y_k+1≈ 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 
π(μ+Δθ|D_n)≫π(θ^*+Δθ|D_n) for any Δθ^2≤ Mlog n with high probability. Therefore, more posterior mass is distributed within the √(Mlog n)-radius ball centered at μ than at the true parameter θ^*.
 
 0.2in
Proof of Theorem <ref>

 The proof of this theorem is quite similar to the proof of Theorem <ref> and <ref>.  We define the same testing function
 as in the proof of theorem <ref> and <ref>, and define the following two sets:
 
    C_n={θ: θ-θ^*≤ M√(n)σ^*ϵ_n, (1-ϵ_n)/(1+ϵ_n)<σ^2/σ^*2<(1+ϵ_n)/(1-ϵ_n)}^c\ B_n,
       B_n={θ: {θ_i-θ_i-1}_i=2^n , δσϵ_n/n}.

Using the same arguments, one can still establish exponential separation results (<ref>) and (<ref>).

To establish (<ref>), we notice that 

    m(y)/f^*(y)=∫_σ^2∫_θσ^*nexp{-θ^*-θ+ε^2/σ^2}
    dθ dσ^2
    /σ^nexp{-ε^2/σ^*2}
    ≥   ∫_σ^2∫_θexp{-θ^*-θ^2/σ^2
    -2(θ^*-θ)^Tε/σ^2+ε^2/σ^*2-ε^2/σ^2-nlog(σ/σ^*)}π(θ,σ^2)dθ dσ^2
    
    
    
    
    ≥   π(max{|θ_i-θ_i^*|}/σ≤log n/n, 0≤σ^2-σ^*2≤σ^*2log n/n)exp{-c_3'log n }

for some constant c_3' and 

    π(max{|θ_i-θ_i^*|}/σ≤log n/n, 0≤σ^2-σ^*2≤σ^*2log n/n)
    ≥   ∑_rπ(max{|θ_r(1)-θ_r(1)^*|, |θ_r(i)-θ_r(i-1)|}/σ≤  |G^*|log n/n^2, 0≤σ^2-σ^*2≤σ^*2log n/n|r)π(r).
 
This ensures (<ref>).

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 ∫_±ϵ_n/n^2π_t,s(x)dx≥ 1-1/n^1+u 
for some sufficiently large u where π_t,s denotes the t density
function with scale parameter s,
then for any ranking r,

    π(B_n|r)≥ 1-π(max{θ_r(i)-θ_r(i-1)}≤σϵ_n/n^2 |r)
    ≥ 1-(1-1/n^1+u)^n ≈ n^-u.

This hence implies that -log(π(B_n)) ≥ ulog n.

 
0.2in
plain


