# Posterior Propriety of an Objective Prior in a 4-Level Normal Hierarchical Model

Authors: Chengyuan Song; Dongchu Sun; Kun Fan; Rongji Mu
Journal: Mathematical Problems in Engineering
 (2020)

---

## Abstract

The use of hierarchical Bayesian models in statistical practice is extensive, yet it is dangerous to implement the Gibbs sampler without checking that the posterior is proper. Formal approaches to objective Bayesian analysis, such as the Jeffreys-rule approach or reference prior approach, are only implementable in simple hierarchical settings. In this paper, we consider a 4-level multivariate normal hierarchical model. We demonstrate the posterior using our recommended prior which is proper in the 4-level normal hierarchical models. A primary advantage of the recommended prior over other proposed objective priors is that it can be used at any level of a hierarchical model.


---

## Body

## 1. Introduction

Bayesian hierarchical models have a wide range of modern applications including engineering [1], astrophysics [2], economics [3], environmental sciences [4], climatology [5], survival analysis [6], and genetics [7]. It is wonderful here to stay, but hyperparameter priors are often chosen in a casual fashion. Statisticians often use improper priors to express ignorance or to provide good frequency properties. However, it is dangerous to implement the Gibbs sampler without checking that the posterior is proper. As Hobert and Casella [8] pointed, without proper precaution, simple noninformative priors can be misused, sometimes unknowingly, and lead to other difficulties, such as the nonconvergence of the Gibbs sampler. Therefore, it is hazardous to skip the demonstration at the risk of making the inference from an improper posterior distribution. There are many examples of this in the statistical and other literatures. See especially, Wang et al. [9]; Hobert and Casella [8, 10]; Berger and Strawderman [11]; Berger et al. [12]; Speckman et al. [13]; Roy and Dey [14]; Michalak and Morris [15]; Ramos et al. [16]; Tomazella et al. [6]; etc. The importance of the posterior propriety motivates us to explore it.The normal hierarchical distribution has received enormous attention and is also of substantial importance in contemporary statistical theory and application. Ning [17] proposed a 2-level multivariate normal hierarchical model for the degradation data of multiple units with change point. Wang and Coit [18] handled the reliability prediction problem using a multivariate normal distribution model, considering multiple degradation measures. Heuristically, one could also build the multivariate normal hierarchical model with adding a normal distribution for unknown mean vectors to deal with this reliability prediction problem. The unknown parameters because of either convenience or a lack of prior information are often modeled with improper objective priors. Some references of objective priors can be found in Berger and Strawderman [11]; Berger [19]; Pollo et al. [20]; and Ferreira et al. [21]. However, formal approaches to objective Bayesian analysis, such as the Jeffreys-rule approach or reference prior approach, are only implementable in simple hierarchical settings. It is thus common to use less formal approaches, such as utilizing formal priors from nonhierarchical models in hierarchical settings. However, this can be fraught with danger. For instance, nonhierarchical Jeffreys-rule priors for variances or covariance matrices result in improper posterior distributions if they are used at higher levels of a hierarchical model (see [12]).Berger et al. [12] addressed the question of choice of hyperpriors in the following 2-level normal hierarchical model:(1)y∼Npθ,Ip,θ∼Npb,Q,where(2)yp×1=y1y2⋮ym,θp×1=θ1θ2⋮θm,bp×1=ββ⋮β,Qp×p=Im⊗V,for i=1,2,…,m, in which the yi are the k×1 observation vectors and the θi are the k×1 unknown mean vectors (thus p=mk), β is a k×1 unknown “hypermean” vector, and V is an unknown k×k “hypercovariance matrix.” Some commonly priors for the hyperparameter β can be considered. For example, the constant prior and conjugate normal prior. For the prior of the covariance matrix V, a large literature can be found. For example, the articles by Dey and Srinivasan [22]; Lin and Perlman [23]; Haff [24]; Yang and Berger [25]; Daniels and Kass [26, 27]; Ledoit and Wolf [28]; Sun and Berger [29]; Rajaratnam et al. [30]; Hoff [31]; And Berger et al. [12] explored posterior propriety and admissibility of various priors. But no overall conclusion was reached as to a specific prior in their paper. Following the seminal work of Berger et al. [12]; Berger et al. [32] recommended a particular objective prior for choice of πβ and πV:(3)πβ∝11+β2k−1/2,πV=1V1−1/2k∏i<jvi−vj,where v1>v2⋯>vk are the ordered eigenvalues of V. Their recommendation is based on considerations of the posterior propriety, admissibility, ease of implementation (including computational considerations), and performance. A primary advantage of the recommended prior over other proposed objective priors is that it can be used at any level of a hierarchical model. But the property of the resulting posterior is not clear, and it represents a very difficult question to answer. Michalak and Morris [15] pointed out that except in the simplest models when improper priors are used, it can be daunting and time-consuming to verify that the resulting posterior distribution is proper. While Berger et al. [32] suspected that the propriety holds in any level normal hierarchical models, they were not strictly able to prove it. They only showed that the posterior is proper in a 3-level hierarchical model. In this paper, we follow the story and consider the posterior propriety of the recommended prior in a 4-level normal hierarchical model. We demonstrate that the posterior using the recommended prior is still proper in the 4-level normal hierarchical model.The structure of this paper is as follows. Section2 presents a 4-level normal hierarchical model. We propose the recommended hyperpriors for “hypermean” vectors and “hypercovariance matrix.” In Section 3, we demonstrate that the posterior using the recommended prior is proper in the 4-level normal hierarchical model. In Section 4, we consider MCMC computation from the posterior. Section 5 gives the performance of this prior, in comparison with other objective priors that were studied in Berger et al. [12], presenting strong numerical evidence of the superiority of (3). Some concluding remarks are provided in Section 6.


