Predictive density estimation under the Wasserstein loss

TL;DR

This paper studies predictive density estimation using Wasserstein loss, showing that Bayesian predictive densities based on posterior means can outperform equivariant methods in normal models.

Contribution

It establishes the optimality of plug-in densities with posterior means under Wasserstein loss for location and scale families, and introduces Bayesian predictive densities that dominate existing methods.

Findings

Bayesian predictive densities based on posterior means are optimal under Wasserstein loss.

Plug-in densities form a complete class for the considered families.

Bayesian methods can outperform equivariant predictors in normal models.

Abstract

We investigate predictive density estimation under the $L^2$ Wasserstein loss for location families and location-scale families. We show that plug-in densities form a complete class and that the Bayesian predictive density is given by the plug-in density with the posterior mean of the location and scale parameters. We provide Bayesian predictive densities that dominate the best equivariant one in normal models.

Equations

$$R(\theta ,\hat{p})=\int L(\tilde{p}(y\mid \theta ),\hat{p}(y\mid x^n))p(x^n\mid \theta )\mathrm{d}{x}^n.$$

$$\underset{\theta }{\sup }R(\theta ,{\hat{p}}_{\ast })=\underset{\hat{p}}{\min }\underset{\theta }{\sup }R(\theta ,\hat{p}).$$

$$R(\theta ,{\hat{p}}_1)\le R(\theta ,{\hat{p}}_2)$$

$$R(\theta ,{\hat{p}}_1)<R(\theta ,{\hat{p}}_2)$$

$$\int R(\theta ,{\hat{p}}_{\pi })\pi (\theta )\mathrm{d}\theta =\underset{\hat{p}}{\min }\int R(\theta ,\hat{p})\pi (\theta )\mathrm{d}\theta .$$

$${\hat{p}}_{\pi }(y\mid x^n)=\int \tilde{p}(y\mid \theta )\pi (\theta \mid x^n)\mathrm{d}\theta =\frac{\int \tilde{p}(y\mid \theta )p(x^n\mid \theta )\pi (\theta )\mathrm{d}\theta }{\int p(x^n\mid \theta )\pi (\theta )\mathrm{d}\theta }.$$

$$W_p(p_1,p_2)=\underset{X,Y}{\inf }\mathrm{E}[d(X,Y{)}^p{]}^{1/p},$$

$$W_2(p_1,p_2)=\underset{X,Y}{\inf }\mathrm{E}[\Vert X-Y{\Vert }^2{]}^{1/2}.$$

$$W_2({\mathrm{N}}_d({\mu }_1,{\Sigma }_1),{\mathrm{N}}_d({\mu }_2,{\Sigma }_2))={\left(\Vert {\mu }_1-{\mu }_2{\Vert }^2+\mathrm{tr}\left({\Sigma }_1+{\Sigma }_2-2({\Sigma }_1^{1/2}{\Sigma }_2{\Sigma }_1^{1/2}{)}^{1/2}\right)\right)}^{1/2}.$$

$$\mathrm{tr}\left({\Sigma }_1+{\Sigma }_2-2({\Sigma }_1^{1/2}{\Sigma }_2{\Sigma }_1^{1/2}{)}^{1/2}\right)=\Vert {\Sigma }_1^{1/2}-{\Sigma }_2^{1/2}{\Vert }_F^2,$$

$$W_2(p_1,p_2)={\left(\Vert {\mu }_1-{\mu }_2{\Vert }^2+W_2(q_1,q_2{)}^2\right)}^{1/2},$$

$$W_2(p_1,p_2)\ge {\left(\Vert {\mu }_1-{\mu }_2{\Vert }^2+\mathrm{tr}\left({\Sigma }_1+{\Sigma }_2-2({\Sigma }_1^{1/2}{\Sigma }_2{\Sigma }_1^{1/2}{)}^{1/2}\right)\right)}^{1/2}.$$

$$p_2(x-{\mu }_2)=p_1(\tilde{\Sigma }(x-{\mu }_1)),$$

