Multivariate Estimation of Poisson Parameters

TL;DR

This paper introduces a new loss function for multivariate Poisson parameter estimation, develops minimax estimators that outperform MLE, and explores Bayesian methods, leading to improved estimation and new multivariate models.

Contribution

It proposes a novel loss function and a class of minimax estimators that dominate MLE, with Bayesian derivations and applications to multivariate Poisson models.

Findings

Minimax estimators outperform MLE under the new loss function.

Bayesian and empirical Bayesian methods derive the estimators.

New multivariate models for rate parameters and counts are constructed.

Abstract

This paper is devoted to the multivariate estimation of a vector of Poisson means. A novel loss function that penalises bad estimates of each of the parameters and the sum (or equivalently the mean) of the parameters is introduced. Under this loss function, a class of minimax estimators that uniformly dominate the maximum likelihood estimator is derived. Crucially, these methods have the property that for estimating a given component parameter, the full data vector is utilised. Estimators in this class can be fine-tuned to limit shrinkage away from the maximum likelihood estimator, thereby avoiding implausible estimates of the sum of the parameters. Further light is shed on this new class of estimators by showing that it can be derived by Bayesian and empirical Bayesian methods. In particular, we exhibit a generalisation of the Clevenson-Zidek estimator, and prove its admissibility. Moreover, a class of prior distributions for which the Bayes estimators uniformly dominate the maximum likelihood estimator under the new loss function is derived. A section is included involving weighted loss functions, notably also leading to a procedure improving uniformly on the maximum likelihood method in an infinite-dimensional setup. Importantly, some of our methods lead to constructions of new multivariate models for both rate parameters and count observations. Finally, estimators that shrink the usual estimators towards a data based point in the parameter space are derived and compared.

Figures

Equations

$$L_0^{\ast }(\theta ,\delta )=\sum _{i=1}^p({\delta }_i-{\theta }_i{)}^2\mathrm{and}{L}_1^{\ast }(\theta ,\delta )=\sum _{i=1}^p({\delta }_i-{\theta }_i{)}^2/{\theta }_i,$$

$$\displaystyle\delta_{\text{CZ},i}(Y)=\Bigl{(}1-\frac{p-1}{p-1+Z}\Bigr{)}Y_{i},\quad{\rm where\ }Z=\sum_{i=1}^{p}Y_{i},$$

$$\displaystyle L(\theta,\delta)=\sum_{i=1}^{p}(\delta_{i}-\theta_{i})^{2}+c\Bigl{(}\sum_{i=1}^{p}\delta_{i}-\gamma\Bigr{)}^{2}=L_{0}^{*}(\theta,\delta)+cp^{2}(\bar{\delta}-\bar{\theta})^{2},$$

$$\begin{array}{r}{L}_c(\theta ,\delta )=\sum _{i=1}^p\frac{({\delta }_i-{\theta }_i{)}^2}{{\theta }_i}+c\frac{({\sum }_{i=1}^p{\delta }_i-{\sum }_{i=1}^p{\theta }_i{)}^2}{{\sum }_{i=1}^p{\theta }_i}=L_1^{\ast }(\theta ,\delta )+cp\frac{(\bar{\delta }-\bar{\theta }{)}^2}{\bar{\theta }},\end{array}$$

$$\displaystyle\delta_{i}^{c}(Y)=\Bigl{\{}1-\frac{p-1}{p-1+(1+c)Z}\Bigr{\}}\,Y_{i}\quad{\rm for\ }i=1,\ldots,p,$$

$$L(\delta ,\theta )=(\delta -\theta {)}^{\mathrm{t}}A(\delta -\theta ),$$

$$\begin{array}{rl}R(\delta ,\theta )-R({\delta }_0,\theta )& ={\mathrm{E}}_{\theta }\{(\delta -\theta {)}^{\mathrm{t}}A(\delta -\theta )-(Y-\theta {)}^{\mathrm{t}}A(Y-\theta )\}\\ & ={\mathrm{E}}_{\theta }\{-2(Y-\theta )A\varphi (Y)+\varphi (Y{)}^{\mathrm{t}}A\varphi (Y)\}\\ & ={\mathrm{E}}_{\theta }\{-2(Y-\theta )\psi (Y)+\psi (Y{)}^{\mathrm{t}}{A}^{-1}\psi (Y)\},\end{array}$$

$${\mathrm{E}}_{\theta }(Y_i-{\theta }_i)\psi (Y)={\mathrm{E}}_{\theta }{Y}_i\{\psi (Y)-\psi (Y-e_i)\},$$

$$D(y,\varphi )=-\sum _{i=1}^{p}2y_i\{\psi (y)-\psi (y-e_i)\}+\psi (y{)}^{\mathrm{t}}{A}^{-1}\psi (y).$$

$${\psi }_i(y)=\frac{d(y)}{B(y)}T(y_i),\text{for }i=1,\dots p,$$

$$0\le d(y)\le (2/M)\{N(y)-2{\}}_{+},$$

$$\displaystyle D(\phi,y)\leq-2\sum_{i=1}^{p}y_{i}\big{\{}\psi_{i}(y)-\psi_{i}(y-e_{i})\big{\}}+\sum_{i=1}^{p}M\psi_{i}(y)^{2}.$$

$$\displaystyle D(\phi,y)\leq-2d(y)\frac{\big{\{}N(y)-2-Md(y)/2\big{\}}_{+}}{B(y)},$$

$$d_0(y)=(1/M)\{N(y)-2{\}}_{+}.$$

$$\delta (y)=y-A^{-1}{\psi }_0(y),\text{where}{\psi }_{0,i}(y)=\frac{1}{M}\frac{T(y_i)}{B(x)}\{N(y)-2{\}}_{+},$$

$$R(\delta,\theta)-R(Y,\theta)\leq-M^{-1}{\rm E}_{\theta}\,\big{\{}(N(Y)-2)_{+}\big{\}}^{2}/B(Y)<0,$$

$${\varphi }_0(y)={\psi }_{i,0}(y)-\frac{c}{1+cp}\sum _{j=1}^p{\psi }_{0,j}(y)=\{N(y)-2{\}}_{+}\frac{T(y_i)-\frac{cp}{1+cp}\bar{T}(y)}{B(y)},$$

$$\delta_{i}(y)\doteq y_{i}-{p-2\over\sum_{j=1}^{p}(\log y_{j})^{2}}\Bigl{\{}\log y_{i}-{cp\over 1+cp}{1\over p}\sum_{j=1}^{p}\log y_{j}\Bigr{\}}.$$

$${\delta }_1^{\ast }(y)$$

$${\delta }_i^{\ast }(y)$$

$${\delta }_p^{\ast }(y)$$

$${\psi }_{0,i}(y)=\frac{1}{4}\frac{T(y_i)}{B(y)}\{N(y)-2{\}}_{+}.$$

$${\hat{\lambda }}_i=\sum _{j=1}^i{Y}_j-{\psi }_{0,i}(Y)+{\psi }_{0,i+1}(Y)\mathrm{for}i\le p-1,$$

$$a_i=\{\mathrm{E}({\theta }_i^{-1}|y){\}}^{-1},a=\sum _{i=1}^p{a}_i,\text{and}b=\{\mathrm{E}({\gamma }^{-1}|y){\}}^{-1}.$$

$$\displaystyle{\rm E}\,\{L_{c}(\theta,\delta)\,|\,y\}=\sum_{i=1}^{p}\big{\{}\delta_{i}^{2}/a_{i}-2\delta_{i}+{\rm E}(\theta_{i}\,|\,y)\big{\}}+c\Bigl{\{}\Bigl{(}\sum_{i=1}^{p}\delta_{i}\Bigr{)}^{2}/b-2\sum_{i=1}^{p}\delta_{i}+{\rm E}(\gamma\,|\,y)\Bigr{\}},$$

$${\delta }_i^B(y)=\frac{1+c}{1+ca/b}{a}_i=\frac{1+c}{1+ca/b}\frac{1}{\mathrm{E}\{(1/{\theta }_i)|y\}}\mathrm{for}i=1,\dots ,p.$$

$$\displaystyle\begin{array}[]{rcl}\delta_{i}^{B}&=&\displaystyle{1+c\over 1+c(p\bar{\alpha}+z-p)/(p\bar{\alpha}+z-1)}{\alpha_{i}+y_{i}-1\over\beta+1}\\ &=&\displaystyle{(1+c)(p\bar{\alpha}+z-1)\over(1+c)(p\bar{\alpha}+z-1)-c(p-1)}{\alpha_{i}+y_{i}-1\over\beta+1}\end{array}$$

$${\delta }_i^B\doteq \frac{p\bar{\alpha }+z-1}{\beta +1}\frac{{\alpha }_i+y_i-1}{p\bar{\alpha }+z-p},\mathrm{with}\mathrm{sum}\frac{p\bar{\alpha }+z-1}{\beta +1},$$

$${\delta }_i^B=\frac{1+c}{1+cz/(p-1+z)}\frac{y_i}{1+\beta }=\frac{(1+c)(p-1+z)}{p-1+(1+c)z}\frac{y_i}{1+\beta }.$$

$$\displaystyle\begin{split}R_{c}(\delta_{i}^{B},\theta)&={\rm E}_{\theta}\,{\rm E}_{\theta}\,\{L(\theta,\delta_{i}^{B})\,|\,Z\}\\ &={\rm E}_{\gamma}\Bigl{\{}\frac{(1+c)^{2}(p-1+Z)^{2}}{p-1+(1+c)Z}\frac{Z}{\gamma\,(1+\beta)^{2}}\\ &\hskip 56.9055pt-2\frac{(1+c)^{2}(p-1+Z)}{p-1+(1+c)Z}\frac{Z}{1+\beta}+(1+c)\gamma\Bigr{\}}.\end{split}$$

$$\text{MBR}(1,\beta )={\mathrm{E}}_{\gamma }{\mathrm{E}}_{\gamma }\{L(\theta ,{\delta }^B)|Z\}=\frac{1+c}{1+\beta }{\mathrm{E}}_{\gamma }\frac{p(p-1)+(p+c)Z}{p-1+(1+c)Z},$$

Full text

Key words: admissibility; Bayes and empirical Bayes; minimax; Poisson; shrinkage; multivariate estimation

1   Introduction

