Bayes Extended Estimators for Curved Exponential Families

TL;DR

This paper introduces a simple approximation for the Bayesian predictive density in curved exponential families using projection onto finite-dimensional exponential families, with theoretical analysis and practical benefits.

Contribution

It proposes the Bayes extended estimator as a practical approximation of the Bayesian predictive density for curved exponential families, with proven asymptotic properties.

Findings

The projection asymptotically matches the plugin density with the posterior mean.

Numerical simulations show the approximation approaches the true density as dimension increases.

The method offers a good trade-off between risk performance and computational efficiency.

Abstract

The Bayesian predictive density has complex representation and does not belong to any finite-dimensional statistical model except for in limited situations. In this paper, we introduce its simple approximate representation employing its projection onto a finite-dimensional exponential family. Its theoretical properties are established parallelly to those of the Bayesian predictive density when the model belongs to curved exponential families. It is also demonstrated that the projection asymptotically coincides with the plugin density with the posterior mean of the expectation parameter of the exponential family, which we refer to as the Bayes extended estimator. Information-geometric correspondence indicates that the Bayesian predictive density can be represented as the posterior mean of the infinite-dimensional exponential family. The Kullback--Leibler risk performance of the approximation is demonstrated by numerical simulations and it indicates that the posterior mean of the expectation parameter approaches the Bayesian predictive density as the dimension of the exponential family increases. It also suggests that approximation by projection onto an exponential family of reasonable size is practically advantageous with respect to risk performance and computational cost.

Tables

Table 1. TABLE I : All probability densities and a full exponential family

$ℱ$: all probability densities	$ℰ$: a full exponential family
$p\left(x\right)$: m-representation	$\eta$: m-affine parameter
$\log p\left(x\right)$: e-representation	$\theta$: e-affine parameter
infinite dimensional model	finite dimensional model

Table 2. TABLE II : Computational time required to evaluate predictive density for 1000 1000 1000 new samples and memory size of each predictive density ( l = 80 , λ = 1 formulae-sequence 𝑙 80 𝜆 1 l=80,\lambda=1 ).

	Bayes	E-plugin
computational time	$1.24\times {10}^2$ s	$1.6\times {10}^{-2}$ s
predictive density size	$1.29\times {10}^6$ bytes	$5.14\times {10}^4$ bytes

Equations

$$\mathcal{P}=\left\{p(x;\omega )\mid \omega \in \Omega \right\}.$$

$$D\{p(y;\omega );\hat{p}(y;\omega )\}=\int p(y;\omega )\log \frac{p(y;\omega )}{\hat{p}(y;\omega )}dy$$

$$E[D\{p(y;\omega );\hat{p}(y;\omega )\}]$$

$$=\int p(x^n;\omega )D\{p(y;\omega );\hat{p}(y;\omega )\}d{x}^n,$$

$$\int \pi (\omega )\int p(x^n;\omega )D\{p(y;\omega );\hat{p}(y;\omega )\}d{x}^{n}d\omega ,$$

$${\hat{p}}_{\pi }(y\mid x^n)=\int p(y;\omega )p_{\pi }(\omega \mid x^n)d\omega$$

$$p_{\pi }(\omega \mid x^n)=\frac{p(x^n;\omega )\pi (\omega )}{\int p(x^n;\omega )\pi (\omega )d\omega }$$

$$\mathcal{P}=$$

$$\omega =({\omega }^a)\in \Omega ,a=1,\dots ,d,i=1,\dots ,m\},$$

$$\mathcal{E}=$$

$$i=1,\dots ,m\}(\Theta \subseteq {\mathbb{R}}^m)$$

$${\eta }_i=E(x_i)=\frac{\partial }{\partial {\theta }^i}\Psi (\theta )(i=1,\dots ,m)$$

$$⟨{\partial }_{a}p(x;\omega ),{\partial }_{b}p(x;\omega )⟩=\int \frac{{\partial }_{a}p(x;\omega ){\partial }_{b}p(x;\omega )}{p(x;\omega )}dx.$$

$$g_{ab}(\omega )=⟨{\partial }_{a}p(x;\omega ),{\partial }_{b}p(x;\omega )⟩.$$

$$\text{raise}1.5pt{{\Gamma }^{\mathrm{e}}}_{abc}^{}(\omega )$$

$$\text{raise}1.5pt{{\Gamma }^{\mathrm{m}}}_{abc}^{}(\omega )$$

$$\text{raise}1.5pt{{\Gamma }^{\mathrm{e}}}_{ab}^c=\text{raise}1.5pt{{\Gamma }^{\mathrm{e}}}_{abd}^{}{g}^{dc},\text{raise}1.5pt{{\Gamma }^{\mathrm{m}}}_{ab}^c=\text{raise}1.5pt{{\Gamma }^{\mathrm{m}}}_{abd}^{}{g}^{dc},$$

$$T_{abc}=$$

$$=$$

$$\times \{{\partial }_c\log p(x;\omega )\}dx.$$

$${\pi }_J(\omega )=\sqrt{|g(\omega )|},$$

$$⟨\frac{\partial }{\partial {\theta }^i}p(x;\theta ),\frac{\partial }{\partial {\eta }_j}p(x;\theta )⟩={\delta }_j^i$$

$$\text{raise}1.5pt{{\Gamma }^{\mathrm{e}}}_{abc}^{}=({\partial }_a{\partial }_b{\theta }^i)({\partial }_c{\eta }_i)\text{mbox}and\text{raise}1.5pt{{\Gamma }^{\mathrm{m}}}_{abc}^{}=({\partial }_a{\partial }_b{\eta }_i)({\partial }_c{\theta }^i),$$

$$\int {\mathrm{o}}_{\mathrm{p}}(1)\mathrm{dy}=\mathrm{o}(1),\int {\mathrm{O}}_{\mathrm{p}}({\mathrm{n}}^{-1})\mathrm{dy}=\mathrm{O}({\mathrm{n}}^{-1}),$$