$${\hat{\mu }}^{\pi }(x)=\int \mu \pi (\mu \mid x)\mathrm{d}\mu =\frac{\int \mu p(x\mid \mu )\pi (\mu )\mathrm{d}\mu }{\int p(x\mid \mu )\pi (\mu )\mathrm{d}\mu }.$$

$${\hat{\mu }}^{\pi }(x,s)=\frac{\int {\sigma }^{-2}\mu \pi (\mu ,\sigma \mid x,s)\mathrm{d}\mu \mathrm{d}\sigma }{\int {\sigma }^{-2}\pi (\mu ,\sigma \mid x,s)\mathrm{d}\mu \mathrm{d}\sigma }=\frac{\int {\sigma }^{-2}\mu p(x,s\mid \mu ,\sigma )\pi (\mu ,\sigma )\mathrm{d}\mu \mathrm{d}\sigma }{\int {\sigma }^{-2}p(x,s\mid \mu ,\sigma )\pi (\mu ,\sigma )\mathrm{d}\mu \mathrm{d}\sigma }.$$

$$\theta \mid \lambda ,\tau \sim {\mathrm{N}}_d(0,{\tau }^{-1}{\lambda }^{-1}(1-\lambda )I_d),$$

$$\lambda \sim {\lambda }^{-1}{1}_{[0,1]},\tau \sim {\tau }^{-1}{1}_{(0,\infty )},$$

$${\hat{\mu }}^{\pi }(x,s)=\int \mu \pi (\mu ,\sigma \mid x,s)\mathrm{d}\mu \mathrm{d}\sigma =\frac{\int \mu p(x,s\mid \mu ,\sigma )\pi (\mu ,\sigma )\mathrm{d}\mu \mathrm{d}\sigma }{\int p(x,s\mid \mu ,\sigma )\pi (\mu ,\sigma )\mathrm{d}\mu \mathrm{d}\sigma }.$$

$$\theta \mid \lambda ,\tau \sim {\mathrm{N}}_d(0,{\tau }^{-1}{\lambda }^{-1}(1-\lambda )I_d),$$

$$\lambda \sim {\lambda }^{-1}{1}_{[0,1]},\tau \sim 1_{(0,\infty )},$$

$${\hat{\sigma }}_0(x,s)=c\sqrt{s},c=\frac{1}{\sqrt{2}}\frac{\Gamma \left(\frac{n+1}{2}\right)}{\Gamma \left(\frac{n+2}{2}\right)}.$$

$${\varphi }_0(w)=\frac{1}{\sqrt{2}}\frac{\Gamma \left(\frac{n+d+1}{2}\right)}{\Gamma \left(\frac{n+d+2}{2}\right)}\frac{{\int }_0^1{\lambda }^{d/2-1}(1+\lambda w{)}^{-(n+d+1)/2}\mathrm{d}\lambda }{{\int }_0^1{\lambda }^{d/2-1}(1+\lambda w{)}^{-(n+d+2)/2}\mathrm{d}\lambda }.$$

$$\theta \mid \lambda ,\tau \sim {\mathrm{N}}_d(0,{\tau }^{-1}{\lambda }^{-1}(1-\lambda )I_d),$$

$$\lambda \sim {\lambda }^{-1}{1}_{[0,1]},\tau \sim 1_{(0,\infty )}.$$

$$g(\tau ,x,s)$$

$$\propto$$

$${\hat{\sigma }}^{\pi }(x,s)$$

$$=\sqrt{\frac{s}{2}}\frac{\Gamma \left(\frac{n+d+1}{2}\right)}{\Gamma \left(\frac{n+d+2}{2}\right)}\frac{{\int }_0^1{\lambda }^{d/2-1}(1+\lambda w{)}^{-(n+d+1)/2}\mathrm{d}\lambda }{{\int }_0^1{\lambda }^{d/2-1}(1+\lambda w{)}^{-(n+d+2)/2}\mathrm{d}\lambda }$$