Let $Y=(Y_{1},\ldots,Y_{p})$ be a vector of independent Poisson random variables with mean vector $\theta=(\theta_{1},\ldots,\theta_{p})$. In this article we consider the problem of estimating the vector $\theta$. The obvious estimator is $\delta_{0}(Y)=Y$, that is, using $\delta_{0,i}=Y_{i}$ for each of the $p$ individual parameters. It is well known that $\delta_{0}$ is the maximum likelihood solution, that it has components with uniformly minimal variance among the unbiased estimators, and that it is admissible in the one-dimensional problem under squared error loss as well as under its weighted version, see e.g. Lehmann (1983, p. 277). In the simultaneous or multivariate problem, however, Peng (1975) and Clevenson & Zidek (1975) were the first to show that $\delta_{0}$ can be improved upon under the loss functions

[TABLE]

if $p\geq 3$ and $p\geq 2$, respectively. In particular, for the $L_{1}^{*}$ loss, Clevenson & Zidek (1975) derived the estimator

[TABLE]

demonstrating that it improves uniformly on the standard procedure $\delta_{0}$. This Stein-type phenomenon has also been observed for other loss functions. Hwang (1982) obtained results for $\sum_{i=1}^{p}(\delta_{i}-\theta)^{2}/\theta_{i}^{m_{i}}$, for fixed integers $m_{i}$. Ghosh & Yang (1988) considered a loss function based on the entropy or Kullback–Leibler distance $L_{e}(\theta,\delta)=\sum_{i=1}^{p}\theta_{i}\{\delta_{i}/\theta_{i}-\log(\delta_{i}/\theta_{i})-1\}$, A good overview article for admissibility issues, for multivariate Poisson means and for other models for discrete data, is Ghosh, Hwang & Tsui (1983), followed by discussion contributions by Berger (1983), Morris (1983), Hudson (1983). A more recent article on these issues is Brown, Greenshtein & Ritov (2013).

More broadly, the books Efron (2013) and Fourdrinier, Strawderman & Wells (2018) contribute to seeing multivariate estimation, using shrinkage and empirical Bayes methods, as as a global phenomenon. Central themes are variations on ‘borrowing strength’, either via direct constructions or in empirical Bayes setups. Our article is a contribution in these general directions, showing that methods developed for the multinormal and more generally spherically symmetric distributions have certain parallels in the world of multivariate Poisson estimation. Importantly, when estimating a particular Poisson parameter $\theta_{i}$, the full multivariate vector of data is being used, as an integral part of our methods.

The estimators that have been found to be better than $\delta_{0}(Y)=Y$ in these earlier publications are essentially all of the shrinking type, pushing the maximum likelihood estimator closer to the bottom corner of the parameter space. A good example of the merits of this type of shrinkage is provided by Clevenson & Zidek (1975), wherein 36 small $\theta_{i}$ – intensities of an oil-well discovery process – are estimated with the estimator in (1.2). Clevenson & Zidek (1975) had access to ‘known’ $\theta_{i}$ and could check the actual loss incurred. The procedure in (1.2) did indeed give much smaller loss $\sum_{i=1}^{36}(\delta_{\text{CZ},i}-\theta_{i})^{2}/\theta_{i}$ than did $\delta_{0}$; the losses are 14.33 and 39.26, respectively. However, another and perhaps not so pleasant feature of their estimation procedure is conspicuous here, namely that the sum $\gamma=\sum_{i=1}^{36}\theta_{i}$ (or equivalently the mean $\bar{\theta}=\gamma/36$) is seriously underestimated. In their oil-well discovery example the true sum is $\gamma=26.98$, the usual $\delta_{0}$ has $\sum_{i=1}^{36}y_{i}=29$, while their $\sum_{i=1}^{36}\delta_{\text{CZ},i}=12.97$ is much too low.

In some situations the sum (or mean) is unimportant and all that matters is to estimate each $\theta_{i}$. In many multiparameter cases on wants to keep track of the sum (or the mean) of the $\theta_{i}$ as well, however, as one surely would in the oil-well discovery process above. Other multiparameter cases where the sum (or the mean) are deemed as important as the individual $\theta_{i}$ abound; think for example of a decision maker having to make budgetary decisions concerning each of the boroughs of a city and the city as a whole, the resources allocated to each of $p$ hospitals and the whole health budget, etc.

These considerations motivate studying loss functions that take into account the need for good individual estimates, while at the same time guarding against the underestimation of the sum. One example of such a loss function is

[TABLE]

Since $\sum_{i=1}^{p}\delta_{0,i}=\sum_{i=1}^{p}Y_{i}$ is the admissible minimax solution under the loss function \big{(}\sum_{i=1}^{p}\delta_{i}-\gamma\big{)}^{2} one might wonder whether the extra penalisation above would secure admissibility of $\delta_{0}$ in the simultaneous problem. In Section 2 we show that this is not the case; $\delta_{0}$ is again inadmissible when $p\geq 3$, under any given quadratic form loss function $(\delta-\theta)^{\rm t}A(\delta-\theta)$. Another loss function that takes the guarding against underestimation of the sum (or the mean) of the parameters into account is the weighted version of (1.3), that is, the loss functions that generalises the one used by Clevenson & Zidek (1975), namely

[TABLE]

where $c$ is a user-defined constant. Under (1.4), the loss in Clevenson and Zidek’s example for $\delta_{0}$ is $39.26+0.35\,c$ while the loss for $\delta_{\text{CZ}}$ is $14.33+6.35\,c$. Again, one might wonder whether $\delta_{0}$ is admissible for certain values of $c$. In Section 3 we show that this is not the case, and derive an estimator which is the natural generalisation of the one in (1.2),

[TABLE]

which is shown to belong to a class of estimators that dominate $\delta_{0}$ uniformly over the parameter space. In the somewhat similar multivariate normal setting, investigations on how to limit the shrinkage of the James–Stein estimator, in order to account for objectives related to estimation the means, have been conducted by Efron & Morris (1971, 1972).

Our paper proceeds as follows. In Section 2 we study the loss function (1.3) and derive a class of estimators that uniformly dominate the maximum likelihood estimator. Section 3 concerns the weighted version of this loss function, the $L_{c}$ loss function of (1.4), and we derive the already mentioned generalisation $\delta^{c}$ of $\delta_{\rm CZ}$. In Section 4 further light is shed on this new class of estimators by showing that it can be derived by Bayesian and empirical Bayesian methods. We are also able to prove that the estimator $\delta^{c}$ of (1.5) is admissible. Classes of alternative estimators are then studied in Sections 5 and 6, involving Bayes and empirical Bayes strategies which shrink the raw data towards data-driven structures for the $\theta_{i}$, such as the data mean, or a submodel. Some of these procedures succeed in having risks well below the minimax benchmark, in large regions of the parameter space, though without achieving uniform dominance. In Section 7 weighted loss functions are studied, which admit relative weights of importance; notably, an infinite-dimensional setup is also included. Then in Section 8 we demonstrate how some of our Bayesian constructions also lead to new multivariate models for rate parameters and for count observations, of interest in their own right, pointing to models for spatially dependent count data. Finally Section 9 offers a list of concluding remarks.

2   General quadratic loss function

We may write the loss function of (1.3) as $L(\delta,\theta)=(\delta-\theta)^{\rm t}A_{0}(\delta-\theta)$, where $A_{0}$ is the matrix with $1+c$ down its diagonal and $c$ elsewhere. The natural generalisation is

[TABLE]

where $A$ is symmetric and positive definite. Below we obtain some results for general $A$ and apply these to two examples. Note that in situations where there is no ordering of the individual $\theta_{i}$ and no reason to estimate some of them more precisely than the others, the loss function in (1.3), that is, using $A_{0}$ with an appropriate choice of $c$, is the natural choice. Our method of proving inadmissibility of $\delta_{0}$ resembles that of Tsui & Press (1982) and Hwang (1982), where $A$ is diagonal.

Let $\delta(y)=y-\phi(y)$ be a competitor. The difference in risk between these two estimators is then

[TABLE]

writing $\psi(y)=A\phi(y)$. Furthermore, since ${\rm E}_{\theta}\,\theta_{i}g(Y)={\rm E}_{\theta}\,Y_{i}g(Y-e_{i})$ for any $g$ with ${\rm E}_{\theta}\,|g(Y)|<\infty$, using $e_{i}$ to denote the unit vector with the $i$th element equal to one, we have that

[TABLE]

and the risk difference (2.2) can be written ${\rm E}_{\theta}\,D(Y,\phi)$, in which

[TABLE]

If a function $\psi(y)$ can be found such that $D(y,\phi)$ is non-positive for all $y$, with strict inequality for at least one datum $y$, then $\delta_{0}=Y$ is inadmissible, being outperformed by $\delta(y)=y-A^{-1}\psi(y)$.

Theorem 1**.**

Let $A$ be symmetric and positive definite. Then $\delta_{0}=Y$ is inadmissible under loss function (2.1), if $p\geq 3$. It is dominated by $\delta(Y)=Y-A^{-1}\psi(Y)$, where $\psi(\cdot)$ is any member of the following class

[TABLE]

where $T(0)=0$ and $T(y)=\sum_{j\leq y}1/j$ for $y\geq 1$, where $B(y)=\sum_{i=1}^{p}T(y_{i})T(y_{i}+1)$, and $d(y)$ is nondecreasing in each argument and obeying

[TABLE]

writing $a_{+}=\max(a,0)$ for truncating to zero. Here $N(y)=\sum_{i=1}^{p}I\{y_{i}\leq 1\}$ and $M$ is the inverse of the smallest eigenvalue of $A$.

Proof.

The choice of $M$ entails $\psi^{\rm t}A^{-1}\psi\leq M\sum_{i=1}^{p}\psi_{i}^{2}$ so that

[TABLE]

The general Theorem 2.1 in Hwang (1982), with accompanying corollaries, can then be used to find $\psi_{i}(\cdot)$ functions that make $D(\phi,y)$ non-positive for all $y$. We skip details but record that Hwang’s method gives

[TABLE]

which is non-positive for each $\psi(\cdot)$ described in the theorem. ∎

A natural choice for $d(y)$ is the following, minimising the upper bound in (2.3),

[TABLE]

This means using $\phi_{0}(y)=A^{-1}\psi_{0}(y)$ with $D(\phi_{0},y)\leq-M^{-1}\{N(y)-2\}^{2}/B(y)$, which shows that