$${\bar{\eta }}_{\pi }=\frac{\int \eta (\omega )p(x^n;\eta (\omega ))\pi (\omega )d\omega }{\int p(x^n;\eta (\omega ))\pi (\omega )d\omega }.$$

$${\bar{\omega }}_{\pi }=\frac{\int \omega p(x^n;\omega )\pi (\omega )d\omega }{\int p(x^n;\omega )\pi (\omega )d\omega }.$$

$$D\{p(y;\theta (\omega ));p(y;\hat{\theta })\}$$

$$=\int p(y;\theta )\log \left(\frac{\exp ({\theta }^i{y}_i-\Psi (\theta ))}{\exp ({\hat{\theta }}^i{y}_i-\Psi (\hat{\theta }))}\right)dy$$

$$=({\theta }^i-{\hat{\theta }}^i){\eta }_i-(\Psi (\theta )-\Psi (\hat{\theta })).$$

$$\int p_{\pi }(\omega \mid x^n)D\{p(y;\theta (\eta ));p(y;\hat{\theta }(\eta ))\}d\omega$$

$$=(\overline{{\theta }^i{\eta }_i}-{\hat{\theta }}^i\overline{{\eta }_i})-(\overline{\Psi (\theta )}-\Psi (\hat{\theta }))$$

$$=({\theta }^i(\overline{\eta })-{\hat{\theta }}^i){\overline{\eta }}^i-(\Psi (\theta (\overline{\eta }))-\Psi (\hat{\theta }))$$

$$+\left(-{\theta }^i(\overline{\eta })\overline{{\eta }_i}+\overline{{\theta }^i{\eta }_i}+\Psi (\theta (\overline{\eta }))-\overline{\Psi (\theta )}\right)$$

$$=D\{p(y;\overline{\eta });p(y;\hat{\theta })\}$$

$$+\left(-{\theta }^i(\overline{\eta }){\overline{\eta }}^i+\overline{{\theta }^i{\eta }_i}+\Psi (\theta (\overline{\eta }))-\overline{\Psi (\theta )}\right),$$

$$\overline{f(\eta )}=\int p_{\pi }(\omega \mid x^n)f(\eta )d\omega .$$

Full text

I Introduction

Constructing predictive densities is a fundamental problem in statistical analysis that aims at predicting the behavior of future samples using past observations. Let us suppose that we have observations $x^{n}=\{x(1),x(2),\dots,x(n)\}$ that are independently distributed according to a probability distribution with density function $p(x;\omega)$ that belongs to a statistical model

[TABLE]

The objective is to provide the predictive density of $y=x(n+1)$ that is independently distributed according to the same density $p(y;\omega)$. We adopt the Kullback–Leibler divergence

[TABLE]

as a loss function of a predictive density $\hat{p}(y;{\omega})$. Then, the risk function and the Bayes risk with respect to a prior $\pi(\omega)$ can be written as

[TABLE]

and

[TABLE]

respectively.

Except for in limited cases, the Bayesian predictive distribution does not belong to any finite-dimensional model, which makes it intractable to obtain the full density although it is optimal with respect to the Bayes risk. The Bayesian predictive density is defined by

[TABLE]

where $p_{\pi}(\omega\mid x^{n})$ is the posterior density

[TABLE]

of $\omega$. It is shown in [1] that the Bayesian predictive density is optimal with respect to the Bayes risk in terms of the Kullback–Leibler divergence in the family of all probability densities, which we denote as $\mathcal{F}$. However, the full Bayesian predictive density is intractable in most problems due to the complex representation that involves averaging plugin densities about the model parameters. It is not included in the model $\mathcal{P}$ or even in any finite-dimensional model in most problems, while plugin densities are always included in $\mathcal{P}$ as they are constructed by plugging-in an estimator $\hat{\omega}(x^{n})$ to the model.

In the present paper, we represent the Bayesian predictive density as the infinite-dimensional limit of a parameterized distribution of an exponential family. We demonstrate that the Bayesian predictive density can be considered as an infinite-dimensional extension of the plugin density with the posterior mean of the expectation parameter of an exponential family. It is shown that theoretical properties including optimality with respect to the Bayes risk are appropriately retained by this extension. The plugin density with the posterior mean of the expectation parameter coincides with the projection of the Bayesian predictive density onto the exponential family with respect to the Fisher metric. It is shown that it approaches the Bayesian predictive density closer with respect to the risk as the projected exponential family increases. There is also information-geometrical correspondence between the Bayesian predictive density and the projected density that comes from the correspondence between $\mathcal{F}$ and the exponential family. In practice, the Bayesian predictive density can be computationally approximated, for example, by taking the mean of plugin densities using Markov chain Monte Carlo simulations, or by performing an approximation of the posterior density using methods like the Laplace method, although the objective of this research is not to develop a complex approximation based on computational methods. Apart from computational approximations, a class of empirical Bayes predictive densities is proposed for multivariate normal models in [2] to avoid the intractable implementations of Bayesian predictive densities. Rather than constructing an approximation of the Bayesian predictive density that has good performance in terms of risk, we aim to formulate a simple interpretation of the Bayesian predictive density that maintains its theoretical properties on a finite-dimensional model.

The outline of the construction of the approximate predictive densities is explained below. We consider a model of a subspace of an exponential family as $\mathcal{P}$, namely we consider a statistical model of a curved exponential family

[TABLE]

where $\Omega\subset\mathbb{R}^{d}$ and $1\leq d\leq m$. The model $\mathcal{P}$ parametrized by $\omega$ is embedded in an exponential family parametrized by $\theta$, thus here we represent $\theta$ as $\theta(\omega)$. Summation over a repeated index is automatically taken according to Einstein’s summation convention: if an index occurs as an upper and lower index in one term, then the summation is implied. Curved exponential families embedded in exponential families can express a variety of models including network models (e.g., [3]) and time series models (e.g., [4]). They can also be applied to stochastic processes [5]. We consider predictive densities in a finite-dimensional full exponential family