## 2. A 4-Level Normal Hierarchical Model

Consider the following 4-level normal hierarchical model:(4)Level1:yi=θi+Nk0,Ik,i=1,2,…,m,Level2:θi=Ziβ+Nk0,V,βt=β1t,…,βst,Level3:βj=Tjη+Np0,W,j=1,2,…,s,ηt=η1t,…,ηnt,Level4:ηl=ξ+Nq0,Σ,l=1,2,…,n,where the yi are the k×1 observation vectors; the θi are the k×1 unknown mean vectors; the βj are the p×1 unknown vectors; the ηl are the q×1 unknown vectors; ξ is an unknown q-dimensional ”hypermean” vectors; the Zi are the k×sp known matrices; the Tj are the p×nq known matrices; and V, W, and Σ are the unknown k×k, p×p, and q×q “hypercovariance matrix,” respectively. Following the assumption of Berger et al. [32], let k=ps and p=nq. We also assume that m≥2, s≥2, n≥2, and q≥1.There are four unknown parameters:ξ,V,W,andΣ. From Berger et al. [32], the recommended (independent) priors for ξ,V,W,andΣ become(5)πξ∝11+ξ2q−1/2,ξ∈IRq,πV∝1V1−1/2k∏1≤i<j≤kvi−vj,V>0,πW∝1W1−1/2p∏1≤i<j≤pwi−wj,W>0,πΣ∝1Σ1−1/2q∏1≤i<j≤qσi−σj,Σ>0,where v1>⋯>vk>0 are the ordered eigenvalues of V, w1>⋯>wp>0 are the ordered eigenvalues of W, and σ1>⋯>σq>0 are the ordered eigenvalues of Σ.As discussed by Berger et al. [12], the prior πξ can be represented by the following hierarchical structure:(6)ξλ∼Nq0,λIq,λ∼λ−1/2e−1/2λ.For both intuitive and technical reasons, it is convenient to writeV=O1V∗O1′, W=O2W∗O2′, and Σ=O3Σ∗O3′, where V∗=diagv1,…,vk,W∗=diagw1,…,wp, Σ∗=diagσ1,…,σq, and Oi being the matrix of eigenvectors corresponding to V∗,W∗, and Σ∗, i=1,2,3, respectively. Consider the one-to-one transformation from V to V∗,O1 and rewrite the prior as(7)πVdV=πV∗,O1dV∗dO1,where dV=∏i≤jdvij,dV∗=∏i=1kdvi, and dO1 denotes the invariant Haar measure over the space of orthonormal matrices O=O1:O1′O1=Ik (see Anderson [33] for definition). From Farrell [34], the functional relationship between πV and πV∗,O1 is(8)πV∗,O1=πO1V∗O1′∏i<jvi−vj.Therefore, the prior forV becomes the prior of V∗,O1:(9)πV∗,O1∝1V∗1−1/2k.Use of the invariant prior onO1 (essentially a uniform prior over rotations) is natural and noncontroversial. This transformation reveals a significant difficulty of any prior that can be written as a function of V; in the V∗,O1 space, such priors contain the term ∏i<jvi−vj, which gives low mass to close eigenvalues and hence effectively force the eigenvalues apart. It is a criticism of inverse Wishart and Jeffreys priors. This is contrary to the common intuition, in that one often chooses a prior that pushes the eigenvalues closer together.Similarly, one can obtain the prior ofW∗,O2 and Σ∗,O3, given by(10)πW∗,O2∝1W∗1−1/2p,πΣ∗,O3∝1Σ∗1−1/2q.


## 3. Posterior Propriety