[TABLE]

achieves

[TABLE]

for all $\theta$. Note in particular that the same $\psi_{0}(\cdot)$ function works for a large class of loss function (2.1).

Example 2.1**.**

Let $A$ be the square matrix with $1+c$ down its diagonal and $c$ elsewhere. This matrix might be written $A=I+cee^{\rm t}$, with $I$ the identity matrix and $e=(1,\ldots,1)^{\rm t}$ the $p\times 1$ vector of ones. This choice of $A$ gives the loss function (1.3). Here $M=1$ and $A^{-1}=I-\{c/(1+cp)\}ee^{\rm t}$, and $\psi^{\rm t}A^{-1}\psi\leq\sum_{i=1}^{p}\psi_{i}^{2}$ follows. The natural estimator is then $\delta(y)=y-\phi_{0}(y)$, where

[TABLE]

in which $\bar{T}(y)=p^{-1}\sum_{i=1}^{p}T(y_{i})$. Note that $\delta$ does not shrink the $y$ in any particular direction, but rather pushes the components $y_{i}$ in different directions according to the sign of $T(y_{i})-\{cp/(1+cp)\}\bar{T}(y)$. Note further that if the $y_{i}$ are moderate or large, then

[TABLE]

Example 2.2**.**

Samples of independent Poisson variables arise naturally when one or more Poisson processes are observed over time. Dividing the time interval into $p$ parts gives counts $y_{1},\ldots,y_{p}$ with certain means $\theta_{1},\ldots,\theta_{p}$. In the nonparametric setting, where the intensity of the process is unknown, these parameters are also completely unknown. If one wishes to estimate the cumulative intensity of the process, then $(\lambda_{1},\ldots,\lambda_{p})$ are more important than $(\theta_{1},\ldots,\theta_{p})$, where $\lambda_{i}=\sum_{j=1}^{i}\theta_{j}$. This suggests using the loss function $L(\theta,\delta)=\sum_{i=1}^{p}(\widehat{\lambda}_{i}-\lambda_{i})^{2}$, where $\widehat{\lambda}_{i}=\sum_{j=1}^{i}\delta_{j}$. But this is seen to be a special case of (2.1), with elements $a_{i,j}=p+1-\max(i,j)$ filling the $A$ matrix. Its inverse $A^{-1}$ has first row $(1,-1,0,\ldots,0)$, last row $(0,\ldots,0,-1,2)$, and in between we find $(0,\ldots,0,-1,2,-1,0,\ldots,0)$. One has $\psi^{\rm t}A^{-1}\psi\leq 4\psi^{\rm t}\psi$, and can use $M=4$ when applying the theorem. Hence the following estimator improves on $\delta_{0}=Y$:

[TABLE]

where this time

[TABLE]

Notice finally that the corresponding improved estimators for the cumulative $\lambda_{i}$ become

[TABLE]

while $\widehat{\lambda}_{p}=\sum_{j=1}^{p}Y_{j}-\psi_{0,p}(Y)$.**

3   The $L_{c}$ loss function

The main consideration leading to the loss function $(\delta_{i}-\theta_{i})^{2}/\theta_{i}$ is that the statistician seeks precise estimates of small values of $\theta_{i}$. A loss function that penalises heavily for bad estimates of small parameters is then a natural choice. Related to this is the obvious fact that when the parameters are small, they can only be badly overestimated, zero being the boundary of the parameter space. The corresponding multiparameter version of this is $\sum_{i=1}^{p}(\delta_{i}-\theta_{i})^{2}/\theta_{i}$, and is the one treated by Ghosh, Hwang & Tsui (1983), Hwang (1982), Tsui & Press (1982), Clevenson & Zidek (1975), and others. Note that $\delta_{0}=Y$ has constant risk $p$ with this loss, and it is not difficult to establish that it is minimax.

The above mentioned authors obtain classes of estimators that perform uniformly better than $\delta_{0}$ if only $p\geq 2$. As discussed in Section 1 and illustrated by the oil-well example, these shrinkage estimators do not take into account the additional desideratum, namely a precise estimate of the sum $\gamma=\sum_{i=1}^{p}\theta_{i}$. We will now study the $L_{c}$ loss function of (1.4), $L_{c}(\theta,\delta)=L_{1}^{*}(\theta,\delta)+c(\sum_{i=1}^{p}\delta_{i}-\gamma)^{2}/\gamma$. If we consider the second term in (1.4) by itself, we recognise the one-dimensional loss function $(\delta-\gamma)^{2}/\gamma$. It is well known that $\widehat{\gamma}=Z$ is admissible and the unique minimax solution under this loss function, and can therefore not be uniformly improved upon (Lehmann, 1983, p. 277). Consequently, since $z$ is the sum of the $y_{i}$, higher values of $c$ will result in estimators that lie closer the $\delta_{0}$. On the other hand, we know that for $c=0$ the estimator in (1.2) uniformly dominates $\delta_{0}$. Hence, the user defined constant $c$ determines how to compromise between $\delta_{0}$ and $\delta_{\text{CZ}}$.

Before we derive a class of estimators that dominate $\delta_{0}$ under $L_{c}$ in Section 3.2, we derive formulae for the Bayes solution and show that $\delta_{0}$ is minimax.

3.1   The $\delta_{0}=Y$ estimator is minimax.

The maximum likelihood estimator $\delta_{0}=Y$ has constant risk $p+c$ under $L_{c}$ and is a natural candidate for being minimax. We demonstrate minimaxity by exhibiting a sequence of priors with minimum Bayes risks $\text{BR}(\delta,\pi)={\rm E}\,R_{c}(\delta,\theta)$ which converge towards $p+c$; that this is sufficient follows from well-known arguments, as exposited e.g. in Robert (2007, Ch. 2.4). Some analysis is required to characterise the Bayes solution. We first find the values $\delta_{1},\ldots,\delta_{p}$ that minimise the posterior expected loss, i.e. given a dataset $y=(y_{1},\ldots,y_{p})$, with respect to some distribution over the parameter space. With $\gamma=\sum_{i=1}^{p}\theta_{i}$, introduce

[TABLE]

Then

[TABLE]

assuming the moments to exist. Some analysis shows that the minimum takes place for

[TABLE]

This generalises the familiar result that the Bayes solution is $\{{\rm E}(\theta_{i}^{-1}\,|\,Y)\}^{-1}=a_{i}$ under $L_{1}^{*}(\theta,\delta)$, that is, when $c=0$. Note also that if ${\rm E}\,\theta_{i}^{-1}=\infty$ for some $i$, and ${\rm E}\,\theta_{i}$ is finite, then only $\delta_{i}=0$ gives a finite risk, which means that (3.1) is correct even in such cases.

To illustrate this, suppose $\theta_{i}$ has a Gamma prior with parameters $(\alpha_{i},\beta)$, which we write as ${\rm Gamma}(\alpha_{i},\beta)$, i.e. with prior mean $\alpha_{i}/\beta$, and that these are independent. Then $\theta_{i}\,|\,y\sim{\rm Gamma}(\alpha_{i}+y_{i},\beta+1)$, and some calculations lead to the Bayes estimators

[TABLE]

for $i=1,\ldots,p$, writing $\bar{\alpha}=(1/p)\sum_{i=1}^{p}\alpha_{i}$ and again $z=\sum_{i=1}^{p}y_{i}$. If $\alpha_{i}\leq 1$, then the Bayes estimate is zero if $y_{i}=0$, by the comment above about (3.1). We note that for large $c$, corresponding to the loss being essentially related to estimating the sum $\gamma=\sum_{i=1}^{p}\theta_{i}$ well under $(\widehat{\gamma}-\gamma)^{2}/\gamma$ loss,

[TABLE]

tying in with $\gamma\,|\,y$ being a Gamma with parameters $(p\bar{\alpha}+z,\beta+1)$.

Now consider the special case where the $\theta_{i}$ are independent Gammas with parameters $(1,\beta)$. The Bayes solution then takes the form

[TABLE]

It now remains to show that the minimum Bayes risk for this prior, say $\text{MBR}(1,\beta)$, tends to $p+c$ as $\beta\to 0$. Using that $Y_{i}$ given $Z$ is binomial with mean $Z\theta_{i}/\gamma$, provided $Z\geq 1$, and that the $Y$ vector and hence the $\delta^{B}$ estimator are equal to zero when $Z=0$, the risk of $\delta_{i}^{B}$ can be expressed as

[TABLE]

Since the risk depends on the $\theta_{i}$ only through the sum $\gamma$, the minimum Bayes risk may be written

[TABLE]

in which the expectation on the right is with respect to the marginal distribution of $Z$. Since $Z$ tends in probability to infinity as $\beta\to 0$, the function above converges in probability,

[TABLE]

when $\beta\to 0$. Furthermore, this function is bounded by $p$, so by the bounded convergence theorem $\text{MBR}(1,\beta)$ tends to $p+c$ as $\beta$ goes to zero, as was to be shown.

3.2   A dominating class of estimators.

We will now develop a class of estimators with uniformly smaller risk than $p+c$ under the $L_{c}$ loss function, that is, estimators that uniformly dominate the maximum likelihood estimator. Consider estimators of the form $\delta_{i}(Y)=\{1-\phi(Z)\}Y_{i}$. Write $D(\phi,y)=L_{c}(\theta,\delta)-L_{c}(\theta,\delta_{0})$ for the difference in loss. Then

[TABLE]

Now, use the fact that for any real valued function $h$ with finite mean ${\rm E}_{\theta}\,h(Y)$, and with the property that $h(y)=0$ whenever $y_{i}=0$, the following identity holds:

[TABLE]

Using this identity we obtain an expression for the difference in risk ${\rm E}_{\theta}\,D(\phi,Y)=R(\delta^{*},\theta)-R(Y,\theta)$, namely

[TABLE]

This can hence be expressed as ${\rm E}_{\gamma}\,D^{*}(\phi,Z)$, with the $D^{*}(\phi,z)$ function not depending on the parameters; in particular, the risk function $R_{c}(\delta,\theta)=p+c+{\rm E}_{\gamma}\,D^{*}(\phi,Z)$ depends on the parameter vector only via $\gamma=\sum_{i=1}^{p}\theta_{i}$. Also, any function $\phi(\cdot)$ that ensures that $D^{*}(\phi,z)\leq 0$ for all $z$, with strict inequality for at least one datum $z$, yields an estimator that uniformly dominates the $\delta_{0}$. This leads to the following result.