[TABLE]

that includes the original curved exponential model $\mathcal{P}$. We refer to plugin densities in $\mathcal{E}$ as extended plugin densities. The inclusion relation is $\mathcal{P}\subseteq\mathcal{E}\subseteq\mathcal{F}$ and we consider the middle layer of the three-layer structure. The coordinate system $\theta=(\theta^{i})\ (i=1,\dots,m)$ is called the natural parameter of exponential families. Another coordinate system $\eta=(\eta_{i})$ defined by

[TABLE]

is called the expectation parameter. The posterior mean of $\eta$ is closely related to the Bayesian predictive density, and the extended plugin density with the posterior mean of $\eta$ is considered in this paper. Based on the idea of covering $\mathcal{F}$ by extending exponential families, we specify the models that can be embedded in exponential families, and the theoretical properties described in the following sections are based on this embedding. The policy of expressing a probability density by extending exponential families has been investigated, for example, in [6], in which a log-density is approximated using series of polynomials and the rate of convergences is obtained. It should be noted that the practical advantage illustrated in the numerical experiments in Section IV can be attained for other models if it is possible to find an appropriate exponential family onto which we project the Bayesian predictive density.

From the viewpoint of information geometry, the posterior mean of $\eta$ can be considered as the correspondence to the Bayesian predictive density in $\mathcal{E}$. Table I represents the infinite-finite correspondence of m and e representations. Here, “m” and “e” are short notations denoting “mixture” and “exponential,” respectively. Exponential families and mixture families are important dual families in information geometry, and their typical representations are denoted as m-representation and e-representation, respectively. Concerning exponential families, the e-representation is $\log p(x)=\theta^{i}x_{i}-\Psi(\theta)+\log s(x),$ and $\theta$ is called the e-affine parameter, as the basis vector fields $(\partial/\partial\theta_{i})p(x)~{}(i=1,\dots,m)$ are parallel vector fields with respect to the e-connection as defined in Section II of this paper. Concerning mixture families, the m-representation can be written as $p(x)=\eta_{i}q^{i}(x)+c(x)$, and $\eta$ is the affine coordinate system about the m-connection (also defined in Section II). We can also set the m-affine coordinate system in exponential families, and it is defined by (1). Here, the m-affine (or e-affine) parameters are the finite-dimensional typical representations of m-representation (or e-representation, respectively). The Bayesian predictive density is the posterior mean about the m-representation (that is, density functions), and its finite-dimensional correspondence is the posterior mean about $(\eta_{i})$. Since the Bayesian predictive density is optimal in the infinite-dimensional exponential family, we might expect that the posterior mean of $(\eta_{i})$ exhibit the same properties as the Bayesian predictive density, such as optimality with respect to the Bayes risk.

The properties of the posterior mean of ${\eta}$ are investigated in following sections as follows. In Section III-A, we show that the extended plugin density with the posterior mean of $\eta$ is optimal with respect to the Bayes risk in the finite-dimensional exponential family. We denote the posterior mean of $\eta$ as the Bayes extended estimator. In Section III-B, the extended plugin density with the Bayes extended estimator is proved to be the projection of the Bayesian predictive density onto $\mathcal{E}$ in terms of the Fisher metric. In Section III-C, its optimality with respect to the risk along orthogonal shift from the model is shown to be common with the property of optimality in the case of the Bayesian predictive density. The relation between the projection angle under the Fisher metric and the risk difference between the Bayesian predictive density and the extended plugin density with the Bayes extended estimator is investigated. In Section IV, we compare the risk performance of the Bayes estimator, the Bayes extended estimator, and the Bayesian predictive density by conducting numerical simulations on the Gaussian spiked covariance models. We confirm that the projection angle converges to zero as the dimension of $\mathcal{E}$ increases. The simulation results also suggest that the projection of the Bayesian predictive density is practically effective in approximating it with respect to the Kullback–Leibler risk and the computational cost.

II Preliminaries

In this section, we prepare some information-geometric notions. For details of the notions and notation concerning the differential geometry of curved exponential families, refer to [7].

Let $a,b,\dots$ be indices for $\omega$. Let $T_{\omega}\mathcal{P}$ be the tangent space of $\mathcal{P}$ at a point $\omega$. The tangent space $T_{\omega}\mathcal{P}$ is identified with the vector space spanned by $\partial_{a}p(x;\omega)~{}(a=1,\dots,d)$, where $\partial_{a}$ denotes $\partial/\partial\omega^{a}$. We define inner products in the tangent space by

[TABLE]

In a statistical model $\mathcal{P},$ each component of the Fisher information matrix is defined by

[TABLE]

Let $g^{ab}$ be a component of the inverse matrix of $(g_{ab})$. Then, e-connection and the m-connection coefficients are defined as:

[TABLE]

and

[TABLE]

respectively. We define

[TABLE]

and

[TABLE]

The Jeffreys prior density is given by

[TABLE]

where $|g(\omega)|$ is the determinant of the matrix $(g_{ab}(\omega))$.

The coordinate systems $(\theta^{i})$ and $(\eta_{i})$ of the exponential family $\mathcal{E}$ are dual to each other in the sense that

[TABLE]

where $\delta^{i}_{j}$ is the Kronecker delta. In a curved exponential family, e-connection and m-connection coefficients are expressed as:

[TABLE]

respectively.

In the rest of the paper, we assume regularity conditions to ensure that equalities such as

[TABLE]

hold. For the details of the regularity conditions, see [8].

III Main results

III-A Optimality with respect to Bayes risk

The posterior mean of $\eta$ (we denote as $\bar{\eta}_{\pi}$) of $\mathcal{E}$ is evaluated as follows:

[TABLE]

Note that $\bar{\eta}_{\pi}(x^{n})\neq\eta(\bar{\omega}_{\pi}(x^{n}))$ in general where $\bar{\omega}_{\pi}(x^{n})$ is the posterior mean of $\omega$:

[TABLE]

We demonstrate that $p(y;\bar{\eta}_{\pi})$ is optimal in $\mathcal{E}$ with respect to the Bayes risk based on a prior $\pi$.

Proposition III.1**.**

The Bayes risk of $p(y;\hat{\eta})$, where $\hat{\eta}$ is an estimator of $\eta$, is minimized when $\hat{\eta}=\bar{\eta}_{\pi}.$

Proof.

Let $\hat{\theta}$ be an estimator of $\theta$. Note that $\theta$ and $\eta$ are functions of $\omega$. The Kullback–Leibler loss of ${p}(y;\hat{\theta})$ is

[TABLE]

Hence

[TABLE]

where, for a function $f(\eta)$,

[TABLE]

It is minimized when $\hat{\theta}=\theta(\overline{\eta})=\theta(\overline{\eta}_{\pi})$.

By multiplying (III-A) with $\int p(x^{n};\omega)\pi(\omega)d\omega$ and then integrating with respect to $x^{n}$, it is shown that $p(y;\bar{\eta}_{\pi})$ is optimal with respect to the Bayes risk in $\mathcal{E}$. ∎

We refer to $\bar{\eta}_{\pi}$ as the Bayes extended estimator. Hereinafter, we denote the Bayes extended estimator and the Bayes estimator of $\omega$ about a prior $\pi$ as $\hat{\eta}_{\pi}(=\bar{\eta}_{\pi})$ and $\hat{\omega}_{\pi}$, respectively.

From Proposition III.1, the extended plugin density with $\hat{\eta}_{\pi}$ is the projection of the Bayesian predictive density onto $\mathcal{E}$ about the Bayes risk. It is nearest to the Bayesian predictive density in $\mathcal{E}$ regarding the Bayes risk, because the Bayesian predictive density is optimal about the Bayes risk in $\mathcal{F}$. In fact, $p(y,\hat{\eta}_{\pi})$ coincides with the projection of the Bayesian predictive density onto $\mathcal{E}$ regarding the Fisher metric asymptotically, as shown in Section III-B.

The choice of $\mathcal{E}$ does not require to be fixed, and we can consider situations in which the size of the extended model $\mathcal{E}$ can be increased, for example, by employing sequences of exponential families as in [6] and [9]. In those situations, the extended plugin density with the Bayes extended estimator $\hat{\eta}_{\pi}$ approaches the Bayesian predictive density as $\mathcal{E}$ grows, as $\mathcal{E}$ approaches the set of all probability distributions $\mathcal{F}$.

Here, we use a simple example to illustrate the difference of the plugin density with $\hat{\omega}_{\pi}$, the extended plugin density with $\hat{\eta}_{\pi}$, and the Bayesian predictive density.

Example  (Fisher circle model) We consider two dimensional Gaussian distribution $\mathrm{N}(\mu,I_{2})$ with unknown mean vector $\mu$ and the identity covariance matrix $I_{2}$. The density function is

[TABLE]

When the mean vector $\mu$ is expressed as

[TABLE]

the one-dimensional submodel is called the Fisher circle model. Here, the following holds:

[TABLE]

Then, we derive the Bayes estimator $\hat{\omega}_{\pi}$, the Bayes extended estimator $\hat{\eta}_{\pi}$, and the Bayesian predictive density. For $x^{n}=\{x(1),x(2),\dots,x(n)\}$,

[TABLE]

where $\bar{x}=\sum^{n}_{t=1}x(t)/n$. Let $\bar{x}=(\|\bar{x}\|\cos\phi,\|\bar{x}\|\sin\phi)^{\top}$. Then by the law of cosines,

[TABLE]

and

[TABLE]

where $I_{0}(\cdot)$ is the modified Bessel function of the first kind. See [10] (pp. 138–140) for the details. When the uniform prior

[TABLE]

is adopted, the posterior density is

[TABLE]

It follows that the plugin of the Bayes estimator is $\mathrm{N}\left(\frac{\bar{x}}{\|\bar{x}\|},I_{2}\right)$:

[TABLE]

The extended plugin with the Bayes extended estimator is $\mathrm{N}\left(\frac{I_{1}(n\|\bar{x}\|)}{I_{0}(n\|\bar{x}\|)}\frac{\bar{x}}{\|\bar{x}\|},I_{2}\right)$:

[TABLE]

where $\hat{\eta}_{\pi}$ is not included in the circle parametrized by $\omega$ and $\hat{\eta}_{\pi}\neq\eta(\hat{\omega}_{\pi})$. Here $I_{1}(\cdot)$ is the modified Bessel function of the first kind. On the other hand, the Bayesian predictive density is given by

[TABLE]

Therefore, $p_{\pi}(y\mid x^{n})$ is not included in $\mathcal{P}$ or $\mathcal{E}$ because it is not a two-dimensional Gaussian with a covariance matrix $I_{2}$.

III-B Projection of Bayesian predictive densities in terms of Fisher metric

Here, we demonstrate that $p(y;\hat{\eta}_{\pi})$ is the projection of the Bayesian predictive density regarding the Fisher metric. It is shown via asymptotic expansion of $p(y;\hat{\eta}_{\pi})$, which is represented as a point in $\mathcal{E}$ that is parallelly and orthogonally shifted from the plugin density of the maximum likelihood estimator of $\mathcal{P}$ as shown in Figure 1.

Here, we proceed to obtain the asymptotic expansion of $\hat{\eta}_{\pi}$ around $\eta(\hat{\omega}_{\rm MLE})$.

Theorem III.1**.**

The Bayes extended estimator based on a prior $\pi(\omega)$ is expanded as

[TABLE]

where ${\pi}_{J}$ is the density of the Jeffreys prior and $T_{a}=T_{abc}g^{bc}$.

Proof.

See Appendix A. ∎

