arXiv:1712.00346v3 [] 13 Mar 2018

Bayes Minimax Competitors of Preliminary Test Estimators in k
Sample Problems
Ryo Imai∗, Tatsuya Kubokawa†and Malay Ghosh‡
March 14, 2018

Abstract
In this paper, we consider the estimation of a mean vector of a multivariate normal
population where the mean vector is suspected to be nearly equal to mean vectors of k − 1
other populations. As an alternative to the preliminary test estimator based on the test
statistic for testing hypothesis of equal means, we derive empirical and hierarchical Bayes
estimators which shrink the sample mean vector toward a pooled mean estimator given
under the hypothesis. The minimaxity of those Bayesian estimators are shown, and their
performances are investigated by simulation.
Key words and phrases: Admissibility, decision theory, empirical Bayes, hierarchical
Bayes, k sample problem, minimaxity, pooled estimator, preliminary test estimator, quadratic
loss, shrinkage estimator, uniform prior.

1

Introduction

Suppose that there are k laboratories, say Laboratory L1 , . . ., Lk , and a certain instrument is
designed to measure several characteristics at each laboratory and several vector-valued measurements are recorded. Also, suppose that we want to estimate the population mean of Laboratory
L1 . When similar instruments are used at k laboratories, it is suspected that k population means
are nearly equal, in which case, the sample means of the other laboratories L2 , . . ., Lk are used
to produce more efficient estimators than the sample mean based on data from only L1 . This
problem was studies by Ghosh and Sinha (1988) and recently revisited by Imai, Kubokawa and
Ghosh (2017) in the framework of simultaneous estimation of k population means.
The k sample problem described above is expressed in the following canonical form: p-variate
random vectors X 1 , . . . , X k and positive scalar random variable S are mutually independently
distributed as
X i ∼Np (µi , σ 2 V i ),
S/σ

2

∼χ2n ,

for i = 1, . . . , k,

(1.1)

∗
Graduate School of Economics, University of Tokyo, 7-3-1 Hongo, Bunkyo-ku, Tokyo 113-0033, JAPAN.
E-Mail: imronuronnar13@g.ecc.u-tokyo.ac.jp
†
Faculty of Economics, University of Tokyo, 7-3-1 Hongo, Bunkyo-ku, Tokyo 113-0033, JAPAN.
E-Mail: tatsuya@e.u-tokyo.ac.jp
‡
Department of Statistics, University of Florida, 102 Griffin-Floyd Hall, Gainesville, Florida 32611.
E-Mail: ghoshm@stat.ufl.edu

1

where the p-variate means µ1 , . . . , µk and the scale parameter σ 2 are unknown, and V 1 , . . . , V k
are p×p known and positive definite matrices. In this model, we consider to estimate µ1 relative
to the quadratic loss function
L(δ 1 , µ1 , σ 2 ) = kδ 1 − µ1 k2Q /σ 2 = (δ 1 − µ1 )⊤ Q(δ 1 − µ1 )/σ 2 ,

(1.2)

where kak2A = a⊤ Aa for the transpose a⊤ of a, Q is a positive definite and known matrix, and δ 1
is an estimator of µ1 . Estimator δ 1 is evaluated by the risk function R(ω, δ 1 ) = E[L(δ 1 , µ1 , σ 2 )]
for ω = (µ1 , . . . , µk , σ 2 ), a set of unknown parameters.
In this paper, we consider the case that the means µi ’s are suspected to be nearly equal,
namely close to the hypothesis
H0 : µ1 = · · · = µk .
A classical approach towards solution to this problem is the development of a preliminary test
estimator which uses the pooled mean estimator upon acceptance of the null hypothesis H0 , and
uses separate mean estimators upon rejection of H0 . The test statistic for H0 is
F =

k
X
i=1

b k2V −1 /S,
kX i − ν

b is the pooled estimator defined as
where ν
b=A
ν

k
X

V

(1.3)

i

−1
i X i,

A=

k
X
i=1

i=1

V −1
i

−1

.

The preliminary test estimator of µ1 is

X 1 if F > (p(k − 1)/n)Fp(k−1),n,α
PT
b1 =
µ
b
ν
otherwise

(1.4)

(1.5)

where Fp(k−1),n,α is the upper α point of the F distribution with (p(k − 1), n) degrees of freedom,
b is the pooled estimator given in (1.4). However, the preliminary test estimator is not
and ν
smooth and does not necessarily improve on X 1 .

As an alternative approach to the preliminary test estimator, we consider a Bayesian method
under the prior distributions of µi having a common mean ν as a hyper-parameter. Empirical
and hierarchical Bayes estimators are derived when the uniform prior distribution is assumed
for ν. We also provide empirical Bayes estimator under the assumption of the normal prior
distribution for ν. It is shown that these Bayesian estimators improve on X 1 , namely, are
minimax.

The topic treated in this paper is related to the so-called Stein problem in simultaneous
estimation of multivariate normal means. This problem has long been of interest in the literature.
See, for example, Stein (1956, 1981), James and Stein (1961), Strawderman (1971, 1973), Efron
and Morris (1973, 1976) and Berger (1985). For recent articles exrending to prediction and
high-dimensional problems, see Komaki (2001), Brown, George and Xu (2008) and Tsukuma
and Kubokawa (2015). As articles related to this paper, Sclove, Morris and Radhakrishnan
(1972) showed the inadmissibility of the preliminary test estimator, Smith (1973) provided
2