Theorem 2**.**

For each function $\phi(\cdot)$ such that

[TABLE]

for all $z$, the estimator $\delta_{c}=\{1-\phi(Z)\}Y$ uniformly dominates $\delta_{0}=Y$. These conditions are met for functions of the type $\phi(z)=\psi(z)/\{p-1+(1+c)z\}$ where $\psi(\cdot)$ is nondecreasing with $0<\psi(z)<2(p-1)$. In particular, $\delta_{0}$ is inadmissible if $p\geq 2$.

Proof.

Using the expression for $D^{*}(\phi,Z)$ derived above, the following holds:

[TABLE]

This is valid for all $\gamma$ since $\phi(z)z$ is a strictly increasing function of $z$. ∎

We denote by $\mathcal{D}_{c}$ the class of estimators

[TABLE]

and with $\psi(\cdot)$ satisfying the conditions of Theorem 2. The optimal choice of $\psi(\cdot)$ in terms of minimising risk, based on the simple upper bound of ${\rm E}_{\theta}\,D^{*}(\phi,Z)$, is $\psi(z)=p-1$, leading to the estimator

[TABLE]

Note that $\delta^{c}$ appropriately generalises the Clevenson–Zidek estimator of (1.2). Importantly, it is clearly seen how fine-tuning of $c$ determines the amount of shrinkage away from the $\delta_{0}$. We can use the expression for $D^{*}(\phi,z)$ derived above to find the risk function for the estimator in (3.5),

[TABLE]

Note that the risk depends on $\theta_{i}$ only through the sum $\gamma$, and that numerical evaluation is easy because $Z$ is Poisson with mean $\gamma$. The risk function starts at

[TABLE]

for $\theta=0$, and then increases continuously towards the minimax risk $p+c$. As illustrated in Figure 1, the improvement over $\delta_{0}=Y$ is substantial for small to moderate values of $\gamma$, and always lies below the risk of the usual estimator $\delta_{0}$.

3.3   Loss function robustness.

Robustness of performance statements with respect to the specific loss function used is often overlooked in the literature, as if the loss function worked with had been handed down from above with absolute precision. The matter is clearly of importance, however, as e.g. briefly pointed to in comments by both Berger (1983) and Morris (1983). If an estimator performs well with respect to one loss function ${\rm Lo}_{a}$, but not for another loss function ${\rm Lo}_{b}$, say, even when these two are close, it is a cause for concern. We argue here, however, that our estimator (3.5), derived under loss function $L_{c}$ of (1.3), is somewhat robust regarding the precise value of $c$.

First consider the estimator $\delta_{\text{CZ}}$ of (1.2), which uses $\phi_{0}(z)=(p-1)/\{p-1+z\}$. It satisfies the second requirement of Theorem 2, i.e. $\phi_{0}(z)z<\phi_{0}(z+1)(z+1)$ for all $z$, but it does not necessarily satisfy the first requirement, namely $0<\phi(z)\leq 2(p-1)/\{p-1+(1+c)z\}$ for all $z$. It is easy to verify that if $0\leq c\leq 1$ then it does satisfy this requirement, however, showing that $\delta_{\text{CZ}}$ has certain robustness properties with respect to the $L_{c}$ loss function: It is minimax and uniformly dominates $\delta_{0}$ under $L_{c}$, provided $0\leq c\leq 1$.

One can similarly study how the $\delta^{c_{0}}(Y)=[(1+c_{0})Z/\{p-1+(1+c_{0})Z\}]\,Y$ of (3.5) fares when the loss function is not quite the $L_{c_{0}}$ under which it was derived, but rather $L_{c}$, with another penalty value of $c$, that is, with a somewhat different penalty paid to incorrect estimation of the sum $\gamma$. From the first condition of Theorem 2 we see that $\delta^{c_{0}}$ is still minimax and uniformly better than $\delta_{0}=y$ under $L_{c}$, provided that $0\leq c\leq 2c_{0}+1$. An immediate implication of this is that all estimators in the class $\mathcal{D}_{c}$ are minimax and uniformly better than $\delta_{0}$ under the $L_{1}^{*}$ loss function, showing that the more prudent estimation strategy $\delta^{c}$, in the sense that it shrinks less that $\delta_{\text{CZ}}$, is robust to $c=0$.

4   Bayes, empirical Bayes, and admissibility

In this section a certain class of priors is studied along with Bayes and empirical Bayes consequences. The estimator $\delta^{c}$ of (3.5) will be shown to be both a natural generalised Bayes estimator with respect to a certain noninformative prior, and a natural empirical Bayes estimator with respect to independent Gamma priors. In addition, a class of proper Bayes estimators belonging to the class $\mathcal{D}_{c}$ is derived from another prior construction. Finally, we are also able to show that for each given $c$, the $\delta^{c}$ of (3.5) is an admissible estimator, under the $L_{c}$ loss function. In particular, the Clevenson–Zidek estimator (1.2) is admissible under the $L_{0}$ function for which it was derived. Our result hence generalises that of Johnstone (1984, Theorem 4.1).

4.1   Priors with sum independent of proportions.

In the following we model the means $(\theta_{1},\ldots,\theta_{p})$ in terms of the sum $\gamma=\sum_{i=1}^{p}\theta_{i}$ and proportions $\pi_{i}=\theta_{i}/\gamma$ for $i=1,\ldots,p$.

Lemma 1**.**

Suppose that $\gamma$ and $\pi=(\pi_{1},\ldots,\pi_{p})$ are independent with simultaneous prior density $q(\gamma)r(\pi_{1},\ldots,\pi_{p-1})$. This corresponds to a density

[TABLE]

for $(\theta_{1},\ldots,\theta_{p})$. The posterior distribution for $\theta$ given the data is on the same form: $\gamma$ and $(\pi_{1},\ldots,\pi_{p})$ are still independent, and

[TABLE]

Proof.

The first part in an exercise in the transformation of random variables, involving calculating the determinant of Jacobi matrix |\partial(\gamma,\pi_{1},\ldots,\pi_{p-1})/\partial(\theta_{1},\ldots,\theta_{p})|=\big{(}\sum_{i=1}^{p}\theta_{i}\big{)}^{-(p-1)}. The second part follows because the combined Poisson likelihood $\prod_{i=1}^{p}\theta_{i}^{y_{i}}e^{-\theta_{i}}/y_{i}!$ is proportional to $\gamma^{z}e^{-\gamma}\pi_{1}^{y_{1}}\cdots\pi_{p-1}^{y_{p-1}}(1-\pi_{1}-\cdots-\pi_{p-1})^{y_{p}}$. ∎

If in particular $(\pi_{1},\ldots,\pi_{p})$ has a Dirichlet prior distribution with parameters $(\alpha_{1},\ldots,\alpha_{p})$, then the posterior is another Dirichlet with updated parameters $(\alpha_{1}+y_{1},\ldots,\alpha_{p}+y_{p})$, and this holds regardless of the prior used for $\gamma$. The important case of independent $\theta_{i}$ from a Gamma $(\alpha,\beta)$ corresponds to a Gamma prior $(p\alpha,\beta)$ for $\gamma$ and an independent Dirichlet $(\alpha,\ldots,\alpha)$ for the proportions.

Suppose $(\pi_{1},\ldots,\pi_{p})$ comes from a symmetric Dirichlet $(\alpha,\ldots,\alpha)$ independent of $\gamma$, the latter coming from a suitable prior $q(\cdot)$. The Bayes estimator under the $L_{c}$ loss function takes the form (3.1), i.e. $\delta_{i}=\{(1+c)/(1+ca/b)\}a_{i}$, with

[TABLE]

and b=\big{\{}{\rm E}\,(\gamma^{-1}\,|\,Y)\big{\}}^{-1}=K(z)/K(z-1), writing $K(z)=\int_{0}^{\infty}\gamma^{z}e^{-\gamma}q(\gamma)\,{\rm d}\gamma$. Letting in particular $\pi$ be uniform over the simplex we obtain the estimator

[TABLE]

We are now in a position to give three pleasing interpretations of the estimator $\delta^{c}$ of (3.5).

First, it is a generalised Bayes estimator. For if $\gamma$ is given a flat prior $q(\gamma)=1$ on the halfline, then $K(z)=\Gamma(z+1)$, which inserted in (4.1) gives $\delta_{B}=\delta^{c}$.

Second, it is a limit of Bayes estimators. Let the proportions $(\pi_{1},\ldots,\pi_{p})$ have a flat Dirichlet $(1,\ldots,1)$ prior, with an independent $\gamma$ from a Gamma $(1,\beta)$. The Bayes solution is then

[TABLE]

and the limit as $\beta\to 0$ is again $\delta_{i}^{c}$ of (3.5).

Thirdly, it is a natural empirical Bayes estimator. To see one of several such constructions let the $\theta_{i}$ be independent ${\rm Gamma}(1,\beta)$. The corresponding exact Bayes solution is given in (3.3). Now $\beta$ is to be estimated from the data. The marginal distribution of $Z$ is found from the facts that $Z$ given $\gamma$ is Poisson with mean $\gamma$, and $\gamma$ is Gamma distributed with parameters $(p,\beta)$, so

[TABLE]

The sum $Z$ is sufficient for $\beta$ and $Z/(p-1+Z)$ is found to be an unbiased estimator for $1/(1+\beta)$, for any $p\geq 2$. Inserting this data-based value in (3.3) produces once again $\delta^{c}$.

Finally, the Gamma $(1,\beta)$ prior construction may be extended to an hierarchical setup where a prior is put on the hyperparameter $\beta$. This extends the approach of Ghosh & Parsian (1981) from the $L_{1}^{*}$ of (1.1) to the $L_{c}$ setting. Let the parameter $\beta$ be distributed according to

[TABLE]

Utilising that

[TABLE]

the Bayes solution is then

[TABLE]

since $Z$ is sufficient for $\beta$. Also, $\beta$ given $Z=z$ is distributed as

[TABLE]