$$={\varphi }_0(w)\sqrt{s}.$$

$$p(z\mid \mu )=f(z-\mu )$$

$$\int f(z)\mathrm{d}z=1,\int zf(z)\mathrm{d}z=0.$$

Full text

1 Introduction

Suppose that we have $n$ independent observations $x^{n}=(x_{1},\ldots,x_{n})$ from a probability distribution $p(x\mid\theta)$ and predict the future observation $y$ from a probability distribution $\tilde{p}(y\mid\theta)$ by using a predictive density $\hat{p}(y\mid x^{n})$, where $\theta$ is an unknown parameter. Let $L(\tilde{p}(y\mid\theta),\hat{p}(y\mid x^{n}))$ be a loss function that measures the closeness of a predictive density $\hat{p}(y\mid x^{n})$ to the true distribution $\tilde{p}(y\mid\theta)$. Then, the risk function of the predicive density $\hat{p}(y\mid x^{n})$ is defined as the expected loss:

[TABLE]

A predictive density $\hat{p}_{*}$ is called minimax if it minimizes the maximum risk:

[TABLE]

A predictive density $\hat{p}_{1}$ is said to dominate another predictive density $\hat{p}_{2}$ if

[TABLE]

holds for every $\theta$ and also

[TABLE]

holds for some $\theta$. A predictive density $\hat{p}_{\pi}$ is called the Bayesian predictive density with respect to a prior $\pi(\theta)$ if it minimizes the average risk:

[TABLE]

Predictive density estimation under the Kullback-Leibler loss has been well studied. Aitchison (1975) showed that the Bayesian predictive density with respect to a prior $\pi(\theta)$ is given by

[TABLE]

For the normal model with known variance, the Bayesian predictive density based on the uniform prior is the best equivariant and also minimax (Liang and Barron, 2004). Komaki (2001) showed that the Bayesian predictive density based on the Stein prior dominates that based on the uniform prior in dimension larger than or equal to three. George et al. (2006) generalized this result and proved that Bayesian predictive densities based on superharmonic priors dominate that based on the uniform prior. Brown et al. (2008) gave a characterization of admissible predictive densities. For the normal model with unknown variance, Kato (2009) found a Bayesian predictive density that dominates the best equivariant one. For general location-scale families, Komaki (2007) showed that the Bayesian predictive densities based on superharmonic priors asymptotically dominate the best equivariant one.

Predictive density estimation has been investigated for other loss functions as well. Corcuera and Giummole (1999) studied the class of alpha divergence losses, which includes the Kullback-Leibler loss, and derived an explicit form of the Bayesian predictive density. Kubokawa et al. (2015) and Kubokawa et al. (2017) considered $L^{2}$ loss and $L^{1}$ loss, respectively. These studies showed that improved predictive density estimation depends on the loss function.

In this study, we investigate predictive density estimation when the loss function is defined by the Wasserstein distance, which is a distance function between probability measures based on the metric of the underlying space (Olkin and Pukelsheim, 1982; Villani, 2010) and has been widely used in machine learning and computer vision (Peyré and Cuturi, 2018). Since the Wasserstein distance is derived as the optimal transportation cost between distributions, predictive density estimation under the Wasserstein loss is suitable for resource allocation. Namely, suppose that we predict the spatial distribution of demand for some resource and then distribute it accordingly. To reduce the cost of re-transportation, predictive densities with smaller Wasserstein loss are preferable. Here, we consider $L^{2}$ Wassestein distance and focus on location families and location-scale families. For both families, plug-in densities are shown to form a complete class (Berger, 1980). For location families, we prove that the Bayesian predictive density is given by the plug-in density with the posterior mean of the location parameter. For location-scale families, we prove that the Bayesian predictive density is given by the plug-in density with the posterior mean of the location and scale parameters. We give Bayesian predictive densities that dominate the best equivariant one in normal models.

