Bayesian prediction regions and density estimation with type-2 censored data
111.

Akbar Asgharzadeha, Éric Marchandb & Ali Saadati Nika

a Department of Statistics, University of Mazandaran, P.O. Box 47146-1407,
Balbosar, IRAN
b Département de mathématiques, Université de Sherbrooke, Sherbrooke Qc,
CANADA, J1K 2R1
(e-mails<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
Summary

For exponentially distributed lifetimes, we consider the prediction of future
order statistics based on having observed the first $m$ order statistics. We
focus on the previously less explored aspects of predicting: (i) an arbitrary
pair of future order statistics such as the next and last ones, as well as
(ii) the next $N$ future order statistics. We provide explicit and exact
Bayesian credible regions associated with Gamma priors, and constructed by
identifying a region with a given credibility $1-\lambda$ under the Bayesian
predictive density. For (ii), the HPD region is obtained, while a two-step
algorithm is given for (i). The predictive distributions are represented as
mixtures of bivariate Pareto distributions, as well as multivariate Pareto
distributions. For the non-informative prior density choice, we demonstrate
that a resulting Bayesian credible region has matching frequentist coverage
probability, and that the resulting predictive density possesses the
optimality properties of best invariance and minimaxity.

AMS 2020 subject classifications:  62F15, 62N01, 62N05, 62C10, 62C20.

Keywords and phrases: Bayesian predictive density; Credibility; Coverage
probability; Mixtures; Multivariate Pareto; Prediction region; Type-2
censoring.

##  1 Introduction

Predictive analysis based on censored data in life testing experiments is
fundamental and leads to interesting challenges. In this work, we are
concerned with prediction regions for future order statistics based on the
first $m$ order statistics generated by exponentially distributed data. There
has been some previous work on such problems, but we focus here on the less
explored: (a) multivariate aspects, and (b) use of Bayesian predictive
densities to generate prediction regions and their related theoretical
properties.

We consider an i.i.d. sample of size $n$ from an exponential distribution with
density $\theta e^{-\theta t}\mathbb{I}_{(0,\infty)}(t)$, but we are only able
to observe the first $m$ order statistics $X_{1:n},\ldots,X_{m:n}$, commonly
referred to as type-II censoring. Our objectives include:

  1. (I)

the joint prediction of two future order statistics $X_{r:n}$ and $X_{s:n}$,
with $m<r<s\leq n$, as well as

  2. (II)

the joint prediction of the next $N$ order statistics with $2\leq N\leq n-m$.

Scenario (II) seems a natural one to consider and includes the case of all
future order statistics with $N=n-m$. Scenario (I) in more relevant in
situations where the focus is for instance on the next and last order
statistics (i.e., $r=m+1,s=n$), or the last two order statistics (i.e.,
$r=n-1,s=n$). Obviously, both scenarios overlap for $N=2$. Surprisingly
perhaps, for Gamma priors which we consider, the specification of Bayesian
predictive densities and regions leads to more complex representations in the
bivariate case. The bivariate case (I) was considered recently by Bagheri et
al.  (2022) and we refer to this work for further motivation and historical
aspects of the problem. Their prediction regions are non-Bayesian however and
based on a pair pivotal quantities which will not arise as a Bayesian
solution.

It seems natural to generate prediction regions for future order statistics
via a Bayesian predictive density, but such problems seem to have been
relatively unexplored. An exception is given by Dunsmore  (1974), but his work
concerns a single order statistic. We proceed in doing so for Gamma prior
densities as well as the non-informative density
$\pi_{0}(\theta)=\frac{1}{\theta}\,\mathbb{I}_{(0,\infty)}(\theta)$. The
obtained expressions are quite tractable and, interestingly, bring into play
finite mixtures of bivariate Pareto distributions, and multivariate Pareto
distributions. The obtained mixtures are “non-probabilistic” in the sense that
the coefficients $a_{i}$ in $\sum_{i}a_{i}f_{i}(t)$; the $f_{i}$’s being
densities; take on both positive and negative values. This aspect is less
appreciated than “probabilistic mixtures”, but it does not hinder the
usefulness of the representation for computational purposes of moments and
cumulative distributions functions.

For a given credibility $1-\lambda$, we fully describe the highest posterior
density (HPD) prediction region for scenario (II), while we propose an
algorithm to generate exact prediction regions for the bivariate scenario (I),
based on a natural decomposition of the joint predictive density into its
marginal and conditional parts. The analysis is carried out and much
facilitated by considering the prediction of spacings between future order
statistics, which can be converted back to the prediction of the order
statistics themselves.

Many of the recent studies on predictive density estimation focus on
theoretical properties of the density itself in a decision-theoretic context.
For scenarios (I) and (II), we report on the best invariant and minimax
properties for Kullback-Leibler divergence loss of the Bayesian predictive
density associated with the usual non-informative prior density
$\pi(\theta)=\frac{1}{\theta}$, attributable to the work of Liang & Barron
(2004). Moreover, we show that the frequentist probabilty of coverage for a
prediction region generated by such a Bayesian predictive density matches its
exact credibility. Therefore, such a Bayesian density not only possesses
optimality properties in a decision-theoretic framework, but also provides a
satisfactory frequentist option that compares favourably with previous
solutions (e.g., Bagheri et al., 2022).

The organization of the manuscript is as follows. The preliminary results of
Chapter 2 cover model densities for scenarios (I) and (II), some general
aspects on Bayesian predictive densities, and multivariate Pareto
distributions. Bayesian predictive densities and regions are derived and
illustrated in Section 3. In Section 4, we demonstrate in a more general
context including ours that a Bayesian credible region with credibility
$1-\lambda$ associated with the non-informative prior density $\pi_{0}$ yields
matching frequentist coverage probability $1-\lambda$ for all $\theta>0$.
Finally in Section 5, we concisely review optimality properties of the
Bayesian predictive densities again associated with density $\pi_{0}$ and
pertaining to the best invariant and minimaxity properties.

##  2 Preliminary results

###  2.1 Model densities

For the first $m$ order statistics $X_{1:n},\ldots,X_{m:n}$ among $n$
generated from i.i.d. $Exp(\theta)$ data, it is well known (e.g., Lawless
(1971)) $X=\sum_{i=1}^{m}X_{i:n}+(n-m)X_{m:n}$ is a sufficient statistic and
Gamma distributed $\mathcal{G}(m,\theta)$. Hereafter, we therefore consider
such a summary and corresponding density, that is

$X\sim\frac{\theta^{m}}{\Gamma(m)}\,x^{m-1}\,e^{-\theta
x}\,\mathbb{I}_{(0,\infty)}(x).$ (2.1)

For the prediction of $X_{r:n}$ and $X_{s:n}$ based on $X$, it is convenient
and equivalent to consider

$Y=(Y_{1},Y_{2})^{\top}\hbox{ with
}Y_{1}\,=\,X_{r:n}-X_{m:n},\,Y_{2}\,=\,X_{s:n}-X_{r:n}.$ (2.2)

The equivalence stems from the correspondence between a prediction region $R$
for $Y$ and the inverse mapping
$\\{(y_{1}+x_{m:n},y_{1}+y_{2}+x_{m:n}):(y_{1},y_{2})\in R\\}$ as a prediction
region for $(X_{r:n},X_{s:n})$.

For the multivariate version where the objective is to predict jointly the $N$
next future order statistics, it is analogously useful to consider

$Z=(Z_{1},\ldots,Z_{N})\,,\hbox{ with }\,Z_{i}=X_{m+i:n}-X_{m+i-1:n}\hbox{ for
}i\in\\{1,\ldots,N\\},2\leq N\leq n-m\,,$ (2.3)

as the objects of prediction.

As described with the next result, both transformations (2.2) and (2.3) lead
to convenient underlying model densities for $Y$ or $Z$.

######  Lemma 2.1.

Given set-up (2.1, 2.2, and 2.3) and fixed $\theta$:

  1. (a)

$Y$ and $X$ are independently distributed;

  2. (b)

$Z$ and $X$ are independently distributed;

  3. (c)

$Y_{1}$ and $Y_{2}$ are independently distributed with joint density on
$\mathbb{R}_{+}^{2}$ given by $q_{\theta}(y)\,=\,\theta^{2}\,q_{1}(\theta
y_{1},\theta y_{2})$ and with

$q_{1}(y)\,=\,\frac{(n-m)!}{(r-m-1)!\,(s-r-1)!\,(n-s)!}\,\big{(}1-e^{-y_{1}}\big{)}^{r-m-1}\,\big{(}1-e^{-y_{2}}\big{)}^{s-r-1}\,e^{-(n-r+1)y_{1}-(n-s+1)y_{2}}\,;$
(2.4)

  4. (d)

$Z_{1},\ldots,Z_{N}$ are independently distributed with $Z_{i}\sim
Exp\big{(}(n-m-i+1)\,\theta\big{)}$.

Proof. Part (b) follows from a familiar renewal property of the exponential
distribution, and furthermore implies (a) since $Y$ is a function $h(Z)$ of
$Z$. For (c), since $(U=X_{r:n}-X_{m:n},V=X_{s:n}-X_{m:n})$ and $X_{m:n}$ are
independently distributed, the distribution of $(U,V)$ matches that of its
conditional distribution given $X_{m:n}$, and the latter joint distribution
can by seen to be equivalent, with the above-mentioned renewal property, to
that of the $(r-m)^{\hbox{th}}$ and $(s-m)^{\hbox{th}}$ order statistics from
a sample of size $n-m$ from an $Exp(\theta)$ distribution. Such a joint
density is given by

$f(u,v)\,=\,\frac{(n-m)!}{(r-m-1)!\,(s-r-1)!\,(n-s)!}\,\big{(}1-e^{-\theta
u}\big{)}^{r-m-1}\,\big{(}e^{-\theta u}-e^{-\theta
v}\big{)}^{s-r-1}\,(e^{-\theta v})^{n-s}\,\theta^{2}\,e^{-\theta(u+v)},\,$
(2.5)

for $0<u<v$. The result follows by transforming $(U,V)$ to $Y$. Finally for
(d), since $Z$ and $X_{m:n}$ are independently distributed, the distribution
of $Z$ matches that of its conditional distribution given $X_{m:n}$, and as
above the latter joint distribution matches that of the distribution of the
first $N$ order statistics spacings from a sample of size $N$ from an
$Exp(\theta)$ distribution. It is well known that such order statistics
spacings are independently distributed as $Exp\big{(}(n-m-i+1)\theta\big{)}$
(e.g., Lehmann & Casella, 1998, problem 6.18, page 71) which yields the
result. ∎

######  Remark 2.1.

The marginal densities of $Y_{1}$ and $Y_{2}$ in (2.4) can be extracted from
part (c) and are of the form
$\mathcal{B}^{-1}(c_{1},c_{2})\,(1-e^{-u})^{c_{1}-1}\,e^{-c_{2}u}\,\mathbb{I}_{(0,\infty)}(u)$,
with $\mathcal{B}(c_{1},c_{2})$ the Beta function. An alternative and readily
verified representation for the distribution of the above $Y_{1}$ and $Y_{2}$
is:

$e^{-\theta Y_{1}}\sim\hbox{Beta}(n-r+1,r-m)\,,\,e^{-\theta
Y_{2}}\sim\hbox{Beta}(n-s+1,s-r)\,\hbox{ independent}\,.$

###  2.2 Bayesian predictive densities

A general set-up for predictive density estimation relates to the following
model density:

$(X,Y)|\theta\sim
p_{\theta}(x)\,q_{\theta}(y|x)\,,x\in\mathbb{R}^{d_{1}},y\in\mathbb{R}^{d_{2}}.$

We observe $X$ according to $p_{\theta}$ and we wish to estimate the density
$q_{\theta}(\cdot|x)$. Except for $\theta$, the densities are known. The
observed $X$ provides information about $\theta$ and determines the
conditional density $q_{\theta}(\cdot|x)$ to estimate. Much, but not all
(e.g., Fourdrinier et al. 2019), previous work on properties of predictive
densities focusses on models for which $X$ and $Y$ are conditionally
independent, i.e., $q_{\theta}(\cdot|x)\equiv q_{\theta}(\cdot)$ for all $x$,
but it is natural in general to estimate the density $q_{\theta}(\cdot|x)$
when there is dependence.

Assume that we have a prior density $\pi$ for $\theta$, and a resulting
posterior density $\pi(\cdot|x)$ taken to be absolutely continuous with
respect to measure $\nu$. A natural estimator of the density
$q_{\theta}(\cdot|x)$ is the Bayes predictive density given by the conditional
density $q(\cdot|x)$ of $Y|X=x$. Integrating out $\theta$, we obtain the
estimator

$\hat{q}_{\pi}(y|x)\,=\,\int_{\Theta}q_{\theta}(y|x)\,\pi(\theta|x)\,d\nu(\theta).$
(2.6)

The above is a fully Bayesian procedure that can be used for obtaining
predictive point estimates or prediction regions for $Y$. The estimator or
density $\hat{q}_{\pi}(\cdot;X)$ will naturally depend on $X$ through a
sufficient statistics $T(X)$.

###  2.3 Multivariate Pareto densities

The predictive densities that we elicit below bring into play univariate,
bivariate, and multivariate type II Pareto distributions, as well as non-
probabilistic mixtures of such distributions. Such distributions have been
extensively studied (see for instance Section 52.4 of Kotz et al. 2000, or
Arnold 2014, and the references therein). With convenient expressions for the
moments and cumulative distribution functions (c.d.f.’s) relative to the
$f_{i}$’s, the mixture representations clearly facilitate expressions for
moments and c.d.f.’s for the full distribution.

Multivariate Pareto distributions possess Pareto marginals with parameters
$\ell,h>0$, densities and survival functions, which we denote and write as

$f_{\ell,h}(t)\,=\,\frac{\ell\,h}{\;\;\;(1+ht)^{\ell+1}}\,\mathbb{I}_{(0,\infty)}(t)\,,\hbox{
and }\bar{F}_{l,h}(t)\,=\,(1+ht)^{-\ell}\,,$ (2.7)

for $t>0$, respectively. Here is a multivariate Pareto definition; which
includes the bivariate case; and some useful properties.

######  Definition 2.1.

A random vector $Z=(Z_{1},\ldots,Z_{N})^{\top}$ has a multivariate Pareto type
II distribution, denoted $Z\sim\mathcal{P}2(m,h_{1},\ldots,h_{N})$ when it has
density

$g_{m,h_{1},\ldots,h_{N}}(z)\,=\,\frac{(m)_{N}\prod_{i=1}^{N}h_{i}}{(1+\sum_{i=1}^{N}h_{i}z_{i})^{m+N}}\,\mathbb{I}_{\mathbb{R}_{+}^{N}}(z),$
(2.8)

for parameters $h_{1},\ldots h_{N},m>0$, and
$(m)_{N}=\frac{\Gamma(m+N)}{\Gamma(m)}$.

Such distributions form a $N+1$ parameter family with scale parameters
$1/h_{i},i=1,\ldots,N$, and shape parameter $m$. As recorded with the
following lemma containing a selection of properties that are known and
readily verified, such multivariate Pareto distributions possess univariate
Pareto marginals and conditionals, subvectors that are also distributed as
multivariate Pareto with densities, and a joint survival function having a
rather simple form.

######  Lemma 2.2.

Consider $Z\sim\mathcal{P}2(m,h_{1},\ldots,h_{N})$. Let $i\neq j$. Then, we
have: (i) $Z_{i}\sim f_{m,h_{i}}$, (ii) $Z_{j}|Z_{i}=z_{i}\sim
f_{m+1,h_{j}/(1+h_{i}z_{1})}$, (iii) $(h_{1}Z_{i},\ldots,h_{N}Z_{N})$ has
density $g_{m,1,\ldots,1}$, (iv) the joint survival function is given by
$\mathbb{P}(\cap_{i=1}^{N}\\{Z_{i}\geq
z_{i}\\})\,=\,\\{(1+\sum_{i=1}^{N}h_{i}z_{i})\\}^{-m}$ for $z_{i}\geq
0,i=1,\ldots,N$, (v) the simple regressions are linear with (vi)
$\mathbb{E}(Z_{j}|Z_{i})\,=\,\frac{1+h_{i}Z_{i}}{h_{j}m}\,$, and (vi) the
correlation between $Z_{i}$ and $Z_{j}$ is given by
$\rho(Z_{i},Z_{j})=\frac{1}{m}$ for $m>2$.

Furthermore, we will require the following.

######  Lemma 2.3.

Let $Z\sim\mathcal{P}2(m,h_{1},\ldots,h_{N})$ and
$W=\sum_{i=1}^{N}h_{i}Z_{i}$. Then $W$ is Beta type-II distributed (denoted
$W\sim B2(N,m)$) with p.d.f.
$f_{W}(w)\,=\,\frac{\Gamma(N+m)}{\Gamma(N)\Gamma(m)}\frac{w^{N-1}}{(1+w)^{N+m}}$,
for $w\in(0,\infty)$.

Proof. The result is known but we provide a proof for completeness. Given the
multivariate Pareto representation
$(h_{1}Z_{1}\,\ldots,h_{N}Z_{N})=^{d}(\frac{E_{1}}{G},\ldots,\frac{E_{N}}{G})$
with $E_{1},\ldots E_{N},G$ independently distributed, and with $E_{i}\sim
Exp(1)$ and $G\sim\mathcal{G}(m,1)$, we see that $W$ is distributed as the
ratio of two independent $\mathcal{G}(N,1)$ and $\mathcal{G}(m,1)$ variables,
hence Beta type II with the given parameters. ∎

##  3 Predictive densities and regions

Based on $X$ as in (2.1), we provide in this section Bayesian predictive
densities and regions for $Y$ and $Z$ as defined in (2.2) and (2.3). We
consider Gamma $\mathcal{G}(\alpha,\beta)$ prior densities
$\pi_{\alpha,\beta}(\theta)\propto\theta^{\alpha-1}\,e^{-\beta\theta}\,\mathbb{I}_{(0,\infty)}(\theta)$,
including the usual non-informative case
$\pi_{0}(\theta)\,=\,\frac{1}{\theta}\,\mathbb{I}_{(0,\infty)}(\theta)$ for
$\alpha=\beta=0$.

###  3.1 Predictive densities

We begin with the future next $N$ order statistics.

######  Theorem 3.1.

The Bayes predictive density of $Z$ in (2.3), based on
$X\sim\mathcal{G}(m,\theta)$ and prior density $\pi_{\alpha,\beta}$ for
$\theta$, is that of $\mathcal{P}2(m+\alpha,h_{1},\ldots,h_{N})$ distribution
with $h_{i}=\frac{n-m-i+1}{x+\beta}$ for $i=1,\ldots,N$.

Proof. With
$q_{\theta}(z|x)\,=\,\frac{\theta^{N}\,(n-m)!}{(n-m-N)!}\,e^{-\theta\sum_{i=1}^{N}(n-m-i+1)z_{i}}$
and $\theta|x\sim\mathcal{G}(\alpha+m,x+\beta)$, we obtain from (2.6):

$\displaystyle\hat{q}_{\pi_{\alpha,\beta}}(z|x)\,$ $\displaystyle=$
$\displaystyle\,\frac{(n-m)!\,(\beta+x)^{\alpha+m}}{(n-m-N)!\,\Gamma(\alpha+m)\,}\int_{0}^{\infty}\theta^{N+\alpha+m-1}\,e^{-\theta\big{(}\beta+x+\sum_{i=1}^{N}(n-m-i+1)z_{i}\big{)}}\,d\theta\,$
$\displaystyle=$
$\displaystyle\,\frac{(n-m)!\,(\beta+x)^{\alpha+m}}{(n-m-N)!\,\Gamma(\alpha+m)\,}\,\Gamma(N+\alpha+m)\,\big{\\{}\beta+x+\sum_{i=1}^{N}(n-m-i+1)z_{i}\big{\\}}^{-(N+\alpha+m)},$

which is indeed a $\mathcal{P}2(m+\alpha,h_{1},\ldots,h_{N})$ density. ∎

Observe that the density has univariate Pareto marginals and multivariate
Pareto distributed subvectors in accordance to Lemma 2.2.

We now turn to the Bayesian predictive density for two future order statistics
and demonstrate a non-probabilistic mixture of bivariate Pareto densities as
given in Definition 2.1 for $N=2$.

######  Theorem 3.2.

The Bayes predictive density of $Y$ in (2.2), based on
$X\sim\mathcal{G}(m,\theta)$ and prior density $\pi_{\alpha,\beta}$ for
$\theta$, is given by

$\hat{q}_{\pi_{\alpha,\beta}}(y|x)\,=\,\sum_{i=0}^{r-m-1}\sum_{j=0}^{s-r-1}\omega_{i,j}\,g_{m+\alpha,\frac{a_{i}}{x+\beta},\frac{b_{j}}{x+\beta}}\big{(}y\big{)}\,,$
(3.9)

with $g$ a bivariate Pareto density given in (2.8) for $N=2$,

$a_{i}=n-r+i+1,b_{j}=n-s+j+1,\hbox{and
}\omega_{i,j}\,=\,\frac{(n-m)!}{(n-s)!}\,\frac{(-1)^{i+j}}{i!\,j!}\,\,\frac{1}{(r-m-i-1)!\,(s-r-j-1)!}\frac{1}{a_{i}b_{j}}.$

Proof. From (2.4) and (2.6), with $\theta|x\sim\mathcal{G}(m+\alpha,x+\beta)$,
we obtain

$\hat{q}_{\pi_{\alpha,\beta}}(y|x)\,=\,\frac{(n-m)!}{(n-s)!}\frac{\,(x+\beta)^{m+\alpha}}{\Gamma(m+\alpha)}\,\int_{0}^{\infty}\theta^{m+\alpha+1}\,\frac{(1-e^{-\theta
y_{1}})^{r-m-1}}{(r-m-1)!}\,\frac{(1-e^{-\theta
y_{2}})^{s-r-1}}{(s-r-1)!}e^{-\theta(L(y)+x+\beta)}\,d\theta,$

with $L(y)\,=\,(n-r+1)y_{1}+(n-s+1)y_{2}$. Binomial expansions and an
interchange of sum and integral yield

$\displaystyle\hat{q}_{\pi_{\alpha,\beta}}(y|x)\,$ $\displaystyle=$
$\displaystyle\,\frac{\,(x+\beta)^{m+\alpha}}{\Gamma(m+\alpha)}\,\sum_{i=0}^{r-m-1}\sum_{j=0}^{s-r-1}a_{i}b_{j}\,\omega_{i,j}\,\int_{0}^{\infty}\theta^{m+\alpha+1}\,e^{-\theta(L(y)+iy_{1}+jy_{2}+x+\beta)}\,d\theta\,$
$\displaystyle=$
$\displaystyle\,(m+\alpha)\,(m+\alpha+1)\,(x+\beta)^{m+\alpha}\,\sum_{i=0}^{r-m-1}\sum_{j=0}^{s-r-1}a_{i}b_{j}\,\omega_{i,j}\,\big{(}L(y)+iy_{1}+jy_{2}+x+\beta\big{)}^{-(m+\alpha+2)},$

which leads to (3.9). ∎

We point out the following result for the next two order statistics, which
follows immediately from either Theorem 3.1 or Theorem 3.2.

######  Corollary 3.1.

For the particular case where $r=m+1,s=m+2$, the predictive density of $Y$ as
in (2.2) is that of a $\mathcal{P}2(m+\alpha,h_{1},h_{2})$ distribution (3.9)
with $h_{1}=\frac{n-m}{x+\beta}$ and $h_{2}=\frac{n-m-1}{x+\beta}$.

######  Remark 3.2.

Interestingly, the above weights $\omega_{i,j}$ arise through a series
expansion of the Beta function. Indeed, by the the binomial expansion of
$(1-t)^{d-1}$ for $d\in\mathbb{N}_{+}$, we obtain for $c>0$

$\displaystyle\int_{0}^{1}t^{c-1}\,(1-t)^{d-1}\,dt\,$ $\displaystyle=$
$\displaystyle\,\frac{\Gamma(c)\Gamma(d)}{\Gamma(c+d)}$
$\displaystyle\Longrightarrow$
$\displaystyle\sum_{k=0}^{d-1}\gamma_{c,d,k}=1\,,$

with
$\gamma_{c,d,k}\,=\,\frac{\Gamma(c+d)}{\Gamma(c)}\,\frac{(-1)^{k}}{k!\,(d-1-k)!}\,\frac{1}{c+k}$.
Furthermore, observe that $\omega_{i,j}=\omega_{1,i}\,\omega_{2,j}$ with
$\omega_{1,i}=\frac{(n-r)!}{(n-s)!}\,\gamma_{n-r+1,r-m,i}$ and
$\omega_{2,j}=\frac{(n-s)!}{(n-r)!}\,\gamma_{n-s+1,s-r,j}$, demonstrating
alternatively that $\sum_{i,j}\omega_{i,j}=1$.

The mixture representation of the predictive density in Theorem 3.2 coupled
with the bivariate Pareto properties of Lemma 2.2 facilitate the evaluation of
the marginal and conditional distributions associated with (3.9). As shown
below, mixture representations of univariate Pareto distributions arise.

######  Corollary 3.2.

For $m<r\leq n-1$, the marginal density of $Y_{1}$ associated with the Bayes
predictive density (3.9) is given by

$\hat{q}_{\pi}(y_{1}|x)\,=\,\sum_{i=0}^{r-m-1}\gamma_{n-r+1,r-m,i}\,f_{m+\alpha,\frac{a_{i}}{x+\beta}}(y_{1})\,,$
(3.10)

where $\gamma_{c,d,k}$ is given in Remark 3.2 and $f_{l,h}$ is a univariate
Pareto density as given in (2.7). In the particular case where $r=m+1$, we
have that $\hat{q}_{\pi}(y_{1};x)=f_{m+\alpha,\frac{n-m}{x+\beta}}(y_{1})$ for
$y_{1}>0$.

Proof. We have from Theorem 3.2

$\displaystyle\hat{q}_{\pi}(y_{1}|x)\,$ $\displaystyle=$
$\displaystyle\,\int_{0}^{\infty}\hat{q}_{\pi}(y|x)\,dy_{2}$ $\displaystyle=$
$\displaystyle\,\sum_{i=0}^{r-m-1}\gamma_{n-r+1,r-m,i}\sum_{j=0}^{s-r-1}\gamma_{n-s+1,s-r,j}\;\int_{0}^{\infty}g_{m+\alpha,\frac{a_{i}}{x+\beta},\frac{b_{j}}{x+\beta}}\big{(}y\big{)}\,dy_{2},$

which leads to (3.10) since the joint density
$g_{m+\alpha,\frac{a_{i}}{x+\beta},\frac{b_{j}}{x+\beta}}$ for $Y$ has
marginal $f_{m+\alpha,\frac{a_{i}}{x+\beta}}$ for $Y_{1}$, for all $i,j$
(Lemma 2.2) and since $\sum_{j=0}^{s-r-1}\gamma_{n-s+1,s-r,j}=1$. ∎

The result in itself is not new and was obtained with the univariate analysis
carried out by Dunsmore  (1974). A similar development establishes that (3.10)
holds for the marginal predictive density of $X_{n:n}-X_{m:n}$. For the
conditional distributions, we have the following.

######  Corollary 3.3.

For the Bayes predictive density (3.9):

  1. (a)

The conditional density of $Y_{2}$ given $Y_{1}=y_{1}$ is given by the mixture
representation

$\hat{q}_{\pi}(y_{2}|y_{1};x)\,=\,\sum_{i=0}^{r-m-1}\sum_{j=0}^{s-r-1}\alpha_{i}(y_{1})\,\beta_{j}\,f_{m+\alpha+1,\frac{b_{j}}{x\,+\,\beta\,+\,a_{i}\,y_{1}}}(y_{2}),$

with
$\alpha_{i}(y_{1})\propto\gamma_{n-r+1,r-m,i}\,f_{m+\alpha,\frac{a_{i}}{x+\beta}}(y_{1})$
such that $\sum_{i=0}^{r-m-1}\alpha_{i}(y_{1})=1$, and
$\beta_{j}=\gamma_{n-s+1,s-r,j}\,$;

  2. (b)

The conditional density of $Y_{1}$ given $Y_{2}=y_{2}$ is given by the mixture
representation

$\hat{q}_{\pi}(y_{1}|y_{2};x)\,=\,\sum_{i=0}^{r-m-1}\sum_{j=0}^{s-r-1}\xi_{j}(y_{2})\,\nu_{i}\,f_{m+\alpha+1,\frac{a_{i}}{x\,+\,\beta\,+\,\beta_{j}\,y_{2}}}(y_{1}),$

with
$\xi_{j}(y_{2})\propto\gamma_{n-s+1,s-r,j}\,f_{m+\alpha,\frac{b_{j}}{x+\beta}}(y_{2})$
such that $\sum_{j=0}^{s-r-1}\xi_{j}(y_{2})=1$, and
$\nu_{i}=\gamma_{n-r+1,r-m,i}$.

Proof. Part (a) follows directly with the properties of Lemma 2.2 upon writing
$g_{m+\alpha,\frac{a_{i}}{x+\beta},\frac{b_{j}}{x+\beta}}\big{(}y\big{)}\,=f_{m+\alpha,\frac{a_{i}}{x\,+\,\beta}}(y_{1})\,f_{m+\alpha+1,\frac{b_{j}}{x\,+\,\beta\,+\,a_{i}y_{1}}}(y_{2})$.
The same approach leads to (b). ∎

######  Remark 3.3.

With the above marginal and conditional distributions expressible as mixtures
of univariate Pareto distributions, corresponding moments are readily
available. For instance, one obtains

$\displaystyle\mu_{1}=\mathbb{E}_{\pi}(Y_{1}|x)$ $\displaystyle=$
$\displaystyle\frac{x+\beta}{m+\alpha-1}\,\sum_{i=0}^{r-m-1}\gamma_{n-r+1,r-m,i}\,\frac{1}{a_{i}}\,$
$\displaystyle\hbox{ and
}\mu_{2}(y_{1})\,=\,\mathbb{E}_{\pi}(Y_{2}|Y_{1};x)\,$ $\displaystyle=$
$\displaystyle\,\sum_{i=0}^{r-m-1}\sum_{j=0}^{s-r-1}\alpha_{i}(y_{1})\,\beta_{j}\,\frac{x+\beta+a_{i}y_{1}}{(m+\alpha)b_{j}}$
$\displaystyle=$
$\displaystyle\,y_{1}\,\sum_{i=0}^{r-m-1}\sum_{j=0}^{s-r-1}\alpha_{i}(y_{1})\,\beta_{j}\,\frac{a_{i}}{(m+\alpha)b_{j}}\,+\,\frac{x+\beta}{m+\alpha}\sum_{j=0}^{s-r-1}\,\,\frac{\beta_{j}}{b_{j}}.$

Observe that the conditional expectation $\mathbb{E}_{\pi}(Y_{2}|y_{1};x)$ is
affine linear of the form $C+Dy_{1}$.

###  3.2 HPD Bayesian prediction regions

A Bayesian prediction region with credibility $1-\lambda$ for $Z|\theta$ based
on an observed value $x$ of $X|\theta$ and for a given prior density $\pi$ is
such that

$\int_{R(x)}\hat{q}_{\pi}(z|x)\,dz\,=\,1-\lambda,$ (3.11)

with $\hat{q}_{\pi}$ given in (2.6). One such choice which minimizes volume is
the ubiquitous HPD region which is of the form

$R_{HPD}(x)\,=\,\big{\\{}z\in\mathbb{R}_{+}^{N}\,:\,\hat{q}_{\pi}(z|x)\,\geq
k\big{\\}},$ (3.12)

where $k$ is chosen so that (3.11) is satisfied. Here is an explicit form of
the HPD credible region for the future order statistics spacings
$Z_{1},\dots,Z_{N}$.

######  Theorem 3.3.

Based on $X_{i:n},i=1,\ldots,m$, the first $m$ order statistics among $n$ from
i.i.d. Exp$(\theta)$ data, setting $X$ as in (2.1) and $Z$ as in (2.3), the
HPD region of credibility $1-\lambda$ for $Z$ associated with Gamma prior
$\pi_{\alpha,\beta}$ for $\theta$ is given by

$R_{HPD}(x)\,=\,\big{\\{}z\in\mathbb{R}_{+}^{N}\,:\,\sum_{i=1}^{N}(n-m-i+1)\,z_{i}\,\leq
c_{0}(x+\beta)\big{\\}},$ (3.13)

with $c_{0}$ the quantile of order $1-\lambda$ of a $B2(N,m+\alpha)$
distribution.

Proof. It follows from Theorem 3.1 that $R_{HPD}$ defined generally in (3.12)
is of form (3.13). From Lemma 2.3, the posterior predictive distribution of
$W=\sum_{i=1}^{N}\frac{(n-m-i+1)\,z_{i}}{x+\beta}$ is $B2(N,m+\alpha)$
distributed and the result follows by setting $c_{0}$ such that
$\mathbb{P}(W\leq c_{0})=1-\lambda$. ∎

######  Remark 3.4.

With simple integral or finite sum forms for the $B2(N,m+\alpha)$ c.d.f., the
evaluation of the above quantiles is rather straightforward. For instance in
the bivariate case, we have $\mathbb{P}(W\leq
c)=1-\frac{1+c\,(m+\alpha+1)}{(1+c)^{m+\alpha+1}}$, so that $c_{0}$ is
expressible as the solution in $c>0$ of

$\frac{1+c\,(m+\alpha+1)}{(1+c)^{m+\alpha+1}}\,=\,\lambda.$

######  Remark 3.5.

HPD credible regions are in general not invariant with respect to
transformations, but they are whenever the transformation is affine linear.
Since the original order statistics $T=(X_{m+1:n},\ldots,X_{m+N;n})^{\top}$
are related to the spacings $Z=(Z_{1},\ldots,Z_{N})^{\top}$ by the affine
linear transformation $T=b+AZ$, with $b=(X_{m:n},\ldots,X_{m:n})^{\top}$ and
$A$ a lower triangular matrix with non-zero elements $a_{i,j}=1$ for $i\leq
j$, the transformation of $R_{HPD}(x)$ given by (3.13) to the order statistics
$X_{m:n},\ldots,X_{m+N;n}$ is also HPD.

###  3.3 Prediction regions: an algorithm for the bivariate case

We present here a general bivariate case solution to obtain a prediction
region with a given credibility based on Theorem 3.2’s predictive density for
$Y$ as in (2.2). There are several options available, but we opt for a rather
explicit strategy to construct a Bayesian credible region. It can be applied
for any choice of $(r,s,m,n)$, prior $\pi_{\alpha,\beta}$ such that
$m+\alpha>1$, and the targeted credibility $1-\lambda$. We refer to notation
used throughout this paper, namely $a_{i}$ and $b_{j}$ as in Theorem 3.2,
$\gamma_{c,d,k}$ as in Remark 3.2, $f_{l,h}$ and $\bar{F}_{l,h}$ as the
univariate Pareto density and survival functions in (2.7), $\beta_{j}$ and
$\alpha_{i}(y_{1})$ as in Corollary 3.3, and $\mu_{1}$ and $\mu_{2}(y_{1})$ as
in Remark 3.3.

  1. Step 1. Select a prediction region $A$ for $Y_{1}$ of credibility $\sqrt{1-\lambda}$ based on the marginal density in Corollary 3.2. Such a choice would desirably be an interval with relatively high levels of the predictive density. A suitable choice is:

$A\,=\,\big{[}\mu_{1}-\Delta_{1},\mu_{1}+\Delta_{1}\big{]}\,\cap[0,\infty)\,,$

with $\Delta_{1}>0$ uniquely chosen such that

$\displaystyle\int_{\mu_{1}-\Delta_{1}}^{\mu_{1}+\Delta_{1}}\hat{q}_{\pi_{\alpha,\beta}}(y_{1}|x)\,dy_{1}\,=\,\sqrt{1-\lambda}$
$\displaystyle\Longleftrightarrow$
$\displaystyle\sum_{i=0}^{r-m-1}\gamma_{n-r+1,r-m,i}\,\Big{(}\bar{F}_{m+\alpha,\frac{a_{i}}{x+\beta}}(\mu_{1}-\Delta_{1})_{+}\,-\,\bar{F}_{m+\alpha,\frac{a_{i}}{x+\beta}}(\mu_{1}+\Delta_{1})\Big{)}\,=\,\sqrt{1-\lambda},$

with $z_{+}=\max\\{0,z\\}$. Since this last expression is strictly increasing
in $\Delta_{1}$, it is rather straightforward to approach $\Delta_{1}$
numerically.

  2. Step 2. For each $y_{1}\in A$, select a prediction region $B(y_{1})$ of conditional credibility $\sqrt{1-\lambda}$ for $y_{2}$ based on the conditional density of $Y_{2}|Y_{1}=y_{1}$ given in part (a) of Corollary 3.3. The challenge here is similar to the one in Step 1 but to be repeated for all $y_{1}$. Analogously to Step 1, $B(y_{1})$ can reasonably be constructed pivoting around the mean as

$B(y_{1})\,=\,\big{[}\mu_{2}(y_{1})-\Delta_{2}(y_{1}),\mu_{2}(y_{1})+\Delta_{2}(y_{1})\big{]}\,\cap[0,\infty),$

with $\Delta_{2}(y_{1})>0$ chosen such that

$\displaystyle\\!\sum_{i=0}^{r-m-1}\sum_{j=0}^{s-r-1}\alpha_{i}(y_{1})\,\beta_{j}\,\int_{\mu_{2}(y_{1})-\Delta_{2}(y_{1})}^{\mu_{2}(y_{1})+\Delta_{2}(y_{1})}\,f_{m+\alpha+1,\frac{b_{j}}{x+\beta+a_{i}y_{1}}}(y_{2})\,dy_{2}\,=\,\sqrt{1-\lambda}$
$\displaystyle\Longleftrightarrow\\!$
$\displaystyle\sum_{i=0}^{r-m-1}\sum_{j=0}^{s-r-1}\alpha_{i}(y_{1})\,\beta_{j}\,\,\Big{(}\bar{F}_{m+\alpha+1,\frac{b_{j}}{x+\beta+a_{i}y_{1}}}(\mu_{2}(y_{1})-\Delta_{2}(y_{1}))_{+}\,-\,\bar{F}_{m+\alpha+1,\frac{b_{j}}{x+\beta+a_{i}y_{1}}}(\mu_{2}(y_{1})+\Delta_{2}(y_{1}))\Big{)}\,\,$
$\displaystyle=\sqrt{1-\lambda}.$

The resulting prediction region $R=\\{(y_{1},y_{2}):y_{1}\in A,y_{2}\in
B(y_{1})\\}$ has credibility $1-\lambda$ indeed since

$\displaystyle\mathbb{P}_{\pi}\big{(}Y_{1}\in A,Y_{2}\in B(y_{1})|x\big{)}$
$\displaystyle=$
$\displaystyle\int_{A}\hat{q}_{\pi}(y_{1}|x)\big{\\{}\int_{B(y_{1})}\hat{q}_{\pi}(y_{2}|y_{1};x)\,dy_{2}\big{\\}}\,dy_{1}$
$\displaystyle=$
$\displaystyle\int_{A}\hat{q}_{\pi}(y_{1}|x)\,(\sqrt{1-\lambda})\,dy_{1}$
$\displaystyle=$ $\displaystyle 1-\lambda.$

######  Remark 3.6.

Other regions can also be selected. For instance, one-sided choices with $A$
of the form $[0,\bar{\Delta}_{1}]$ or $[\underline{\Delta}_{1},\infty)$, and
with $B$ of the form $[0,\bar{\Delta}_{2}(y_{1})]$ or
$[\underline{\Delta}_{2}(y_{1}),\infty)$. Another alternative would be to aim
for different credibilities $1-\lambda_{1}$ and $1-\lambda_{2}$ in Steps 1 and
2, respectively, such that $(1-\lambda_{1})(1-\lambda_{2})=1-\lambda$.

###  3.4 Example

The following dataset from Murthy et al.  (2004) shows $n=30$ ordered failure
times for repairable items:

0.11 | 0.30 | 0.40 | 0.45 | 0.59 | 0.63 | 0.70 | 0.71 | 0.74 | 0.77
---|---|---|---|---|---|---|---|---|---
0.94 | 1.06 | 1.17 | 1.23 | 1.23 | 1.24 | 1.43 | 1.46 | 1.49 | 1.74
1.82 | 1.86 | 1.97 | 2.23 | 2.37 | 2.46 | 2.63 | 3.46 | 4.36 | 4.73

For the purpose of illustration, suppose that only the $m=20$ first values are
observed and that we wish to predict either:

  1. (i)

the next two order statistics spacings
$(Z_{1},Z_{2})\,=\,\big{(}X_{21:30}-X_{20:30},X_{22:30}-X_{21:30}\big{)}$,

or

  1. (ii)

the next and last order statistics $\big{(}X_{21:30},X_{30:30}\big{)}$.

In both cases, we use the non-informative prior density
$\pi_{(}\theta)\,=\,\frac{1}{\theta}\,\mathbb{I}_{(0,\infty)}(\theta)$, and
consider credibility $1-\lambda=0.95$. As expanded upon in the next section,
the frequentist coverage of such a prediction region matches the credibility
for all $\theta>0$. The data yields
$x\,=\,\sum_{i=1}^{20}x_{i:30}\,+\,10\,x_{20:30}\,=\,35.79$.

For case (i), Theorem 3.3 applies and the HPD credible region for
$(Z_{1},Z_{2})$ is given by

$\displaystyle
R_{\textrm{HPD}}(35.79)\,=\,\bigl{\\{}(z_{1},z_{2})\in\mathbb{R}_{+}^{2}\,:\,0.2794\,z_{1}\,+\,0.2515z_{2}\,\leq
0.2606\big{\\}},$

$c_{0}\,=\,0.2606$ being the quantile of order $0.95$ of a $B2(2,20)$
distribution.

For case (ii), we illustrate the use of Section 3.3’s algorithm which is
prescribed for $(Y_{1},Y_{2})$ with $Y_{1}\,=\,X_{21:30}-X_{20:30}$ and
$Y_{2}\,=X_{30:30}-X_{21:30}$. Figure 1 (a) presents the resulting Bayesian
prediction region along with the regression function or conditional
expectation $\mathbb{E}(Y_{2}|y_{1},x=35.79)$ as a function of $y_{1}$. The
first step yields $A=[0,0.722]$, while the $B(y_{1})$ intervals are shown in
the figure. For example, $B(0.50)=\,[0,13.059]$. Observe that a left
neighbourhood of $0$ is included in $B(y_{1})$ for all $y_{1}\in A$ as the
second step credibility, equal to $(0.95)^{1/2}\approx 0.9747$, is quite
large. In contrast, the Bayesian prediction region with credibility $0.80$ is
displayed in Figure 1 (b) resulting in $A=[0,0.4258]$ and intervals $B(y_{1})$
centered at $\mathbb{E}(Y_{2}|y_{1},x=35.79)$ which exclude values close to
$0$ for all $y_{1}\in A$. Finally, the corresponding prediction regions for
$\big{(}X_{21:30},X_{30:30}\big{)}\,$ and credibilities $0.95$ and $0.80$ are
obtained as $\\{(y_{1}+x_{20:30},y_{1}+y_{2}+x_{20:30}):(y_{1},y_{2})\in R\\}$
and displayed in Figure 2.

(a) $1-\lambda=0.95$ (b) $1-\lambda=0.80$

Figure 1: Bayesian prediction regions and $\mathbb{E}(Y_{2}|y_{1},x=35.79)$
(dashed)

(a) $1-\lambda=0.95$ (b) $1-\lambda=0.80$

Figure 2:  Bayesian prediction regions for $(X_{21:30},X_{30:30})$ in case
(ii)

##  4 Frequentist coverage and credibility

Through Theorem 3.2’s Bayesian predictive densities, the prediction regions
for future order statistics described in the previous section are constructed
in order to attain exact Bayesian credibility. This includes the non-
informative prior prediction densities $\hat{q}_{\pi_{0}}(\cdot;X)$ given by
Theorems 3.1 and 3.2 for $\alpha=\beta=0$. Moreover, prediction regions based
on $\hat{q}_{\pi_{0}}(\cdot;X)$ also lead to exact frequentist coverage as
expanded on below. In fact, we cast the result in a more general scale
parameter family model setting with scale parameter densities

$p_{\sigma}(x)=\frac{1}{\sigma}\,p_{1}(\frac{x}{\sigma})\hbox{ and
}q_{\sigma}(y^{\prime})=\frac{1}{\sigma^{d_{2}}}\,q_{1}(\frac{y_{1}^{\prime}}{\sigma},\ldots,\frac{y_{d_{2}}^{\prime}}{\sigma}).$
(4.14)

We will make use of the following intermediate result, a univariate version of
which was first given by L’Moudden et al.  (2017) (i.e., $d_{2}=1$).

######  Lemma 4.4.

Under model (4.14) and prior density
$\pi_{0}(\sigma)\,=\,\frac{1}{\sigma}\,\mathbb{I}_{(0,\infty)}(\sigma)$, the
Bayesian predictive density $\hat{q}_{\pi_{0}}$ is given by

$\hat{q}_{\pi_{0}}(y^{\prime}|x)\,=\,\frac{1}{x^{d_{2}}}\,h(\frac{y^{\prime}}{x}),$
(4.15)

where $h$ is the frequentist density of
$R=\frac{Y^{\prime}}{X}=(R_{1},\ldots,R_{d_{2}})$, which is free of $\sigma$,
and given by

$h(r)\,=\,\int_{0}^{\infty}u^{d/2}\,q_{1}(r_{1}u,\ldots,r_{d_{2}}u)\,p(u)\,du.$
(4.16)

Furthermore, it is the case that

$\big{(}\frac{Y_{1}^{\prime}}{x},\ldots,\frac{Y^{\prime}_{d_{2}}}{x}\big{)}\big{|}x\,=^{d}\,\big{(}\frac{Y_{1}^{\prime}}{X},\ldots,\frac{Y^{\prime}_{d_{2}}}{X}\big{)}\big{|}\sigma\,\hbox{
for all }x,\sigma,$ (4.17)

i.e., the posterior predictive and frequentist distributions of $R$ match and
are furthermore independent of the observed value of $x$ and the parameter
value of $\sigma$.

Proof. Identity (4.17) follows from the first part of the lemma. A direct
evaluation yields (4.16) as the density of $R|\sigma$. Finally, the posterior
density of $\sigma$ is given by
$\pi(\sigma|x)\,=\,\frac{x}{\sigma^{2}}\,p_{1}(\frac{x}{\sigma})$, which leads
to the predictive density

$\displaystyle\hat{q}_{\pi_{0}}(y^{\prime}|x)\,$ $\displaystyle=$
$\displaystyle\,\int_{0}^{\infty}\frac{1}{\sigma^{d_{2}}}\,q_{1}(\frac{y_{1}^{\prime}}{\sigma},\ldots,\frac{y_{d_{2}}^{\prime}}{\sigma})\,\frac{x}{\sigma^{2}}\,p_{1}(\frac{x}{\sigma})\,d\sigma$
$\displaystyle=$
$\displaystyle\frac{1}{x^{2}}\int_{0}^{\infty}u^{d_{2}}q_{1}\big{(}\frac{y_{1}^{\prime}u}{x},\ldots,\frac{y_{d_{2}}^{\prime}u}{x}\big{)}\,p(u)\,du$
$\displaystyle=$
$\displaystyle\frac{1}{x^{d_{2}}}\,h(\frac{y^{\prime}}{x})\,.\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\qed$

######  Remark 4.7.

Identity (4.17) is quite general and clarifies why Bayesian analysis with
respect to the prior density $\pi_{0}$ matches pivotal based analysis that
stems from the right-hand side. In the univariate case and for the prediction
of a single future order statistic, as remarked upon by Dunsmore (1974), the
identity explains leads to his Bayesian solutions matching the pivotal based
solution of Lawless (1971).

######  Theorem 4.4.

Consider model (4.14) and a Bayesian prediction region $R(X)$ with credibility
$1-\lambda$ associated with the prior density
$\pi_{0}(\theta)=\frac{1}{\theta}\,\mathbb{I}_{(0,\infty)}(\theta)$. Then,
$R(X)$ has exact frequentist coverage probability, i.e., $\mathbb{P}(R(X)\ni
Y^{\prime}|\theta)\,=\,1-\lambda$ for all $\theta>0$.

Proof. Since $R(X)$ has credibility $1-\lambda$, we have

$\mathbb{P}\big{(}Y^{\prime}\in
R(x))|x\big{)}\,=\,\mathbb{P}\big{(}\frac{Y^{\prime}}{x}\in
R^{*}(x)\big{|}x\big{)}\,=\,1-\lambda,$