Bayes estimators in one-way and two-way models, Ghosh and Sinha (1988) derived hierarchical
and empirical Bayes estimators with minimaxity, and Sun (1996) provided Bayesian minimax
estimators in multivariate two-way random effects models.
In Section 2, we treat two classes of shrinkage estimators, say Class 1 and Class 2. Estimators
b , and estimators in Class 2 incorporate a
in Class 1 shrink X 1 toward the pooled estimator ν
b in addition to estimators in Class 1. Ghosh and Sinha (1988) obtained
part of shrinking ν
conditions for minimaxity of estimators in Class 1 in the two sample problem. For k = 2, the
test statistic F can be simply
(X 1 − X 2 )⊤ (V 1 + V 2 )−1 (X 1 − X 2 )/S, while F is
Pk written as⊤F =
b ) V −1
b )/S in the k sample case. This means that
more complicated as F = i=1 (X i − ν
i (X i − ν
the minimaxity is harder to establish in the k sample problem. The key tool in the extension is
the inequality given in Lemma 2.2. Using the inequality, we obtain conditions for minimaxity
of shrinkage estimators in Class 1 and Class 2.
In Section 3, we suggest three kinds of Bayesian methods for estimation of µ1 . Empirical and
hierarchical Bayes estimators are derived under the assumption of the uniform prior distribution
for ν. Because these estimators belong to Class 1, we can provide conditions for their minimaxity
using the result in Section 2. Also empirical Bayes estimator is obtained under the normal prior
distribution for ν and belongs to Class 2. Conditions for the minimaxity are given from the result
in Section 2. The performances of those Bayesian procedures are investigated by simulation in
Section 5. The results in Section 2 are extended in Section 4 to the problem of estimating a
linear combination of µ1 , . . . , µk under the quadratic loss. Concluding remarks are given in
Section 6.

2

Classes of Minimax Estimators

Motivated from the preliminary test estimator (1.5), we consider a class of estimators shrinking
b , given by
X 1 toward the pooled estimator ν
b 1 (φ) = X 1 −
µ

φ(F, S)
b ),
(X 1 − ν
F

(2.1)

b
where φ(F, S) is an absolutely continuous function, F is the test statistic given in (1.3) and ν
b given in (1.4).
is the pooled estimator ν
b 1 (φ) in Theorem 2.1, which is proved
We can derive conditions on φ(F, S) for minimaxity of µ
in the end of this section.
Theorem 2.1 Assume the conditions
Chmax {(V 1 − A)Q} =
6 0

and

tr {(V 1 − A)Q}/Chmax {(V 1 − A)Q} > 2.

(2.2)

b 1 (φ) is minimax relative to the quadratic loss (1.2) if φ(F, S) satisfies the
Then, the estimator µ
following conditions:
(a) φ(F, S) is non-decreasing in F and non-increasing in S.
(b) 0 < φ(F, S) ≤ 2[tr {(V 1 − A)Q}/Chmax {(V 1 − A)Q} − 2]/(n + 2), where Chmax (C)
denotes the maximum characteristic value of matrix C.
This theorem provides an extension of Ghosh and Sinha (1988) to the k-sample problem. In
the case of V 1 = · · · = V k = Q−1 , the conditions in (2.2) are expressed as k 6= 1 and p > 2,
and the latter condition is well known in the Stein problem.
3

We next consider the class of double shrinkage estimators
b 1 (φ, ψ) = X 1 −
µ

φ(F, S)
ψ(G, S)
b) −
b,
(X 1 − ν
ν
F
G

(2.3)

where φ(F, S) and ψ(G, S) are absolutely continuous functions, and
G = kb
ν k2A−1 /S.
Theorem 2.2 Assume condition (2.2) and
Chmax (AQ) 6= 0

and

tr (AQ)/Chmax (AQ) > 2.

(2.4)

b 1 (φ, ψ) in (2.3) is minimax relative to the quadratic loss
Then, the double shrinkage estimator µ
(1.2) if φ(F, S) and ψ(G, S) satisfy the following conditions:
(a) φ(F, S) is non-decreasing in F and non-increasing in S.
(b) 0 < φ(F, S) ≤ [tr {(V 1 − A)Q}/Chmax {(V 1 − A)Q} − 2]/(n + 2).
(c) ψ(G, S) is non-decreasing in G and non-increasing in S.
(d) 0 < ψ(G, S) ≤ {tr (AQ)/Chmax (AQ) − 2}/(n + 2).
For the proofs, the Stein identity due to Stein (1981) and the chi-square identity due to Efron
and Morris (1976) are useful. See also Bilodeau and Kariya (1989) for a multivariate version of
the Stein identity.
Lemma 2.1 (1) Assume that Y = (Y1 , . . . , Yp )⊤ is a p-variate random vector having Np (µ, Σ)
and that h(·) is an absolutely continuous function from Rp to Rp . Then, the Stein identity is
given by
 

E[(Y − µ)⊤ h(Y )] = E tr Σ∇Y h(Y )⊤ ,
(2.5)

provided the expectations in both sides exist, where ∇Y = (∂/∂Y1 , . . . , ∂/∂Yp )⊤ .

(2) Assume that S is a random variable such that S/σ 2 ∼ χ2n and that g(·) is an absolutely
continuous function from R to R. Then, the chi-square identity is given by
E[Sg(S)] = σ 2 E[ng(S) + 2Sg ′ (S)],

(2.6)

provided the expectations in both sides exist.
Proof of Theorem 2.1 The risk function is decomposed as
i
h
1 h φ2
φ i
b k2Q
b 1 (φ)) =E[kX 1 − µ1 k2Q /σ 2 ] − 2E (X 1 − µ1 )⊤ Q(X 1 − ν
b) 2
+ 2 E 2 kX 1 − ν
R(ω, µ
σ F
σ
F
=I1 − 2I2 + I3 ,
(say)
(2.7)
for φ = φ(F, S) and ψ = ψ(G, S).
P
−1
It is easy to see I1 = tr (V 1 Q). Note that F is expressed as F = ( ki=1 X ⊤
i V i Xi −
⊤ −1
−1
b A ν
b )/S. Letting ∇1 = ∂/∂X 1 , we have ∇1 F = 2V 1 (X 1 − ν
b )/S and
ν
n
n
φ
φ
φF o 1
φ(F, S) o
−1
⊤
b
b
b )⊤
= (I − V −1
A)
−
+
2V
(X
−
ν
)(X
−
ν
)
+
,
∇1 (X 1 − ν
1
1
1
1
F
F
F2
F S
4

where φF = (∂/∂F )φ(F, S). Then from (2.5), it is seen that
n
h h
φ(F, S) oii
b )⊤
I2 =E tr QV 1 ∇1 (X 1 − ν
F
h
b )(X 1 − ν
b )⊤ } n
φ
tr {QV 1 V −1
φF oi
φ
1 (X 1 − ν
−
+
2
+
=E tr {QV 1 (I − V −1
A)}
1
F
S
F2
F
h
i
φ
φ
=E tr {(V 1 − A)Q} − 2B + 2BφF ,
(2.8)
F
F
where

k
X
b )⊤ V −1
b ).
b ) Q(X 1 − ν
b )/
(X j − ν
B = (X 1 − ν
j (X j − ν
⊤

(2.9)

j=1

Also, from (2.6), it is observed that

i
1 hS
2
Bφ
E
σ2
F
hn
 F n φ2
φφS i
φφF o
=E
− 2 +2
+ 2SB
Bφ2 + 2SB −
F
S
F
F
F
h n+2
φφS i
2
=E B
.
φ − 4BφφF + 4SB
F
F

I3 =

(2.10)

Thus, the risk function is expressed as
h
φ
b 1 (φ)) =E tr (V 1 Q) + [(n + 2)Bφ − 2tr {(V 1 − A)Q} + 4B]
R(ω, µ
F
φφS i
.
− 4BφF − 4BφφF + 4SB
F

(2.11)

Using Lemma 2.2 given below, we can see that
B≤
which implies that

b )⊤ Q(X 1 − ν
b)
(X 1 − ν
≤ Chmax ((V 1 − A)Q),
⊤
−1
b ) (V 1 − A) (X 1 − ν
b)
(X 1 − ν

b 1 (φ)) − tr (V 1 Q)
R(ω, µ
o
h
tr {(V 1 − A)Q}
φn
(n + 2)φ − 2
+4
≤E Chmax ((V 1 − A)Q)
F
Chmax ((V 1 − A)Q)
i
φφS
.
− 4BφF − 4BφφF + 4SB
F
b 1 (φ)) ≤ tr (V 1 Q) under the conditions in Theorem 2.1.
Hence, R(ω, µ

(2.12)


Lemma 2.2 It holds that
k
X
j=1

b )⊤ V −1
b ) ≥ (x1 − ν
b )⊤ (V 1 − A)−1 (x1 − ν
b ).
(xj − ν
j (xj − ν
5

(2.13)

P
b = AV −1
Proof. Let C = (V 1 − A)−1 and x∗ = kj=2 V −1
1 x1 + Ax∗ and
j xj . Then, ν
−1 −1
b = C V 1 x1 − Ax∗ . The RHS of (2.13) is rewritten as
x1 − ν
−1 −1 −1
−1
⊤
b )⊤ C(x1 − ν
b ) = x⊤
(x1 − ν
V 1 x1 − 2x⊤
1V1 C
1 V 1 Ax∗ + x∗ ACAx∗ .

(2.14)

On the other hand, it can be observed that

−1 −1 −1 −1 −1
−1 −1 −1
−1
b )⊤ V −1
b ) =x⊤
(x1 − ν
V 1 C V 1 x1 − 2x⊤
V 1 Ax∗ + x⊤
1V1 C
1V1 C
∗ AV 1 Ax∗ ,
1 (x1 − ν

k
k
X
X
−1
−1
⊤
⊤
⊤ −1
b ) V j (xj − ν
b) =
x⊤
(xj − ν
j V j xj − 2x∗ AV 1 x1 − 2x∗ Ax∗
j=2

j=2

−1
−1
−1
−1
−1
+ x⊤
− V −1
− V −1
1 V 1 A(A
1 )AV 1 x1 + 2x1 V 1 A(A
1 )Ax∗

−1
+ x⊤
− V −1
∗ A(A
1 )Ax∗ ,

which gives
k
X
b)
b )⊤ V −1
(xj − ν
j (xj − ν
j=1

−1
−1
−1 −1 −1
V 1 C + A(A−1 − V −1
=x⊤
1 V 1 {C
1 )A}V 1 x1
−1
−1 −1
V 1 + I − A(A−1 − V −1
− 2x⊤
1 V 1 {C
1 )}Ax∗

− x⊤
∗ Ax∗ +

k
X

−1
x⊤
j V j xj .

(2.15)

j=2

−1
−1 − V −1 ) = I
It is noted that C −1 V −1
+ A(A−1 − V −1
A, C −1 V −1
1 C
1 )A = V 1 −P
1 + I − A(|A
1
k
−1 −1
−1 −1
−1
−1
and ACA + A = A(V 1 − A) V 1 = (A − V 1 ) = ( j=2 V j ) . From (2.14) and (2.15),
the inequality (2.13) is equivalent to
k
X

−1
x⊤
j V j xj

j=2

≥

k
X

V

−1
j xj

k
⊤  X
j=2

j=2

V

−1
j

k
−1  X
j=2


V −1
j xj .

(2.16)

P
Pk
−1
−1
To show the inequality (2.16), let ν ∗ = ( kj=2 V −1
j ) ( j=2 V j xj ). Then, it can be seen that
k
k
k
k

−1  X
⊤  X
X
X
−1
−1
−1
−1
x
V
V
x
V
x
−
x⊤
V
j
j
j
j
j
j
j
j

=

k
X
j=2

=

j=2

j=2

j=2

−1
⊤
x⊤
j V j xj − ν ∗

k
X
j=2

j=2


ν∗
V −1
j

k
X

(xj − ν ∗ )⊤ V −1
j (xj − ν ∗ ),

j=2

which is nonnegative, and the proof of Lemma 2.2 is complete.

6



b 1 (φ, ψ) is
Proof of Theorem 2.2 The risk function of µ

1 h
ψi
b 1 (φ, ψ)) =R(ω, µ
b 1 (φ)) − 2 2 E (X 1 − µ1 )⊤ Qb
R(ω, µ
ν
σ
G
i
h ψ2
i
1 h φψ
b )⊤ Qν + E 2 2 ν
b ⊤ Qb
(X 1 − ν
ν ,
+ 2 2E
σ
FG
σ G

for φ = φ(F, S) and ψ = ψ(G, S). Because of (X 1 −b
ν )⊤ Qb
ν ≤ {(X 1 −b
ν )⊤ Q(X 1 −b
ν )}1/2 {b
ν ⊤ Qb
ν }1/2 ,
we have
b )⊤ Q(X 1 − ν
b )}1/2 on {b
2S 2 n {(X 1 − ν
2 φψ
ν ⊤ Qb
ν }1/2 o
⊤
b
(X
−
ν
)
Qb
ν
≤
φ
ψ
P
1
k
σ2 F G
σ2
b ⊤ A−1 ν
b
ν
b k2 −1
j=1 kX j − ν
Vj

≤

=
for B defined in (2.9). Thus,

2S 2 n

φ2

σ2

φ2

o
b)
b ⊤ Qb
(X 1 −
ν
ν
1−ν
+ ψ2
Pk
b )2
2(b
ν ⊤ A−1 ν
b k2 −1 }2
2{ j=1 kX j − ν
V
b )⊤ Q(X
ν

ψ2

j

⊤

b Qb
S
S
ν
ν
,
B+ 2
⊤
2
−1
σ F
σ Gν
b A ν
b

1 h
ψi
b 1 (φ, ψ)) =R(ω, µ
b 1 (φ)) − 2 2 E (X 1 − µ1 )⊤ Qb
R(ω, µ
ν
σ
G
h S φ2 i
h S ψ2 i
W +E 2 B
+ 2E 2
σ G
σ F
b 1 (φ)) − 2J1 + 2J2 + J3 , (say)
=R(ω, µ

b ⊤ Qb
b . The Stein identity is applied to rewrite J1 as
where W = ν
ν /b
ν ⊤ A−1 ν
J1 = E

hψ

G

tr (AQ) − 2W

i
ψ
+ 2W ψG .
G

(2.17)

(2.18)

The chi-square identity is used to rewrite J2 as
h
ψ2
ψψS i
J2 = E (n + 2)W
.
− 4W ψψG + 4W S
G
G

(2.19)

The calculation of J3 is given in (2.10), and from (2.11), (2.17), (2.18) and (2.19), it follows that
b 1 (φ, ψ)) − tr (V 1 Q)
R(ω, µ
hφ
ψ
2(n + 2)Bφ − 2tr {(V 1 − A)Q} + 4B +
2(n + 2)W ψ − 2tr (AQ) + 4W
=E
F
G
ψψS i
φφS
.
(2.20)
− 4W ψG − 8W ψψG + 8W S
− 4BφF − 8BφφF + 8BS
F
G
b 1 (φ, ψ)) ≤ tr (V 1 Q) under the condiBecause W ≤ Chmax (AQ), it can be verified that R(ω, µ
tions (a)-(d) in Theorem 2.2, which is proved.


7

3

Hierarchical and Empirical Bayes Minimax Estimators

3.1

Empirical Bayes estimator under the uniform prior for ν

We begin with assuming the prior distribution
µi | ν, τ 2 ∼Np (ν, τ 2 V i ),

for i = 1, . . . , k,

(3.1)

p

ν ∼Uniform(R ),

where Uniform(Rp ) denotes the improper uniform distribution over Rp , and τ 2 is an unknown
parameter. The posterior distribution of µi given X i and ν, the posterior distribution of ν
given X 1 , . . . , X k and the marginal distribution of X 1 , . . . , X k are


b ∗i (σ 2 , τ 2 , ν), (σ −2 + τ −2 )−1 , i = 1, . . . , k,
µi | X i , ν, τ 2 , σ 2 ∼Np µ
ν | X 1 , . . . , X k , τ 2 , σ 2 ∼Np (b
ν , (τ 2 + σ 2 )A),
fπ (x1 , . . . , xk | τ 2 , σ 2 ) ∝

1
exp
(τ 2 + σ 2 )p(k−1)/2

n

−

(3.2)

Pk

b k2V −1 o
i=1 kxi − ν
i
,
2(τ 2 + σ 2 )

b ∗i (σ 2 , τ 2 , ν) = X i − {σ 2 /(τ 2 + σ 2 )}(X i − ν). Then, the Bayes estimator of µ1 is
where µ

σ2
b ).
(X 1 − ν
(3.3)
τ 2 + σ2
P
b k2V −1 /{p(k−1)−2}
Because τ 2 +σ 2 and σ 2 are unknown, we estimate τ 2 +σ 2 by ki=1 kX i − ν
2 2
bB
µ
µ∗1 (σ 2 , τ 2 , ν) | X 1 , . . . , X k ] = X 1 −
1 (σ , τ ) = E[b

i

from the marginal likelihood in (3.2). When σ 2 is estimated by σ̂ 2 = S/(n + 2), the resulting
empirical Bayes estimator is
a 
0
b ),
b EB
, 1 (X 1 − ν
(3.4)
µ
=
X
−
min
1
1
F

for a0 = {p(k − 1) − 2}/(n + 2). It follows from Theorem 2.1 that the empirical Bayes estimator
b EB
µ
is minimax for 0 < a0 ≤ 2[tr {(V 1 − A)Q}/Chmax {(V 1 − A)Q} − 2]/(n + 2).
1

3.2

Hierarchical Bayes minimax estimator under the uniform prior for ν

We consider the prior distribution for σ 2 in (1.1) and τ 2 in (3.1), namely, in addition of (3.1),
we assume that
 σ 2 a+1
π(τ 2 | σ 2 ) ∝ 2
,
τ + σ2
(3.5)
π(σ 2 ) ∝(σ 2 )c−2 , for σ 2 ≤ 1/L,
where a and c are constants, and L is a positive constant. From (3.2), the posterior distribution
of (τ 2 , σ 2 ) given X 1 , . . . , X k , S is
π(τ 2 , σ 2 | x1 , . . . , xk , S)
n
σ 2 p(k−1)/2+a+1  1 {n+p(k−1)}/2+2−c
exp
−
∝
τ 2 + σ2
σ2


8

Pk

i=1 kxi

b k2V −1
−ν

2(τ 2 + σ 2 )

i

−

S o
.
2σ 2

Then, the hierarchical Bayes estimator of µ1 relative to the quadratic loss (1.2) is written as
b HB
µ
=E[µ1 /σ 2 | X 1 , . . . , X k , S]/E[1/σ 2 | X 1 , . . . , X k , S]
1
=X 1 −

E[(τ 2 + σ 2 )−1 | X 1 , . . . , X k , S]
b ).
(X 1 − ν
E[(σ 2 )−1 | X 1 , . . . , X k , S]

(3.6)

Making the transformation λ = σ 2 /(τ 2 +σ 2 ) and η = 1/σ 2 with the Jacobian |∂(τ 2 , σ 2 )/∂(λ, η)| =
1/(λ2 η 3 ) gives
π(λ, η | x1 , . . . , xk , S) ∝ λp(k−1)/2+a−1 η {n+p(k−1)}/2−c−1 exp

n

−

k
λη X
η o
b k2V −1 − S ,
kxi − ν
i
2
2
i=1

where 0 < λ < 1 and η ≥ L. Thus, we have
R 1 R ∞ p(k−1)/2+a {n+p(k−1)}/2−c
η
exp{− ηS
E[(τ 2 + σ 2 )−1 | X 1 , . . . , X k , S]
2 (λF + 1)}dηdλ
0 L λ
=R 1 R ∞
p(k−1)/2+a−1
{n+p(k−1)}/2−c
E[(σ 2 )−1 | X 1 , . . . , X k , S]
λ
η
exp{− ηS (λF + 1)}dηdλ
0

2

L

=φHB (F, S)/F ,

for
HB

φ

RF R∞

xp(k−1)/2+a v {n+p(k−1)}/2−c exp{−v(x + 1)/2}dvdx
,
(F, S) = R F0 R ∞LS
p(k−1)/2+a−1 v {n+p(k−1)}/2−c exp{−v(x + 1)/2}dvdx
x
0 LS

(3.7)

where the transformations x = F λ and v = Sη are used.

It is noted that the hierarchical Bayes estimator belongs to the class (2.1), and the minimaxity can be shown by checking the conditions (a) and (b) in Theorem 2.1. By differentiating
φHB (F, S) with respect to F and S, the condition (a) can be easily verified. The condition (a)
implies that
R ∞ p(k−1)/2+a
x
/(1 + x){n+p(k−1)}/2+1−c dx
HB
HB
φ (F, S) ≤ lim lim φ (F, S) = R ∞0 p(k−1)/2+a−1
F →∞ S→0
/(1 + x){n+p(k−1)}/2+1−c dx
0 x
B(p(k − 1)/2 + a + 1, n/2 − a − c)
p(k − 1) + 2a
=
,
=
B(p(k − 1)/2 + a, n/2 − a − c + 1)
n − 2(a + c)
for the beta function B(a, b) if a > −p(k − 1)/2 and a + c < n/2. Thus, the condition (b) is
satisfied if
a > −p(k − 1)/2, a + c < n/2,
p(k − 1) + 2a
≤ 2[tr {(V 1 − A)Q}/Chmax {(V 1 − A)Q} − 2]/(n + 2),
n − 2(a + c)

(3.8)

which provides conditions on a and c for the minimaxity of the hierarchical Bayes estimator.

3.3

Hierarchical empirical Bayes minimax estimator under the normal prior
for ν

b EB
b HB
The empirical Bayes estimator µ
and the hierarchical Bayes estimator µ
are derived under
1
1
the uniform prior distribution for ν. Instead of the uniform prior, we here assume the normal
9

prior distribution for ν, namely,
µi | ν, τ 2 ∼Np (ν, τ 2 V i ),

for i = 1, . . . , k,

ν | γ 2 , ∼Np (0, γ 2 A),
−1
P
k
−1
. Then the posterior distributions are given by
for γ > 0 and A =
i=1 V i


b ∗i (σ 2 , τ 2 , ν), (σ −2 + τ −2 )−1 ,
i = 1, . . . , k,
µi | X i , ν, τ 2 , σ 2 ∼Np µ

γ 2 (τ 2 + σ 2 ) 
γ2
b
A ,
ν
,
ν | X 1 , . . . , X k , τ 2 , γ 2 , σ 2 ∼Np 2
γ + τ 2 + σ2 γ 2 + τ 2 + σ2

(3.9)

(3.10)

b ∗i (σ 2 , τ 2 , ν) is defined below (3.2). Thus, the Bayes estimator is
where µ

σ2
(X 1 − E[ν | X 1 , . . . , X k ])
τ 2 + σ2
σ2
σ2
b
b.
=X 1 − 2
(X
−
ν
)
−
ν
1
τ + σ2
γ 2 + τ 2 + σ2

2 2 2
bB
µ
1 (σ , τ , γ ) =X 1 −

(3.11)

Because the marginal density of X 1 , . . . , X k is

1
fπ (x1 , . . . , xk | τ , γ , σ ) ∝ 2
exp
(τ + σ 2 )p(k−1)/2
2

2

2

n

−

Pk

i=1 kxi

b k2V −1 o
−ν
i

2(τ 2 + σ 2 )
kb
ν k2A−1

(3.12)
n
o
1
exp −
,
× 2
2(γ 2 + τ 2 + σ 2 )
(γ + τ 2 + σ 2 )p/2
P
b k2V −1 /{p(k − 1) − 2} and kb
ν k2A−1 /(p − 2),
we can estimate τ 2 + σ 2 , γ 2 + τ 2 + σ 2 by ki=1 kX i − ν
i

respectively, from the marginal likelihood. When σ 2 is estimated by σ̂ 2 = S/(n+2), the resulting
hierarchical empirical Bayes estimator is
a 
b 
0
0
b HEB
b
b,
µ
=
X
−
min
,
1
(X
−
ν
)
−
min
,1 ν
(3.13)
1
1
1
F
G

for a0 = {p(k − 1) − 2}/(n + 2) and b0 = (p − 2)/(n + 2). This estimator belongs to the
b HEB
class (2.3). It follows from Theorem 2.2 that the hierarchical empirical Bayes estimator µ
1
is minimax if 0 < a0 ≤ [tr {(V 1 − A)Q}/Chmax {(V 1 − A)Q} − 2]/(n + 2) and 0 < b0 ≤
{tr (AQ)/Chmax (AQ) − 2}/(n + 2).

4

Extension to Estimation of Linear Combinations

We here extend the result in Section 2 to estimation of the linear combination of µ1 , . . . , µk ,
namely,
k
X
di µi ,
θ=
i=1

where d1 , . . . , dk are constants. Based on (2.1), we consider the class of the estimators
b
θ(φ)
=

k
o
n
X
φ(F, S)
b) .
(X i − ν
di X i −
F
i=1

10

(4.1)

b
When the estimator is evaluated in light of risk relative to the loss function kθ(φ)
− θk2Q /σ 2 , we
obtain conditions for minimaxity of the estimators (4.1).

P
P
P
Theorem 4.1 Assume the conditions Chmax (( ki=1 d2i V i −( ki=1 di )2 A)Q) 6= 0 and tr {( ki=1 d2i V i −
P
P
P
b
( ki=1 di )2 A)Q}/Chmax (( ki=1 d2i V i − ( ki=1 di )2 A)Q) > 2. Then, the estimator θ(φ)
is minimax if φ(F, S) satisfies the following conditions:
(a) φ(F, S) is non-decreasing in F and non-increasing in S.
(b) φ(F, S) satisfies the inequality
P
P
i
2 h tr {( ki=1 d2i V i − ( ki=1 di )2 A)Q}
−
2
.
0 < φ(F, S) ≤
Pk
P
n + 2 Chmax (( i=1 d2i V i − ( ki=1 di )2 A)Q)

(4.2)

b
Proof. The risk function of the estimator θ(φ)
is
b
R(ω, θ(φ))
=

k
X
2 i
d2i h
φ(F, S)
b
(X
−
ν
)
E
X
−
µ
−
i
i
i
Q
σ2
F
i=1
o⊤ n
oi
X di dj hn
φ(F, S)
φ(F, S)
b
b
(X
−
ν
)
(X
−
ν
)
E
X
−
µ
−
Q
X
−
µ
−
+
i
j
i
j
i
j
σ2
F
F
i6=j

=K1 + K2 .

(say)

Concerning K2 , the same arguments as in (2.8) and (2.10) are used to get
n
h h
oii
1 h
φi
⊤ φ(F, S)
⊤
b
b
(X
−
ν
)
=E
tr
V
Q∇
E
(X
−
µ
)
Q(X
−
ν
)
j
i
i
i
j
i
σ2
F
F
h
i
φ
φ
=E − tr (AQ) − 2Bij + 2Bij φF ,
F
F
and
Bij E
where

h
h S φ2 i
φ2
φφS i
=B
E
(n
+
2)
,
−
4φφ
+
4S
ij
F
σ2 F
F
F

k
X
b )⊤ V −1
b ) Q(X j − ν
b )/
(X a − ν
Bij = (X i − ν
a (X a − ν).
⊤

(4.3)

a=1

Thus, K2 is written as
K2 =

X
i6=j

h
n φ
φ
φ2
φφS oi
.
di dj 2tr (AQ) + Bij 4 − 4φF + (n + 2) − 4φφF + 4S
F
F
F
F

(4.4)

Concerning K1 , from (2.11), it follows that
K1 =

h
φ
d2i E tr (V i Q) − 2tr {(V i − A)Q}
F
i=1
n
φ2
φ
φφS oi
+ Bii (n + 2) + 4 − 4φF − 4φφF + 4S
.
F
F
F

k
X

11

(4.5)

It is here observed that

k
X

d2i Bii +

i=1

X

di dj Bij = B(X),

i6=j

where X = (X 1 , . . . , X k ) and
P
P
b )}⊤ Q{ ki=1 di (xi − ν
b )}
{ ki=1 di (xi − ν
.
B(x1 , . . . , xk ) =
Pk
−1
⊤
b ) V j (xj − ν)
j=1 (xj − ν

(4.6)

Combining (4.4) and (4.5), one gets
b
R(ω, θ(φ))
=

k
X

d2i tr (V i Q)

i=1

+E

hφn
F

−2

k
X

tr (d2i V i Q) +

i=1

k
2
X
o
2
di tr (AQ) + 4B(X) + (n + 2)φB(X)
i=1

n
φφS oi
+ B(X) − 4φF − 4φφF + 4S
,
F
P
which is smaller than R(ω, ki=1 di X i ) under the conditions (a) and (b) in Theorem 4.1, because
P
P
B(x) ≤ Chmax (( ki=1 d2i V i − ( ki=1 di )2 A)Q) from Lemma 4.1. Hence, the proof of Theorem
4.1 is complete.

Lemma 4.1 For B(x1 , . . . , xk ) given in (2.9), it holds that
B(x1 , . . . , xk ) ≤ Chmax

k
 X

d2i V i

k
 
X
di )2 A Q .
−(
i=1

i=1

b and x = (x1 , . . . , xk ). Then, it is noted that
Proof. Let y i = xi − ν
Pk
−1
b ) = 0, which means that
i=1 V i (xi − ν
   
y1
y1
 ..   .. 
W  .  =  . ,
yk

for

yk

 
A
 ..  −1
W = block diag(I, . . . , I) −  .  (V 1 , . . . , V −1
k ).
A

12

Pk

−1
i=1 V i y i

=

Then, the numerator and the denominator of B(x) are


 
d1 I
y1
k
k
X
X
 .. 
⊤
⊤  .. 
⊤
b )} = (y 1 , . . . , y k )  .  Q(d1 I, . . . , dk I)  . 
b )} Q{
di (xi − ν
di (xi − ν
{
i=1
i=1
dk I
yk


 
d1 I
y1
 .. 
⊤  .. 
⊤
⊤
= (y 1 , . . . , y k )W  .  Q(d1 I, . . . , dk I)W  .  ,
dk I

yk


y1
−1
−1  .. 
⊤
⊤
b )⊤ V −1
(xj − ν
j (xj − ν) = (y 1 , . . . , y k )block diag(V 1 , . . . , V k )  .  .
j=1
yk

k
X



Thus, we get an upper bound given by


d1 I




B(x) ≤Chmax block diag(V 1 , . . . , V k )W ⊤  ...  Q(d1 I, . . . , dk I)W
dk I


d1 I 



=Chmax Q(d1 I, . . . , dk I)W block diag(V 1 , . . . , V k )W ⊤  ... 
dk I
k
k
 
 X
X
di )2 A Q ,
d2i V i − (
=Chmax
i=1

i=1

which shows Lemma 4.1.



We here give some examples of the
P condition (4.2)
P in specific cases. For example, in the
case of k = 2, it is observed that tr {( 2i=1 d2i V i − ( 2i=1 di )2 A)Q} = tr {(d1 V 1 − d2 V 2 )(V 1 +
V 2 )−1 (d1 V 1 − d2 V 2 )Q} = tr (H), where
H = (V 1 + V 2 )−1/2 (d1 V 1 − d2 V 2 )Q(d1 V 1 − d2 V 2 )(V 1 + V 2 )−1/2 .
P
P
Also, it can be seen that Chmax (( ki=1 d2i V i − ( ki=1 di )2 A)Q) = Chmax (H). Hence, the condition (b) in Theorem 4.1 is expressed as
i
2 h tr (H)
−2 .
0 < φ(F, S) ≤
n + 2 Chmax (H)
For example, in the case that V 1 = · · · = V k = Q = I, we have
k
k
k
X
X
X
2
2
(di − d)2 ,
di ) A)Q} = p
di V i − (
tr {(
i=1

i=1

i=1

Pk
−1

Pk
Pk
Pk
2
2
2
for d = k
i=1 di . Similarly, Chmax (( i=1 di V i − ( i=1 di ) A)Q) =
i=1 (di − d) , which
implies that the condition (b) is expressed as
0 < φ(F, S) ≤ 2(p − 2)/(n + 2).
13

5

Simulation Studies

We investigate the numerical performances of the risk functions of the preliminary-test estimator
and several empirical and hierarchical Bayes estimators through simulation. We employ the
quadratic loss function L(δ 1 , µ1 , σ 2 ) in (1.2) for Q = V −1
1 .
The estimators which we compare are the following five:
PT: the preliminary-test estimator given in (1.5)

X 1 if F > (p(k − 1)/n)Fp(k−1),n,α ,
PT
b
µ1 =
b
ν
otherwise,

JS: the James-Stein estimator

b JS
µ
1 = X1 −

S
p−2
n + 2 kX 1 k2

X 1,

V −1
1

EB: the empirical Bayes estimator given in (3.4)
a 
0
b ),
b EB
, 1 (X 1 − ν
µ
=
X
−
min
1
1
F
−1
for a0 = [tr {(V 1 − A)V −1
1 }/Chmax {(V 1 − A)V 1 } − 2]/(n + 2) (one can see that this constant
choice is optimal with respect to the upper bound of the risk difference),
HB: the hierarchical Bayes estimator given in (3.6) and (3.7),
φHB (F, S)
b ),
(X 1 − ν
F
HEB: the hierarchical empirical Bayes estimator given in (3.13)
b 
a 
0
0
b
b,
b HEB
,
1
(X
−
ν
)
−
min
,1 ν
µ
=
X
−
min
1
1
1
F
G
−1
−1
−1
for a0 = [tr {(V 1 −A)V −1
1 }/Chmax {(V 1 −A)V 1 }−2]/{2(n+2)} and b0 = {tr (AV 1 )/Chmax (AV 1 )−
2}/{2(n + 2)} (these constants are also optimal choices).
b HB
µ
= X1 −
1

It is noted that the James-Stein estimator does not use X 2 , . . . , X k , but is minimax. Conb HB
cerning the hierarchical Bayes estimator µ
1 , the constants c and L are c = 1 and L = 0, and
a is the solution of the equation
p(k − 1) + 2a
−1
(n + 2) = tr {(V 1 − A)V −1
1 }/Chmax {(V 1 − A)V 1 } − 2,
n − 2(a + 1)

which guarantees the minimaxity from the condition (3.8).
In this simulation, we generate random numbers of X 1 , . . . , X k and S based on the model
(1.1) for p = k = 5, n = 20, σ 2 = 2 and V i = (0.1 × i)I p , i = 1, . . . , k. For the mean vectors µi ,
we treat the 12 cases:
(µ1 , . . . , µ5 )
=(0, 0, 0, 0, 0), (1j 5 , 1j 5 , 1j 5 , 1j 5 , 1j 5 ), (2j 5 , 2j 5 , 2j 5 , 2j 5 , 2j 5 ), (3j 5 , 3j 5 , 3j 5 , 3j 5 , 3j 5 ),
(−0.4j 5 , −0.2j 5 , 0, 0.2j 5 , 0.4j 5 ), (2j 5 , −0.5j 5 , −0.5j 5 , −0.5j 5 , −0.5j 5 ),
(4j 5 , −1j 5 , −1j 5 , −1j 5 , −1j 5 ),
(1.2j 5 , 1.4j 5 , 1.6j 5 , 1.8j 5 , 2j 5 ), (0, 2j 5 , 2j 5 , 2j 5 , 2j 5 ), (0, 4j 5 , 4j 5 , 4j 5 , 4j 5 ), (2j 5 , 0, 0, 0, 0, 0),

14

where j p = (1, . . . , 1)⊤ ∈ Rp . The first four are the cases of equal means, the next three are the
P
cases with 5i=1 µi = 0 and the last four are various unbalanced cases.

b 1 , based on 5,000 replication of simulation, we obtain an approximated
For each estimator µ
b 1 ) = E[L(b
value of the risk function R(ω, µ
µ1 , µ1 , σ 2 )]. Table 1 reports the percentage relative
b 1 over X 1 , defined by
improvement in average loss (PRIAL) of each estimator µ
b 1 )}/R(ω, X 1 ).
PRIAL = 100 × {R(ω, X 1 ) − R(ω, µ

Table 1: Values of PRIAL of estimators PT, JS, EB, HB and HEB
(µ1 , . . . , µ5 )
(0, 0, 0, 0, 0) ⊗ j 5
(1, 1, 1, 1, 1) ⊗ j 5
(2, 2, 2, 2, 2) ⊗ j 5
(3, 3, 3, 3, 3) ⊗ j 5
(−0.4, −0.2, 0, 0.2, 0.4) ⊗ j 5
(2, −0.5, −0.5, −0.5, −0.5) ⊗ j 5
(4, −1, −1, −1, −1) ⊗ j 5
(1.2, 1.4, 1.6, 1.8, 2) ⊗ j 5
(0.2, 2, 2, 2, 2) ⊗ j 5
(0.4, 4, 4, 4, 4) ⊗ j 5
(2, 0, 0, 0, 0) ⊗ j 5

PT
52.15317
52.15317
52.15317
52.15317
37.34717
−56.8291
0.7375904
37.34717
−98.94453
−2.492994
−94.45962

JS
53.97469
12.89115
4.066823
2.268442
37.20396
4.066823
1.620614
9.463467
49.73141
37.56052
4.066823

EB
14.66425
14.66425
14.66425
14.66425
13.01833
3.213333
1.358956
13.01833
4.947591
1.805795
4.439434

HB
14.57437
14.57437
14.57437
14.57437
12.97352
3.213459
1.358821
12.97352
4.949466
1.80584
4.440511

HEB
26.38606
9.891098
8.516356
8.249028
22.64692
6.031053
2.098222
8.397694
5.183324
2.071347
4.479298

It is revealed from Table 1 that the performance of the preliminary test estimator PT strongly
depends on the setup of parameters, namely it is good under the hypothesis of equal means, but
not good for parameters close to the hypothesis. The James-Stein estimator JS is good for small
µ1 , but not good for large µ1 . The empirical Bayes estimator EB and the hierarchical Bayes
estimator HB perform similarly and they remain good even for large µ1 as long as µ1 = · · · µ5 .
The performance
of HEB depends on parameters and is good for smaller means. For means
P
with 5i=1 µi = 0, HEB is better than EB and HB. Thus, EB, HB and HEB are used as an
alternative to PT.

6

Concluding Remarks

An interesting query is to find an admissible and minimax estimator of µ1 . In the framework of
simultaneous estimation of (µ1 , . . . , µk ), Imai, et al. (2017) demonstrated that all the estimators
within the class (2.1) are improved on by an estimator belonging to the class (2.3), which means
that the hierarchical Bayes estimators against the uniform prior of ν is inadmissible. However,
we could not show the same story in the single estimation of µ1 . Because the hierarchical Bayes
b HB
estimator µ
is derived under the uniform prior for ν, it could be supposedly inadmissible.
1
Thus, whether it is admissible or not is an open question.

15

An approach to the admissible and minimax estimation is the proper prior distribution
σ 2 a+1
,
τ 2 + σ2

b+1
σ2
π(γ 2 | τ 2 , σ 2 ) ∝ 2
,
γ + τ 2 + σ2
π(σ 2 ) ∝(σ 2 )c−3 , for σ 2 ≤ 1/L,

π(τ 2 | σ 2 ) ∝

(6.1)

where a, b and c are constants and L is a positive constant. As seen from Imai, et al. (2017),
the resulting Bayes estimator has the form
b F1 B = X 1 −
µ

ψ F B (F, G, S)
φF B (F, G, S)
b) −
b,
(X 1 − ν
ν
F
G

(6.2)

where φF B (F, G, S) and ψ F B (F, G, S) are functions of F , G and S. Unfortunately, we could not
establish the minimaxity for this type of estimators, which is another interesting query.
In this paper, we investigated the minimaxity of Bayesian alternatives to the preliminary test
estimator. Beyond this framework, we consider the estimation of µ1 based on X 1 and S without
X 2 , . . . , X k . Using the argument as in Strawderman (1973), we can derive an admissible and
minimax estimator of µ1 . Taking the hypothesis H0 : µ1 = · · · = µk into account, we can use
the same argument as in Strawderman (1973) under the assumption that X 2 , . . . , X k are given
and fixed.
we consider the prior distributions (3.1) and (3.5) where ν is replaced with
Pk Namely,
−1 −1 Pk
−1
∗
ν = ( i=2 V i )
i=2 V i X i . Then, the Bayes estimator is
X1 −

φS (kX 1 − ν ∗ k2 , S)
(X 1 − ν ∗ ),
kX 1 − ν ∗ k2 /S

where φS is a function derived from the Bayes estimator against the Strawderman type prior.
This shrinks X 1 towards ν ∗ and is minimax under some condition on φS . Thus, it is admissible
and minimax in the framework that X 2 , . . . , X k are given and fixed.

Acknowledgments
We would like to thank the Editor, the Associate Editor and the reviewer for valuable comments
and helpful suggestions which led to an improved version of this paper. Research of the second
author was supported in part by Grant-in-Aid for Scientific Research (15H01943 and 26330036)
from Japan Society for the Promotion of Science.

References
[1] J. Berger, Statistical Decision Theory and Bayesian Analysis, 2nd. ed. (1985) Springer.
[2] M. Bilodeau, T. Kariya, Minimax estimators in the normal MANOVA model, J. Multivariate Anal. 28 (1989) 260–270.
[3] L.D. Brown, E.I. George, X. Xu, Admissible predictive density estimation, Ann. Statist. 36
(2008) 1156–1170.
16

[4] B. Efron, C.N. Morris, Steinfs estimation rule and its competitors : an empirical Bayes
approach, J. Amer. Statist. Assoc. 68 (1973) 117–130.
[5] B. Efron, C. Morris, Families of minimax estimators of the mean of a multivariate normal
distribution, Ann. Statist. 4 (1976) 11–21.
[6] M. Ghosh, B.K. Sinha, Empirical and hierarchical Bayes competitors of preliminary test
estimators in two sample problems, J. Multivariate Anal. 27 (1988) 206–227.
[7] R. Imai, T. Kubokawa, M. Ghosh, Bayesian simultaneous estimation for means in k sample
problems, (2017) arXiv:1711.10822.
[8] W. James, C. Stein, Estimation with quadratic loss, in: Proc. Fourth Berkeley Sympos.
Math. Statist. and Prob., Vol.I, Univ. California Press, Berkeley, 1961, pp.361–379.
[9] F. Komaki, A shrinkage predictive distribution for multivariate normal observables.
Biometrika 88 (2001) 859–864.
[10] S.L. Sclove, C. Morris, R. Radhakrishnan, Nonoptimality of preliminary test estimators for
the multinormal mean, Ann. Math. Statist. 43 (1972) 1481–1490.
[11] A.F.M.N. Smith, Bayes estimates in one-way and two-way models, Biometrika 60 (1973)
319–329.
[12] C. Stein, Inadmissibility of the usual estimator for the mean of a multivariate normal
distribution, Proc. Third Berkeley Symp. Math. Statist. Probab. 1 (1956) 197–206, Univ.
California Press, Berkeley.
[13] C. Stein, Estimation of the mean of a multivariate normal distribution, Ann. Statist. 9
(1981) 1135–1151.
[14] W.E. Strawderman, Proper Bayes minimax estimators of the multivariate normal mean,
Ann. Math. Statist. 42 (1971) 385–388.
[15] W.E. Strawderman, Proper Bayes minimax estimators of the multivariate normal mean
vector for the case of common unknown variances, Ann. Statist. 1 (1973) 1189–1194.
[16] H. Tsukuma, T. Kubokawa, A unified approach to estimating a normal mean matrix in
high and low dimensions, J. Multivariate Anal. 139 (2015) 312–328.

17