In this section, we consider the posterior propriety of the recommended prior in the 4-level normal hierarchical model (4).


### 3.1. The Caseq≥2

Theorem 1.
Consider the hierarchical Bayes model (4) with q≥2. Assume that Z has rank ps and ZT has rank nq. Then, the posterior distribution is always proper.Proof.
For the technical reason, we define(11)y=y1y2⋮ymmk×1,θ=θ1θ2⋮θmmk×1,1n=11⋮1n×1,Z=Z1Z2⋮Zmmk×sp,T=T1T2⋮Tssp×nq.
Then, (4) is equivalent to(12)Level1:yθ∼Nmkθ,Imk,Level2:θβ,V∼NmkZβ,Im⊗V,Level3:βη,W∼NspTη,Is⊗W,Level4:ηξ,Σ∼Nnq1n⊗ξ,In⊗Σ.
We marginalize out overθ, β, and η, one by one. It yields the marginal distribution of y given ξ,V,W,andΣ is(13)yξ,V,W,Σ∼NmkGξ,Δ,where Gmk×q=ZT1n⊗Iq and(14)Δ=Imk+Im⊗V+ZIs⊗WZ′+ZTIn⊗ΣT′Z′.
Note thatξλ∼Np0,λIq, and marginalizing out ξ yields the marginal distribution of y given V,W,Σ,andλ:(15)yV,W,Σ,λ∼Nmk0,Δ+λGG′.
Therefore, the likelihood ofy given V,W,Σ,andλ is(16)LV,W,Σ,λ∝1Δ+λGG′1/2exp−12y′Δ+λGG′−1y≤1Imk+Im⊗Vh1Imk+ZIs⊗WZ′h2Imk+ZTIn⊗ΣT′Z′h3Imk+λGG′h4,where h1=1/2mk+ϵ,h2=1/2sp+ϵ,h3=1/2−1/2mk−1/2sp−1/2q−3ϵ,h4=1−h1−h2−h3=1/2q+ϵ, and ϵ is a positive constant, assuring h3≥0. Note that h1+h2+h3+h4=1/2. Let my be the marginal density of y. Clearly,(17)my=∫∫∫∫LV,W,Σ,λπVπWπΣπλdVdWdΣdλ≤∫πVImk+Im⊗Vh1dV×∫πWImk+ZIs⊗WZ′h2dW×∫πλImk+λGG′h4dλ×∫πΣImk+ZTIn⊗ΣT′Z′h3dΣ.
Note that(18)∫πVImk+Im⊗Vh1dV=∫1v1>v2>⋯>vkIk+V1/2k+mϵV1−1/2k∏1≤i<j≤kvi−vjdV.
From relationship (8) between V and V∗,O1,(19)∫πVImk+Im⊗Vh4dV=∬1v1>v2>⋯>vkIk+V∗1/2k+mϵV∗1−1/2kdV∗dO1≤∏i=1k∫0∞11+vi1/2k+mϵvi1−1/2kdvi<∞.
Next, consider the integration overW. From the matrix identity A+XBX′=ABB−1+X′A−1X, we obtain(20)Imk+ZIs⊗WZ′=WsIs⊗W−1+Z′Z.
Letz0>0 be the smallest eigenvalue of Z′Z, then(21)Is⊗W−1+Z′Z≥W−1+z0Isps.
Thus,(22)Imk+ZIs⊗WZ′≥Ip+z0Ws.
Using the above inequality, it yields(23)∫πWImk+ZIs⊗WZ′h2dW≤∫1Ik+z0W∗1/2p+sϵW∗1−1/2pdW∗≤∏j=1p∫0∞11+z0wj1/2p+sϵwj1−1/2pdwi<∞.
Letu1≥⋯uq>0 be the eigenvalues of G′G. Using the matrix identity A+XBX′=ABB−1+X′A−1X again, it follows that(24)Imk+λGG′=λIq1λIq+G′G≥1+upλq.
Therefore,(25)∫1Imk+λGG′h4πλdλ≤∫0∞1λ1/21+upλ1/2+qϵe−2/λdλ<∞.
Last, to consider the integration overΣ. Let t0 be the smallest eigenvalue of T′Z′ZT. Then,(26)Imk+ZTIn⊗ΣT′Z′≥Iq+t0Σ∗n.
Thus,(27)∫πΣImk+ZTIn⊗ΣT′Z′h3dΣ≤∫1Ik+t0Σ∗n1/2−1/2mk−1/2sp−1/2q−3ϵΣ∗1−1/2qdΣ∗≤∏l=1q∫0∞11+t0σln1/2−1/2mk−1/2sp−1/2q−3ϵσl1−1/2qdσl,which is finite if(28)n12−12mk−12sp−12q−3ϵ+1−12q>1.
Sinceϵ can be chosen arbitrarily small, the integration over Σ is finite if(29)nq−1>1+nqmk+nqsp=1+1ms+1s.
Sinceq,n,m,ands≥2, (29) is true. The theorem is proved.