where $R^{*}(x)\,=\,\\{r\in\mathbb{R}_{+}^{d_{2}}:\frac{r}{x}\in R(x)\\}$.
Now, since the predictive distribution of $\frac{Y^{\prime}}{x}$ is free of
$x$ with p.d.f. $h$ (Lemma 4.4), it follows that $R^{*}(x)$ is free of $x$
such that $\int_{R^{*}(x)}h(z)\,dz\,=\,1-\lambda$. On the other hand, the
frequentist coverage of $R(X)$ is given by

$\mathbb{P}(Y^{\prime}\in
R(X)|\,\theta\,)\,=\,\mathbb{P}\big{(}\frac{Y^{\prime}}{X}\in
R^{*}(X)\big{|}\theta\big{)}\,=\,\int_{R^{*}(x)}h(r)\,dr\,=\,1-\lambda,$

since $h$ is also the density of $\frac{Y^{\prime}}{X}|\theta$ (Lemma 4.4). ∎

To conclude, the above (with $Y^{\prime}=Y$ or $Y^{\prime}=Z$) applies to our
set-ups as follows.

######  Corollary 4.4.

Based on $X$ as in (2.1) and the non-informative prior density
$\pi_{0}(\theta)=\frac{1}{\theta}\,\mathbb{I}_{(0,\infty)}(\theta)$, Bayesian
predictive regions for $Z$ or $Y$ with credibility $1-\lambda$, based on the
predictive densities given in Theorems 3.1 and 3.2, have matching frequentist
coverage probability $1-\lambda$ for all $\theta>0$.