We can obtain the asymptotic expansion of the extended plugin density with $\hat{\eta}_{\pi}$.

Theorem III.2**.**

The extended plugin density with $\hat{\eta}_{\pi}$ is expanded as

[TABLE]

where $\partial^{i}={\partial}/{\partial\eta_{i}}$.

Proof.

Symbols such as $\eta(\hat{\omega}_{\rm MLE}),\partial_{a}\eta(\hat{\omega}_{\rm MLE})$, and $\partial_{a}\partial_{b}\eta(\hat{\omega}_{\rm MLE})$ are abbreviated to $\hat{\eta},\partial_{a}\hat{\eta}$, and $\partial_{ab}\hat{\eta}$, respectively. Considering the asymptotic expansion introduced in Theorem III.1, we obtain the following:

[TABLE]

∎

The shift from $p(y;\eta(\hat{\omega}_{\rm MLE}))$ to $p(y;\hat{\eta}_{\pi})$ in Theorem III.2 is composed of two components, one “parallel” and the other “orthogonal” to the model $\mathcal{P}$. That is, the term

[TABLE]

is included in the tangent space spanned by $\partial_{a}p(y;\eta)~{}(a=1,\dots,d)$ and the term

[TABLE]

is orthogonal to $\partial_{a}p(x;\eta)~{}(a=1,\dots,d)$ with respect to the inner product (2), because

[TABLE]

by using (3) and (4).

We can compare the orthogonal shifts (III-B) to the orthogonal shifts from $p(y;\eta(\hat{\omega}_{\rm MLE}))$ to the Bayesian predictive density, and we show (III-B) is the projection of the orthogonal shifts to the Bayesian predictive density onto $\mathcal{E}$. In [11], the Bayesian predictive density $p_{\pi}(y\mid x)$ is asymptotically expanded as

[TABLE]

The parallel shift is identical to that of $p(y;\hat{\eta}_{\pi})$. On the other hand, the orthogonal shift

[TABLE]

is different from that of $p(y;\hat{\eta}_{\pi})$ and it is not included in the tangent space of $\mathcal{E}$. Therefore, the shifted density $p_{\pi}(y\mid x)$ is not included in $\mathcal{E}$, while $p(y;\hat{\eta}_{\pi})$ is in $\mathcal{E}$.

To cope with the shifts that are orthogonal to $\mathcal{P}$, we introduce a coordinate system to the subspace of $\mathcal{E}$ that is orthogonal to $\mathcal{P}$. We divide the tangent vectors of $\mathcal{E}$ at $\eta$ into two parts, namely, into those parallel to $\mathcal{P}$ and those orthogonal to $\mathcal{P}$. For each point $\eta\in\mathcal{E}$, the tangent space $T_{\eta}\mathcal{E}$ is identified with the vector space spanned by

[TABLE]

The tangent space $T_{\omega}\mathcal{P}$ is a subspace of $T_{\eta}\mathcal{E}$. Let $A(\omega)$ be an $(m-d)$-dimensional smooth submanifold of $\mathcal{E}$ attached to each point $\omega\in\mathcal{P}$ and assume that $A(\omega)$ orthogonally transverses $\mathcal{P}$ at $\eta(\omega)$. Such a family of submanifolds $\{A(\omega)\mid\omega\in\Omega\}$ is called an ancillary family. We introduce an adequate coordinate system $\xi=(\xi^{\kappa})~{}(\kappa=d+1,\dots,m)$ to $A(\omega)$ so that a pair $(\omega,\xi)$ uniquely specifies a point of $\mathcal{E}$ in the neighborhood of $\eta(\omega)$. We adopt a coordinate system $\xi$ on $A(\omega)$ so that $\eta(\omega,\xi)\in\mathcal{P}$ if $\xi=0$. Then, we have

[TABLE]

where $\partial_{\kappa}p(x;\eta)=({\partial}/{\partial\xi^{\kappa}})p(x;\eta)$. Since $A(\omega)$ orthogonally transverses $\mathcal{E}$, we have

[TABLE]

Now we show that the extended plugin density with $\hat{\eta}_{\pi}$ is the projection of the Bayesian predictive density onto $\mathcal{E}$ asymptotically, as shown in Figure 1. The projection of (7) onto the tangent space of $\mathcal{E}$ is

[TABLE]

Because $\langle\partial_{c}p(x;\eta),\partial_{\lambda}p(x;\eta)\rangle=0~{}(c=1,\dots,d)$, (III-B) is

[TABLE]

where

[TABLE]

$(a,b=1,\dots,d,~{}\kappa=d+1,\dots,m)$ is the mixture embedding curvature of $\mathcal{P}$ in $\mathcal{E}$ and \raise 1.5pt\hbox{\displaystyle\mathop{H}^{\rm m}}_{ab}^{\,~{}\kappa}=\raise 1.5pt\hbox{\displaystyle\mathop{H}^{\rm m}}_{ab\lambda}^{\,}g^{\kappa\lambda}.

Here we represent \raise 1.5pt\hbox{\displaystyle\mathop{\Gamma}^{\rm m}}_{ab\kappa}^{\,} as \raise 1.5pt\hbox{\displaystyle\mathop{H}^{\rm m}}_{ab\kappa}^{\,} to clarify its geometrical interpretation as the mixture embedding curvature, distinguishing from \raise 1.5pt\hbox{\displaystyle\mathop{\Gamma}^{\rm m}}_{abc}^{\,},\raise 1.5pt\hbox{\displaystyle\mathop{\Gamma}^{\rm m}}_{a\kappa b}^{\,} and \raise 1.5pt\hbox{\displaystyle\mathop{\Gamma}^{\rm m}}_{\kappa ab}^{\,}.

We confirm that it coincides with the orthogonal shift to $p(y;\hat{\eta}_{\pi})$ asymptotically as follows. Let

[TABLE]

then the orthogonal component in Theorem III.2 is

[TABLE]