### 3.2. The Caseq=1

Theorem 2.
Consider the hierarchical Bayes model (4) with q=1. Assume that Z has rank ps and ZT has rank n. Then, the posterior distribution is proper if(30)1m+1s+2n+1sn+1msn<1.Proof.
In the caseq=1, note that Σ is just a variance, so that(31)Δ=Imk+Im⊗V+ZIs⊗WZ′+ΣZTT′Z′.
The prior ofξ becomes(32)πξ∝1.
Integrating outξ in (22) with its constant prior and dropping all exponential terms of the likelihood (as they are less than one), the marginal likelihood for V and W satisfies(33)LV,W<1Δ1/2G′Δ−1G1/2.
Thus, for the posterior propriety, we need to verify(34)∫∫1Δ1/2G′Δ−1G1/2πVπWπΣdΣdVW<∞.
Denotet1 be the largest eigenvalue of T′Z′ZT. Let λ0A and λ1A denote the minimum and maximum eigenvalue of A, respectively. Then,(35)λ1Δ≤λ11+v1Imk+w1ZZ′+ΣZTT′Z′≤1+v1+z1w1+t1Σ≤1+v11+z1w11+t1Σ.
Noting thatG is an mk×1 vector, it follows that(36)G′Δ−1G>λ0Δ−1G′G=λ1−1ΔG′G>11+v11+z1w11+t1ΣG′G.
Letb1=1/2mk+1/2m+ϵ and b2=1/2sp+1/2s+ϵ, where ϵ is a small positive constant, satisfying 1/2−b1−b2≥0. Then,(37)Δ1/2≥Imk+Im⊗Vb1Imk+ZIs⊗WZ′b2Imk+ΣZTT′Z′1/2−b1−b2≥Ik+Vmb1Ip+z0Wsb21+t0Σn1/2−b1−b2.
Clearly,(38)my≤C∭1+v11+z1w11+t1Σ1/2πVπWπΣIk+Vmb1Ip+z0Wsb21+t0Σn1/2−b1−b2dVdWdΣ,where C=G′G−1/2. Note that(39)∫1+v11/2Ik+Vmb1πVdV≤∫1+v11/2Ik+Vmb1V1−1/2kdV∗≤∫0∞11+v1mb1−1/2v11−1/2kdv1×∏i=2k∫0∞11+vimb1vi1−1/2kdvi,which is finite since mb1−1/2+1−1/2k=1+mϵ>1. Similarly,(40)∫1+z1w11/2Ik+Wsb2πWdW≤∫0∞1+z1w11/21+w1sb2w11−1/2pdw1×∏j=2p∫0∞11+wjsb2wj1−1/2pdwj,which is finite since sb2+1−1/2p−1/2=1+sϵ>1. Next, consider the integration over Σ. It follows that(41)∫0∞1+t1Σ1/21+t0Σn1/2−b1−b2Σ1/2dΣ<∞,if n1/2−b1−b2=n1/2−1/2mk−1/2m−1/2sp−1/2s−2ϵ>1. Since ϵ can be chosen arbitrarily small, p=n and k=sn; this is true if(42)1m+1s+2n+1sn+1msn<1.
This completes the proof.


## 3.1. The Caseq≥2