##  5 Best invariance and minimaxity

As seen in Section 3, Bayesian predictive densities such as those given in
Theorem 3.2 facilitate the construction of a prediction region for future
values of $Y$ or $Z$ with a given credibility. Non-Bayesian predictive
densities are also available, such as plug-in densities and those obtained by
likelihood or pivotal-based methods. Kullback-Leibler (KL) divergence loss and
accompanying risk can be used to evaluate the frequentist performance of
density estimators $\hat{q}(\cdot;X)$. For our problem, these are given by

$L_{KL}(\theta,\hat{q}(\cdot;x))\,=\,\int
q_{\theta}(t)\,\log\big{(}\frac{q_{\theta}(t)}{\hat{q}(t;x)}\big{)}\,dt\,$

(assuming Lebesgue densities), and

$R_{KL}(\theta,\hat{q})\,=\,\mathbb{E}_{\theta}\big{\\{}L_{KL}(\theta,\hat{q}(\cdot;X))\big{\\}}.$

It is of interest to assess the efficiency of predictive densities for KL
risk, and decision-theoretic properties of invariance and minimaxity
applicable to our contexts are reviewed in this section. An early contribution
to the determination of a best invariant density is due to Murray  (1977),
while a more exhaustive treatment of best invariant densities, as well as
minimaxity in predictive density estimation, appears in Liang & Barron
(2004).

