Harmonic Bayesian prediction under
α-divergence

arXiv:1605.05899v4 [math.ST] 28 Jul 2017

Yuzo Maruyama∗ and Toshio Ohnishi†
University of Tokyo∗ and Kyushu University†
e-mail: maruyama@csis.u-tokyo.ac.jp; ohnishi@econ.kyushu-u.ac.jp
Abstract: We investigate Bayesian shrinkage methods for constructing
predictive distributions. We consider the multivariate normal model with
a known covariance matrix and show that the Bayesian predictive density with respect to Stein’s harmonic prior dominates the best invariant
Bayesian predictive density, when the dimension is greater than three.
Alpha-divergence from the true distribution to a predictive distribution
is adopted as a loss function.
AMS 2000 subject classifications: Primary 62C20; secondary 62J07.
Keywords and phrases: harmonic prior, minimaxity, Bayesian predictive
density.

1. Introduction
Let X ∼ Nd (µ, vx I) and Y ∼ Nd (µ, vy I) be independent d-dimensional multivariate normal vectors with common unknown mean µ. We assume that d ≥
3 and that vx and vy are known. Let φ(·, σ 2 ) be the probability density of
Nd (0, σ 2 I). Then the probability density of X and that of Y are φ(x − µ, vx )
and φ(y − µ, vy ), respectively.
Based on only observing X = x, we consider the problem of obtaining a
predictive density p̂(y | x) for Y that is close to the true density φ(y − µ, vy ). In
most earlier papers on such prediction problems, a predictive density p̂(y | x) is
often evaluated by
Z
φ(y − µ, vy )
dy,
(1.1)
DKL {φ(y − µ, vy ) || p̂(y | x)} =
φ(y − µ, vy ) log
p̂(y | x)
Rd
which is called the Kullback-Leibler divergence loss (KL-div loss) from φ(y −
µ, vy ) to p̂(y | x). The overall quality of the procedure p̂(y | x) for each µ is then
summarized by the Kullback-Leibler divergence risk
Z
RKL {φ(y − µ, vy ) || p̂(y | x)} =
DKL {φ(y − µ, vy ) || p̂(y | x)} φ(x − µ, vx )dx.
Rd

(1.2)
Aitchison (1975) showed that the Bayesian solution with respect to a prior
π(µ) under KL-div loss given by (1.1) is the Bayesian predictive density
Z
p̂π (y | x) =
φ(y − µ, vy )π(µ | x)dµ,
(1.3)
Rd

1

Y. Maruyama and T. Ohnishi/A Bayesian prediction under α-divergence

2

where π(µ | x) = φ(x − µ, vx )π(µ)/mπ (x, vx ) is the posterior density corresponding to π(µ) and
Z
φ(x − µ, v)π(µ)dµ

mπ (x, v) =

(1.4)

Rd

is the marginal density of X ∼ Nd (µ, vI) under the prior π(µ).
For the prediction problems in general, many studies suggest the use of the
Bayesian predictive density rather than plug-in densities of the form
φ(y − µ̂(x), vy ),
where µ̂(x) is an estimated value of µ. Liang and Barron (2004) showed that
the Bayesian predictive density with respect to the uniform prior
πU (µ) = 1,

(1.5)

which is given by
Z
p̂U (y | x) =

φ(y − µ, vy )πU (µ | x)dµ = φ(y − x, vx + vy )

(1.6)

Rd

is best invariant and minimax. Although the best invariant Bayesian predictive
density is generally a good default procedure, it has been shown to be inadmissible in some cases. Specifically, Komaki (2001) showed that the Bayesian
predictive density with respect to Stein’s (1974) harmonic prior
πH (µ) = kµk−(d−2)

(1.7)

dominates the best invariant Bayesian predictive density p̂U (y | x). George, Liang
and Xu (2006) extended Komaki’s (2001) result to general shrinkage priors
including Strawderman’s (1971) prior.
From a more general viewpoint, the KL-div loss given by (1.1) is in the class
of α-divergence loss (α-div loss) introduced by Csiszár (1967) and defined by


Z
p̂(y | x)
φ(y − µ, vy )dy, (1.8)
Dα {φ(y − µ, vy ) || p̂(y | x)} =
fα
φ(y − µ, vy )
Rd
where



2
(1+α)/2

,
 4/(1 − α ) 1 − z
fα (z) = z log z,


− log z,

|α| < 1,
α = 1,
α = −1.

When α = −1, we have
D−1 {φ(y − µ, vy ) || p̂(y | x)} = DKL {φ(y − µ, vy ) || p̂(y | x)} ,
where DKL is given by (1.1). When α = 0, we have
f0 (z) = 4(1 − z 1/2 ),

Y. Maruyama and T. Ohnishi/A Bayesian prediction under α-divergence

Z
D0 {φ(y − µ, vy ) || p̂(y | x)} = 2

3

n
o2
p̂1/2 (y | x) − φ1/2 (y − µ, vy ) dy,

Rd

p
where D0 {φ(y − µ, vy ) || p̂(y | x)}/2 is the Hellinger distance between p̂(y | x)
and φ(y − µ, vy ). As in the Kullback-Leibler divergence risk given by (1.2),
the overall quality of the procedure p̂(y | x) for each µ is summarized by the
α-divergence risk
Z
Rα {φ(y − µ, vy ) || p̂(y | x)} =
Dα {φ(y − µ, vy ) || p̂(y | x)} φ(x − µ, vx )dx.
Rd

Corcuera and Giummolè (1999) showed that a Bayesian predictive density under
α-div loss is
Z
2
 1−α
1−α



,
−1 ≤ α < 1,
φ 2 (y − µ, vy )φ(x − µ, vx )π(µ)dµ
d
RZ

p̂π (y | x; α) ∝


exp
{log φ(y − µ, vy )} φ(x − µ, vx )π(µ)dµ , α = 1.
Rd

(1.9)
By (1.9), in the prediction problem under α-div loss with α = 1 from the
Bayesian point of view, the Bayesian solution is the normal density
p̂π (y | x; 1) = φ(y − µ̂π (x), vy ),
where µ̂π (x) is the posterior mean given by
Z
µ̂π (x) =
µπ(µ | x)dµ = x + vx ∇x log m(x, vx ).

(1.10)

(1.11)

Rd

In general, the Bayesian prediction problem under α = 1 reduces to the estimation problem under the KL-div loss in the case of the exponential family density.
This is because the exponential family density is closed under the calculation in
(1.9) with α = 1, as pointed out in Yanagimoto and Ohnishi (2009).
As demonstrated in Maruyama and Strawderman (2012), the α-div loss in
the case of α = 1 is written as
D1 {φ(y − µ, vy ) || φ(y − µ̂π (x), vy )} =

kµ̂π (x) − µk2
,
2vy

and hence the prediction problem under α = 1 reduces to the estimation problem
of µ under the quadratic loss. Stein (1981) showed that
"
#
1/2


∆x mπ (X, vx )
2
2
,
(1.12)
EX kµ̂π (X) − µk = dvx + 4vx EX
1/2
mπ (X, vx )
which implies that the risk difference under α = 1 is expressed as
R1 {φ(y − µ, vy ) || p̂U (y | x; 1)} − R1 {φ(y − µ, vy ) || p̂π (y | x; 1)}
"
#
1/2
2vx2
∆x mπ (X, vx )
=
EX −
.
1/2
vy
mπ (X, vx )

(1.13)

Y. Maruyama and T. Ohnishi/A Bayesian prediction under α-divergence

4

Under the KL-div loss or α-div loss with α = −1, George, Liang and Xu
(2006) showed that the risk difference is given by
R−1 {φ(y − µ, vy ) || p̂U (y | x; −1)} − R−1 {φ(y − µ, vy ) || p̂π (y | x; −1)}
"
#
Z vx
1/2
∆z mπ (Z, v)
EZ −
=2
dv,
1/2
v∗
mπ (Z, v)

(1.14)