where $B(a,b)=\Gamma(a)\Gamma(b)/\Gamma(a+b)$ is the Beta function. Provided $\beta$ comes from the $s(\beta)$ of (4.2), and $0<\eta\leq(p-2-c)/(1+c)$, then the Bayes solution in (4.3), by virtue of belonging to the class $\mathcal{D}_{c}$, is minimax and uniformly dominates $\delta_{0}$. To see this, insert

[TABLE]

in (4.3), yielding

[TABLE]

By some algebra we obtain that for the Bayes solution we here consider

[TABLE]

This function is non-decreasing for all $z\geq 0$. Moreover, we see that it is bounded above by

[TABLE]

since $\eta\leq(p-2-c)/(1+c)$. This means that the class of Bayes solutions in (4.4), with $\eta$ obeying the constraint mentioned, satisfy both conditions of Theorem 2.

4.2   Admissibility.

So far we have studied estimators that dominate the maximum likelihood estimator under the $L_{c}$ loss function. In this section we will prove that the estimator $\delta^{c}$ of (3.5) cannot be uniformly improved upon, that is, $\delta^{c}$ is admissible. The basic ingredient in this proof is the characterisation of admissibility given by Brown & Farrel (1988, Theorem 2.6). According to this theorem an estimator $\delta$ is admissible for $\theta$ under $(\delta-\theta)^{2}/\theta$ if and only if there exists a sequence of finite measures $\nu_{n}$ such that the Bayes solution with respect to $\nu_{n}$, say $\delta_{n}$, converges to $\delta$ as $n\to\infty$, and

[TABLE]

This prior sequence has to satisfy certain conditions, the details of which are stated in Brown & Farrel (1988) and Johnstone (1984, pages 237-238). For our purpose, it is sufficient to know that such a sequence exists if $\delta$ is admissible.

Consider the class of estimators given by

[TABLE]

The Bayes estimators of (4.1) are on this form with $\kappa(z)=K(z)/K(z-1)$; in particular the estimator $\delta^{c}$ of (3.5) is on this form, with $\kappa(z)=z$. As in Johnstone (1984), it turns out that estimators of the class (4.5) are admissible provided that $\kappa(Z)$ is admissible for $\gamma$ under the loss function $(\kappa-\gamma)^{2}/\gamma$. The theorem below is in part a restatement of his Theorem 4.1.

Theorem 3**.**

Estimators of the form (4.5) are admissible for $(\theta_{1},\ldots,\theta_{p})$ under $L_{c}$ if and only if $\kappa(Z)$ is admissible for $\gamma$ under $L(\gamma,\delta)=(\delta-\gamma)^{2}/\gamma$.

Proof.

The risk function of (4.5) under $L_{c}$ can be written

[TABLE]

with $L(\gamma,\delta)=(\kappa-\gamma)^{2}/\gamma$; see Appendix A for the derivation of this claim. Introduce $u(z)=(1+c)^{2}z/\{p-1+(1+c)z\}$ and write $\bar{L}(\delta,\gamma,z)=u(z)L(\gamma,\delta)$. Let $\delta$ and $\delta^{\prime}$ be of the form (4.5) with $\kappa(Z)$ and $\kappa^{\prime}(Z)$ respectively. If {\rm E}_{\gamma}\big{\{}\bar{L}(\kappa^{\prime}(Z),\gamma,Z)\big{\}}-{\rm E}_{\gamma}\big{\{}\bar{L}(\kappa(Z),\gamma,Z)\big{\}}\leq 0, with strict inequality for some $\gamma$, then $R_{c}(\delta,\theta)-R_{c}(\delta^{\prime},\theta)\leq 0$ with strict inequality for some $\theta$. This shows that if $\kappa(Z)$ is inadmissible under $L^{*}$, then $\delta$ is inadmissible under $L_{c}$.

The contrapositive statement is more enlightening: If $\delta$ is admissible under $L_{c}$ then $\kappa(z)$ is admissible under $\bar{L}$. Conversely, assume that $\kappa(Z)$ is admissible for $\gamma$ under $L$, let $\{\nu_{n}({\rm d}\gamma)\}_{n\geq 1}$ be a sequence of prior measures satisfying the conditions of Brown & Farrel (1988, Theorem 2.6), and let $\rho_{n}(A)=\Gamma(p)\int_{A}\nu_{n}({\rm d}\gamma)\prod_{i=1}^{p}{\rm d}\pi_{i}$ be the prior measure over $(0,\infty)\times S$, where $S$ is the $(p-1)$-dimensional simplex. Then, since $0<u(z)<1+c$ for all $z\geq 1$,

[TABLE]

The right hand side is non-negative for all $n$ since $\delta_{n}$ is the Bayes solution, and it also tends to zero by Theorem 2.6 in Brown & Farrel (1988) because $\kappa(z)$ is admissible. This implies admissibility of $\delta$ under $L_{c}$. ∎

The immediate corollary to Theorem 3 is that the estimator $\delta^{c}$ of (3.5) is admissible under $L_{c}$, because $\kappa(Z)=Z$ is admissible under the weighted squared error loss function $L(\gamma,\delta)=(\delta-\gamma)^{2}/\gamma$. As a special case, the Clevenson–Zidek estimator (1.2) is admissible under the $\sum_{i=1}^{p}(\delta_{i}-\theta_{i})^{2}/\theta_{i}$ loss function.

5   Smoothing towards the mean

In the following we consider different strategies for smoothing the maximum likelihood estimator towards the mean of the observations. Pushing the maximum likelihood estimator towards a data-based point should in many cases yield more reduction in risk than pushing $\delta_{0}$ towards the origin, particularly when the $\theta_{i}$ are not small and not too spread out. This is clearly visible in Figure 1 where the improvement in risk of $\delta^{c}$ compared to $\delta_{0}$ becomes smaller as $\gamma=\sum_{i=1}^{p}\theta_{i}$ grows; the improvement in risk is, in other words, biggest near the ‘point of attraction’.

Ghosh, Hwang & Tsui (1983) considered a modification of $\delta_{\text{CZ}}$ that shrinks the maximum likelihood estimator towards the minimum of the observations, and were able to prove uniform dominance under the weighted loss $L_{1}^{*}$ of (1.1), for their modified estimator. Albert (1981) took the Bayes estimator under $L_{0}(\theta,\delta)=\sum_{i=1}^{p}(\delta_{i}-\theta_{i})^{2}$ as his point of departure, and developed an estimator that pushes the observations towards the mean $\bar{y}$ and performs better than the usual estimator in large parts of the parameter space.

A complication when working with the $L_{c}$ loss function of (1.4) is that we compete with the maximum likelihood estimator on two turfs, so to speak, namely under $\sum_{i=1}^{p}(\delta_{i}-\theta_{i})^{2}/\theta_{i}$ and under \big{(}\sum_{i=1}^{p}\delta_{i}-\gamma\big{)}^{2}/\gamma. One reason for choosing the mean as our ‘point of attraction’ is that several of the estimators we construct preserve the mean, that is, $\sum_{i=1}^{p}\delta_{i}(Y)=Z$. In this way such new estimators will always ‘match’ the risk performance of $\delta_{0}$ when it comes estimating $\gamma$, and the penalty parameter $c$ becomes immaterial.

5.1   A restricted minimax estimator.

Consider estimators of the form

[TABLE]

where $\bar{Y}=Z/p$. Notice that $\sum_{i=1}^{p}\delta_{i}=Z=\sum_{i=1}^{p}\delta_{0,i}$, which means that in calculations of the difference in risk between estimators of the form (5.1) and $\delta_{0}=Y$, the $c$-term in the $L_{c}$ loss function disappears. The value of $g(0)$ is immaterial, since $\delta$ is then equal to zero, and we may take $g(0)=0$ for convenience.

Using that $Y$ given $Z=z$ is multinomial with cell probabilities $\pi_{i}$, for each $z\geq 1$, the risk difference $R_{c}(\delta,\theta)-R_{c}(\delta_{0},\theta)$ can be expressed as

[TABLE]

writing

[TABLE]

for the mean of the inverse proportions. It is not possible to find a function $g(\cdot)$ such that the estimator in (5.1) dominates $\delta_{0}=Y$ over the entire parameter space $\Theta=(0,\infty)^{p}$. This is the ‘tyranny of the boundary’ phenomenon; as a single $\theta_{i}\to 0$ the sum $B(\pi)\to\infty$ and the risk blows up for non-null choices of $g$.

Estimators can be found, however, that dominate $\delta_{0}$ over large proportions of the parameter space. Let $\Theta(b_{0})\subset\Theta$ be the subset of the parameter space where $B(\pi)\leq pb_{0}$. The minimum value of $B(\pi)$ is $p$, so $b_{0}>1$. In some situations one might have prior grounds for believing that the $\theta_{i}$ are somewhat similar, that is, not too far from the mean $\bar{\theta}$. One might for example have that each $\theta_{i}\geq\varepsilon\bar{\theta}=\varepsilon\gamma/p$, or equivalently $\pi_{i}\geq\varepsilon/p$, for some small $\varepsilon>0$. This implies that $B(\pi)\leq p/\varepsilon$, so $b_{0}=1/\varepsilon$ may be used, so the risk difference can be bounded:

[TABLE]

Based on this upper bound we derive the estimator with

[TABLE]

This estimator succeeds in having

[TABLE]

for all $\theta\in\Theta(b_{0})$, with the last inequality following from Jensen’s.

5.2   A restricted Bayes estimator.

The estimator in (5.4) was derived with the aim of risk function dominance in a given large region $\Theta(b_{0})$ of the parameter space. We may also derive the restricted Bayes solution, that is, the best estimator among those of the form (5.1), under a prior of the type discussed in Section 4. Since $\gamma$ and $\pi$ are independent, the Bayes risk of such an estimator is

[TABLE]

where ${\rm E}^{q}\,(\cdot)$ and ${\rm E}^{r}\,(\cdot)$ are the expectations of $\gamma$ and $\pi$ with respect to their prior distributions, cf. Lemma 1. If we now let ${\rm E}^{r}\,B(\pi)=pb_{0}$, this reproduces the estimator in (5.4), but with a differently interpreted $b_{0}$. The risk function is

[TABLE]