A density $\hat{q}_{m}$ is minimax whenever $\hat{q}_{m}$ minimizes among all
densities the supremum frequentist risk, i.e., in our cases when

$\sup_{\theta>0}R_{KL}(\theta,\hat{q}_{m})\,=\,\inf_{\hat{q}}\sup_{\theta>0}R_{KL}(\theta,\hat{q})\,.$

For our problems, and more generally model (4.14) with $\sigma=1/\theta$,
Kullback-Leibler divergence loss, a predictive density $\hat{q}$ is invariant
under changes of scale whenever it satisfies the scale parameter family
property

$\hat{q}(y^{\prime};cx)\,=\,\frac{1}{c^{d_{2}}}\,\hat{q}(\frac{y^{\prime}}{c};x)\,,y^{\prime}\in\mathbb{R}_{+}^{d_{2}}\,,$
(5.18)

for all $c,x>0$. The class of invariant densities here includes
$\hat{q}_{\pi_{0}}$; as can be verified directly from the expressions given in
Theorems 3.1 and 3.2 with $\alpha=\beta=0$, or even by (4.15); as well as
plug-in densities $q_{\hat{\theta}}$ with $\hat{\sigma}(x)=kx$, i.e., a scale
invariant point estimator $\hat{\sigma}$ of $\theta$ satisfying
$\hat{\sigma}(cx)\,=\,c\,\hat{\theta}(x)$ for $c,x>0$; such as the maximum
likelihood choice with $c=1/m$.