Theorem 1.
Consider the hierarchical Bayes model (4) with q≥2. Assume that Z has rank ps and ZT has rank nq. Then, the posterior distribution is always proper.Proof.
For the technical reason, we define(11)y=y1y2⋮ymmk×1,θ=θ1θ2⋮θmmk×1,1n=11⋮1n×1,Z=Z1Z2⋮Zmmk×sp,T=T1T2⋮Tssp×nq.
Then, (4) is equivalent to(12)Level1:yθ∼Nmkθ,Imk,Level2:θβ,V∼NmkZβ,Im⊗V,Level3:βη,W∼NspTη,Is⊗W,Level4:ηξ,Σ∼Nnq1n⊗ξ,In⊗Σ.
We marginalize out overθ, β, and η, one by one. It yields the marginal distribution of y given ξ,V,W,andΣ is(13)yξ,V,W,Σ∼NmkGξ,Δ,where Gmk×q=ZT1n⊗Iq and(14)Δ=Imk+Im⊗V+ZIs⊗WZ′+ZTIn⊗ΣT′Z′.
Note thatξλ∼Np0,λIq, and marginalizing out ξ yields the marginal distribution of y given V,W,Σ,andλ:(15)yV,W,Σ,λ∼Nmk0,Δ+λGG′.
Therefore, the likelihood ofy given V,W,Σ,andλ is(16)LV,W,Σ,λ∝1Δ+λGG′1/2exp−12y′Δ+λGG′−1y≤1Imk+Im⊗Vh1Imk+ZIs⊗WZ′h2Imk+ZTIn⊗ΣT′Z′h3Imk+λGG′h4,where h1=1/2mk+ϵ,h2=1/2sp+ϵ,h3=1/2−1/2mk−1/2sp−1/2q−3ϵ,h4=1−h1−h2−h3=1/2q+ϵ, and ϵ is a positive constant, assuring h3≥0. Note that h1+h2+h3+h4=1/2. Let my be the marginal density of y. Clearly,(17)my=∫∫∫∫LV,W,Σ,λπVπWπΣπλdVdWdΣdλ≤∫πVImk+Im⊗Vh1dV×∫πWImk+ZIs⊗WZ′h2dW×∫πλImk+λGG′h4dλ×∫πΣImk+ZTIn⊗ΣT′Z′h3dΣ.
Note that(18)∫πVImk+Im⊗Vh1dV=∫1v1>v2>⋯>vkIk+V1/2k+mϵV1−1/2k∏1≤i<j≤kvi−vjdV.
From relationship (8) between V and V∗,O1,(19)∫πVImk+Im⊗Vh4dV=∬1v1>v2>⋯>vkIk+V∗1/2k+mϵV∗1−1/2kdV∗dO1≤∏i=1k∫0∞11+vi1/2k+mϵvi1−1/2kdvi<∞.
Next, consider the integration overW. From the matrix identity A+XBX′=ABB−1+X′A−1X, we obtain(20)Imk+ZIs⊗WZ′=WsIs⊗W−1+Z′Z.
Letz0>0 be the smallest eigenvalue of Z′Z, then(21)Is⊗W−1+Z′Z≥W−1+z0Isps.
Thus,(22)Imk+ZIs⊗WZ′≥Ip+z0Ws.
Using the above inequality, it yields(23)∫πWImk+ZIs⊗WZ′h2dW≤∫1Ik+z0W∗1/2p+sϵW∗1−1/2pdW∗≤∏j=1p∫0∞11+z0wj1/2p+sϵwj1−1/2pdwi<∞.
Letu1≥⋯uq>0 be the eigenvalues of G′G. Using the matrix identity A+XBX′=ABB−1+X′A−1X again, it follows that(24)Imk+λGG′=λIq1λIq+G′G≥1+upλq.
Therefore,(25)∫1Imk+λGG′h4πλdλ≤∫0∞1λ1/21+upλ1/2+qϵe−2/λdλ<∞.
Last, to consider the integration overΣ. Let t0 be the smallest eigenvalue of T′Z′ZT. Then,(26)Imk+ZTIn⊗ΣT′Z′≥Iq+t0Σ∗n.
Thus,(27)∫πΣImk+ZTIn⊗ΣT′Z′h3dΣ≤∫1Ik+t0Σ∗n1/2−1/2mk−1/2sp−1/2q−3ϵΣ∗1−1/2qdΣ∗≤∏l=1q∫0∞11+t0σln1/2−1/2mk−1/2sp−1/2q−3ϵσl1−1/2qdσl,which is finite if(28)n12−12mk−12sp−12q−3ϵ+1−12q>1.
Sinceϵ can be chosen arbitrarily small, the integration over Σ is finite if(29)nq−1>1+nqmk+nqsp=1+1ms+1s.
Sinceq,n,m,ands≥2, (29) is true. The theorem is proved.


## 3.2. The Caseq=1