with $g(z)=(p-1)/\{p-1+(b_{0}-1)z\}$. Consider again the prior construction of Section 4. If $(\pi_{1},\ldots,\pi_{p})$ comes from a Dirichlet distribution with parameters $(\alpha,\ldots,\alpha)$, then $b_{0}=(\alpha-1/p)/(\alpha-1)$, provided that $\alpha>1$. Recall that the inequality in 5.3 only holds if $b_{0}>1$, which means that the subspace $\Theta(b_{0})$ is empty if our prior knowledge dictates $0<\alpha\leq 1$. On the other side of the spectrum, as $\alpha\to\infty$ this estimator would assign the estimate $1/p$ to each of the proportions, and since $b_{0}$ goes to one as $\alpha$ goes to infinity, the estimator (5.4) tends to $\bar{y}=z/p$. In other words, having smaller values of $\alpha>1$ expands the space $\Theta(b_{0})$ but results in less gain in risk.

The (5.4) estimator may also use a data-based value for $b_{0}$. With $\pi$ and $\gamma$ independent, the marginal distribution of the data is

[TABLE]

This likelihood can be maximised to obtain an estimate $\widehat{\alpha}$ which is plugged into $b_{0}=(\alpha-1/p)/(\alpha-1)$, again provided that $\widehat{\alpha}>1$.

5.3   More careful smoothing towards the mean.

We now consider the construction

[TABLE]

where $h(y)$ is a function such that if one or more $y_{i}=0$, then $h(y)=h(y_{1},\ldots,y_{p})=0$. The intention is that of more careful smoothing towards the mean than with (5.1), to achieve risk improvement in potentially larger parts of the parameter space. Note that $\sum_{i=1}^{p}\widehat{\theta}_{i}=\sum_{i=1}^{p}y_{i}$, so any risk difference $R_{c}(\widehat{\theta},\theta)-R_{c}(Y,\theta)$ does not depend on the $c$ term of the loss function. Also, the value of $g(z)$ at $z=0$ is immaterial, so we may take $g(0)=0$, for convenience.

To work with the risk functions, we start from

[TABLE]

The risk difference $R_{c}(\widehat{\theta},\theta)-R_{c}(Y,\theta)$ may hence be expressed as

[TABLE]

with

[TABLE]

The property that $y_{i}=0$ implies $h(y)=0$ is actively used here; without such a constraint, more complicated terms need to be added to the ${\rm E}_{\theta}\,D(Y)$ here.

Several choices may be considered and worked with for further fine-tuning, regarding the $h(y)$ function. For the present report we limit attention to the special case of $h_{0}(y)=I\{{\rm each\ }y_{i}\geq 1\}$. We need to study and bound the $D(y)$ function with this choice of $h(y)$. Note that $\sum_{i=1}^{p}h_{0}(y+e_{i})$ is $p$, if all $y_{i}\geq 1$; is 1, if one and only one of $y_{i}=0$; and is 0, otherwise. Similarly, $\sum_{i=1}^{p}y_{i}h_{0}(y+e_{i})=z$, if all $y_{i}\geq 1$; and 0, otherwise. Furthermore,

[TABLE]

(i) Assume first that all $y_{i}\geq 1$. Then

[TABLE]

The function p-1-x+\hbox{1\over 2}x^{2}/p for $x\geq 1$ has its minimum value at position $x=p$, with minimum value \hbox{1\over 2}p-1, which is positive as long as $p\geq 3$. (ii) When there is only one $y_{i}=0$, the rest $y_{j}\geq 1$, we find

[TABLE]

(iii) Otherwise, which means that the number of $y_{i}=0$ is between 2 and $p-1$, we find that $D(y)=0$.

Our best choice for $g(z)$, based on the upper risk bound for the case of all $y_{i}\geq 1$, is

[TABLE]

The estimator

[TABLE]

therefore, has risk function $R_{c}(\widehat{\theta})=p+c+{\rm E}_{\theta}\,D_{0}(Y)$, where an exact expression for $D_{0}(Y)$ is found via the above. We also know that

[TABLE]

that

[TABLE]

for the cases of $y$ where precisely one $y_{i}=0$, the other $y_{j}\geq 1$; and finally that $D_{0}(y)=0$ for those $y$ for which the number of $y_{i}=0$ is between 2 and $p-1$. As long as all $\theta_{i}$ are at least moderate, so that $h_{0}(Y)=1$ with high probability, there is clear risk improvement on the minimax risk $p+c$. The corner cases, however, where one $\theta_{i}$ is small and the others not, are the ‘bad cases’ for the (5.7) estimator, where the risk might become larger than $p+c$. Since $D_{0}(y)$ flattens to zero when $z=\sum_{i=1}^{p}y_{i}$ increases, the risk converges to $p+c$ for all $\theta=\gamma(\pi_{1},\ldots,\pi_{p})$ where $\gamma$ tends to infinity.

6   Bayes and empirical Bayes with more structure

In the previous sections we have constructed and analysed estimators essentially symmetric in the observations. Sometimes some structure is anticipated in the parameters, however, as with setting up regressions or log-linear models for Poisson tables, see e.g. Agresti (2019), or with classes of priors. The present section briefly complements our earlier efforts by examining risk function consequences for estimators that favour asymmetric structures.

6.1   Risk functions with Gamma priors.

A natural class of priors takes independent Gamma priors, with parameters $(\alpha_{i},\beta)$, for the $\theta_{i}$. As was seen in Section 3, the Bayes estimator then takes the form

[TABLE]

writing $\bar{\alpha}$ for $(1/p)\sum_{i=1}^{p}\alpha_{i}$; also, for cases where $\alpha_{i}<1$ and $y_{i}=0$, the estimator is zero.

The present task is to study the consequent risk functions, for such estimators, outside the special and simplest case where each $\alpha_{i}=1$. The point will be that estimators then lose the minimax property, with the risk exceeding the $p+c$ benchmark level when one or more of the $\theta_{i}$ come close to zero, but that the risk otherwise can be lower than $p+c$ in big and reasonable parts of the parameter space. For simplicity of presentation we restrict attention here to the case of $c=0$. Similar results and insights may be reached for the general loss function $L_{c}$, using more laborious methods we develop and exploit for somewhat different purposes in Section 6.2. For the ensuing estimator $(\alpha_{i}+Y_{i}-1)/(\beta+1)$, some calculations yield

[TABLE]

This is smaller than or equal to $R(Y,\theta)=p$ when

[TABLE]

This defines a fairly large parameter region, $\{\theta\colon A(\theta)\leq 2+\beta\}$, with the $1/\theta_{i}$ not being too far away from the prior mean values $\beta/(\alpha_{i}-1)$, and where using the Bayes estimator hence is better than with the $\delta_{0}=Y$ method. Under the prior itself, the random $A(\theta)$ has mean

[TABLE]

and variance of order $O(1/p)$, showing that only rather unlikely values of $\theta$ will have risk above the benchmark $p$.

6.2   Risks for a class of empirical Bayes estimators.

With the independent Gamma priors used above we next note that the marginal distribution of $y_{i}$ becomes

[TABLE]

for $y_{i}=0,1,2,\ldots$. In the setup where the $\alpha_{i}$ are taken known but $\beta$ an unknown parameter, we see that $Z=\sum_{i=1}^{p}Y_{i}$ is sufficient. Since $Z$ given the parameters is Poisson with mean $\gamma=\sum_{i=1}^{p}\theta_{i}$, and $\gamma$ is Gamma $(p\bar{\alpha},\beta)$, writing $\bar{\alpha}=(1/p)\sum_{i=1}^{p}\alpha_{i}$, its distribution may be written

[TABLE]

From this we can derive

[TABLE]

provided only that $p\geq 2$. Hence $Z/(p\bar{\alpha}-1+Z)$ can be used as an estimator for the quantity $1/(\beta+1)$. With pre-specified $\alpha_{i}$, then, a natural empirical Bayes version of estimator (6.1) emerges:

[TABLE]

with

[TABLE]

In particular, under $L_{0}$ loss, with $c=0$, the natural empirical Bayes estimator is

[TABLE]

generalising the earlier symmetric case with all $\alpha_{i}=1$, which yields the already studied minimax and admissible estimator (1.5).

Expressions for the risk function $R_{c}(\widehat{\theta},\theta)$ can now be worked out, using the fact that $y_{i}\,|\,z$ is binomial $(z,\pi_{i})$, with $\pi_{i}=\theta_{i}/\gamma$. Consider the general class of estimators

[TABLE]

The $h(z)$ functions we are encountering all have $h(0)=0$, and will in fact have the form $q(z)z$, for suitable $q(z)$. Also, they will be nondecreasing with $h(z)\rightarrow 1$ as $z$ increases. The task now is to develop formulae for their risk functions, through suitable representations of the form

[TABLE]

and then showing, for relevant choices of $h(z)$, that this is less than $p+c$ for large and relevant parameter regions.

We start from

[TABLE]

For the case of $c=0$ this leads to

[TABLE]

For the ‘$c$ term’ part of the risk, note that $\widehat{\theta}_{i}=h(z)(p\bar{\alpha}+z-p)$, leading to

[TABLE]

This may be exploited further using the identity ${\rm E}_{\gamma}\,r(Z)/\gamma={\rm E}_{\gamma}\,r(Z+1)/(Z+1)$ for functions $r(z)$ with $r(0)=0$.

For brevity we limit attention here to the case of $c=0$; extensions can be worked out using the same methods. We use the identity pointed to for

[TABLE]

and find $R(\widehat{\theta},\theta)={\rm E}_{\gamma}\,Q(Z,\theta)$, with

[TABLE]

The risk function may hence be expressed as

[TABLE]

in which $C(\pi)=(1/p)\sum_{i=1}^{p}(\alpha_{i}-1)^{2}/\pi_{i}$, and

[TABLE]

Here $m_{1}(\gamma)$ and $m_{2}(\gamma)$ are inside $(0,1)$, and increase to 1; the $R_{1}(\gamma)$ and $R_{2}(\gamma)$ are bounded and converge to respectively one and zero as $\gamma$ increases. For larger $\gamma$, therefore, the risk goes to the minimax risk $p$. The risk function (6.8) may exceed the minimax threshold level $p$ if one or more of the $\pi_{i}=\theta_{i}/\gamma$ are small, but even for small $\pi_{i}$ the risk decreases with increasing $\gamma$. Otherwise the situation is that the risk may become significantly smaller than $p$ in parts of the parameter space not disagreeing much from what is judged likely under the prior, and that it can be smaller than $p$ also in other larger parameter regions. An upper bound is