The present invariance structure implies that an invariant density has
constant risk as a function of $\theta$ as long as it is finite, from which it
follows that there exists an optimal choice among invariant densities. The
risk constancy follows for vastly more general settings (e.g., Berger, 1985),
but can be also derived directly as follows:

$\displaystyle R_{KL}(c\sigma,\hat{q})\,$ $\displaystyle=$
$\displaystyle\int_{(0,\infty)}p_{c\sigma}(x)\,\int_{(0,\infty)^{d_{2}}}q_{c\sigma}(y^{\prime})\,\log\big{(}\frac{q_{c\sigma}(y^{\prime})}{\hat{q}(y;x)}\big{)}\,dy^{\prime}\,dx$
$\displaystyle=$
$\displaystyle\int_{(0,\infty)}\frac{1}{c}\,p_{\sigma}(\frac{x}{c})\,\int_{(0,\infty)^{d_{2}}}q_{\sigma}(\frac{y^{\prime}}{c})\frac{1}{c^{d_{2}}}\,\log\big{(}\,\frac{q_{\sigma}(\frac{y^{\prime}}{c})}{\hat{q}(\frac{y^{\prime}}{c};\frac{x}{c})}\big{)}\,dy^{\prime}\,dx$
$\displaystyle=$
$\displaystyle\int_{(0,\infty)}p_{\sigma}(t)\,\int_{(0,\infty)^{d_{2}}}q_{\sigma}(u)\,\log\big{(}\frac{q_{\sigma}(u)}{\hat{q}(u;t)}\big{)}\,du\,dt$
$\displaystyle=$ $\displaystyle\,R_{KL}(\sigma,\hat{q}),$