Since $h_{ab}~{}(a,b=1,\dots,d)$ are included in the space spanned by $\partial_{\kappa}p(x;\omega)~{}(\kappa=d+1,\dots,m)$, the following holds:

[TABLE]

As $\langle\partial_{c}p(x;\eta),\partial_{\lambda}p(x;\eta)\rangle=0$

$~{}{(c=1,\dots,d,~{}\lambda=d+1,\dots,m)}$,

[TABLE]

III-C Projection angle and risk difference

In this section, we demonstrate that the extended plugin density with $\hat{\eta}_{\pi}$ is optimal with respect to the risk along the orthogonal shift from the model $\mathcal{P}$. This property is parallel to those of Bayesian predictive densities investigated in [11], as explained later in this section. Orthogonal shifts from $\mathcal{P}$ can asymptotically improve the Kullback–Leibler risk from plugin densities of $\mathcal{P}$. By evaluating those risk improvements, the risk difference between $p(y;\hat{\eta}_{\pi})$ and Bayesian predictive densities can be corresponded to the projection angle (as represented in Figure 2).

In the following discussion, we consider extended plugin densities $p(y;\hat{\eta})$ with estimators $\hat{\eta}=\eta(\hat{\omega},\hat{\xi})$ where $\hat{\omega}$ and $\hat{\xi}$ can be expressed in the following forms, respectively:

[TABLE]

Here, $\alpha^{a}(\omega),\beta^{\kappa}(\omega)$ are smooth functions of $\mathrm{O}_{p}(1)$. The density can be expanded as

[TABLE]

The following results hold also for other asymptotically efficient estimators of $\omega$ other than $\hat{\omega}_{\rm MLE}$, although here we consider only $\hat{\omega}_{\rm MLE}$ for simplicity. This class of extended plugin densities include $p(x;\hat{\omega}_{\rm MLE})$ and $p(x;\hat{\eta}_{\pi})$. For $\hat{\omega}=\hat{\omega}_{\rm MLE}$ and $\hat{\xi}=0$, the density is the plugin density with the maximum likelihood estimator $\hat{\omega}_{\rm MLE}$. The extended plugin density $p(y;\hat{\eta}_{\pi})$ in Theorem III.2 is given by, using (9) and (III-B),

[TABLE]

We derive the Kullback–Leibler risk of extended plugin densities in this class.

Proposition III.2**.**

The Kullback–Leibler risk of an extended plugin density $p(y;\hat{\eta})$ where $\hat{\eta}=(\hat{\omega},\hat{\xi})$ is expressed as (III-C) is asymptotically expanded as follows:

[TABLE]

where \raise 1.5pt\hbox{\displaystyle\mathop{\nabla}^{\rm e}}_{a}^{\,}\alpha^{b}=\partial_{a}\alpha^{b}+\raise 1.5pt\hbox{\displaystyle\mathop{\Gamma}^{\rm e}}_{ac}^{\,~{}b}\alpha^{c}.

Proof.

See Appendix B. ∎

The risk expansion (III.2) confirms that the risk can be improved from a predictive density in $\mathcal{P}$ by selecting an appropriate orthogonal shift $\beta$. We obtain the optimal orthogonal shift.

Theorem III.3**.**

The optimal $\beta^{\kappa}$ in (III-C) is given by

[TABLE]

Proof.

The risk in Proposition III.2 is

[TABLE]

Therefore, $\beta$ is optimal when

[TABLE]

∎

Accordingly, the orthogonal component of the shift in Theorem III.2 is optimal, and the extended plugin density with $\hat{\eta}_{\pi}$ has the optimal shift (as illustrated in Figure 3). The risk improvement by the optimal shift is evaluated by the inner product of the optimal shifts $\beta^{\kappa}_{\rm opt}$ and given by

[TABLE]

which does not depend on the parallel shift $\alpha$. Here, \raise 1.5pt\hbox{\displaystyle\mathop{H}^{\rm m}}_{ab}^{\,~{}\lambda}\raise 1.5pt\hbox{\displaystyle\mathop{H}^{\rm m}}_{cd}^{\,~{}\kappa}g^{ab}g^{cd}g_{\kappa\lambda} is the mixture mean curvature of $\mathcal{P}$ embedded in $\mathcal{E}$ at $\omega$, thus this risk improvement has a geometrical interpretation.

These results are related to the properties of Bayesian predictive densities. In [11], it is shown that Bayesian predictive densities are optimal along orthogonal shift from $\mathcal{P}$. The orthogonal shift is not included in the tangent space of $\mathcal{E}$, as explained in Section III-B. The risk improvement achieved by the orthogonal shift is the mixture mean curvature of $\mathcal{P}$ embedded in $\mathcal{F}$, while the risk improvement of $p(y;\hat{\eta}_{\pi})$ is the mixture mean curvature of $\mathcal{P}$ embedded in $\mathcal{E}$. The risk improvement is evaluated by inner products of the optimal shifts, and as a result the cosine of the angle between the two orthogonal shifts, as shown in Figure 2, can be found as the square root of the ratio of the two risk improvements.

Example  (Fisher circle model, continued) We have $g_{\omega\omega}=1$ and \raise 1.5pt\hbox{\displaystyle\mathop{\Gamma}^{\rm m}}_{\omega\omega}^{\,~{}\omega}=0. Thus, the optimal orthogonal shift is

[TABLE]

and the risk improvement obtained by the optimal orthogonal shift is

[TABLE]

The risk improvement corresponding to the optimal shift (7) is ${3}/({8n^{2}})+o(n^{-2})$.

If the variance of $x_{1},x_{2}$ is $\sigma^{2}$, the risk improvements corresponding to the optimal orthogonal shift and to the shift (7) can be obtained as follows, respectively:

[TABLE]

Therefore, when $\sigma^{2}$ is large, the risk improvement obtained by $p(y;\hat{\eta}_{\pi})$ becomes relatively significant as well, and the performance of the Bayes extended estimator is close to that of the Bayesian predictive density. The cosine of the angle between the two shift vectors is