[TABLE]

with $B(\pi)$ as in (5.2). This may in particular be investigated further, with the choice $h_{0}(z)$, corresponding to the estimator (6.6). The risk function is bounded; converges to the minimax value $p$ when $\gamma$ increases, regardless of proportions $(\pi_{1},\ldots,\pi_{p})$; may offer substantial improvement for sizeable portions of the parameter space; and its maximum value is often not much bigger than $p$.

An illustration of the empirical Bayes strategy (6.6) is provided in Figure 2, in a situation with simulated regression data $(x_{i},y_{i})$. The prior takes the $\theta_{i}$ to stem from Gamma $(\alpha_{i},\beta)$, with $\alpha_{i}=\exp(\gamma_{0}+\gamma_{1}x_{i})$, for suitable prior guess values $(\gamma_{0},\gamma_{1})$, and then estimates $\beta$ from data, as per (6.2). Similar Bayes and empirical Bayes methods can be developed for priors of the type $\theta_{i}\sim{\rm Gamma}(d\theta_{0,i},d)$, with either hyperpriors on the prior parameters $\theta_{0,i}$ and $d$, or involving estimators for these from the data.

7   Estimation with a weighted loss function

Above we have worked with our loss function $L_{c}$, a natural extension of the Clevenson–Zidek loss function $L_{1}^{*}(\theta,\delta)=\sum_{i=1}^{p}(\delta_{i}-\theta_{i})^{2}/\theta_{i}$, to account for not shrinking the mean too much. Another useful extension is to allow for weighting, with

[TABLE]

where $w_{1},\ldots,w_{p}$ are fixed, positive, and context driven, reflecting relative importance. This is e.g. important when the $y_{i}$ result from different levels of exposure, as with $y_{i}$ stemming from a Poisson with parameter $w_{i}\theta_{i}$. The $\delta_{0}=Y$ estimator is again the natural benchmark; it has constant risk $w_{0}=\sum_{i=1}^{p}w_{i}$, and it is minimax. To prove this second claim, one may work with the prior where the $\theta_{i}$ are independent and Gamma distributed, with parameters $(\alpha,\beta)$, and with $\alpha\geq 1$ to avoid a certain technical issue. Some work shows that the Bayes estimator becomes $\widehat{\theta}_{i}=1/{\rm E}\,(\theta_{i}\,|\,y)=(\alpha+y_{i}-1)/(\beta+1)$, with associated minimum Bayes risk as simple as $w_{0}/(\beta+1)$. Letting $\beta\rightarrow 0$ we have convergence to the constant risk of $\delta_{0}$.

This estimator may be uniformly improved upon, however, as we now demonstrate, yielding another generalisation of the Clevenson–Zidek estimator. Consider estimators of the form

[TABLE]

Using the identity (3.4) we may express the risk difference $r_{w}$ between $\widehat{\theta}$ and $\delta_{0}$, i.e. ${\rm E}_{\theta}\,\{L_{w}(\theta,\widehat{\theta})-L_{w}(\theta,\delta_{0})\}$, as

[TABLE]

say, where

[TABLE]

As we have argued on previous occasions in our paper, if we succeed in finding a function $\phi(v)$ such that $D_{w}(y)\leq 0$ for all $y$, with strict inequality for at least one $y$, we have demonstrated inadmissibility of $\delta_{0}=Y$.

To work with this we set up some mild requirements regarding the $w_{i}$ weights. We take all $w_{i}$ to be inside some $[a,b]$ interval, situated inside $(0,1)$, and also stipulate that $w_{0}=\sum_{i=1}^{p}w_{i}>1$. For estimators of the form (7.1), consider

[TABLE]

with (i) $\psi(v)$ nondecreasing, (ii) with $0<\psi(v)<2(w_{0}-1)$ for all $v=v(y)$, and (iii) such that $v\phi(v)$ is increasing. We then find

[TABLE]

which is demonstrably negative. Our preferred generalisation of the Clevenson–Zidek estimator, to the present case of weighted loss, becomes

[TABLE]

The special case of equal weights $w_{i}=1$ leads back to the Clevenson–Zidek estimator (1.2).

Remarkably, the apparatus above allows extension to the case of infinitely many Poisson parameters. Suppose $Y_{1},Y_{2},\ldots$ are independent Poisson counts with means $\theta_{1},\theta_{2},\ldots$, and that loss incurred by estimators $\delta_{1},\delta_{2},\ldots$ is taken to be $L(\theta,\delta)={\sum_{i=1}^{\infty}}w_{i}(\delta_{i}-\theta_{i})^{2}/\theta_{i}$. Here the sequence of weights is such that $w_{0}={\sum_{i=1}^{\infty}}w_{i}$ is finite, and the parameter space to be considered is $\Omega$, the set of sequences of $\theta_{i}$ for which ${\sum_{i=1}^{\infty}}w_{i}\theta_{i}$ is finite (including in particular each bounded sequence).

The benchmark procedure is again $\delta_{0}$, with components $\delta_{0,i}=y_{i}$. It has constant risk $w_{0}$, and our previous arguments may be extended to demonstrate that this procedure is minimax. Also, crucially, the estimator

[TABLE]

offers uniform risk improvement over $\delta_{0}$, where now $v(y)={\sum_{i=1}^{\infty}}w_{i}y_{i}$. This follows from arguments used to reach the corresponding statement for the finitely-many procedure (7.2), in the light of two necessary remarks. The first is that the identity (3.4) continues to hold, for all $h(Y_{1},Y_{2},\ldots)$ with finite mean, with the property that $h(y)=0$ as long as $y_{i}=0$. The second is that the upper bound reached above for $D_{w}(y)$, under the condition that $a\leq\theta_{i}\leq b$ for all $i$, applies here too, but now we need $a=0$, since the ${\sum_{i=1}^{\infty}}w_{i}$ is finite. In other words,

[TABLE]

remains correct.

Above we stipulated a scaling of the importance weights $w_{i}$ so that their sum $w_{0}={\sum_{i=1}^{\infty}}w_{i}$ is above 1. This is partly in order for the estimator (7.3) to be a natural generalisation of the Clevenson–Zidek estimator. Similar reasoning goes through for estimators $[\psi(v(Y))/\{w_{0}-\varepsilon+v(Y)\}]\,Y$, if we instead stipulate $w_{0}>\varepsilon$.

8   Multivariate models for count variables

The methods developed in Section 4 utilised certain constructions which also involve multivariate models for rates and for count variables, of interest in their own right. Models can be built with both positive and negative correlations betwen rates $\theta_{i}$ and between count observations $Y_{i}$. These modelling ideas also point to Bayesian nonparametrics, cf. our Remark B in Section 9.

8.1   Sum and proportions models.

Suppose first in general terms that given $\theta=(\theta_{1},\ldots,\theta_{p})$, the observations $Y_{1},\ldots,Y_{p}$ have independent Poisson distributions with these parameters, and that the $\theta$ has a background distribution, which we for simplicity of presentation here take to be symmetric with finite variances. Let us write ${\rm E}\,\theta_{i}=\theta_{0}$, ${\rm Var}\,\theta_{i}=\sigma_{0}^{2}$, $\text{cov}(\theta_{i},\theta_{j})=\rho\sigma_{0}^{2}$ for $i\not=j$. We then deduce

[TABLE]

A class of multivariate models for the $\theta_{i}$, and by implication also for the $Y_{i}$, emerges from the construction of Lemma 1, with a prior $q(\gamma)$ for $\gamma=\sum_{i=1}^{p}\theta_{i}$ and a symmetric prior ${\rm Dir}(\alpha,\ldots,\alpha)$ for the proportions $\pi_{i}=\theta_{i}/\gamma$. Write ${\rm E}\,\gamma=\gamma_{0}=p\theta_{0}$ and ${\rm Var}\,\gamma=p\tau_{0}^{2}$. The $\pi_{i}$ have means $1/p$, variances $(1/p)(1-1/p)/(p\alpha+1)$, and covariances $-(1/p^{2})/(p\alpha+1)$. From these facts we first find ${\rm E}\,\theta_{i}={\rm E}\,(\gamma\pi_{i})=\theta_{0}$ and then

[TABLE]

Similarly, some calculations lead to $\text{cov}(\theta_{i},\theta_{j})=(\tau_{0}^{2}\alpha-\theta_{0}^{2})/(p\alpha+1)$, so that the correlation parameter $\rho$ of (8.1) may be expressed as

[TABLE]

For the special case of $\gamma\sim{\rm Gamma}(p\alpha,\beta)$, we have $\theta_{0}=\alpha/\beta$ and $\tau_{0}^{2}=\alpha/\beta^{2}$, the covariance is zero, and the formula for ${\rm Var}\,\theta_{i}$ gives $\alpha/\beta^{2}$. This is indeed the familiar case of independent $\theta_{i}\sim{\rm Gamma}(\alpha,\beta)$.

For other models for $\gamma$, however, the construction above leads to useful multivariate models for $(\theta_{1},\ldots,\theta_{p})$ and for $(Y_{1},\ldots,Y_{p})$. In general terms, if $\gamma$ has density $q(\gamma)$, the distribution of the data vector can be written as

[TABLE]

with $K(z)=\int_{0}^{\infty}\exp(-\gamma)\gamma^{z}q(\gamma)\,{\rm d}\gamma$ as in Section 4. If in particular $\gamma\sim{\rm Gamma}(\alpha_{0},\beta_{0})$, then

[TABLE]

The point is that this has a multiplicative independence structure only if $\alpha_{0}=p\alpha$, then with negative binomial marginal distributions. In this light, the (8.4) construction amounts to an extended class of models for count data, allowing both positive and negative correlations. The generalisation to the nonsymmetric case takes $\pi\sim{\rm Dir}(\alpha_{1},\ldots,\alpha_{p})$, and leads to

[TABLE]

with the particular choice $\gamma\sim{\rm Gamma}(\sum_{i=1}^{p}\alpha_{i},\beta)$ yielding independent negative binomials.

8.2   Poisson processes with dependence in time and space.