where p̂U (y | x; −1) is given by (1.6), Z ∼ Nd (µ, vI) and v∗ = vx vy /(vx + vy ).
From this viewpoint, George, Liang and Xu (2006) and Brown, George and Xu
(2008) considered the prediction problem under α-div loss with two extreme
cases α = ±1 and found a beautiful relationship of risk differences for two cases
via ∆z {mπ (z, v)}1/2 for some v. Under both risks R1 and R−1 , any shrinkage
prior of the satisfier of the superharmonicity
(
∀v ∈ (v∗ , vx ) for α = −1,
1/2
∆z mπ (z, v) ≤ 0 for
(1.15)
v = vx for α = 1,
implies the improvement over the best invariant Bayesian procedure. It is wellknown that the superharmonicity of π(µ), ∆µ π(µ) ≤ 0, implies the superharmonicity of mπ (z, v), ∆z mπ (z, v) ≤ 0. Further the superharmonicity of
mπ (z, v) implies the superharmonicity of {mπ (z, v)}1/2 . Hence the harmonic
prior πH (µ) = kµk−(d−2) is one of the satisfiers of the superharmonicity of
{mπ (z, v)}1/2 .
Because of the relationship given by (1.13), (1.14) and (1.15), it is of great
interest to find the corresponding link via ∆z {mπ (z, v)}1/2 for α-div loss with
general α ∈ (−1, 1). To our knowledge, decision-theoretic properties seem to
depend on the general structure of the problem (the general type of problem
(location, scale), and the dimension of the parameter space) and on the prior in
a Bayesian-setup, but not on the loss function, as Brown (1979) pointed out in
the estimation problem.
In this paper, we investigate the risk difference, diffRα,U,π , in the case of
α-div loss, defined by
diffRα,U,π = Rα {φ(y − µ, vy ) || p̂U (y | x; α)}
− Rα {φ(y − µ, vy ) || p̂π (y | x; α)} .

(1.16)

In (1.16), p̂π (y | x; α) is given by (1.9) and p̂U (y | x; α) is the Bayesian predictive
density under the uniform prior (1.5), the form of which will be derived in (2.5)
of Section 2. As a generalization of Liang and Barron’s (2004) result, p̂U (y | x; α)
for general α ∈ (−1, 1) is best invariant and minimax, as shown in Appendix A.
Further, analyzing diffRα,U,π , we provide some asymptotic results and a nonasymptotic decision-theoretic result.
Asymptotic results We show not only somewhat expected relationship
lim diffRα,U,π = diffR1,U,π ,

α→1−0

lim

α→−1+0

diffRα,U,π = diffR−1,U,π ,

(1.17)

Y. Maruyama and T. Ohnishi/A Bayesian prediction under α-divergence

5

where diffR1,U,π and diffR−1,U,π are given in (1.13) and (1.14) respectively,
but also the asymptotic relationship for general α ∈ (−1, 1),
lim

vx /vy →+0

diffRα,U,π = diffR1,U,π .

(1.18)

Hence, the asymptotic situation vx /vy → 0 corresponds to the case α → 1 and
∆z {mπ (z, v)}1/2 plays an important role for general α ∈ (−1, 1).
Non-asymptotic result We particularly investigate a decision-theoretic property of the Bayesian predictive density with respect to πH (µ) = kµk−(d−2)
under α-div loss with general α ∈ (−1, 1). We show that, the Bayesian predictive density with respect to πH (µ) = kµk−(d−2) dominates the best invariant
Bayesian predictive density with respect to πU (µ) = 1 if

d+2
2


if
is a positive integer,

vx
d(1
+
α)
1
−
α
2
≤ 
2
d + 2 1 − {κ − 2/(1 − α)}

vy


otherwise,
1−α
d
2κ(κ − 1)
where κ is the smallest integer larger than 2/(1 − α).
The organization of this paper is as follows. In Section 2, we derive the exact
form of p̂π (y | x; α), propose a general sufficient condition for diffRα,U,π ≥ 0,
where diffRα,U,π is given by (1.16), and demonstrate the asymptotic relationship
described in (1.17) and (1.18). In Section 3, we propose the non-asymptotic
result under the harmonic prior πH (µ) = kµk−(d−2) described in the above.
Some technical proofs are given in Sections A and B of Appendix.
2. Bayesian predictive density under α-divergence loss
As in (1.9), the Bayes predictive density under α-div loss is
Z
p̂π (y | x; α) ∝

1/β
φ (x − µ, vx ) φβ (y − µ, vy )π(µ)dµ
,

(2.1)

Rd

where

1−α
.
2
Clearly, it follows from α ∈ (−1, 1) that 0 < β < 1. Let
β=

γ=

1
.
1 + βvx /vy

(2.2)

(2.3)

Since the relation of completing squares with respect to µ, for φ (x − µ, vx ) φβ (y−
µ, vy ), is given by
1
β
kx − µk2 + ky − µk2
vx
vy

Y. Maruyama and T. Ohnishi/A Bayesian prediction under α-divergence

=
=
=
=

6



1
1−γ
2
2
kx − µk +
ky − µk
vx
γ


kγx + (1 − γ)yk2
1−γ
1 1
2
2
2
kµ − {γx + (1 − γ)y}k −
+ kxk +
kyk
vx γ
γ
γ


1 1
kµ − {γx + (1 − γ)y}k2 + (1 − γ)ky − xk2
vx γ
γ
1
kµ − {γx + (1 − γ)y}k2 + β ky − xk2 ,
vx γ
vy

we have the identity,
φ (x − µ, vx ) φβ (y − µ, vy )
= γ (1−β)d/2 φ(γx + (1 − γ)y − µ, vx γ)φβ (y − x, vy /γ).

(2.4)

Under the uniform prior πU (µ) = 1, we have, from (2.4),
Z
φ (x − µ, vx ) φβ (y − µ, vy )πU (µ)dµ = γ (1−β)d/2 φβ (y − x, vy /γ)
Rd

in (2.1). Therefore the Bayesian predictive density under the uniform prior is
p̂U (y | x; α) = φ(y − x, vy /γ) = φ(y − x, vy + βvx ),

(2.5)

which is the target predictive density so that the risk difference
diffRα,U,π = Rα {φ(y − µ, vy ) || p̂U (y | x; α)} − Rα {φ(y − µ, vy ) || p̂π (y | x; α)}
is going to be investigated in this paper. As shown in Appendix A, p̂U (y | x; α)
for general α ∈ (−1, 1) is best invariant and minimax, which is regarded as a
generalization of Liang and Barron’s (2004) minimaxity result. Hence p̂π (y | x; α)
with diffRα,U,π ≥ 0 for all µ ∈ Rd is minimax.
The exact form of Bayes predictive density p̂π (y | x; α) for (2.1) with normalizing constant, which is regarded as a generalization of Theorem 1 of Komaki
(2001) as well as Lemma 2 of George, Liang and Xu (2006), is provided as
follows.
Theorem 2.1. The Bayes predictive density under π(µ) is
1/β

p̂π (y | x; α) =

mπ (γx + (1 − γ)y, vx γ)
h
i p̂U (y | x; α),
1/β
EZ1 mπ (x + ξZ1 , vx γ)

(2.6)

where Z1 ∼ Nd (0, I) and
ξ = (1 − γ)(vy /γ)1/2 .

(2.7)

Proof. By (2.1), (2.4) and (2.5), we have
p̂π (y | x; α) ∝ φ(y − x, vy /γ)m1/β
π (γx + (1 − γ)y, vx γ).

(2.8)

Y. Maruyama and T. Ohnishi/A Bayesian prediction under α-divergence

7

The normalizing constant of (2.8) is
Z
φ(y − x, vy /γ)m1/β
π (γx + (1 − γ)y, vx γ)dy
Rd
Z


x + (1 − γ)(vy /γ)1/2 z1 , vx γ dz1
=
φ(z1 , 1)m1/β
π
Rd
h
i
= EZ1 m1/β
π (x + ξZ1 , vx γ) ,
where the first equality is from the transformation, z1 = (γ/vy )1/2 (y − x).
In the following, as a generalization of the Bayes predictive density, we consider
f (γx + (1 − γ)y)
p̂f (y | x; α) =
p̂U (y | x; α)
(2.9)
EZ1 [f (x + ξZ1 )]
R
where f : Rd → R+ is general. As in the proof of Theorem 2.1, p̂f (y | x; α)dy =
1 follows. Also p̂f (y | x; α) is nonnegative for any y ∈ Rd and hence p̂f (y | x; α)
is regarded as a predictive density.
By the definition of the α-div loss given by (1.8), the risk difference between
p̂U and p̂f is written as
diffRα,U,f
= Rα {φ(y − µ, vy ) || p̂U (y | x; α)} − Rα {φ(y − µ, vy ) || p̂f (y | x; α)}
1−β 
1−β )
Z (
p̂f (y | x; α)
p̂U (y | x; α)
1
−
=
β(1 − β) R2d
φ(y − µ, vy )
φ(y − µ, vy )

(2.10)

× φ(x − µ, vx )φ(y − µ, vy )dxdy.
Then we have a following result.
Theorem 2.2.
1. The risk difference diffRα,U,f given by (2.10) is written
by E[ρ(W, Z)] where W ∼ Nd (µ, vx γ), Z ∼ Nd (0, I) and
Z
4γ (1−β)d/2 ξ −∆w %(w + tz; t; f )
ρ(w, z) = 2 β−1
t
dt
(2.11)
β f
(w) 0 %2/β−1 (w + tz; t; f )
where
β/2

%(u; t; f ) = {EZ1 [f (tZ1 + u)]}

, for Z1 ∼ Nd (0, I).

(2.12)

2. A sufficient condition for diffRα,U,f ≥ 0 for ∀µ ∈ Rd is
∆u %(u; t; f ) ≤ 0

∀u ∈ Rd , 0 ≤ ∀t ≤ ξ.

(2.13)

Proof. Part 2 easily follows from Part 1 and, in the following, we show Part 1.
By (2.4), (2.5), and (2.9), the integrand of (2.10) is rewritten as
(
β−1 
β−1 )
φ(y − µ, vy )
φ(y − µ, vy )
−
φ(y − µ, vy )φ(x − µ, vx )
p̂f (y | x; α)
p̂U (y | x; α)

Y. Maruyama and T. Ohnishi/A Bayesian prediction under α-divergence

=γ

(1−β)d/2

(

EZ1 [f (x + ξZ1 )]
f (γx + (1 − γ)y)

β−1

8

)
−1

× φ(γx + (1 − γ)y − µ, vx γ)φ(y − x, vy /γ).
By the change of variables, w = γx + (1 − γ)y and z = −(γ/vy )1/2 (y − x), where
Jacobian of the matrix below is (γ/vy )d/2 ,
  
 
γId
(1 − γ)Id
w
x
=
,
(2.14)
z
y
(γ/vy )1/2 Id −(γ/vy )1/2 Id
the risk difference is expressed as
#
"
β−1

f (W + ξ(Z1 + Z))
γ (1−β)d/2
−1
EW,Z
EZ1
β(1 − β)
f (W )


γ (1−β)d/2
EW f (W )1−β {g(ξ; W ) − g(0; W )}
β(1 − β)
#
"
Z ξ
γ (1−β)d/2
∂
=
EW f (W )1−β
g(t; W )dt ,
β(1 − β)
0 ∂t

=

(2.15)

where ξ = (1 − γ)(vy /γ)1/2 as in (2.7), W ∼ Nd (µ, vx γI), Z1 ∼ Nd (0, I),
Z ∼ Nd (0, I) and
h
i
β−1
g(t; w) = EZ EZ1 [f (w + t{Z1 + Z})]
.
(2.16)
In the following, EZ1 [f ] = EZ1 [f (w + t{Z1 + z})] for notational simplicity.
Then we have


∂
∂
β−1
g(t; w) = EZ
{EZ1 [f ]}
∂t
∂t
h
i
β−2
(2.17)
= (β − 1)EZ {EZ1 [f ]}
EZ1 [(Z1 + Z)T ∇w f ]
h
i
β−2
= (β − 1)EZ {EZ1 [f ]}
(EZ1 [Z1T ∇w f ] + Z T EZ1 [∇w f ]) .
In (2.17), we have


1
1
EZ1 [Z1T ∇w f ] = EZ1 Z1T ∇z1 f = EZ1 [∆z1 f ]
t
t

(2.18)

= tEZ1 [∆w f ] = t∆w EZ1 [f ]
where the second equality follows from the Gauss divergence theorem. Similarly

Y. Maruyama and T. Ohnishi/A Bayesian prediction under α-divergence

9

we have
h
i
β−2 T
(β − 1)EZ {EZ1 [f ]}
Z EZ1 [∇w f ]


β−2 T 1
= (β − 1)EZ {EZ1 [f ]}
Z EZ1 [∇z f ]
t
h
i
1
β−2 T
= (β − 1)EZ {EZ1 [f ]}
Z ∇z EZ1 [f ]
t
h
i
1
β−1
= EZ Z T ∇z {EZ1 [f ]}
t
h
i
1
β−1
= EZ ∆z {EZ1 [f ]}
t h
i
= tEZ ∆w {EZ1 [f ]}

β−1

(2.19)

,

where the fourth equality follows from the Gauss divergence theorem. By (2.17),
(2.18) and (2.19), we have
h
i
∂
β−1
β−2
g(t; w) = tEZ ∆w {EZ1 [f ]}
+ (β − 1) {EZ1 [f ]}
∆w EZ1 [f ] . (2.20)
∂t
Recall the formula of Laplacian for a function h(u),


∆u h(u)
a
a
2
∆u h (u) = ah (u)
+ (a − 1)k∇u log h(u)k ,
h(u)

(2.21)

for a 6= 0. Then, in (2.20), we have
∆w {EZ1 [f ]}
=
=
=

β−1

(β − 1)
1−β

{EZ1 [f ]}
2(β − 1)
{EZ1 [f ]}

1−β

2(β − 1)
{EZ1 [f ]}

1−β

β−2

+ (β − 1) {EZ1 [f ]}
∆w EZ1 [f ]


∆w EZ1 [f ]
+ (β − 2)k∇w log EZ1 [f ] k2
2
EZ1 [f ]


∆w EZ1 [f ]
2
+ (β/2 − 1) k∇w log EZ1 [f ] k
EZ1 [f ]
∆w {EZ1 [f ]}

(2.22)

β/2
β/2

(β/2) {EZ1 [f ]}
β/2

=

4(β − 1) ∆w {EZ1 [f ]}
.
1−β/2
β
{EZ1 [f ]}

By (2.15), (2.20) and (2.22), we completes the proof.
Remark 2.1. In the previous version of this article as well as George, Liang and
Xu (2006), not only the Stein identity but also the heat equation
∂
1
φ(u, v) = ∆u φ(u, v),
∂v
2
was efficiently applicable for deriving a nice expression of the risk difference, like
Part 1 of Theorem 2.2. It seemed to us that the heat equation was an additional

Y. Maruyama and T. Ohnishi/A Bayesian prediction under α-divergence

10

necessary tool for investigating the Stein phenomenon of predictive density. But
it is not true, the heat equation is no longer necessary. As seen in the proof of
Theorem 2.2, only the Stein identity or the Gauss divergence theorem is the
key, as in Stein “estimation” problem.
The superharmonicity of f implies the superharmonicity of EZ1 [f (tZ1 + u)].
Furthermore, using the relationship (2.21), we see that the superharmonicity of
EZ1 [f (tZ1 + u)] implies the superharmonicity of
%(u; t; f ) = {EZ1 [f (tZ1 + u)]}β/2
for β ∈ (0, 1). Hence, for Part 2 of Theorem 2.2, we have a following corollary.
Corollary 2.1. Suppose f : Rd → R+ is superharmonic. Then the predictive
density p̂f (y | x; α) given by (2.9) as
p̂f (y | x; α) =

f (γx + (1 − γ)y)
p̂U (y | x; α),
EZ1 [f (x + ξZ1 )]

dominates p̂U (y | x; α).
In Section 3, we will investigate the properties of the Bayesian predictive
density p̂π (y | x; α) where
f (u) = {mπ (u, vx γ)}1/β
is assumed in Theorem 2.2 and Corollary 2.1. Actually in this case, Corollary 2.1
is not useful since the superharmonicity of {mπ (u, vx γ)}1/β for β ∈ (0, 1) is very
restrictive. Recall the relationship given by (2.21). For example, the superharmonicity of mπ (u, vx γ) does not imply the superharmonicity of {mπ (u, vx γ)}1/β .
Hence, in Section 3, we will seriously consider the superharmonicity of
n
h
ioβ/2
1/β
%(u; t; m1/β
.
π ) = EZ1 {mπ (tZ1 + u, vx γ)}


Further, When 1/β = 2/(1−α) is not an integer, EZ1 {mπ (tZ1 + u, vx γ)}1/β
in Part 2 of Theorem 2.2 is not tractable for our current methodology in Section
3. Thus we propose a variant of Theorem 2.2 with f (u) = {mπ (u, vx γ)}1/β , for
a non-integer 1/β as follows. Let κ be the smallest integer among integers which
is strictly greater than 1/β,
κ = min{n ∈ Z | n > 1/β}.

(2.23)

Then κ − 1 < 1/β < κ. As in (2.15), the risk difference is expressed as
Rα {φ(y − µ, vy ) || p̂U (y | x; α)} − Rα {φ(y − µ, vy ) || p̂π (y | x; α)}


"
1/β #β−1
m
(W
+
ξ(Z
+
Z),
v
γ)
γ (1−β)d/2
π
1
x
EW,Z EZ1
− 1
=
β(1 − β)
mπ (W, vx γ)

Y. Maruyama and T. Ohnishi/A Bayesian prediction under α-divergence

11

where W ∼ Nd (µ, vx γI), Z1 ∼ Nd (0, I) and Z ∼ Nd (0, I). From Jensen’s inequality, we have
h
i
EZ1 m1/β
π (w + ξ(Z1 + Z), vx γ)
i
h
(2.24)
= EZ1 {mκπ (w + ξ(Z1 + Z), vx γ)}1/(βκ)
1/(βκ)

≤ {EZ1 [mκπ (w + ξ(Z1 + Z), vx γ)]}

,

since 0 < 1/(βκ) < 1 and hence
Rα {φ(y − µ, vy ) || p̂U (y | x; α)} − Rα {φ(y − µ, vy ) || p̂π (y | x; α)}
#
"
(β−1)/(βκ)
 κ
γ (1−β)d/2
mπ (W + ξ(Z1 + Z), vx γ)
−1 .
EW,Z EZ1
≥
β(1 − β)
mκπ (W, vx γ)
Applying the same technique starting (2.15) through (2.22) to the lower bound
above, we have a variant of Part 2 of Theorem 2.2.
Theorem 2.3. Assume 1/β is not a positive integer. Let κ be the smallest
integer greater than 1/β. A sufficient condition for diffRα,U,π ≥ 0 is
c(β)/κ

∆u {EZ1 [mκπ (tZ1 + u, vx γ)]}

≤ 0,

∀u ∈ Rd , 0 ≤ ∀t ≤ ξ

(2.25)

where Z1 ∼ Nd (0, I) and
c(β) =

κ − 1/β + 1
∈ (1/2, 1).
2

(2.26)

2.1. Asymptotics
1/β

In this subsection, using Theorem 2.2 with f = mπ , we investigate asymptotics
of the risk difference
diffRα,U,π = Rα {φ(y − µ, vy ) || p̂U (y | x; α)} − Rα {φ(y − µ, vy ) || p̂π (y | x; α)}
where p̂U (y | x; α) and p̂π (y | x; α) are given by (2.5) and (2.6), respectively.
2.1.1. α → −1
Let v∗ = vx vy /(vx + vy ). When α → −1 or equivalently β → 1, we have
γ→
and hence

1
v∗
vx2
=
and ξ 2 →
= vx − v∗
1 + vx /vy
vx
vx + vy
2γ (1−β)d/2
{mπ (w, vx γ)}1/β−1 → 2,
β2

(2.27)

Y. Maruyama and T. Ohnishi/A Bayesian prediction under α-divergence

12

which are parts of ρ(w, z) given by (2.11). Further, in %(t; u) given by (2.12),
we have
EZ1 [mπ (tZ1 + u, vx γ)] = mπ (u, vx γ + t2 ) → mπ (u, v∗ + t2 ).

(2.28)

By (2.27) and (2.28), we have
2
%(t; u) → m1/2
π (u, v∗ + t ),
√
Z vx −v∗ Z
1/2
−∆u mπ (u, v∗ + t2 )
t
EZ [ρ(w, Z)] → 4
φ(u − w, t2 )dudt
1/2
2
d
(2.29)
0
R
mπ (u, v∗ + t )
Z vx −v∗ Z
1/2
−∆u mπ (u, v∗ + t)
=2
φ(u − w, t)dudt.
1/2
0
Rd
mπ (u, v∗ + t)

By (2.29), we have
Z

Z

vx −v∗

Z

EW,Z [ρ(W, Z)] → 2
Rd

Z
=2

Rd

0

!

1/2

−∆u mπ (u, v∗ + t)
1/2

φ(u − w, t)dudt

mπ (u, v∗ + t)

× φ(w − µ, v∗ )dw
Z
1/2
−∆u mπ (u, v∗ + t)

!

vx −v∗

1/2

mπ (u, v∗ + t)
"
#
Z vx
1/2
∆z mπ (Z, v)
=2
EZ −
dv
1/2
v∗
mπ (Z, v)
Rd

0

φ(u − µ, v∗ + t)du dt

= R−1 {φ(y − µ, vy ) || p̂U (y | x; −1)} − R−1 {φ(y − µ, vy ) || p̂π (y | x; −1)},
where Z ∼ Nd (µ, vI) and v∗ = vx vy /(vx + vy ). The last equality follows from
George, Liang and Xu’s (2006) result which was already explained in (1.14) of
Section 1. Hence we have
lim

α→−1+0

diffRα,U,π = diffR−1,U,π .

2.1.2. (1 − α)vx /vy → 0
Consider the asymptotic situation where
(1 − α)vx /vy → 0 ⇔ β(vx /vy ) → 0 ⇔ γ → 1.
Note that EZ [ρ(w, Z)] is rewritten as the product ρ1 (w)ρ2 (w) where
2γ (1−β)d/2
{mπ (w, vx γ)}1/β−1 ξ 2 ,
β2

Z ξ Z
2
−∆u %(t; u)
2
ρ2 (w) = 2
t
φ(u − w, t )du dt
2/β−1 (t; u)
ξ 0
Rd %
ρ1 (w) =

(2.30)

Y. Maruyama and T. Ohnishi/A Bayesian prediction under α-divergence

=

Z

1
ξ2

ξ2

Z
Rd

0

13

√

−∆u %( t; u)
√
φ(u
−
w,
t)du
dt.
%2/β−1 ( t; u)

Since ξ 2 is rewritten as
ξ2 =

(1 − γ)2 vy
=
γ



1−γ
γ

2
vy γ =

vx2 2
β γ,
vy

(2.31)

we have
ρ1 (w) = 2

vx2 (1−β)d/2+1
γ
{mπ (w, vx γ)}1/β−1
vy

and
lim ρ1 (w) = 2

γ→1

vx2
{mπ (w, vx )}1/β−1 .
vy

(2.32)

When γ → 1, we have ξ 2 → 0 by (2.31) and hence
√
Z

−∆u %( t; u)
√
lim ρ2 (w) = lim
φ(u − w, t)du
2/β−1 ( t; u)
γ→1
t→0
Rd %
√


Z
−∆u %( t; u)
√
=
lim
δ(u − w)du,
2/β−1 ( t; u)
Rd t→0 %

(2.33)

where δ(·) is the Dirac delta function. By (2.33) and
√

Z

lim %( t; u) =

t→0
γ→1

Rd

m1/β
π (u1

β/2
+ u, vx γ)δ(u1 )du1

= m1/2
π (u, vx ),

we have


1/2−1/β
lim ρ2 (w) = −∆w m1/2
(w, vx ).
π (w, vx ) mπ

γ→1

(2.34)

By (2.32) and (2.34), we have
1/2

lim EZ [ρ(w, Z)] = lim ρ1 (w)ρ2 (w) = 2

γ→1

γ→1

vx2 −∆w mπ (w, vx )
,
vy m1/2
π (w, vx )

which implies that
lim diffRα,U,π

α→1

lim

vx /vy →0

vy
diffRα,U,π
vx

"
#
1/2
vx2
−∆w mπ (W, vx )
= diffR1,U,π = 2 E
,
1/2
vy
mπ (W, vx )
"
#
1/2
vy
−∆w mπ (W, vx )
=
diffR1,U,π = 2vx E
.
1/2
vx
mπ (W, vx )

Therefore the asymptotic situation vx /vy → 0 corresponds to the case α → 1
and ∆z {mπ (z, v)}1/2 plays an important role for general α ∈ (−1, 1).

Y. Maruyama and T. Ohnishi/A Bayesian prediction under α-divergence

14

3. Improvement under the harmonic prior
Under the harmonic prior πH (µ) = kµk−(d−2) , let
Z
mH (w, v) =
φ(w − µ, v)πH (µ)dµ.

(3.1)

Rd

Let ν be an integer larger than or equal to 2. The superharmonicity related to
EZ1 [mνH (tZ1 + u, v)] with Z1 ∼ Nd (0, I) is as follows.
Theorem 3.1. Let c ∈ (0, 1) and Z1 ∼ Nd (0, I). Let ν be an integer larger than
or equal to 2. Then, we have
c/ν

∆u {EZ1 [mνH (tZ1 + u, v)]}

≤ 0,

when

0≤t≤

(d + 2)(1 − c)v
dν(ν − 1)

∀u ∈ Rd ,

1/2
.

(3.2)

Proof. Section B of Appendix.
When 1/β is an integer larger than or equal to 2, namely,
α = 0, 1/3, 1/2, 3/5, 2/3, . . . ,
β = 1/2, 1/3, 1/4, 1/5, 1/6, . . . ,

(3.3)

let ν = 1/β, v = vx γ and c = 1/2 in Theorem 3.1 and compare (3.2) in Theorem
3.1 with 0 ≤ t2 ≤ ξ 2 = β 2 vx2 γ/vy in Theorem 2.2. If
β 2 vx
(d + 2)(1 − c)
vx γ
vx γ ≤
vy
dν(ν − 1)
or equivalently
vx
d+2
d+2
=
,
≤
vy
d(1 + α)
2d(1 − β)
mH (w, vx γ) satisfies the sufficient condition of Theorem 2.2 and we have a following result of the Bayesian predictive density with respect to Stein’s harmonic
prior πH (µ) = kµk−(d−2) , which is given by
1/β

p̂H (y | x; α) =

mH (γx + (1 − γ)y, vx γ)
i p̂U (y | x; α).
h
1/β
EZ1 mH (x + ξZ1 , vx γ)

(3.4)

Theorem 3.2. Suppose 2/(1−α) is an positive integer for α ∈ (−1, 1). Suppose
vx
d+2
≤
.
vy
d(1 + α)

(3.5)

Then, under α-div loss, the Bayesian predictive density p̂H (y | x; α) with respect
to the harmonic prior πH (µ) = kµk−(d−2) dominates the best invariant Bayesian
predictive density p̂U (y | x; α) = φ(y − x, vy /γ).

Y. Maruyama and T. Ohnishi/A Bayesian prediction under α-divergence

15

Remark 3.1. For any d ≥ 3 and α ∈ (−1, 1), we have
1
d+2
> .
d(1 + α)
2
Note that, in most typical situations,
vx
1
≤ ,
vy
2
is easily assumed as follows. Suppose that we have a set of observations x1 , . . . , xn
from Nd (µ, σ 2 I). An unobserved set xn+1 , . . . , xn+m from the same distribution
is predicted by using a predictive density as a function of x1 , . . . , xn . From
sufficiency,
Xn
Xm
x = n−1
xi ∼ Nd (µ, σ 2 I/n) and y = m−1
xn+i ∼ Nd (µ, σ 2 I/m)
i=1

i=1

and clearly vx /vy = m/n in this case. Since, m is typically 1 or 2 whereas n is
relatively large, the condition (3.5) is satisfied.
When β = 2/(1 − α) is not an integer, Theorem 2.3 can be applied. Let κ be
the smallest integer greater than 1/β. Suppose
β2

(d + 2){1 − c(β)}vx γ
vx
vx γ ≤
,
vy
dκ(κ − 1)

(3.6)

where c(β) is given by (2.26) as c(β) = c({1 − α}/2) = {κ − 2/(1 − α) + 1}/2,
the left-hand side is the upper bound of t of Theorem 2.3 and the right-hand
side is the upper bound of t of Theorem 3.1. When
vx
≤
vy



2
1−α

2

d + 2 1 − {κ − 2/(1 − α)}
,
d
2κ(κ − 1)

which is equivalent to (3.6), mH (w, vx γ) satisfies the sufficient condition of Theorem 2.3 and we have a following result.
Theorem 3.3. Suppose 2/(1 − α) is not an positive integer for α ∈ (−1, 1). Let
κ be the smallest integer greater than 2/(1 − α). Suppose
vx
≤
vy



2
1−α

2

d + 2 1 − {κ − 2/(1 − α)}
.
d
2κ(κ − 1)

(3.7)

Then the Bayesian predictive density p̂H (y | x; α) with respect to the harmonic
prior πH (µ) = kµk−(d−2) dominates the best invariant Bayesian predictive density p̂U (y | x; α) = φ(y − x, vy /γ).
By the definition of κ,
κ−1<

2
< κ.
1−α

16

0.8
0.6
0.4
0.0

0.2

upper bound of vx /vy

1.0

1.2

Y. Maruyama and T. Ohnishi/A Bayesian prediction under α-divergence

−1.0

−0.5

0.0

0.5

α
Fig 1. The upper bound of vx /vy in Theorems 3.2 and 3.3

As 2/(1 − α) ↑ κ, the upper bound given by (3.7) approaches (d + 2)/{d(1 + α)}
which is exactly the upper bound given by (3.5) of Theorem 3.2. On the other
hand, as 2/(1 − α) ↓ κ − 1, the upper bound given by (3.7) approaches 0.
Figure 1 gives a graph of behavior of the upper bound of vx /vy for improvement in Theorems 3.2 and 3.3. This undesirable discontinuity with respect to
the upper bound of Theorem 3.3 is due to Jensen’s inequality (2.24) which
was not used in the proof of Theorem 2.2. However, we would like to emphasize that, for any α ∈ (−1, 1), there exists a positive upper bound of of
vx /vy for improvement. We can naturally make a conjecture that the lower
bound of vy /vx for improvement, d(1 + α)/(d + 2), of Theorem 3.2 is still
valid even if 2/(1 − α) is not an
 integer. For that purpose,
 the methodology
for appropriately treating EZ1 {mH (tZ1 + u, vx γ)}2/(1−α) or more generally


EZ1 {mπ (tZ1 + u, vx γ)}2/(1−α) for non-integer 2/(1 − α) is needed and it remains an open problem.
Appendix A: Minimaxity of p̂U (y | x; α)
In this section, we show that
p̂U (y | x; α) = φ(y − x, vy /γ) = φ(y − x, vy + βvx )

(A.1)

is minimax, by following Sections II and III of Liang and Barron (2004). We
start with the definition of invariance under location shift.

Y. Maruyama and T. Ohnishi/A Bayesian prediction under α-divergence

17

Definition A.1. A predictive density p̂(y | x) is invariant under location shift,
if for all a ∈ Rd and all x, y, p̂(y + a | x + a) = p̂(y | x).
Hence any invariant predictive density should be of the form
p̂(y | x) = q(y − x)
which satisfies
Z
q(y)dy = 1.
Rd

Clearly p̂U (y | x; α) is invariant under location shift. Note that invariant procedures have constant risk since the risk of the invariant predictive density q(y −x)
is
Rα {φ(y − µ, vy ) || q(y − x)}



Z Z
q(y − x)
φ(y − µ, vy )dy φ(x − µ, vx )dx
=
fα
φ(y − µ, vy )
Rd
Rd



Z Z
q(zy − zx )
=
fα
φ(zy , vy )dzy φ(zx , vx )dzx
φ(zy , vy )
Rd
Rd

(A.2)

where zx = x − µ and zy = y − µ, which does not depend on µ. More specifically,
the risk of the invariant predictive density q(y − x) is as follows.
Lemma A.1. The risk of an invariant predictive density q(y − x) is
Rα {φ(y − µ, vy ) || q(y − x)}
=

1 − γ (1−β)d/2
+ γ (1−β)d/2 Dα {φ(z, vy /γ) || q(z)} .
β(1 − β)

(A.3)

Proof. By (A.2) and the definition of α-div loss,
Rα {φ(y − µ, vy ) || q(y − x)}


Z Z
1
1−
q 1−β (y − x)φβ (y, vy )φ(x, vx )dxdy .
=
β(1 − β)
Rd Rd
By the identity (2.4) with µ = 0, we have
φ (x, vx ) φβ (y, vy ) = γ (1−β)d/2 φ(γx + (1 − γ)y, vx γ)φβ (y − x, vy /γ),
and hence
1
Rα {φ(y − µ, vy ) || q(y − x)} =
β(1 − β)


1−γ

(1−β)d/2

Z
Rd

Z

q 1−β (y − x)

Rd

× φβ (y − x, vy /γ)φ(γx + (1 − γ)y, vx γ)dxdy .
By the change of variables,
  
w
γId
=
z
−Id

(1 − γ)Id
Id

 
x
y

(A.4)

Y. Maruyama and T. Ohnishi/A Bayesian prediction under α-divergence

18

where Jacobian of the matrix is 1, we have
Rα {φ(y − µ, vy ) || q(y − x)}


Z Z
1
=
1 − γ (1−β)d/2
q 1−β (z)φβ (z, vy /γ)φ(w, vx γ)dzdw
β(1 − β)
Rd Rd


Z
1
1 − γ (1−β)d/2
q 1−β (z)φβ (z, vy /γ)dz
=
β(1 − β)
Rd
=

1 − γ (1−β)d/2
+ γ (1−β)d/2 Dα {φ(z, vy /γ) || q(z)} .
β(1 − β)

In (A.3) of Lemma A.1, Dα {φ(z, vy /γ) || q(z)} is non-negative and takes
zero if and only if q(z) = φ(z, vy /γ). Hence the best invariant procedure is
p̂U (y | x; α) = φ(y − x, vy /γ), where the constant risk is
1 − γ (1−β)d/2
.
β(1 − β)
Since the risk is constant for invariant predictive density, the best invariant p̂U (y |
x; α) is the minimax procedure among all invariant procedures. If a constant risk
procedure is shown to have an extended Bayes property defined in the below,
then it is, in fact, minimax over all procedures. See Theorem 5.18 of Berger
(1985) and Theorem 5.1.12 of Lehmann and Casella (1998) for the detail.
Definition A.2. A predictive procedure p̂∗ (y | x) is called extended Bayes, if
there exists a sequence of Bayes procedures p̂πc (y | x; α) with proper prior densities πc (µ) for c = 1, . . . , such chat their Bayes risk differences go to zero, that
is,
Z
lim
Rα {φ(y − µ, vy ) || p̂∗ (y | x)}πc (µ)dµ
c→∞
Rd

Z
−
Rα {φ(y − µ, vy ) || p̂πc (y | x; α)}πc (µ)dµ = 0.
Rd

Recall that
Z
p̂π (y | x; α) ∝

1/β
φβ (y − µ, vy )φ(x − µ, vx )π(µ)dµ

(A.5)

Rd

for β = (1 − α)/2 and α ∈ (−1, 1). Under the prior µ ∼ Nd (0, {cvx γ}I) with the
density πc (µ) = φ(µ, cvx γ), the Bayesian solution is


cγ
1+c
p̂πc (y | x; α) = φ y −
x, vy
1 + cγ
1 + cγ

Y. Maruyama and T. Ohnishi/A Bayesian prediction under α-divergence

19

by the identity
φβ (y − µ, vy )φ (x − µ, vx ) φ(µ, cvx γ)
d(1−β)/2 


γx + (1 − γ)y cvx γ
1 + cγ
φ µ−c
,
=
1+c
1+c
1+c


cγx
1
+
c
× φβ y −
, vy
φ (x, vx (1 + cγ)) .
1 + cγ
1 + cγ

(A.6)

Furthermore, by the identity (A.6), the Bayes risk of p̂πc (y | x; α) is given by

1−β
Z Z Z 
1
p̂πc (y | x; α)
1−
β(1 − β)
φ(y − µ, vy )
Rd Rd Rd

× φ (x − µ, vx ) φ(y − µ, vy )φ(µ, cvx γ)dxdydµ
(

d(1−β)/2 Z Z Z


1 + cγ
γx + (1 − γ)y cvx γ
1
1−
,
φ µ−c
=
β(1 − β)
1+c
1+c
1+c
Rd Rd Rd



cγx
1+c
× φ y−
, vy
φ (x, vx (1 + cγ)) dµdydx
1 + cγ
1 + cγ
(

d(1−β)/2 )
1
1 + cγ
=
1−
,
β(1 − β)
1+c
which approaches (1 − γ (1−β)d/2 )/{β(1 − β)} as c goes to infinity, the constant
risk of p̂U (y | x; α). Hence p̂U (y | x; α) is extended Bayes and hence minimax.
Appendix B: Proof of Theorem 3.1
Recall the identity
kµk−(d−2) = b

Z
0

∞



kµk2
g d/2−2 exp −g
dg
2v

(B.1)

for any v > 0, where b = 1/{Γ(d/2 − 1)2d/2−1 v d/2−1 }. Then we have
Z
mH (w, v) =
φ(w − µ, v)kµk−(d−2) dµ
d
R


Z ∞
Z
kµk2
1
kw − µk2
exp
−
−
g
dµ
=b
g d/2−2 dg
d/2 v d/2
2v
2v
0
Rd (2π)


Z ∞
g d/2−2
gkwk2
=b
exp −
dg
2(g + 1)v
(1 + g)d/2
0


Z 1
λkwk2
=b
λd/2−2 exp −
dλ,
2v
0

Y. Maruyama and T. Ohnishi/A Bayesian prediction under α-divergence

20

where the third equality is from the relation of completing squares with respect
to µ
kw − µk2 + gkµk2 = (g + 1)kµ − w/(g + 1)k2 + {g/(g + 1)}kwk2
and the fourth equality is from the transformation λ = g/(g + 1).
Note that mνH (w, v) for a positive integer ν is expressed as
ν

mH (w, v) = b

Z

ν

ν
Y

Dν i=1

d/2−2
λi

Y
 Pν
2
i=1 λi kwk
dλi ,
exp −
2v

Qν
where Dν is ν-dimensional unit hyper-cube. In the following, dλ denotes Q
i=1 dλi
for
notational
simplicity.
Furthermore
the
subscript
and
superscript
of
and
P
is omitted for simplicity if they are i = 1 and i = ν respectively. Hence
mνH (w, v) in the above is written as
 P

Z Y
λi kwk2
d/2−2
ν
ν
λi
exp −
mH (w, v) = b
dλ.
2v
Dν
For the calculation of
EZ1 [mνH (tZ1 + u, v)] =

Z

mνH (x + u, v)φ(x, t2 )dx

(B.2)

Rd

under Z1 ∼ Nd (0, I), note the relation of completing squares with respect to x,
P
o
kxk2
1 nX
( λi ) kx + uk2
+ 2 =
λi kx + uk2 + skxk2
t
v
)
(v
P
P
(B.3)
2

s λi
1 X
λi
2
=
u +P
kuk ,
λi + s x + P
λi + s
λi + s
v
where s = v/t2 . Then, by (B.3), we have
EZ1 [mνH (tZ1 + u, v)] =

bν v d/2
td

Z
Dν

Q d/2−2
P


λ
s λi kuk2
P
P i
exp
−
dλ.
v( λi + s) 2
( λi + s)d/2

Re-define u := {s/v}1/2 u and let
Z
ψ(u; ν, s) =
Dν

Q d/2−2
P


λi
λi kuk2
P
P
exp −
dλ.
λi + s 2
( λi + s)d/2
c/ν

By (2.21), the super-harmonicity of {EZ1 [mνH (tZ1 + u, v)]}
u ∈ Rd is equivalent to
c

− 1 k∇u ψk2 + ψ∆u ψ ≤ 0, ∀u ∈ Rd .
ν

(B.4)

with respect to

(B.5)

Y. Maruyama and T. Ohnishi/A Bayesian prediction under α-divergence

21

The integrand of ψ given by (B.4) is denoted by
P
Q d/2−2


λi
λ
P
exp
−
z
ζ(λ) = ζ(λ1 , . . . , λν ) = P i
λi + s
( λi + s)d/2
where z = kuk2 /2. Then we have
∂
ψ = −uj
∂uj

Z

P
λi
dλ,
ζ(λ) P
λi + s

for j = 1, . . . , d and
∂2
ψ=
∂u2j

Z

(

P
 P
2 )
λi
λ
i
ζ(λ) − P
+ u2j P
dλ.
λi + s
λi + s

Noting z = kuk2 /2, we have
2

Z

k∇u ψk = 2z

P
2
Z
2
λ1
λi
2
dλ = 2ν z
ζ(λ) P
dλ
ζ(λ) P
λi + s
λi + s

(B.6)

and
P
 P
2
Z
λi
λi
∆u ψ = −d ζ(λ) P
dλ + 2z ζ(λ) P
dλ
λi + s
λi + s
Z
Z
λ1
λ2
dλ
= −dν ζ(λ) P
dλ + 2νz ζ(λ) P 1
( λi + s)2
λi + s
Z
λ1 λ2
+ 2ν(ν − 1)z ζ(λ) P
dλ.
( λi + s)2
Z

(B.7)

In (B.6) and (B.7), the second equalities are from symmetry with respect to
λi ’s.
Let
Z
X
ρ(j1 , j2 , l) =
λj11 λj22 (
λi + s)l ζ(λ)dλ,
Dν
Z

l
X
Y
η(j2 , l) =
λj22 1 +
λi + s ζ(1, λ2 , . . . , λν )
dλi ,
Dν−1

i=2

i=2

where j1 and j2 are nonnegative integers. Then k∇u ψk2 and ∆u ψ given by (B.6)
and (B.7) is rewritten as
k∇u ψk2 = 2ν 2 zρ(1, 0, −1)2 ,
∆u ψ = −dνρ(1, 0, −1) + 2νzρ(2, 0, −2) + 2ν(ν − 1)zρ(1, 1, −2).
Here are some useful relationships and inequalities.

(B.8)

Y. Maruyama and T. Ohnishi/A Bayesian prediction under α-divergence

22

Lemma B.1.
szρ(j1 , j2 , l) = −η(j2 , l + 2) + (j1 + d/2 − 2)ρ(j1 − 1, j2 , l + 2)
+ (l − d/2 + 2)ρ(j1 , j2 , l + 1), for j1 ≥ 1,
ρ(0, 0, l) = νρ(1, 0, l − 1) + sρ(0, 0, l − 1),

(B.9)
(B.10)

ρ(1, 0, l) = ρ(2, 0, l − 1) + (ν − 1)ρ(1, 1, l − 1) + sρ(1, 0, l − 1), (B.11)
η(0, 1) = η(0, 0) + (ν − 1)η(1, 0) + sη(0, 0),
η(0, 1)ρ(0, 0, −1) ≥ η(0, 0)ρ(0, 0, 0),
ρ(1, 0, −1)
1
≥
.
ρ(1, 0, 0)
νd/(d + 2) + s

(B.12)
(B.13)
(B.14)

Proof. See Sub-section B.1.
Applying the identity (B.9) to k∇u ψk2 and ∆u ψ given in (B.8), we have
sk∇u ψk2 = 2ν 2 {szρ(1, 0, −1)}ρ(1, 0, −1)
= ν 2 {−2η(0, 1) + (d − 2)ρ(0, 0, 1) − (d − 2)ρ(1, 0, 0)} ρ(1, 0, −1),
s∆u ψ = −dνsρ(1, 0, −1) + ν{−2η(0, 0) + dρ(1, 0, 0) − dρ(2, 0, −1)}
+ ν(ν − 1) {−2η(1, 0) + (d − 2)ρ(1, 0, 0) − dρ(1, 1, −1)}
= −2ν{η(0, 0) + (ν − 1)η(1, 0)} + ν(ν − 1)(d − 2)ρ(1, 0, 0),
where the second equality of s∆u ψ follows from (B.11). Then we have


s c−ν
2
k∇u ψk + ψ∆u ψ
ν
ν

(B.15)

= (ν − c) [2η(0, 1) − (d − 2){ρ(0, 0, 1) − ρ(1, 0, 0)}] ρ(1, 0, −1)
− 2{η(0, 0) + (ν − 1)η(1, 0)}ρ(0, 0, 0) + (ν − 1)(d − 2)ρ(1, 0, 0)ρ(0, 0, 0).
By applying (B.10), (B.12) and (B.13), the terms of (B.15) including η(·, ·),
divided by 2, is
(ν − c)η(0, 1)ρ(1, 0, −1) − {η(0, 0) + (ν − 1)η(1, 0)}ρ(0, 0, 0)
= (ν − c)η(0, 1)ρ(1, 0, −1) − {η(0, 1) − sη(0, 0)}ρ(0, 0, 0)
= (ν − c)η(0, 1)ρ(1, 0, −1) − η(0, 1) {νρ(1, 0, −1) + sρ(0, 0, −1)}
+ sη(0, 0)ρ(0, 0, 0)

(B.16)

= −cη(0, 1)ρ(1, 0, −1) − s{η(0, 1)ρ(0, 0, −1) − η(0, 0)ρ(0, 0, 0)}
≤ 0,
where the first equality follows from (B.12), the second equality follows from
(B.10) and the inequality follows from (B.13).

Y. Maruyama and T. Ohnishi/A Bayesian prediction under α-divergence

23

The terms of (B.15) not including η(·, ·), divided by (d − 2), are rewritten as
(ν − c) {−ρ(0, 0, 1) + ρ(1, 0, 0)} ρ(1, 0, −1) + (ν − 1)ρ(1, 0, 0)ρ(0, 0, 0)
= −(ν − c)(ν − 1)ρ(1, 0, 0)ρ(1, 0, −1) − (ν − c)sρ(0, 0, 0)ρ(1, 0, −1)
+ (ν − 1)ρ(1, 0, 0)ρ(0, 0, 0)


(ν − c)s
≤−
− (ν − 1) ρ(1, 0, 0)ρ(0, 0, 0)
νd/(d + 2) + s
(1 − c)s − ν(ν − 1)d/(d + 2)
=−
ρ(1, 0, 0)ρ(0, 0, 0),
νd/(d + 2) + s

(B.17)

which is nonpositive for s ≥ ν(ν − 1)d/{(1 − c)(d + 2)}, where the first equality
follows from (B.10) and the inequality follows from (B.14).
By (B.16) and (B.17), we have

c
− 1 k∇u ψk2 + ψ∆u ψ ≤ 0, ∀u ∈ Rd
ν
or equivalently
∆u {EZ1 [mνH (tZ1 + u, v)]}

c/ν

≤ 0,

∀u ∈ Rd ,

when t ≤ {(d + 2)(1 − c)v/{dν(ν − 1)}}1/2 .
B.1. Proof of Lemma B.1
[Part of (B.9)]

Note
P
P




∂
sz
z λi
z λi
P
=− P
exp
−
.
exp − P
λi + s
λi + s
∂λ1
( λi + s)2

(B.18)

Then, by an integration by parts, we have
Z 1
X
sz
λj11 λj22 (
λi + s)l ζ(λ)dλ1
0
d/2−2+j2

= −λ2

Y
i=3

d/2−2

Z

λi

1

d/2−2+j1

λ1

(

X

λi + s)l−d/2+2

0

P


z λi
∂
exp − P
dλ1
×
∂λ1
λi + s
(
P

1
X
Y d/2−2
z λi
d/2−2+j2
d/2−2+j1
= −λ2
λi
λ1
(
λi + s)l−d/2+2 exp − P
λi + s 0
i=3
P


Z 1
X
z λi
d/2−3+j1
l−d/2+2
P
−(d/2 − 2 + j1 )
λ1
(
λi + s)
exp −
dλ1
λi + s
0
P



Z 1
X
z λi
d/2−2+j1
l−d/2+1
dλ1 .
−(l − d/2 + 2)
λ1
(
λi + s)
exp − P
λi + s
0


Y. Maruyama and T. Ohnishi/A Bayesian prediction under α-divergence

24

(B.9) follows from integration with respect to λ2 , . . . , λν in the both hand side
of the above equality.
[Parts of (B.10), (B.11) and (B.12)] The equalities (B.10), (B.11) and (B.12)
easily follows from symmetry with respect to λi ’s.
[Part of (B.13)]

Note that (B.13) is equivalent to

η(0, 0)ρ(0, 0, 0) − η(0, 1)ρ(0, 0, −1)
= {ρ(0, 0, 0) − ρ(0, 0, −1)}η(0, 1) − {η(0, 1) − η(0, 0)}ρ(0, 0, 0)
Z
Z
Y
Y
=
f1 (λ2 , . . . , λν )
dλi
f2 (λ2 , . . . , λν )
dλi
Dν−1

−

Dν−1

i=2

Z
f3 (λ2 , . . . , λν )
Dν−1

Y
i=2

i=2

Z
dλi

f4 (λ2 , . . . , λν )
Dν−1

Y

dλi

i=2

≤ 0,
where
1



1
ζ(λ1 , . . . , λν )dλ1
1− P
λi + s
0
X
f2 (λ2 , . . . , λν ) = (1 +
λi + s)ζ(1, λ2 , . . . , λν )
i=2
X
f3 (λ2 , . . . , λν ) = (
λi + s)ζ(1, λ2 , . . . , λν )
i=2
Z 1
f4 (λ2 , . . . , λν ) =
ζ(λ1 , . . . , λν )dλ1 .
Z

f1 (λ2 , . . . , λν ) =

0

P
P
Since both 1−1/ ( λi + s) and λi +s are increasing in each of its arguments,
we have
o
n
X
X
λi + s (1 +
ξi + s)
1 − 1/
i=2

n
o
X
1
P
(λ1 ∨ 1) +
(λi ∨ ξi ) + s
≤ 1−
i=2
(λ1 ∨ 1) + i=2 (λi ∨ ξi ) + s
X
=
(λi ∨ ξi ) + s,
(B.19)
i=2

where ∨ is the maximum operator, i.e. λi ∨ ξi = max(λi , ξi ). In the following, ∧
denotes the minimum operator, i.e. λi ∧ ξi = min(λi , ξi ). Note that a function
h: Rν → R is said to be multivariate totally positive of order two (MTP2) if it
satisfies
h(x1 , . . . , xν )h(y1 , . . . , yν ) ≤ h(x1 ∨ y1 , . . . , xν ∨ yν )h(x1 ∧ y1 , . . . , xν ∧ yν )
for any x, y ∈ Rν . By Lemma B.2 below, ζ(λ1 , . . . , λν ) is MTP2 as a function
of ν-variate function and hence the inequality
ζ(λ1 , . . . , λν )ζ(1, ξ2 , . . . , ξν )
≤ ζ(λ1 ∨ 1, λ2 ∨ ξ2 , . . . , λν ∨ ξν )ζ(λ1 ∧ 1, λ2 ∧ ξ2 , . . . , λν ∧ ξν )
= ζ(1, λ2 ∨ ξ2 , . . . , λν ∨ ξν )ζ(λ1 , λ2 ∧ ξ2 , . . . , λν ∧ ξν )

(B.20)

Y. Maruyama and T. Ohnishi/A Bayesian prediction under α-divergence

25

follows. By (B.19) and (B.20), we have
f1 (λ2 , . . . λν )f2 (ξ2 , . . . ξν )
Z 1 h nX
o
≤
(λi ∨ ξi ) + s ζ(1, λ2 ∨ ξ2 , . . . , λν ∨ ξν )
i=2
0
i
× ζ(λ1 , λ2 ∧ ξ2 , . . . , λν ∧ ξν ) dλ1

(B.21)

= f3 (λ2 ∨ ξ2 , . . . , λν ∨ ξν )f4 (λ2 ∧ ξ2 , . . . , λν ∧ ξν ).
From Theorem B.1 below, shown by Karlin and Rinott (1980), the theorem
follows.
[Part of (B.14)]

By Jensen’s inequality, we have
Z
1
λ1 ζ(λ)
ρ(1, 0, −1)
Pν
=
dλ
ρ(1, 0, 0)
λ1 + i=2 λi + s ρ(1, 0, 0)
1
.
≥
ρ(2, 0, 0)
ρ(1, 1, 0)
+ (ν − 1)
+s
ρ(1, 0, 0)
ρ(1, 0, 0)

(B.22)

Let f be a probability density given by
ν−1

ν
Y
d d
d/2−2
d/2−1
λi
,
f (λ1 , . . . , λν ) =
−1
λ1
2 2
i=2
which is clearly MTP2. Also let
g1 (λ1 , . . . , λν ) = λ1 ,

g2 (λ1 , . . . , λν ) = −

P
exp (sz/{ λi + s})
P
,
( λi + s)d/2

which are both increasing increasing in each of its arguments. Hence, by so-called
FKG inequality given in Theorem B.2 below,
Z
g1 (λ1 , . . . , λν )g2 (λ1 , . . . , λν )f (λ1 , . . . , λν )dλ
Dν
Z
Z
≥
g1 (λ1 , . . . , λν )f (λ1 , . . . , λν )dλ
g2 (λ1 , . . . , λν )f (λ1 , . . . , λν )dλ
Dν

Dν

or equivalently
R

g1 (λ1 , . . . , λν )g2 (λ1 , . . . , λν )f (λ1 , . . . , λν )dλ
R
g (λ , . . . , λν )f (λ1 , . . . , λν )dλ
Dν 2 1
Z
≤
g1 (λ1 , . . . , λν )f (λ1 , . . . , λν )dλ,
Dν

Dν

since g2 < 0. Since ρ(2, 0, 0)/ρ(1, 0, 0) is expressed as
R
g1 (λ1 , . . . , λν )g2 (λ1 , . . . , λν )f (λ1 , . . . , λν )dλ
ρ(2, 0, 0)
R
= Dν
,
ρ(1, 0, 0)
g (λ , . . . , λν )f (λ1 , . . . , λν )dλ
Dν 2 1

Y. Maruyama and T. Ohnishi/A Bayesian prediction under α-divergence

we have

26

ρ(2, 0, 0)
d
≤
.
ρ(1, 0, 0)
d+2

(B.23)

d−2
d
ρ(1, 1, 0)
≤
≤
.
ρ(1, 0, 0)
d
d+2

(B.24)

Similarly we have

Hence, by (B.22), (B.23) and (B.24), we have
ρ(1, 0, −1)
1
≥
.
ρ(1, 0, 0)
νd/(d + 2) + s
Lemma B.2. Let
Q d/2−2
P


λi
λi
exp − P
z .
ζ(λ1 , . . . , λν ) = P
λi + s
( λi + s)d/2
Then ζ(λ1 , . . . , λν ) is MTP2.
Proof. Note
P




sz
λi
z = exp(−z) exp P
.
exp − P
λi + s
λi + s
From the form of ζ, we have only to check
X
X
X
X
(
λi + s)(
ξi + s) ≥ (
λi ∨ ξi + s)(
λi ∧ ξi + s)
or equivalently
X
X
X
X
(
λi )(
ξi ) ≥ (
λi ∨ ξi )(
λi ∧ ξi ).
We have
X
X
X
X
(
λi )(
ξi ) − (
λi ∨ ξi )(
λi ∧ ξi )
X
=
{λi ξj + λj ξi − (λi ∨ ξi )(λj ∧ ξj ) − (λj ∨ ξj )(λi ∧ ξi )} .
i6=j

Without the loss of generality, assume λi ≥ ξi . Then we have
λi ξj + λj ξi − (λi ∨ ξi )(λj ∧ ξj ) − (λj ∨ ξj )(λi ∧ ξi )
= λi ξj + λj ξi − λi (λj ∧ ξj ) − (λj ∨ ξj )ξi
= λi {ξj − (λj ∧ ξj )} − ξi {(λj ∨ ξj ) − λj }
= (λi − ξi ){ξj − (λj ∧ ξj )}
≥ 0,
which completes the proof.

Y. Maruyama and T. Ohnishi/A Bayesian prediction under α-divergence

27

Theorem B.1 (Theorem 2.1 of Karlin and Rinott (1980)). Let f1 , f2 , f3 and
f4 be nonnegative functions satisfying for all x, y ∈ Rν
f1 (x)f2 (y) ≤ f3 (x ∨ y)f4 (x ∧ y).
Then
Z

Z
f1 (x)dx

Z
f2 (x)dx ≤

Z
f3 (x)dx

f4 (x)dx.

Theorem B.2 (FKG Inequality, e.g. Theorem 2.3 of Karlin and Rinott (1980)).
Let f (x) for x ∈ Rν be a probability density satisfying MTP2. Then for any pair
of increasing functions g1 (x) and g2 (x), we have
Z
Z
Z
g1 (x)g2 (x)f (x)dx ≥ g1 (x)f (x)dx g2 (x)f (x)dx.
References
Aitchison, J. (1975). Goodness of prediction fit. Biometrika 62 547–554.
MR0391353
Berger, J. O. (1985). Statistical decision theory and Bayesian analysis, second
ed. Springer Series in Statistics. Springer-Verlag, New York.
Brown, L. D. (1979). A heuristic method for determining admissibility of
estimators—with applications. Ann. Statist. 7 960–994. MR536501
Brown, L. D., George, E. I. and Xu, X. (2008). Admissible predictive density estimation. Ann. Statist. 36 1156–1170. MR2418653
Corcuera, J. M. and Giummolè, F. (1999). A generalized Bayes rule for
prediction. Scand. J. Statist. 26 265–279. MR1707658
Csiszár, I. (1967). Information-type measures of difference of probability distributions and indirect observations. Studia Sci. Math. Hungar. 2 299–318.
MR0219345
George, E. I., Liang, F. and Xu, X. (2006). Improved minimax predictive
densities under Kullback-Leibler loss. Ann. Statist. 34 78–91. MR2275235
Karlin, S. and Rinott, Y. (1980). Classes of orderings of measures and related correlation inequalities. I. Multivariate totally positive distributions. J.
Multivariate Anal. 10 467–498. MR599685
Komaki, F. (2001). A shrinkage predictive distribution for multivariate normal
observables. Biometrika 88 859–864. MR1859415
Lehmann, E. L. and Casella, G. (1998). Theory of point estimation, second
ed. Springer Texts in Statistics. Springer-Verlag, New York. MR1639875
Liang, F. and Barron, A. (2004). Exact minimax strategies for predictive
density estimation, data compression, and model selection. IEEE Trans. Inform. Theory 50 2708–2726. MR2096988
Maruyama, Y. and Strawderman, W. E. (2012). Bayesian predictive densities for linear regression models under α-divergence loss: some results and
open problems. In Contemporary developments in Bayesian analysis and

Y. Maruyama and T. Ohnishi/A Bayesian prediction under α-divergence

28

statistical decision theory: a Festschrift for William E. Strawderman. Inst.
Math. Stat. (IMS) Collect. 8 42–56. Inst. Math. Statist., Beachwood, OH.
MR3202501
Stein, C. (1974). Estimation of the mean of a multivariate normal distribution.
In Proceedings of the Prague Symposium on Asymptotic Statistics (Charles
Univ., Prague, 1973), Vol. II 345–381. Charles Univ., Prague. MR0381062
Stein, C. (1981). Estimation of the mean of a multivariate normal distribution.
Ann. Statist. 9 1135–1151. MR630098
Strawderman, W. E. (1971). Proper Bayes minimax estimators of the multivariate normal mean. Ann. Math. Statist. 42 385–388. MR0397939
Yanagimoto, T. and Ohnishi, T. (2009). Bayesian prediction of a density
function in terms of e-mixture. J. Statist. Plann. Inference 139 3064–3075.
MR2535183