Theorem 2.
Consider the hierarchical Bayes model (4) with q=1. Assume that Z has rank ps and ZT has rank n. Then, the posterior distribution is proper if(30)1m+1s+2n+1sn+1msn<1.Proof.
In the caseq=1, note that Σ is just a variance, so that(31)Δ=Imk+Im⊗V+ZIs⊗WZ′+ΣZTT′Z′.
The prior ofξ becomes(32)πξ∝1.
Integrating outξ in (22) with its constant prior and dropping all exponential terms of the likelihood (as they are less than one), the marginal likelihood for V and W satisfies(33)LV,W<1Δ1/2G′Δ−1G1/2.
Thus, for the posterior propriety, we need to verify(34)∫∫1Δ1/2G′Δ−1G1/2πVπWπΣdΣdVW<∞.
Denotet1 be the largest eigenvalue of T′Z′ZT. Let λ0A and λ1A denote the minimum and maximum eigenvalue of A, respectively. Then,(35)λ1Δ≤λ11+v1Imk+w1ZZ′+ΣZTT′Z′≤1+v1+z1w1+t1Σ≤1+v11+z1w11+t1Σ.
Noting thatG is an mk×1 vector, it follows that(36)G′Δ−1G>λ0Δ−1G′G=λ1−1ΔG′G>11+v11+z1w11+t1ΣG′G.
Letb1=1/2mk+1/2m+ϵ and b2=1/2sp+1/2s+ϵ, where ϵ is a small positive constant, satisfying 1/2−b1−b2≥0. Then,(37)Δ1/2≥Imk+Im⊗Vb1Imk+ZIs⊗WZ′b2Imk+ΣZTT′Z′1/2−b1−b2≥Ik+Vmb1Ip+z0Wsb21+t0Σn1/2−b1−b2.
Clearly,(38)my≤C∭1+v11+z1w11+t1Σ1/2πVπWπΣIk+Vmb1Ip+z0Wsb21+t0Σn1/2−b1−b2dVdWdΣ,where C=G′G−1/2. Note that(39)∫1+v11/2Ik+Vmb1πVdV≤∫1+v11/2Ik+Vmb1V1−1/2kdV∗≤∫0∞11+v1mb1−1/2v11−1/2kdv1×∏i=2k∫0∞11+vimb1vi1−1/2kdvi,which is finite since mb1−1/2+1−1/2k=1+mϵ>1. Similarly,(40)∫1+z1w11/2Ik+Wsb2πWdW≤∫0∞1+z1w11/21+w1sb2w11−1/2pdw1×∏j=2p∫0∞11+wjsb2wj1−1/2pdwj,which is finite since sb2+1−1/2p−1/2=1+sϵ>1. Next, consider the integration over Σ. It follows that(41)∫0∞1+t1Σ1/21+t0Σn1/2−b1−b2Σ1/2dΣ<∞,if n1/2−b1−b2=n1/2−1/2mk−1/2m−1/2sp−1/2s−2ϵ>1. Since ϵ can be chosen arbitrarily small, p=n and k=sn; this is true if(42)1m+1s+2n+1sn+1msn<1.
This completes the proof.


## 4. Computation