Suppose independent Poisson processes $Y_{1}(t),\ldots,Y_{k}(t)$ are observed over a time period $[0,\tau]$, and divide this period into $p$ cells or windows, say $(t_{j-1},t_{j}]$. We take $Y_{i}(t)$ to have cumulative intensity function $G_{i}(t)$, with consequent Poisson parameters $\theta_{i,j}=G_{i}(t_{j})-G_{i}(t_{j-1})$ for the counts $Y_{i,j}=Y_{i}(t_{j})-Y_{i}(t_{j-1})$, for $j=1,\ldots,p$. There are now different ways of modelling the $k\times p$ matrix of rate parameters, using aspects of the apparatus above. For simplicity of presentation we limit these brief pointers to the neutral cases, where the vectors of fractions involved come from symmetric Dirichlet distributions.

Idea (a) is to allow for dependence over time, for each process:

[TABLE]

and with $(\pi_{i,1},\ldots,\pi_{i,p})$ from a ${\rm Dir}(\alpha_{i},\ldots,\alpha_{i})$. The particular case of $\gamma_{i}\sim{\rm Gamma}(p\alpha_{i},\beta_{i})$ corresponds to independent $\theta_{i,j}\sim{\rm Gamma}(\alpha_{i},\beta_{i})$ for $j=1,\ldots,p$. Idea (b) is to build dependence structure into the sequence of processes:

[TABLE]

and with $(\pi_{1,j},\ldots,\pi_{k,j})$ from a ${\rm Dir}(\alpha_{j},\ldots,\alpha_{j})$. If in particular $\kappa_{j}\sim{\rm Gamma}(k\alpha_{j},\beta_{j})$, then we have independent $\theta_{i,j}\sim{\rm Gamma}(\alpha_{j},\beta_{j})$ for $i=1,\ldots,p$.

Both ideas (a) and (b) have Bayesian counterparts, and motivate extensions of the $L_{c}$ loss function. Suppose we are interested in precise estimates of the full $k\times p$ parameter matrix and in the cumulatives $\gamma_{i}=G_{i}(\tau)$. A natural loss function is then

[TABLE]

where $L_{c}$ is the loss function introduced in Section 3, and $c$ is a positive constant set by the statistician. The Bayes solution is

[TABLE]

where $a_{i,j}=\{{\rm E}\,(\theta_{i,j}^{-1}\mid y)\}^{-1}$, $a_{i}=\sum_{j=1}^{p}a_{i,j}$, and $b_{i}=\{{\rm E}\,(\gamma_{i}^{-1}\mid y)\}^{-1}$. Suppose prior knowledge dictates that intensities of the $k$ processes are functionally somewhat alike, but at different levels; then a natural prior construction takes $(\pi_{i,1},\ldots,\pi_{i,p})\equiv(\pi_{1},\ldots,\pi_{p})\sim{\rm Dir}(\alpha_{1},\ldots,\alpha_{p})$ (i.e., one draw for all $i$) and $\gamma_{i}\sim{\rm Gamma}(\alpha_{0,i},\beta_{0,i})$ for $i=1,\ldots,k$, with these being independent and also independent of $(\pi_{1},\ldots,\pi_{p})$. The posterior is then

[TABLE]

where $Z_{j}=\sum_{i=1}^{k}\{Y_{i}(t_{j})-Y_{i}(t_{j-1})\}$ for $j=1,\ldots,p$. The Bayes estimator that emerges is

[TABLE]

for $j=1,\ldots,p$ and $i=1,\ldots,k$, where $\bar{\alpha}=(1/p)\sum_{j=1}^{p}\alpha_{j}$ and $Z$ is the total sum $\sum_{j=1}^{p}Z_{j}$. This and the accompanying natural frequentist estimator

[TABLE]

are interesting because they borrow information in all directions, so to speak; cross-sectionally through $Z_{j}$; in time through $Y_{i}(\tau)$; and both horizontally and vertically via the total sum $Z$. From Section 3 we know that if the two parts of $L_{c}^{*}$ are viewed separately, the estimators

[TABLE]

uniformly dominate $Y_{i,j}$ and $Y_{i}(\tau)$, under the loss functions $\sum_{i=1}^{k}\sum_{j=1}^{p}(\delta_{i,j}-\theta_{i,j})^{2}/\theta_{i,j}$ and $\sum_{i=1}^{k}(\delta_{i}-\gamma_{i})^{2}/\gamma_{i}$, respectively. The estimators (8.5)–(8.6) provide guidance on how to exploit the multivariate nature of the problem in order to compromise between $\delta_{i,j}^{\prime}$ and $\delta_{i}^{\prime}$, and thereby achieve risk dominance in large parts of the parameter space, or even uniformly.

9   Concluding remarks

We round off our paper by offering a list of concluding remarks, some pointing to further research.

A. Normal approximations and the square-root transformation. When the $\theta_{i}$ are likely to not being small, normal approximations might work well, and multiparameter estimation may proceed via e.g. the approximate model $2y_{i}^{1/2}\sim{\rm N}(2\theta_{i}^{1/2},1)$. The point is that there is a voluminous literature on shrinkage methods for normal setups, with Stein–James estimators etc. This also invites loss functions of the type $\sum_{i=1}^{p}(\delta_{i}^{1/2}-\theta_{i}^{1/2})^{2}$, which may also be motivated via the Hellinger distance between the real and estimated Poisson distributions. Methods of our paper, using exact Poisson calculations as opposed to normal approximations, may indeed be used to demonstrate that estimators of the type $\widehat{\theta}_{i}^{1/2}=\{1-g(z)\}y_{i}^{1/2}$ may be found, which dominate the default method’s $y_{i}^{1/2}$ in large parts of the parameter space. Specifically, if $p\geq 2$, and given a lower positive threshold $\pi_{0}$, methods exist which dominate the default method in the region characterised by having all $\pi_{i}=\theta_{i}/\gamma\geq\pi_{0}$.

B. Bayesian nonparametrics. Suppose a time-inhomogeneous Poisson process $Y=\{Y(t)\colon t\geq 0\}$ is observed over some time window $[0,\tau]$, with cumulative intensity $G=\{G(t)\colon t\geq 0\}$. We wish to estimate this function, with the loss function

[TABLE]

with some fixed weight measure $w$. Note however that this is not an obvious extension of our previous Clevenson–Zidek loss function, since it works via the cumulatives. A natural estimator is $Y$ itself, with constant risk function $r_{0}=\int_{0}^{\tau}w({\rm d}t)$. A natural class of priors takes $G\sim{\rm Gamma}(aG_{0},a)$, say, a Gamma process with independent increments and $G(t)\sim{\rm Gamma}(aG_{0}(t),a)$. Then $G(t)\,|\,{\rm data}\sim{\rm Gamma}(aG_{0}(t)+Y(t),a+1)$, and the Bayes estimator becomes

[TABLE]

Working with expressions for the minimum Bayes risk one may show that this converges as $a\rightarrow 0$ to the value $r_{0}$, proving that the estimator $Y$ is minimax.

A larger class of priors can however also be investigated, inspired by methods and results of our Section 4. Write $G=\gamma F$, with $\gamma=G(\tau)$ the full mass and $F=G/\gamma$ normalised to be a cumulative distribution function. Now construct a prior by having $\gamma$ from some density $q(\gamma)$ independent of a Dirichlet process for $F$, say $F\sim{\rm Dir}(bF_{0})$, i.e. prior mean $F_{0}$ and $b$ the mass or precision parameter. An appropriate extension of our Lemma 1 then shows that (i) $\gamma$ and $F$ are independent, given data; (ii) $\gamma\,|\,{\rm data}$ has a density proportional to $q(\gamma)\gamma^{z}\exp(-\gamma)$, with $z=Y(\tau)$; and (iii) $F\,|\,{\rm data}$ is a Dirichlet $bF_{0}+Y$. The Bayes estimator becomes

[TABLE]

with $K(z)=\int_{0}^{\infty}\gamma^{z}e^{-\gamma}q(\gamma)\,{\rm d}\gamma$ as in Section 4. Note that this construction provides a genuine extension of the Gamma process, and still with easy control of the posterior distribution for $G=\gamma F$. The usual Gamma process, with $G(t)\sim{\rm Gamma}(aG_{0}(t),a)$, corresponds to the special case of $G=\gamma F$, where $\gamma$ has a ${\rm Gamma}(aG_{0}(\tau),a)$ and is independent of $F\sim{\rm Dir}(aG_{0})$, i.e. with prior mean $F_{0}=G_{0}/G_{0}(\tau)$ and mass parameter $aG_{0}(\tau)$.

C. Separating sum and proportions. When estimating Poisson parameters $\theta_{1},\ldots,\allowbreak\theta_{p}$, another type of loss function than those worked with above is to separate the sum $\gamma$ and proportions $\pi_{i}=\theta_{i}/\gamma$, and then work with $(\widehat{\gamma}-\gamma)^{2}/\gamma+z\sum_{i=1}^{p}(\widehat{\pi}_{i}-\pi_{i})^{2}/\pi_{i}$. The maximum likelihood procedure corresponds to $\widehat{\pi}_{i}=Y_{i}/Z$ and $\widehat{\gamma}=Z$, having constant risk function $1+(p-1)=p$. Procedures achieving lower risk in large parameter regions may be constructed via empirical Bayes arguments.

D. Shrinking towards submodels. There is scope for extension of our methods and constructions in several directions for multiparameter Poisson- and Poisson-related inference. It is inherently useful to shrink raw estimates towards meaningful submodels, such as with log-linear setups for analysis of tables of count-data, see e.g. Agresti (2019). There are notable Poisson-related models for small-area estimation, involving also mixed Poisson models, zero-inflated versions, etc. Generally speaking Bayes and empirical Bayes constructions can be brought to such tables, and will tend to work well for sizeable parameter regions, whereas the exact risk calculations worked with in the present article are harder to generalise.

Acknowledgments

The authors are grateful for partial funding from the PharmaTox Strategic Research Initiative at the Faculty of Mathematics and Natural Sciences (for E.Aa.S’s PhD work), and from the Norwegian Research Council for the five-year research project group FocuStat: Focus Driven Statistical Inference with Complex Data, at the Department of Mathematics (led by N.L.H.), both at University of Oslo.

Appendix A Appendix

Here we prove the risk function expression in (4.6), needed in Theorem 3. For the $L_{c}(\theta,\delta)$ loss function, and for estimators of the form

[TABLE]

considered in (4.5), we have

[TABLE]

proving what was required.