This paper is organized as follows. In Section 2, we present useful properties of the Wasserstein distance and also review improved estimation of location and scale parameters. In Sections 3 and 4, we investigate predictive density estimation under the $L^{2}$ Wasserstien loss for location families and location-scale families, respectively. In Section 5, we give concluding remarks.

2 Preliminaries

2.1 Wasserstein distance

The Wasserstein distance is a metric between probability distributions. Let $S$ be a metric space with distance $d$. The $L^{p}$ Wasserstein distance between two probability distributions $p_{1}$ and $p_{2}$ on $S$ is defined by

[TABLE]

where the infimum is taken over all joint distributions of $(X,Y)$ with marginal distributions of $X$ and $Y$ equal to $p_{1}$ and $p_{2}$, respectively. The Wasserstein distance can be interpreted as the optimal transportation cost from $p_{1}$ to $p_{2}$. See Villani (2010) for details.

In this study, we focus on the setting $S=\mathbb{R}^{d}$ and $p=2$. Thus,

[TABLE]

For finiteness of the $L^{2}$ Wasserstein distance, we restrict our attention to distributions with finite mean and covariance in the following.

For normal distributions, the $L^{2}$ Wasserstein distance is obtained in closed form.

Lemma 2.1** (Olkin and Pukelsheim (1982)).**

The $L^{2}$ Wasserstein distance between $d$-dimensional normal distributions ${\rm N}_{d}(\mu_{1},\Sigma_{1})$ and ${\rm N}_{d}(\mu_{2},\Sigma_{2})$ is given by

[TABLE]

Note that if $\Sigma_{1}$ and $\Sigma_{2}$ are commutative ($\Sigma_{1}\Sigma_{2}=\Sigma_{2}\Sigma_{1}$), then

[TABLE]

where $\|A\|_{F}=({\rm tr}(A^{\top}A))^{1/2}$ denotes the Frobenius norm of $A$.

We will also use the following properties of the $L^{2}$ Wasserstein distance.

Lemma 2.2** (Givens and Shortt (1984)).**

Let $p_{1}$ and $p_{2}$ be probability densities on $\mathbb{R}^{d}$ with mean $\mu_{1}$ and $\mu_{2}$, respectively. Then,

[TABLE]

where $q_{1}$ and $q_{2}$ are probability densities with mean zero defined by $q_{1}(x)=p_{1}(x+\mu_{1})$ and $q_{2}(x)=p_{2}(x+\mu_{2})$, respectively.

Lemma 2.3** (Gelbrich (1990)).**

Let $p_{1}$ and $p_{2}$ be probability densities on $\mathbb{R}^{d}$ with mean $\mu_{1}$ and $\mu_{2}$ and covariance matrices $\Sigma_{1}$ and $\Sigma_{2}$, respectively. Then,

[TABLE]

The equality holds if

[TABLE]

where $\tilde{\Sigma}=(\Sigma_{1}^{-1/2}(\Sigma_{1}^{1/2}\Sigma_{2}\Sigma_{1}^{1/2})^{1/2}\Sigma_{1}^{-1/2})^{-1}$.

Unlike other dispersion measures between probability distributions such as the alpha divergence and the $L^{p}$ distance, the Wasserstein distance has the characteristic that it inherits the metric structure of the underlying space. By exploiting this property, the Wasserstein distance has been widely used in machine learning and computer vision (Peyré and Cuturi, 2018).

2.2 Shrinkage estimation of location parameter

We briefly review existing results on Bayes shrinkage estimation of location parameters. For more details, see Fourdrinier et al. (2018).

First, we consider the known scale cases. Suppose that we have an observation $x\sim{\rm N}_{d}(\mu,I_{d})$ and estimate $\mu$ under the quadratic loss $L(\mu,\hat{\mu})=\|\hat{\mu}-\mu\|^{2}$. For this problem, the usual estimator $\hat{\mu}(x)=x$ is minimax and the (generalized) Bayes estimator with respect to a prior $\pi(\mu)$ is given by the posterior mean:

[TABLE]

Many minimax generalized Bayes estimators have been developed. For example, Stein (1973) proved that the generalized Bayes estimator with respect to a superharmonic prior is minimax and dominates the usual estimator $\hat{\mu}(x)=x$. Analogous results for scale mixtures of normal distributions are given in Maruyama (2003) and Fourdrinier et al. (2008). Extensions to more general classes of spherically symmetric distributions are given in Fourdrinier et al. (2008).

Next, we consider the unknown scale cases. Suppose that we have independent observations $x\sim{\rm N}_{d}\left(\mu,\sigma^{2}I_{d}\right)$ and ${s}\sim\sigma^{2}\chi^{2}_{n}$. For the problem of estimating $\mu$ under the scaled quadratic loss $L(\mu,\hat{\mu};\sigma^{2})=\|\hat{\mu}-\mu\|^{2}/\sigma^{2}$, the usual estimator $\hat{\mu}(x,s)=x$ is minimax and the (generalized) Bayes estimator with respect to a prior $\pi(\mu,\sigma)$ is given by

[TABLE]

Maruyama and Strawderman (2005) proposed a class of generalized Bayes minimax estimators in this setting. By putting $a=-1$, $b=0$, $C=D=I_{d}$, $e=-1$ and $\gamma=1$ in their Theorem 2.3, the following is obtained.

Lemma 2.4** (Maruyama and Strawderman (2005)).**

Let $\tau=\sigma^{-2}$. If $(n-2)(d-4)\geq 8$, then the generalized Bayes estimator of $\mu$ under the scaled quadratic loss with respect to the hierarchical prior

[TABLE]

[TABLE]

is minimax and dominates the usual estimator $\hat{\mu}(x,s)=x$.

Later, we will consider estimation of $\mu$ under the unscaled quadratic loss $L(\mu,\hat{\mu})=\|\hat{\mu}-\mu\|^{2}$. For this loss function, the usual estimator $\hat{\mu}(x,s)=x$ is the best equivariant and the (generalized) Bayes estimator with respect to a prior $\pi(\mu,\sigma)$ is given by

[TABLE]

Then, Lemma 2.4 is rewritten as follows.

Lemma 2.5**.**

Let $\tau=\sigma^{-2}$. If $(n-2)(d-4)\geq 8$, then the generalized Bayes estimator of $\mu$ under the unscaled quadratic loss with respect to the hierarchical prior

[TABLE]

[TABLE]

dominates the best equivariant estimator $\hat{\mu}(x,s)=x$.

2.3 Improved estimation of scale parameter under quadratic loss

Suppose that we have independent observations $x\sim{\rm N}_{d}\left(\mu,\sigma^{2}I_{d}\right)$ and ${s}\sim\sigma^{2}\chi^{2}_{n}$ where both $\mu$ and $\sigma^{2}$ are unknown. Stein (1964) studied estimation of the variance $\sigma^{2}$ under the scaled quadratic loss $L(\sigma^{2},\hat{\sigma}^{2})=(\hat{\sigma}^{2}/\sigma^{2}-1)^{2}$ and showed that the estimator $\hat{\sigma}^{2}(x,s)=\min\{s/(n+2),(\|x\|^{2}+s)/(n+d+2)\}$ dominates the best equivariant estimator $\hat{\sigma}_{0}^{2}(x,s)=s/(n+2)$. Note that this improved estimator is not generalized Bayes. Later, Brewster and Zidek (1974), Strawderman (1974) and Maruyama and Strawderman (2006) developed generalized Bayes minimax estimators that dominate $\hat{\sigma}_{0}^{2}(x,s)$.

Here, we consider estimation of the standard deviation $\sigma$ under the unscaled quadratic loss $L(\sigma,\hat{\sigma})=(\hat{\sigma}-\sigma)^{2}$. The best equivariant estimator is given by

[TABLE]