In this section, we consider the MCMC computation from the posterior arising from the model (4). The normal hierarchical models are typically handled today by Gibbs sampling. One difficulty of computation is to sample the covariance matrix efficiently. The main new development discussed in Berger et al. [35] is a new and efficient computational algorithm for dealing with priors on covariance matrices as in (3). It overcomes the computational bottleneck mentioned in the introduction.Fact 1.
Here are the full conditional distributions:(a)
Fori=1,…,m, the conditional posterior of θi given β and V is the usual conjugate posterior density:(43)NkIk+V−1−1yi+V−1Ziβ,Ik+V−1−1.(b)
Define(44)M=Is⊗W−1+Z′Z−1Z′Im⊗VZZ′Z−1−1,u=Is⊗W−1Tη+Z′Z−1Z′Im⊗VZZ′Z−1−1Z′Z−1Z′θ.The conditional posterior ofβ given θ,η,V,W is NspM−1u,M−1.(c)
Definingr1=m/2+1−1/3k, the conditional posterior density of V for given θ,β can be written as(45)πVθ,β∝1Vr1∏i<jvi−vjetr−12V−1H1.Here and in the followingetrA represents exptraceA for a squared matrix A, and(46)H1=H1θ,β=∑i=1mθi−Zi′βθi−Zi′β′.(d)
Define(47)K=In⊗Σ−1+T′T−1T′Is⊗WTT′T−1−1,ν=1n⊗Σ−1ξ+T′T−1T′Is⊗WTT′T−1−1T′T−1T′β.The conditional posterior ofη given β,ξ,W,Σ is NnqK−1ν,K−1.(e)
Definingr2=s/2+1−1/2p, the conditional posterior density of W for given β,η can be written as(48)πWβ,η∝1Wr2∏i<jwi−wjetr−12W−1H2,whereH2=∑j=1sβj−Tj′ηβj−Tj′η′.(f)
To compute with the recommended prior forξ, use the equivalent representation as follows:(49)ξλ∼Nq0,λIq,λ∝λ−1/2e−1/2λ.(50)λξ∝1λq+1/2exp−1+ξ22λ.(51)Nq1λIq+nΣ−1−1Σ−1∑l=1nηi,1λIq+nΣ−1−1.(g)
Definingr3=n/2+1−1/2q, the conditional posterior density of Σ for given β,η can be written as(52)πΣβ,η∝1Σr3∏i<jσi−σjetr−12Σ−1H3,whereH3=∑l=1nηl−ξηl−ξ′.
(i) For sampleλ from its full conditional, the Inverse Gamma q−1/2,1+ξ2/2 density
Ifq=1, this step is not needed as the hyperprior, for ξ is constant.
(ii) Givenλ,Σ,ηi, Gibbs sampling of ξ can be done from its full conditional, which is (when q=1, set λ=∞)
Sampling ofθ, β, η, and ξ can simply be carried out with a Gibbs step, as its full conditional will be a normal distribution. To sample the covariance matrix V, W, and Σ from Fact 1, the Metropolis–Hastings Algorithm [36] and Hit-and-Run Method Chen and Schmeiser [37] could be used, based on proposal distributions that generate full-candidate matrices. But the two well-known methods are both inefficient, especially for high-dimensional data. Fortunately, Berger et al. [35] proposed a new and efficient computational algorithm for sampling the covariance matrix V, W, and Σ from Fact 1, respectively.
Take samplingV from (45) as an example. For the conditional density of V given in (45), we use the eigenvalue-vector decomposition V=O1V∗O1′, where O1 is the orthogonal and V∗ is the diagonal matrix of ordered eigenvalues. For H1 defined in (46), note that trV−1H1=trO1V∗−1O1′H1=trV∗−1O1′H1O1, so the conditional density of V given in (45) can be transformed to(53)πV∗,O1θ,β∝1V∗r1etr−12V∗−1O1′H1O11v1>v2>⋯>vk.
Gibbs sampling ofV∗: following Lemma 1 in Berger et al. [12], i.e.,(54)∬gO1V∗O1′1v1>⋯>vkdV∗dO1=1k!∬gO1V∗O1′dV∗dO1,we can first sample V∗ given O1,H1, from(55)πV∗O1,θ,β∝1∏i=1kvir1etr−12V∗−1O1′H1O1=∏i=1k1vir1exp−civi,where ci is the i,i element of O1′H1O1/2. Thus, we can directly sample vi independently from inverse gamma r1−1,ci distributions, and then simply rearrange the vi, so that v1≥⋯≥vk.
Gibbs sampling ofO1: from (19), the conditional density of O1 given V∗ and H1 is(56)πO1V∗;θ,β∝etr−12H1O1V∗−1O1′.
WriteH1=LUL′, where L is the matrix of normalized eigenvectors LL′=Ik and U=diagu1,…,uk is the diagonal matrix of corresponding eigenvalues with u1≥⋯≥uk. Denote Γ=L′O1. Then, the conditional density of Γ given V∗ and H1 is(57)πΓV∗;θ,β∝etr−12UΓV∗−1Γ′.
Hoff [31] introduced an MCMC algorithm to sample from the posterior of O1 by updating two randomly selected columns of O1. Alternatively, Berger et al. [35] suggested updating two randomly selected rows of Γ (essentially equivalent to Hoff’s method in full rank cases, but considerably faster in situations of less than full rank). For example, to update the first two rows of Γ, we write Γ=diagΩ,Ik−2Γ12′,Γ−12′′, where Γ12 is the first 2 rows of Γ, Γ−12 is the rest k−2 rows of Γ, and(58)Ω=DϵΩ0=ϵ100ϵ2cosω−sinωsinωcosω.
Here,ω∈−π/2,π/2 and ϵi=±1 for i=1,2. Write U1=diagu1,u2 and U2=diagu3,…,uk. Then, the conditional posterior of ω, πωΓ12,Γ−12,V∗−1;H1 is(59)∝etr−12U100U2Ω00Ik−2Γ12Γ−12V∗−1Γ12′,Γ−12′Ω′00Ik−2∝etr−12U1Ω0Γ12V∗−1Γ12′Ω0′.
Write(60)Γ12V∗−1Γ12′=cosϕ−sinϕsinϕcosϕs100s2cosϕsinϕ−sinϕcosϕ,where ϕ∈−π/2,π/2 and s1>s2. From (59), the conditional posterior of θ is(61)πθΓ−12,V∗;θ,β∝expc0cos2ω+ϕ,θ∈−π/2,π/2,where c0=−1/2s1−s2u1−u2≤0. Let α=cos2ω+ϕ. Then, the full conditional density of α is(62)παΓ−12,V∗;θ,β∝ec0αα−1/21−α−1/2,α∈0,1.
Samplingα∈0,1 can be done with a rejection sampler with proposal Beta1/2,1/2 distribution.
From the studies in Berger et al. [35], the new method substantially outperforms the Metropolis and Hit-and-Run algorithms in moderate dimensions and succeeds for k up to 100, whereas the other methods break down in much lower dimensions.


## 5. Simulation