for $c,\sigma>0$, with the change of variables
$(t,u)\,=\,\frac{1}{c}(x,y^{\prime})$, and by making use of (5.18).

Under general conditions for problems that are invariant, which are met here,
a best invariant procedure exists and coincides with the generalized Bayes
estimator associated with a (right) invariant prior density (e.g., Berger,
1985, section 6.6.2). In out set-up, such a prior density is the non-
informative $\pi_{0}$, and it leads to the best invariant property of the
density $\hat{q}_{\pi_{0}}$. Furthermore, $\hat{q}_{\pi_{0}}$ is minimax (see
Liang & Barron  (2004), Theorem 2 and Proposition 3). We conclude this section
by summarizing the above as applicable (with $Y^{\prime}=Y$ or $Y^{\prime}=Z$)
to the problems at hand.

######  Theorem 5.5.

Consider $X_{i:n},i=1,\ldots,m$, the first $m$ order statistics among $n$ from
i.i.d. Exp$(\theta)$ data, and $X$, $Y$, and $Z$ as in (2.1), (2.2) and (2.3).
Then the best invariant densities for estimating the densities of $Y$ and $Z$,
respectively, under KL divergence loss are $\hat{q}_{\pi_{0}}$, as given in
Theorems 3.1 and 3.2 for $\alpha=\beta=0$. Furthermore, its KL risk is
constant as a function of $\theta$ and $\hat{q}_{\pi_{0}}$ is minimax.