Kubokawa (1994) proposed a unified approach to improving on the best equivariant estimators of powers of the scale parameter. Consider a class of estimators $\hat{\sigma}(x,s)=\phi(w)\sqrt{s}$ where $w=\|x\|^{2}/s$ and let

[TABLE]

Then, from Theorem 2.2 of Kubokawa (1994), we obtain the following.

Lemma 2.6**.**

Assume the following conditions:

• •

$\phi(w)$* is non-decreasing.*

• •

$\phi(w)\to c$* as $w\to\infty$, where $c$ is defined in (1).*

• •

$\phi(w)\geq\phi_{0}(w)$* for every $w\geq 0$.*

Then, the estimator $\hat{\sigma}(x,s)=\phi(w)\sqrt{s}$ dominates $\hat{\sigma}_{0}(x,s)$ in (1).

In fact, the lower bound $\phi_{0}(w)$ in Lemma 2.6 corresponds to a generalized Bayes estimator.

Proposition 2.1**.**

Let $\tau=\sigma^{-2}$. The estimator $\hat{\sigma}(x,s)=\phi_{0}(w)\sqrt{s}$ is the generalized Bayes estimator of $\sigma$ under the unscaled quadratic loss with respect to the following hierarchical prior:

[TABLE]

[TABLE]

Proof.

The joint distribution of $(\tau,x,s)$ is

[TABLE]

Therefore, the generalized Bayes estimator of $\sigma=\tau^{-1/2}$ under the unscaled quadratic loss is

[TABLE]

∎

3 Location family

Let

[TABLE]

be a location family with finite second moments, where $z\in\mathbb{R}^{d}$, $f(z)\geq 0$ and

[TABLE]

We consider prediction of $y\sim p(y\mid\mu)$ based on the observation $x\sim p(x\mid\mu)$ under the $L^{2}$ Wasserstein loss.

For this problem, the class of plug-in densities form a complete class (Berger, 1980) as follows. Note that it is sufficient to restrict our attention to predictive densities with finite mean, since the $L^{2}$ Wasserstein loss would diverge otherwise.

Proposition 3.1**.**

Any predictive density $\hat{p}(y\mid x)$ with mean $\hat{\mu}(x)$ is dominated by the plug-in density $p(y\mid\hat{\mu}(x))$.

Proof.

By Lemma 2.2, we obtain

[TABLE]

Therefore, the risk of the plug-in density $p(y\mid\hat{\mu}(x))$ is not larger than that of $\hat{p}(y\mid x)$ for every $\mu$. ∎

Since $W_{2}(p(y\mid\mu),p(y\mid\hat{\mu}(x)))^{2}=\|\mu-\hat{\mu}(x)\|^{2}$, the risk function of the plug-in density $p(y\mid\hat{\mu}(x))$ is equal to the quadratic risk function of the estimator $\hat{\mu}(x)$:

[TABLE]

Thus, predictive density estimation under the $L^{2}$ Wasserstein loss reduces to point estimation of $\mu$ under the quadratic loss. Therefore, the best equivariant predictive density and the Bayesian predictive density is obtained as follows.

Theorem 3.1**.**

The best equivariant predictive density under the $L^{2}$ Wasserstein loss is given by the plug-in density with the usual estimator $\hat{\mu}(x)=x$ and it is minimax.

Theorem 3.2**.**

The Bayesian predictive density with respect to a prior $\pi(\mu)$ under the $L^{2}$ Wasserstein loss is given by the plug-in density with the posterior mean $\hat{\mu}^{\pi}(x)$ of $\mu$:

[TABLE]

*where *

[TABLE]

Its risk function is equal to the quadratic risk function of the Bayes estimator $\hat{\mu}^{\pi}(x)$:

[TABLE]

Proof.

From Proposition 3.1, it is sufficient to restrict our attention to the class of plug-in densities. For this class, the risk function is obtained as

[TABLE]