[TABLE]

which approaches [math] as $\sigma^{2}$ increases. In this way, the Bayesian predictive density and $p(y;\hat{\eta}_{\pi})$ approach each other.

IV Numerical studies

The numerical simulations of the Kullback–Leibler risk of $p(y;\hat{\eta}_{\pi})$ and Bayesian predictive densities are shown in a curved Gaussian model. It confirmed the theoretical results so far and also illustrates the practical importance of the projection of Bayesian predictive densities.

The model $\mathcal{P}$ is the spiked covariance model (for more information about related models, see e.g. [12]), that is $l$-dimensional Gaussian $\mathrm{N}_{l}(0,\Sigma)$ with the covariance matrix $\Sigma$ expressed as

[TABLE]

where the vector \mbox{\boldmathu}\in\mathbb{R}^{l} satisfies \mbox{\boldmathu}^{\top}\mbox{\boldmathu}=1 and $\lambda>0$. The eigenvalues of the matrix $\Sigma$ are $\lambda+1,1,\dots,1$, and $u$ is the first eigenvector.

The model $\mathcal{P}$ parametrized by \omega=(\mbox{\boldmathu},\lambda) is embedded in the larger full exponential family $\mathcal{E}=\{\mathrm{N}_{l}(0,\Sigma)\mid\Sigma\}$, and the expectation parameter $\eta$ comprises the components of $\Sigma$. The extended plugin distribution $p(y;\hat{\eta}_{\pi})$ is $\mathrm{N}_{l}(0,\hat{\Sigma}_{\pi})$, where $\hat{\Sigma}_{\pi}$ is the posterior mean of $\Sigma$. The Bayes estimator (\hat{\lambda}_{\pi},\hat{\mbox{\boldmathu}}_{\pi}) of $\mathcal{P}$ is composed of $\hat{\lambda}_{\pi}$, the posterior mean of $\lambda$, and \hat{\mbox{\boldmathu}}_{\pi}, the first eigen vector of $\hat{\Sigma}_{\pi}$. The plugin distribution of the Bayes estimator is \mathrm{N}_{l}(0,\hat{\lambda}_{\pi}\hat{\mbox{\boldmathu}}_{\pi}\hat{\mbox{\boldmathu}}_{\pi}^{\top}+I_{l}). The settings of $l$ and $n$ are $(l,n)=(5,20)$ and $(80,320)$. The posterior mean of $\lambda$ and $\Sigma$ is computed by the 1000 MCMC samples for $l=5$ and 2000 MCMC samples for $l=80$ produced by Gibbs sampling with $250$ and $500$ burn-in samples, respectively. The Bayesian predictive density is computed by taking the mean of the plugin densities of those MCMC samples of (\lambda,\mbox{\boldmathu}). The Kullback–Leibler risk is derived as the mean of 2000 trials.

The results are illustrated in Figure 4, and it is confirmed that $p(y;\hat{\eta}_{\pi})$ approaches the Bayesian predictive density as the size of $\mathcal{E}$ increases. In Figure 4(a), the three-layer structure of the plugin density, the extended plugin density, and the Bayesian predictive density is seen in the risk comparison. On the other hand, in Figure 4(b), it could be seen that the risk plots of $p(y;\hat{\eta}_{\pi})$ and the Bayesian predictive density are quite close, which means the projection angle between them is close to zero. The dimension of the parameter space of $\mathcal{P}$ is $l$ and that of $\mathcal{E}$ is $l(l+1)/2$; thus, $\mathcal{P}$ is embedded in relatively larger full exponential families when $l$ increases. Therefore, it is natural that the two risk performances approach as $l$ increases because the extended model $\mathcal{E}$ approaches the set of all probability densities $\mathcal{F}$, and $p(y;\hat{\eta}_{\pi})$ approaches the Bayesian predictive density.

Figure 4(b) also illustrates the practical advantage of $p(y;\hat{\eta}_{\pi})$, as it shows that projection of Bayesian predictive densities onto finite-dimensional models of reasonable size is an effective way of approximation. Bayesian predictive densities are typically approximated by the mean of plugin densities, because obtaining the full density function is intractable. In these experiments, the Bayesian predictive density is the mean of 2000 plugin densities. The problem of this approximation is that it requires large space and time to compute density, because storing all MCMC samples is necessary, and we have to take the mean of these plugin densities for each $y$. Figure 4(b) demonstrates that the approximation by $p(y;\hat{\eta}_{\pi})$, a single point in $\mathcal{E}$ is comparable to the mean of 2000 points in $\mathcal{P}$. Approximation by $p(y;\hat{\eta}_{\pi})$ does not require storing MCMC samples, and the full density function is available avoiding taking mean of plugin densities every time. Table II provides the computation time required to evaluate the density of 1000 new samples $y$ and the memory size of the Bayesian predictive density and $p(y;\hat{\eta}_{\pi})$. It is shown that both computational time and memory size are effectively saved when we utilize $p(y;\hat{\eta}_{\pi})$. Therefore, the extended plugin density $p(y;\hat{\eta}_{\pi})$ is an effective approximation of Bayesian predictive densities, with a large advantage in terms of computational costs, which has a natural theoretical interpretation as a projection onto a finite dimensional model.

Appendix A Proof of Theorem III.1

Let $L(\eta)=\frac{1}{n}\sum_{t=1}^{n}\log p(x(t);\eta)$. Then, $\hat{\eta}_{\pi}(x^{n})=\{(\hat{\eta}_{\pi})_{i}(x^{n})\}$ is given by

[TABLE]