##  6 Concluding remarks

We have illustrated the natural usage of Bayesian credible regions for the
prediction of future order statistics under a type-2 censoring scheme with
exponentially distributed data. We have emphasized multivariate aspects and
provided explicit expressions for Bayesian predictive densities and resulting
prediction credible regions. We have also addressed optimality properties
achieved with the non-informative prior density choice, such as the matching
of Bayesian credibility and frequentist probabilty coverage, as well as the
best invariant and minimax properties of the corresponding Bayesian predictive
density.

It would be interesting to extend the analysis more broadly to other types of
probability models and censoring schemes. We have provided a Bayesian HPD
credible region (Theorem 3.3) for the next $N$ future order statistics.
However, for the general bivariate case, such a solution is lacking. We do not
know for instance if the Bayesian predictive density is unimodal which would
facilitate the determination of such a region.

## Acknowledgements

The authors are thankful to the Université de Sherbrooke for its financial
support awarded as part of its visiting researcher’s program. Éric Marchand’s
research is supported in part by the NSERC of Canada.

## References

  * Arnold (2014) Arnold, B.C. (2014). Univariate and multivariate Pareto models. Journal of Statistical Distributions and Applications, 1, article 11.
  * Bagheri et al. (2022) Bagheri, S.F., Asgharzadeh, A., Fernández, A.J., & Pérez-Gonzàlez, C.P. (2022). Joint prediction for future failure times under type-II censoring. IEEE Transactions on Reliability, 71, 100–110.
  * Berger (1985) Berger, J. (1985). Statistical Decision Theory and Related Topics. Springer Texts in Statistics. Springer-Verlag, New York. Second edition.
  * Dunsmore (1974) Dunsmore, I.R. (1974). The Bayesian prediction problem in life testing models. Technometrics, 16, 455–4600.
  * Fourdrinier et al. (2019) Fourdrinier, D., Marchand, É. & Strawderman, W.E. (2019). On efficient prediction and predictive density estimation for spherically symmetric models. Journal of Multivariate Analysis, 173, 18–25.
  * Kotz et al. (2000) Kotz, S., Balakrishnan, N. & Johnson, N.L. (2000). Multivariate continuous distributions: volume 1, Wiley.
  * Lawless (1971) Lawless, J.F. (1971). A prediction problem concerning samples from the exponential distributions, with application in life testing. Technometrics, 13, 725–730.
  * Lehmann & Casella (1998) Lehmann, E.L. & Casella, G. (1998). Theory of Point Estimation. Second edition, Springer, New York.
  * Liang & Barron (2004) Liang, F. & Barron, A. (2004). Exact minimax strategies for predictive density estimation, data compression, and model selection. IEEE Transactions on Information Theory, 62, 2708–2726.
  * L’Moudden et al. (2017) L’Moudden, A., Marchand, É., Kortbi, O. & Strawderman, W.E. (2017). On predictive density estimation for Gamma models with parametric constraints. Journal of Statistical Planning and Inference, 185, 56–68.
  * Murray (1977) Murray, G.D. (1977). A note on the estimation of probability density functions. Biometrika, 64, 150–152.
  * Murthy et al. (2004) Murthy, D.P., Xie, M. & Jiang, R. (2004). Weibull models. John Wiley & Sons.