which is equal to the quadratic risk of $\hat{\mu}$ as an estimator of $\mu$. Therefore, the average risk is minimized by setting $\hat{\mu}$ to the Bayes estimator $\hat{\mu}^{\pi}$ under the quadratic loss, which is given by the posterior mean (Lehmann and Casella, 1998). ∎

Corollary 3.1**.**

If the (generalized) Bayes estimator of $\mu$ with respect to a prior $\pi(\mu)$ is minimax under the quadratic loss, then the Bayesian predictive density with respect to the prior $\pi(\mu)$ is also minimax under the $L^{2}$ Wasserstein loss.

Therefore, Bayesian predictive densities with respect to shrinkage priors are minimax. For example, from Stein (1973), we obtain the following.

Corollary 3.2**.**

For the normal model with known variance, the Bayesian predictive density with respect to a superharmonic prior $\pi(\mu)$ is minimax under the $L^{2}$ Wasserstein loss.

Similarly, improved predictive densities are obtained for scale mixtures of normal distributions by using results in Maruyama (2003) and Fourdrinier et al. (2008) and also for spherically symmetric distributions by using results in Fourdrinier et al. (2008).

4 Location-scale family

Let

[TABLE]

be a location-scale family, where $z\in\mathbb{R}^{d}$, $f(z)\geq 0$ and

[TABLE]

We consider prediction of $y\sim p(y\mid\mu,\sigma)$ based on the independent observations $x_{1},\ldots,x_{n}\sim p(x\mid\mu,\sigma)$ under the $L^{2}$ Wasserstein loss.

For this problem, the class of plug-in densities form a complete class (Berger, 1980) as follows. Note that it is sufficient to restrict our attention to predictive densities with finite mean and covariance, since the $L^{2}$ Wasserstein loss would diverge otherwise.

Proposition 4.1**.**

Any predictive density $\hat{p}(y\mid x^{n})$ with mean $\hat{\mu}(x^{n})$ and covariance $\hat{\Sigma}(x^{n})$ is dominated by the plug-in density $p(y\mid\hat{\mu}(x^{n}),\hat{\sigma}(x^{n}))$ where $\hat{\sigma}(x^{n})=d^{-1}{\rm tr}(\hat{\Sigma}(x^{n})^{1/2})$.

Proof.

By Lemma 2.3, we obtain

[TABLE]

where we used

[TABLE]

in the second inequality. Therefore, the risk of the plug-in density $p(y\mid\hat{\mu}(x^{n}),\hat{\sigma}(x^{n}))$ is not larger than that of $\hat{p}(y\mid x^{n})$ for every $(\mu,\sigma)$. ∎

Since $W_{2}(p(y\mid\mu,\sigma),{p}(y\mid\hat{\mu}(x^{n}),\hat{\sigma}(x^{n})))^{2}=\|\hat{\mu}(x^{n})-\mu\|^{2}+d(\hat{\sigma}(x^{n})-\sigma)^{2}$, the risk function of the plug-in density $p(y\mid\hat{\mu}(x^{n}),\hat{\sigma}(x^{n}))$ is equal to the sum of the (unscaled) quadratic risk functions of the estimators $\hat{\mu}(x^{n})$ and $\hat{\sigma}(x^{n})$:

[TABLE]

Thus, predictive density estimation under the $L^{2}$ Wasserstein loss reduces to joint estimation of $\mu$ and $\sigma$ under unscaled quadratic losses. Therefore, the best equivariant predictive density and the Bayesian predictive density is obtained as follows.

Theorem 4.1**.**

The best equivariant predictive density under the $L^{2}$ Wasserstein loss is given by the plug-in density with the best equivariant estimators of $\mu$ and $\sigma$ under the unscaled quadratic losses.

Theorem 4.2**.**

The Bayesian predictive density with respect to a prior $\pi(\mu,\sigma)$ under the $L^{2}$ Wasserstein loss is given by the plug-in density with the posterior mean $(\hat{\mu}^{\pi}(x^{n}),\hat{\sigma}^{\pi}(x^{n}))$ of $(\mu,\sigma)$:

[TABLE]

*where *

[TABLE]

[TABLE]

Its risk function is equal to the sum of the (unscaled) quadratic risk functions of $\hat{\mu}^{\pi}(x^{n})$ and $\hat{\sigma}^{\pi}(x^{n})$:

[TABLE]

Proof.

From Proposition 4.1, it is sufficient to restrict our attention to the class of plug-in densities $\hat{p}(y\mid x^{n})=p(y\mid\hat{\mu}(x^{n}),\hat{\sigma}(x^{n}))$. For this class, the risk function is obtained as

[TABLE]

which is equal to the sum of the quadratic risks of $\hat{\mu}$ and $\hat{\sigma}$ as estimators of $\mu$ and $\sigma$, respectively. Therefore, the average risk is minimized by setting $\hat{\mu}$ and $\hat{\sigma}$ to the Bayes estimators $\hat{\mu}^{\pi}$ and $\hat{\sigma}^{\pi}$ under the unscaled quadratic loss, which are given by the posterior mean (Lehmann and Casella, 1998). ∎

Corollary 4.1**.**

If both of the Bayes estimators of $\mu$ and $\sigma$ with respect to a prior $\pi(\mu,\sigma)$ dominate the best equivariant estimator under the unscaled quadratic losses, then the Bayesian predictive density with respect to the prior $\pi(\mu,\sigma)$ dominates the best equivariant predictive density under the $L^{2}$ Wasserstein loss.

Finally, we provide a concrete example of the Bayesian predictive density that dominates the best equivariant one in the normal model with unknown variance. Specifically, let $x_{1},\ldots,x_{n}\sim{\rm N}_{d}(\mu,\sigma^{2}I_{d})$ independently and consider prediction of $y\sim{\rm N}_{d}(\mu,\sigma^{2}I_{d})$. The sufficient statistics are

[TABLE]

and

[TABLE]

From Theorem 4.1, the best equivariant predictive density is $\hat{p}(y\mid x^{n})={\rm N}_{d}(\bar{x},c^{2}{s}I_{d})$ where $c$ is defined in (1). By combining Lemma 2.5 and Proposition 2.1 with Corollary 4.1, an improved Bayesian predictive density is obtained as follows.

Theorem 4.3**.**

If $(n-3)(p-4)\geq 8$, then the Bayes solution with respect to the following hierarchical prior dominates the best equivariant predictive density under the $L^{2}$ Wasserstein loss:

[TABLE]

[TABLE]

5 Conclusion

In this study, we investigated predictive density estimation under the $L^{2}$ Wasserstein loss. For both families, plug-in densities form a complete class (Berger, 1980). For location families, the Bayesian predictive density is given by the plug-in density with the posterior mean of the location parameter. For location-scale families, the Bayesian predictive density is given by the plug-in density with the posterior mean of the location and scale parameters. We provided Bayesian predictive densities that dominate the best equivariant one in normal models.

We focused on location families and location-scale families in this study. Extension to the case of an unknown covariance matrix is an interesting future problem. Namely, we consider the model

[TABLE]

which includes the elliptically contoured distributions (Fang and Anderson, 1990) as special cases. In this case, the problem reduces to point estimation of the covairance matrix $\Sigma$ under the loss $L(\Sigma,\hat{\Sigma})=W_{2}(p(x\mid 0,\Sigma),p(x\mid 0,\hat{\Sigma}))^{2}={\rm tr}\left(\Sigma+\hat{\Sigma}-2(\Sigma^{1/2}\hat{\Sigma}\Sigma^{1/2})^{1/2}\right)$, which is different from Stein’s loss or Frobenius loss. In addition, while we considered $L^{2}$ Wasserstein loss in this study, it is an interesting future work to study predictive density estimation under other Wasserstein losses.

Acknowledgements

William Strawderman’s research is partially supported by a grant from the Simons Foundation (#418098).