From the mean square error (MSE) perspective, we compare the performance of 6 objective hyperpriors (considered in Berger et al. [12]), created from the product of three objective hyperpriors for ξ:(1)
Constant prior:π1ξ=1(2)
Conjugate prior:π2ξ∼Nk0,Iq(3)
Recommended prior:π3ξ∝1+ξ2−q−1/2 and two objective hyperpriors for V(1)
Constant prior:πV=1,πW=1 and πΣ=1(2)
Recommended reference prior:(63)πV∝1V1−1/2k∏1≤i<j≤kvi−vj,V>0,πW∝1W1−1/2p∏1≤i<j≤pwi−wj,W>0,πΣ∝1Σ1−1/2q∏1≤i<j≤qσi−σj,Σ>0.Except for the constant prior and recommended reference prior, Berger et al. [12] also studied the nonhierarchical independence Jeffreys prior (πV=V−k+1/2) and the hierarchical independence Jeffreys prior (πV=I+V−k+1/2) in the model (1). From Berger et al. [12], the nonhierarchical independence Jeffreys prior for the covariance matrix cannot be used for the hierarchical models since it results in improper posterior. For the hierarchical independence Jeffreys prior, the common sampling methods are very inefficient to sample V from the posterior distribution. Therefore, we do not consider either of the abovementioned two Jeffreys priors.Setq=n=2,p=s=4,m=k=16. We generate Z=z1,…,zmk′, where zi∼iidNsp0,Isp for i=1,…,mk. Similarly, T=t1,…,tsp′, where tj∼iidNnq0,Inq for j=1,…,sp. We simulate the Bayes risks of the posterior means resulting from the 6 priors in every combination of the following cases:(i)
ξ1=1q or ξ2=301q;(ii)
V1=Ik,W1=Ip,Σ1=Iq or (V2=diag4k−3,…,5,1, W2=diag4p−3,…,5,1, Σ2=diag4q−3,…,5,1).Each choice ofξ,V,W,andΣ specifies a “true” hierarchical model for the simulation and we wish to compute the Bayes risk corresponding to the posterior means for θ that arise from each of the 6 objective priors. To simulate these risks, we generate 2000 random θl∼Nmk1m⊗β,Im⊗V and generate observations ylθl∼Nmkθl,Imk, l=1,…,2000. To obtain the posterior means, δπyl, of θl, under the 6 priors, we run 100,000 MCMC cycles after 3,000 burn-in cycles, using the algorithms described above. Finally, we approximate the Bayes risk as the average observed mean squared error (MSE) as follows:(64)MSE=12000∑l=12000θl−δπyl′θl−δπyl.The results are given in Table1. The recommended prior (in the last line) dominates all the others in terms of risk. Within each three-row segment, comparison of the first row with the third row is a comparison between using the constant prior for ξ and the recommended prior. The gains with the recommended prior are large. Comparison of each of the first three rows with each of the last three rows is a comparison between using the constant prior for V versus the recommended prior. The gains with the recommended prior are also large. The risk results present the strong numerical evidence of the superiority of the recommended priors (5).Table 1
The Bayes risk.
Prior forV,W,andΣ
Prior forξ
ξ
1
,
V
1
ξ
1
,
V
2
ξ
2
,
V
1
ξ
2
,
V
2
Constant prior
π
1
ξ
234.280
301.450
244.154
312.565
π
2
ξ
219.234
284.716
253.791
309.367
π
3
ξ
214.406
267.673
223.430
275.486
Reference prior
π
1
ξ
208.342
280.193
229.257
286.841
π
2
ξ
201.971
264.365
243.239
294.634
π
3
ξ
197.123
255.153
208.372
265.474


## 6. Conclusion

In this paper, we follow Berger et al.’s [12] work and study a 4-level normal hierarchical model. We demonstrate that the posterior using the recommended prior is still proper in the 4-level normal hierarchical model. Herein, we do not demonstrate Berger et al.’s [32] conjecture, i.e., the posterior using the recommended prior is always proper in any level normal hierarchical models. But our method provides one useful guideline for the completion of the story.Besides, the normal hierarchical models are typically handled by Gibbs sampling. One difficulty of computation is to sample the covariance matrix efficiently. The common sampling methods for covariance matrices, for example, Metropolis–Hastings algorithm [36] and Hit-and-Run method [25], are inefficient for higher dimensions. To overcome the computational bottleneck of sampling covariance matrices, we can use a powerful and efficient method proposed by Berger et al. [35] for sampling the conditional density of covariance matrices. Therefore, there is no difficulty in computation for the 4-level normal hierarchical model with the recommended priors (5). In addition, the simulation result presents the strong numerical evidence of the superiority of the recommended priors (5).



---

*Source: 8236934-2020-02-14.xml*