We approximate $\hat{\eta}_{\pi}$ using the Laplace method (e.g., [13], Sec. 3.6 and Sec. 4.6, and [14]). The proposed approximation follows Theorem 4.6.1 in [13]. First we expand $\pi(\omega)\exp(nL(\omega))$ around $\omega=\hat{\omega}_{\rm MLE}$. In the following, symbols such as $\eta(\hat{\omega}_{\rm MLE}),\partial_{a}\eta(\hat{\omega}_{\rm MLE})$, and $\partial_{a}\partial_{b}\eta(\hat{\omega}_{\rm MLE})$ are abbreviated to $\hat{\eta},\partial_{a}\hat{\eta}$, and $\partial_{ab}\hat{\eta}$ respectively. Let $\omega=\hat{\omega}_{\rm MLE}+\phi/\sqrt{n}$. Then,

[TABLE]

[TABLE]

Here, $\hat{J}_{ab}=-\partial_{ab}\hat{L}$ and we denote the terms of $\mathrm{O}_{p}(n^{-1})$ which do not depend on $\pi$ as $C_{1}/n$, which is $\mathrm{O}_{p}(n^{-1})$.

Let $(\hat{J}^{ab})$ be the inverse matrix of $(\hat{J}_{ab})$. Next, we integrate both sides of the above equation with respect to $\omega$. By changing the variables from $\omega$ to $\phi$, and by using the formula of moments of multivariate Gaussian distributions, we obtain

[TABLE]

Here $C_{2}$ does not depend on $\pi$.

Replace $\pi(\omega)$ by $\eta_{i}(\omega)\pi(\omega)$ and we have

[TABLE]

Therefore, the posterior mean of $\eta_{i}$ is expanded as

[TABLE]

Since $\hat{J}_{ab}=\hat{g}_{ab}+o_{p}(1),~{}\partial_{bcd}\hat{L}={E}(\partial_{bcd}\hat{L})+o_{p}(1)$ by the law of large numbers, and

[TABLE]

and

[TABLE]

hold, we obtain

[TABLE]

Appendix B Proof of Proposition III.2

We abbreviate symbols such as $g_{ab}(\omega)$, \raise 1.5pt\hbox{\displaystyle\mathop{H}^{\rm m}}_{ab\kappa}^{\,}(\omega), and $p(y;\omega)$ to $g_{ab}$, \raise 1.5pt\hbox{\displaystyle\mathop{H}^{\rm m}}_{ab\kappa}^{\,}, and $p$, respectively.

The Kullback–Leibler divergence from $p(y;\omega)$ to $p_{\alpha,\beta}(y;\hat{\omega}_{\mathrm{MLE}})$ is expanded as

[TABLE]

and we can expand $\log p(y;\hat{\omega},\hat{\xi})-\log p$ as follows:

[TABLE]

where $\tilde{\omega}=\hat{\omega}-\omega$. Because $E[\partial_{a}\log p]=0,E[-\partial_{ab}\log p]=g_{ab},$ and $g_{a\kappa}=0$ for $a=1,\dots,d,~{}\kappa=d+1,\dots,m$, we have

[TABLE]

where $\tilde{\omega}_{\rm MLE}=\hat{\omega}_{\rm MLE}-\omega$ and

[TABLE]

Therefore, the Kullback–Leibler risk from $p(y;\omega)$ to $p_{\alpha,\beta}(y;\hat{\omega}_{\mathrm{MLE}})$ is expanded as

[TABLE]

Because $\hat{\beta}$ is a smooth function of $\rm{O}_{p}(1)$ and $\hat{\beta}=\beta+{\rm o}_{p}(1)$,

[TABLE]

hold and

[TABLE]

Lastly we derive ${E}(\tilde{\omega}^{a}\tilde{\omega}^{b})$. It can be expanded as (see e.g. (10.19) in [15], as its multivariate version is presented here)

[TABLE]

From the likelihood equation,

[TABLE]

Because $\tilde{\omega}_{\rm MLE}^{b}=g^{ab}\partial_{a}L+{\rm O}_{p}(n^{-1})$ and $\tilde{\omega}_{\rm MLE}=\hat{\omega}_{\rm MLE}-\omega=\hat{\omega}-\alpha/n-\omega+\rm{o}_{p}(n^{-1})$,

[TABLE]

Thus

[TABLE]

and

[TABLE]

Using E[\partial_{d}L\partial_{ce}L]=\raise 1.5pt\hbox{\displaystyle\mathop{\Gamma}^{\rm e}}_{ecd}^{\,}/n and E[\partial_{cde}L]=-\raise 1.5pt\hbox{\displaystyle\mathop{\Gamma}^{\rm e}}_{cde}^{\,}-\raise 1.5pt\hbox{\displaystyle\mathop{\Gamma}^{\rm e}}_{dec}^{\,}-\raise 1.5pt\hbox{\displaystyle\mathop{\Gamma}^{\rm e}}_{ecd}^{\,}-T_{ecd}, we have

[TABLE]

(in the last equation we use g^{be}g^{cd}\raise 1.5pt\hbox{\displaystyle\mathop{\Gamma}^{\rm e}}_{ecd}^{\,}=g^{be}g^{cd}\raise 1.5pt\hbox{\displaystyle\mathop{\Gamma}^{\rm e}}_{dec}^{\,} and g^{ae}g^{cd}\raise 1.5pt\hbox{\displaystyle\mathop{\Gamma}^{\rm e}}_{ecd}^{\,}=g^{ae}g^{cd}\raise 1.5pt\hbox{\displaystyle\mathop{\Gamma}^{\rm e}}_{dec}^{\,}.) Hence we obtain

[TABLE]

Here,

[TABLE]

and

[TABLE]

hold. The first one follows from the definition, and the second one is from the relation

[TABLE]

and $g_{a\kappa}=0$ ($a=1,\dots,d,~{}\kappa=d+1,\dots,m$).

Therefore, the risk is expanded as follows:

[TABLE]

Acknowledgment

The authors greatly appreciate the referees’ comments on an earlier version. The authors are grateful to Tomonari Sei for helpful comments. This paper is based on a part of the first author’s Ph.D. thesis done at the University of Tokyo.