Preface

The contents of chapters 2–8 are summarized in table 1. From the table, we observe the following facts.

1. Each chapter contains a theoretical section, simulations section, and/or a real data section.

2. For the theoretical section, every chapter contains two subsections: Bayes estimators and PESLs, and empirical Bayes estimators of .

3. For the simulations section, every chapter contains four subsections: Two inequalities of Bayes estimators and PESLs, consistencies of moment estimators and Maximum Likelihood Estimators (MLEs), goodness-of-fit of the model, and marginal distributions for various hyperparameters.

4. If a chapter contains a subsection, then we place a . However, if a chapter does not contain a subsection, then we place a .

5. The marginal distributions of the first six hierarchical models (Inverse Gamma-Inverse Gamma (IG-IG), Gamma-Gamma (G-G), Exponential-Inverse Gamma (Exp-IG), Normal-Inverse Gamma (N-IG), Normal-Normal Inverse Gamma (N-NIG), and Uniform-Inverse Gamma (U-IG)) are continuous. Thus they can be used to model continuous data. The Kolmogorov-Smirnov (KS) test is used to perform the goodness-of-fit of the model to the data. The marginal distribution of the last hierarchical model (Poisson-Gamma (P-G)) is discrete, and thus it can be used to model discrete data. The chi-square test is utilized to perform the goodness-of-fit of the model to the data.

TAB. 1 — P: The contents of chapters 2–8.

Subsection	IG-IG	G-G	Exp-IG	N-IG	N-NIG	U-IG	P-G
Theoretical section	Bayes estimators and PESLs
	Empirical Bayes estimators of
	Theoretical comparisons of Bayes estimators and PESLs of three methods
Simulations section	Two inequalities of Bayes estimators and PESLs
	Consistencies of moment estimators and MLEs
	Goodness-of-fit of the model
	Numerical comparisons of Bayes estimators and PESLs of three methods
	Marginal distributions for various hyperparameters
Real data section	A real data example

Chapter 1　Introduction

1.1　Empirical Bayes Method

1.2　The Gamma and Inverse Gamma Distributions

Suppose that and . More specifically, the pdfs of and are respectively given by

It is easy to calculate

and

Suppose that and . The pdfs of and are respectively given by

It is easy to calculate

and

1.3　Hierarchical Models with Positive Parameters

In this section, we will introduce 7 hierarchical models with positive parameters. The hierarchical models are in the following general form:

(1.1)

where are hyperparameters to be estimated, , , or in this book, is the unknown parameter of interest, is the distribution of with parameter , and is the prior distribution of with hyperparameters . It is useful to point out that some hyperparameters may exist in , and thus is more accurate. However, for simplicity, we will use . The hierarchical models that we will consider in this book are IG-IG (2.1), G-G (3.1), Exp-IG (4.1), N-IG (5.1), N-NIG (6.2), U-IG (7.1), and P-G (8.1).

1.4　Estimating the Hyperparameters

The moment method to estimate the hyperparameters is performed by equating the population moments to the sample moments. In general, if there are hyperparameters, then we need to calculate the first origin moments of , . We can use the iterated expectation method to calculate , that is,

where . Assume that can be calculated. Then

and this expectation may be calculated by noting .

The MLE method to estimate the hyperparameters proceeds as follows. First, we calculate the likelihood function of :

(1.2)

where

is the marginal distribution of the hierarchical model (1.1), is the probability density function (pdf) or probability mass function (pmf) of , and is the prior pdf of . Second, we obtain the log-likelihood function of :

Third, taking partial derivatives with respect to and setting them to zeros, we obtain

Fourth, after some algebra, the above equations reduce to

(1.3)

(1.4)

(1.5)

In general, the analytical calculations of the Maximum Likelihood Estimators (MLEs) of by solving the equations (1.3), (1.4), ..., and (1.5) are impossible, and thus we have to resort to numerical solutions. Finally, we can exploit Newton’s method to solve the equations (1.3), (1.4), ..., and (1.5) and to numerically obtain the MLEs of . The iterative scheme of Newton’s method is

where is the Jacobian matrix of , and . Note that the MLEs of are very sensitive to the initial estimators, and the moment estimators are usually proved to be good initial estimators. The Jacobian matrix can be calculated as follows:

where

It is important to note that in calculating (1.2), we have implicitly used the property of independency of . And this independency property is guaranteed by

(1.6)

If (1.6) is true, then we have

that is,

(1.7)

Moreover, if (1.7) is true, then we have (1.6) is true. In other words, (1.6) and (1.7) are equivalent. In addition, (1.7) implies

(1.8)

It is also important to note that

(1.9)

can not guarantee (1.8), and vice versa.

1.5　Stein’s Loss Function

In this section, we will introduce Stein’s loss function and justify why Stein’s loss function is better than the squared error loss function on .

The (weighted) squared error loss function has been used by many authors for the problem of estimating the variance, , based on a random sample from a normal distribution with an unknown mean (see, for example, Maatta and Casella (1990); Stein (1964)). As pointed out by Casella and Berger (2002), the (weighted) squared error loss function penalizes equally for overestimation and underestimation, and it is fine for the unrestricted parameter space . In the positive parameter space where 0 is a natural lower bound, and the estimation problem is not symmetric, we should not select the (weighted) squared error loss function, but select a loss function which penalizes gross overestimation and gross underestimation equally, that is, an action will incur an infinite loss when it tends to 0 or ∞. Stein’s loss function has this property, and hence it is recommended to use for the positive parameter space by many authors (see for instance Li et al. (2025); Shi et al. (2025); Zhang (2025); Sun et al. (2024); Zhang et al. (2024); Sun et al. (2021); Bobotas and Kourouklis (2010); Zhang et al. (2019b); Xie et al. (2018); Zhang et al. (2018); Zhang (2017); Oono and Shinozaki (2006); Petropoulos and Kourouklis (2005); Parsian and Nematollahi (1996); Brown (1990, 1968); James and Stein (1961)).

Now, let us give the justifications of why Stein’s loss function is better than the squared error loss function on . Stein’s loss function is given by

(1.10)

while the squared error loss function is given by

(1.11)

Note that on the positive parameter space , Stein’s loss function penalizes gross overestimation and gross underestimation equally, that is, an action will incur an infinite loss when it tends to 0 or ∞. However, the squared error loss function does not penalize gross overestimation and gross underestimation equally, as an action will incur a finite loss (in fact ) when it tends to 0 and incur an infinite loss when it tends to ∞. Figure 1.1 shows Stein’s loss function and the squared error loss function on when .

FIG. 1.1 — I: Stein’s loss function and the squared error loss function on when .

1.6　The Bayes Estimators and the PESLs

In this section, we will calculate the Bayes estimator of under Stein’s loss function , the Bayes estimator of under the usual squared error loss function , and the Posterior Expected Stein’s Losses (PESLs) at and ( and ) for any hierarchical model such that the posterior expectations exist.

Similar to Zhang (2017), the two Bayes estimators and the two PESLs are respectively given by

(1.12)

(1.13)

(1.14)

(1.15)

To calculate the two Bayes estimators and the two PESLs, it remains to calculate

It has been shown in Zhang (2017) that

(1.16)

by exploiting Jensen’s inequality. Moreover,

(1.17)

which is a direct consequence of the general methodology for finding a Bayes estimator. According to construction, minimizes the Posterior Expected Stein’s Loss (PESL). In the simulations section and the real data section, we will exemplify the two inequalities (1.16) and (1.17).

1.7　Theoretical Comparisons of the Bayes Estimators and the PESLs of Three Methods

Note that the subscripts 0, 1, and 2 below are for the oracle method, the moment method, and the MLE method, respectively. Similar to the derivations in Sun et al. (2021); Zhang (2017), the PESL functions of the three methods are respectively given by

for .

Now we calculate the Bayes estimators and , and the PESLs and following the route of Sun et al. (2021); Zhang (2017). The Bayes estimators of under Stein’s loss function, , minimize the corresponding PESL, that is,

where is the action space, is an action (estimator), given by (1.10) is Stein’s loss function, and is the unknown parameter of interest. Similar to Sun et al. (2021); Zhang (2017), it is easy to obtain

(1.18)

The Bayes estimators of under the squared error loss function are given by

(1.19)

The PESLs evaluated at the Bayes estimators are given by

(1.20)

The PESLs evaluated at the Bayes estimators are given by

(1.21)

Similar to Sun et al. (2021), our primary objectives are to estimate the Bayes estimators and the PESLs

by the oracle method. However, the hyperparameters are unknown. The oracle method knows the hyperparameters only in simulations. But in reality, the hyperparameters are unknown.

The good news is that we can actually obtain the Bayes estimators and the PESLs

by the moment method, and

by the MLE method, once the data are given.

Similar to Sun et al. (2021), the Bayes estimators , the PESL functions , and the PESLs are depicted in figure 1.2. From the figure, we see that the Bayes estimators are the minimizers of the PESL functions , and the PESLs are the PESL functions evaluated at the Bayes estimators . Note that for the Bayes estimators, and maybe larger or smaller than . Similarly, for the PESLs, and maybe larger or smaller than .

FIG. 1.2 — I: The Bayes estimators , the PESL functions , and the PESLs . (a) The Bayes estimators and the PESLs by the oracle method are larger. (b) The Bayes estimators and the PESLs by the oracle method are smaller.

1.8　Simulation Techniques

1.8.1　Consistencies of the Moment Estimators and the MLEs

In this subsection, taking the hyperparameters from Sun et al. (2021) as an example, we will introduce a simulation technique to numerically exemplify that the moment estimators (, , and ) and the MLEs (, , and ) are consistent estimators of the hyperparameters (, , and ) of the hierarchical inverse gamma and inverse gamma model (2.1). Note that only are used in this subsection.

Let denote the hyperparameter , , or . Then the consistency means that

for , where is for the moment estimator, is for the MLE, means convergence in probability, and is the sample size. Alternatively, the consistency means that

(1.22)

for every and every . The probabilities are approximated by the corresponding frequencies:

for , where is the indicator function of , which is equal to 1 if is true and otherwise, and is the number of simulations. Therefore, the frequencies (, , , or , ) tend to means that the estimators are consistent.

1.8.2　Goodness-of-Fit of the Model

In this subsection, we will introduce two simulation techniques to calculate the goodness-of-fit of the hierarchical model to the simulated data (see Ross (2013); Xue and Chen (2007)). The first simulation technique is the chi-square test, which is mainly used for discrete distributions. The second simulation technique is the Kolmogorov-Smirnov (KS) test, which is mainly used for continuous distributions. Note that only are used in this subsection.

Let us take the hierarchical Poisson and gamma model (8.1) as an example to illustrate the process of the chi-square test. Two cases of the goodness-of-fit will be considered. In the first case, the hyperparameters and are assumed known. In the second case, the hyperparameters and are unknown, and this is also the case encountered in real applications.

Case 1. The hyperparameters and are assumed known.

In this case, the hyperparameters and are assumed known. For example, and . Let the null hypothesis be

where is the marginal distribution of the hierarchical Poisson and gamma model (8.1) with the marginal pmf given by (8.3).

The chi-square goodness-of-fit is performed as follows. We first divide the domain of , , into groups:

Let the theoretical probabilities under on these subintervals be

where

and is the probability when is distributed under . Let , denote the number of that lie in the th subinterval . Then the chi-square statistics (Ross (2013); Xue and Chen (2007))

where is convergence in distribution. Moreover, we can compute the p-value, which gives the probability that a value of as large as would have occurred if the null hypothesis were true. Hence,

where pchisq(), which calculates the cumulative distribution function (cdf) of a chi-square random variable, is an R built-in function (R Core Team (2023)). Note that a large p-value (>0.05 in the usual case) indicates that the model specified by fits the (simulated) data well, while a small p-value (≤0.05 in the usual case) indicates that the model specified by does not fit the (simulated) data well. The larger the p-value, the better the model specified by fits the (simulated) data.

Case 2. The hyperparameters and are unknown.

Let the null hypothesis be

where and are unknown. First, the hyperparameters and need to be estimated by the sample. The estimators could be the moment estimators or the MLEs. Let and be given in Case 1. The theoretical probabilities under on the subintervals are calculated by

that is, the unknown hyperparameters and are estimated by their estimators and based on the sample. Then the chi-square statistics (Ross (2013); Xue and Chen (2007))

Note that the degree of freedom is now lost by 2, since two unknown parameters are estimated by the sample. Moreover, the p-value is given by

The Kolmogorov-Smirnov statistic is the distance between the empirical cdf and the population cdf , that is,

(1.23)

The return value of ks.test() is a list containing the components statistic (the value of the test statistic, or the value) and p-value (the p-value of the test). It is well known in the literature that a smaller value or a larger p-value indicates a better fit of the model to the simulated data. Inspired by and based on the value and p-value, Sun et al. (2021) propose five indices (, , , , and ) to compare the three methods, namely, the oracle method, the moment method, and the MLE method in simulations. is the average values (1.23) of simulations (the smaller the better). is the average p-value of simulations (the larger the better). is the percentage of simulations that attain the minimum value in the three methods (the larger the better). The three values should sum to . is the percentage of simulations which attain the maximum p-value in the three methods (the larger the better). The three values should sum to . is the percentage of accepting (defined as p-value >0.05) in simulations for each method (the larger the better). Each value should between 0%–100%.

1.8.3　Numerical Comparisons of the Bayes Estimators and the PESLs of Three Methods

We will only give the mathematical formulas of the averages and proportions of the absolute errors from the oracle method by the moment method and the MLE method for the Bayes estimator under Stein’s loss function . These quantities for the estimators of the hyperparameters, the Bayes estimator , and the PESLs and are similar, and thus they are omitted.

To calculate the averages and proportions of the absolute errors from the oracle method by the moment method, and the MLE method for , we need the absolute error vectors by the moment method and the MLE method for , which are respectively given by

and

The averages of the absolute errors from the oracle method by the moment method, and the MLE method for are the sample averages of the absolute error vectors and , and the averages are respectively given by

and

With the two absolute error vectors and , we can compute a parallel minima vector of the absolute error vectors by

where pmin() is an R built-in function which returns a single vector giving the parallel minima of the argument vectors. Finally, the proportions and of the absolute errors from the oracle method by the moment method, and the MLE method for are computed by

and

where is the indicator function of which is equal to if is true and otherwise. To avoid the case of equal absolute errors from the oracle method by the moment method, and the MLE method for , we could compute

to ensure that the two proportions sum to .

Let be a hyperparameter. Let be an estimator of θ, where is for the moment estimator and is for the MLE. The MSE of the estimator is defined by (see Casella and Berger (2002))

Similarly, the MAE of the estimator is defined by (see Casella and Berger (2002))

Moreover, the MEE of the estimator is defined by

where the entropy loss function is also known as Stein’s loss function, given by (1.10).

1.9　R Codes

Chapter 2　The Empirical Bayes Estimators of the Rate Parameter of the Inverse Gamma Distribution with a Conjugate Inverse Gamma Prior under Stein’s Loss Function

2.1　Introduction

Since the rate parameter of the inverse gamma distribution with a conjugate inverse gamma prior is a positive parameter, , the Bayes estimator of under the squared error loss function, is not appropriate. In contrast, we should select Stein’s loss function because it penalizes gross overestimation and gross underestimation equally. To determine the unknown hyperparameters, we adopt the empirical Bayes method. In this chapter, we calculate the Bayes estimator of under Stein’s loss function and the corresponding PESL. We also obtain the Bayes estimator of under the squared error loss function and the corresponding PESL. Moreover, we obtain the empirical Bayes estimators of the rate parameter of the inverse gamma distribution with a conjugate inverse gamma prior by the moment method and the MLE method under Stein’s loss function.

2.2　Theoretical Results

Suppose that we observe from the hierarchical inverse gamma and inverse gamma model:

(2.1)

where , , and are hyperparameters to be estimated, is the unknown parameter of interest, is the inverse gamma distribution with shape parameter and rate parameter , and is the inverse gamma distribution with shape parameter and rate parameter . The pdf of the inverse gamma distribution can be found in section 1.2. As described in Deely and Lindley (1981), the statistician observes data and wishes to make an inference about . Therefore, provides direct information about the parameter , while supplementary information is also available. The connection between the prime data and the supplementary information is provided by the common distributions and .

Now we give the justifications of why is the only parameter of interest. The pdf of is

for and . We can only handle the case of unknown rate parameter , letting the shape parameter be a hyperparameter to be determined. If the shape parameter is also an unknown parameter of interest, then we have to deal with the part in the posterior distribution, which is very complicated and has no analytical solutions. It seems that the Bayesian community avoids dealing with such a situation by letting be a known constant or assuming be a hyperparameter to be determined.

2.2.1　The Bayes Estimators and the PESLs

For the hierarchical inverse gamma and inverse gamma model (2.1), the posterior density of and the marginal density of are given by the following theorem, whose proof can be found in appendix A.1.

Theorem 2.1. For the hierarchical inverse gamma and inverse gamma model (2.1), the posterior density of is

where

(2.2)

Moreover, the marginal density of is

(2.3)

for and .

From Theorem 2.1, we have

Since is a rate parameter of the inverse gamma distribution, the Bayes estimator of under the squared error loss function is not appropriate. In contrast, we should choose Stein’s loss function (James and Stein (1961); see also Brown (1990)) because it penalizes gross overestimation and gross underestimation equally. The Bayes estimator of under Stein’s loss function is given by (see Zhang (2017))

(2.4)

where and are given by (2.2). Moreover, we also calculate the Bayes estimator of under the usual squared error loss function,

(2.5)

It is easy to show that

(2.6)

which exemplifies the theoretical study of (1.16). Furthermore, from Zhang (2017), the PESLs at and are respectively given by

(2.7)

and

(2.8)

where

is the digamma function. It is easy to show that

(2.9)

which exemplifies the theoretical study of (1.17). The numerical simulations will exemplify (2.6) and (2.9).

It is worth noting that the Bayes estimators ( and ) and the PESLs and in this subsection assume that the hyperparameters , , and are known. In other words, the Bayes estimators and the PESLs in this subsection are calculated by the oracle method, which will be further discussed in subsection 2.2.3.

2.2.2　The Empirical Bayes Estimators of θ_n+1

To obtain the empirical Bayes estimators of , we need to estimate the hyperparameters from the supplementary information . There are two common methods to estimate the hyperparameters: the moment method and the MLE method.

The estimators of the hyperparameters of the model (2.1) by the moment method , , and and their consistencies are summarized in the following theorem, whose proof can be found in appendix A.2.

Theorem 2.2. The estimators of the hyperparameters of the model (2.1) by the moment method are

(2.10)

(2.11)

(2.12)

where , , is the sample kth moment of . Moreover, the moment estimators are consistent estimators of the hyperparameters.

The estimators of the hyperparameters of the model (2.1) by the MLE method , , and and their consistencies are summarized in the following theorem whose proof can be found in appendix A.3.

Theorem 2.3. The estimators of the hyperparameters of the model (2.1) by the MLE method , , and are the solutions to the following equations:

(2.13)

(2.14)

(2.15)

The analytical calculations of the MLEs of , , and by solving the equations (2.13)–(2.15) are impossible, and thus we have to resort to numerical solutions. We can exploit Newton’s method to solve the equations equations (2.13)–(2.15) and to obtain the MLEs of , , and . Note that the MLEs of , , and are very sensitive to the initial estimators, and the moment estimators are usually proved to be good initial estimators.

Theorem 2.4. The empirical Bayes estimator of the parameter of the model (2.1) under Stein’s loss function by the moment method is given by (2.4) with the hyperparameters estimated by in Theorem 2.2. Alternatively, the empirical Bayes estimator of the parameter of the model (2.1) under Stein’s loss function by the MLE method is given by (2.4) with the hyperparameters estimated by numerically determined in Theorem 2.3.

2.2.3　Theoretical Comparisons of the Bayes Estimators and the PESLs of Three Methods

Note that the subscripts 0, 1, and 2 below are for the oracle method, the moment method, and the MLE method, respectively. The PESLs of the three methods are respectively given by

where

, , and are unknown hyperparameters, , , and are the moment estimators of the hyperparameters given in Theorem 2.2, and , , and are the MLEs of the hyperparameters numerically determined in Theorem 2.3.

The Bayes estimators of under Stein’s loss function, are given by

and

The Bayes estimators of under the squared error loss function are given by

for .

The PESLs evaluated at the Bayes estimators are given by

The PESLs evaluated at the Bayes estimators are given by

2.3　Simulations

The simulated data are generated according to model (2.1) with the hyperparameters specified by , , and . The reason why we choose these values is that , , and are required in the moment estimations of the hyperparameters. Other numerical values of the hyperparameters can also be specified.

2.3.1　Two Inequalities of the Bayes Estimators and the PESLs

First, we fix , , and . Then we set a seed number 1 in R software and draw from . After that, we draw from . Figure 2.1 shows the histogram of and the density estimation curve of . It is that we find to minimize the PESL. Numerical results show that

and

which exemplify the theoretical studies of (2.6) and (2.9).

FIG. 2.1 — IG-IG: The histogram of and the density estimation curve of .

In figure 2.2, we fix , , and , but allow to change from 1 to 10. From the figure, we see that the Bayes estimators and PESLs are functions of . The numerical values of the Bayes estimators and the PESLs in the figure are displayed in table 2.1. We see from plot (a) or the first two lines of table 2.1 that the Bayes estimators are decreasing functions of , and are unanimously smaller than , and thus (2.6) is exemplified. Plot (b) or the last two lines of table 2.1 exhibit that the PESLs do not depend on , and are unanimously smaller than , and thus (2.9) is exemplified.

FIG. 2.2 — IG-IG: The Bayes estimators and the PESLs as functions of . (a) Bayes estimators. (b) PESLs.

TAB. 2.1 — IG-IG: The numerical values of the Bayes estimators and the PESLs in figure 2.2: changes.

	1	2	3	4	5	6	7	8	9	10
	0.2143	0.1429	0.1190	0.1071	0.1000	0.0952	0.0918	0.0893	0.0873	0.0857
	0.2500	0.1667	0.1389	0.1250	0.1167	0.1111	0.1071	0.1042	0.1019	0.1000
	0.0731	0.0731	0.0731	0.0731	0.0731	0.0731	0.0731	0.0731	0.0731	0.0731
	0.0856	0.0856	0.0856	0.0856	0.0856	0.0856	0.0856	0.0856	0.0856	0.0856

Now we allow one of the three parameters , , and to change, holding other parameters fixed. Moreover, we also assume that the datum is fixed, as is the case for the real data. Figure 2.3 shows the Bayes estimators and the PESLs as functions of , , and . We see from the left plots of the figure that the Bayes estimators depend on , , and , and (2.6) is exemplified. The right plots of the figure exhibit that the PESLs depend on and , but not on , and (2.9) is exemplified. Furthermore, tables 2.2–2.4 display the numerical values of the Bayes estimators and the PESLs in figure 2.3. In summary, the results of figure 2.3 and tables 2.2–2.4 exemplify the theoretical studies of (2.6) and (2.9).

FIG. 2.3 — IG-IG: The Bayes estimators and the PESLs as functions of , , and . (a) Bayes estimators vs. . (b) PESLs vs. . (c) Bayes estimators vs. . (d) PESLs vs. . (e) Bayes estimators vs. . (f) PESLs vs. .

TAB. 2.2 — IIG-IG: The numerical values of the Bayes estimators and the PESLs in figure 2.3: changes.

	1	2	3	4	5	6	7	8	9	10
	2.8077	2.4066	2.1058	1.8718	1.6846	1.5315	1.4038	1.2958	1.2033	1.1231
	3.3692	2.8077	2.4066	2.1058	1.8718	1.6846	1.5315	1.4038	1.2958	1.2033
	0.0856	0.0731	0.0638	0.0566	0.0508	0.0461	0.0422	0.0390	0.0361	0.0337
	0.1033	0.0856	0.0731	0.0638	0.0566	0.0508	0.0461	0.0422	0.0390	0.0361

TAB. 2.3 — IIG-IG: The numerical values of the Bayes estimators and the PESLs in figure 2.3: changes.

	1	2	3	4	5	6	7	8	9	10
	2.4780	2.4066	2.3828	2.3709	2.3637	2.3590	2.3556	2.3530	2.3510	2.3494
	2.8910	2.8077	2.7799	2.7660	2.7577	2.7521	2.7481	2.7452	2.7429	2.7410
	0.0731	0.0731	0.0731	0.0731	0.0731	0.0731	0.0731	0.0731	0.0731	0.0731
	0.0856	0.0856	0.0856	0.0856	0.0856	0.0856	0.0856	0.0856	0.0856	0.0856

TAB. 2.4 — IIG-IG: The numerical values of the Bayes estimators and the PESLs in figure 2.3: changes.

	5	6	7	8	9	10	11	12	13	14
	6.7943	5.9450	5.2844	4.7560	4.3236	3.9633	3.6585	3.3971	3.1707	2.9725
	7.9267	6.7943	5.9450	5.2844	4.7560	4.3236	3.9633	3.6585	3.3971	3.1707
	0.0731	0.0638	0.0566	0.0508	0.0461	0.0422	0.0390	0.0361	0.0337	0.0316
	0.0856	0.0731	0.0638	0.0566	0.0508	0.0461	0.0422	0.0390	0.0361	0.0337

Since the Bayes estimators and and the PESLs and depend on and , where and , we can plot the surfaces of the Bayes estimators and the PESLs on the domain via the R function persp3d() in the R package rgl (see Adler et al. (2017)). We remark that the R function persp() in the R package graphics can not add another surface to the existing surface, but persp3d() can. Moreover, persp3d() allows one to rotate the perspective plots of the surface according to one’s wishes. Figure 2.4 plots the surfaces of the Bayes estimators and the PESLs, and the surfaces of the differences of the Bayes estimators and the PESLs. The domain for is for all the plots. a is for and b is for in the axes of all the plots. The red surface is for and the blue surface is for in the upper two plots. From the left two plots of the figure, we see that for all on . From the right two plots of the figure, we see that for all on . The results of the figure exemplify the theoretical studies of (2.6) and (2.9).

FIG. 2.4 — IG-IG: (a) The Bayes estimators as functions of and . (b) The PESLs as functions of and . (c) The surface of which is positive for all on . (d) The surface of which is also positive for all on .

2.3.2　Consistencies of the Moment Estimators and the MLEs

In this subsection, similar to subsection 1.8.1, we will numerically exemplify that the moment estimators and the MLEs are consistent estimators of the hyperparameters , , and of the hierarchical inverse gamma and inverse gamma model (2.1). The motivation of this subsection is that in Theorems 2.2 and 2.3, we have theoretically shown that the moment estimators and the MLEs are consistent estimators of the hyperparameters. Note that only are used in this subsection. The simulation design of this subsection is detailed in appendix A.4.

The frequencies of the moment estimators (, , and ) and the MLEs (, , and ) of the hyperparameters (, , and ) as varies for and , 0.5, and 0.1 are reported in table 2.5. From this table, we observe the following facts.

1. Given , 0.5, or 0.1, the frequencies of the estimators (, , , or , ) tend to 0 as increases to infinity, which means that the moment estimators and the MLEs are consistent estimators of the hyperparameters. For , the frequencies of the estimators are still very large. However, we observe the tendencies of declining to 0 as increases to infinity.

2. Comparing the frequencies corresponding to , 0.5, and 0.1, we observe that as gets smaller, the frequencies tend to be larger, since the constraints

are easier to meet.

3. Comparing the moment estimators and the MLEs of the hyperparameters , , and , we see that the frequencies of the MLEs are smaller than those of the moment estimators, which means that the MLEs are better than the moment estimators when estimating the hyperparameters in terms of consistency.

TAB. 2.5 — IIG-IG: The frequencies of the moment estimators and the MLEs of the hyperparameters as varies for and , 0.5, and 0.1.

		Moment estimators			MLEs
	n
	1e4	0.01	0.25	0.33	0.01	0.03	0.04
	2e4	0.00	0.09	0.12	0.00	0.00	0.01
	4e4	0.00	0.01	0.07	0.00	0.00	0.00
	8e4	0.00	0.00	0.00	0.00	0.00	0.00
	1e4	0.07	0.62	0.63	0.05	0.07	0.12
	2e4	0.05	0.55	0.60	0.02	0.05	0.07
	4e4	0.03	0.45	0.48	0.01	0.06	0.09
	8e4	0.01	0.20	0.27	0.00	0.01	0.02
	1e4	0.82	0.94	0.88	0.14	0.56	0.75
	2e4	0.83	0.92	0.93	0.09	0.40	0.65
	4e4	0.71	0.92	0.90	0.10	0.32	0.53
	8e4	0.62	0.87	0.88	0.07	0.19	0.34

2.3.3　Goodness-of-Fit of the Model: KS Test

In this subsection, similar to subsection 1.8.2, we will calculate the goodness-of-fit of the hierarchical inverse gamma and inverse gamma model (2.1) to the simulated data. The motivation of this subsection is that we want to check whether the marginal distribution of the hierarchical inverse gamma and inverse gamma model (2.1) fits the simulated data well. Note that only are used in this subsection.

In our problem, the null hypothesis specifies that

where is the marginal distribution of the hierarchical inverse gamma and inverse gamma model (2.1). The marginal density of the distribution is given by (2.3), which is obviously one-dimensional continuous. Hence, the KS test can be utilized as a measure of the goodness-of-fit.

The results of the KS test goodness-of-fit of the model (2.1) to the simulated data are reported in table 2.6. Note that the data are simulated according to the hierarchical inverse gamma and inverse gamma model (2.1) with , , and . In the table, the hyperparameters , , and are estimated by three methods. The first method is the oracle method, in that we know the hyperparameters , , and . The second method is the moment method, in that the hyperparameters , , and are estimated by their moment estimators (see Theorem 2.2). The third method is the MLE method, in that the hyperparameters , , and are estimated by their MLEs (see Theorem 2.3). In the table, the sample size is , and the number of simulations is .

1. From table 2.6, we observe the following facts.

The values for the three methods are respectively given by 0.2983, 0.0338, and 0.0255, which means that the MLE method is the best method, the moment method is the second-best method, and the oracle method is the worst method. A possible explanation for this phenomenon is that in (1.23), the empirical cdf is based on data, and the population cdfs for the MLE method and the moment method are also based on data, while the population cdf for the oracle method is not based on data.

2. The values for the three methods are respectively given by 0.0677, 0.4077, and 0.6503, which also means that the MLE method ranks first, the moment method ranks second, and the oracle method ranks third. The possible explanation for this phenomenon is described in the previous paragraph.

3. The values for the three methods are respectively given by 0.02, 0.40, and 0.58. The value for the MLE method accounts for over half of the simulations. The order of preference for the three methods is the MLE method, the moment method, and the oracle method.

4. The values for the three methods are respectively given by 0.02, 0.40, and 0.58. A small value corresponds to a large p-value. Hence, the smallest value corresponds to the largest p-value. Therefore, the value and the value for the three methods are the same. The order of preference for the three methods is the MLE method, the moment method, and the oracle method.

5. The values for the three methods are respectively given by 0.14, 0.81, and 0.91. Once again, the order of preference for the three methods is the MLE method, the moment method, and the oracle method. The values for the moment method and the MLE method are over , which means that the two methods have good performances in terms of goodness-of-fit.

6. In summary, for the five indices (, , , , ), the order of preference for the three methods is the MLE method, the moment method, and the oracle method. Comparing the moment method and the MLE method, we find that the MLE method has a better performance than the moment method in terms of all five indices.

TAB. 2.6 — IIG-IG: The results of the KS test goodness-of-fit of the model (2.1) to the simulated data.

	Oracle method	Moment method	MLE method
	0.2983	0.0338	0.0255
	0.0677	0.4077	0.6503
	0.02	0.40	0.58
	0.02	0.40	0.58
	0.14	0.81	0.91

The boxplots of the values and the p-values for the three methods are displayed in figure 2.5. From the figure, we observe the following facts.

1. The values of the oracle method are significantly larger than those of the other two methods. Since for the value, the smaller the better, the order of preference for the three methods is the MLE method, the moment method, and the oracle method.

2. The p-values of the oracle method are significantly smaller than those of the other two methods. Since for the p-value, the larger the better, the order of preference for the three methods is the MLE method, the moment method, and the oracle method.

3. Small values correspond to large p-values, and large values correspond to small p-values.

4. The MLE method has a better performance than the moment method in terms of the values and the p-values.

FIG. 2.5 — IG-IG: The boxplots of the values and the p-values for the three methods. (a) values. (b) p-values.

2.3.4　Numerical Comparisons of the Bayes Estimators and the PESLs of Three Methods

In this subsection, we will numerically compare the Bayes estimators and the PESLs of the three methods (the oracle method, the moment method, and the MLE method). The motivation of this subsection is that the theoretical comparisons of the Bayes estimators and the PESLs of the three methods can be found in subsection 2.2.3. Note that the full data are used in this subsection.

Note that the data are simulated according to the hierarchical inverse gamma and inverse gamma model (2.1) with , , and . Moreover, the oracle method knows the hyperparameters , , and in simulations.

Comparisons of the Bayes estimators and the PESLs of the three methods for sample size and number of simulations are displayed in figure 2.6. From the figure, we observe the following facts.

1. Plot (a): For the Bayes estimators of under Stein’s loss function , the MLE method is slightly closer to the oracle method than the moment method.

2. Plot (b): For the Bayes estimators of under the squared error loss function , the MLE method is also slightly closer to the oracle method than the moment method.

3. Plot (c): For the PESLs , the MLE method is much closer to the oracle method than the moment method.

4. Plot (d): For the PESLs , the MLE method is also much closer to the oracle method than the moment method.

5. All four plots indicate that the MLE method is better than the moment method, as the Bayes estimators and the PESLs of the MLE method are closer to those of the oracle method than those of the moment method.

The boxplots of the absolute errors from the oracle method by the moment method, and the MLE method for sample size and number of simulations are displayed in figure 2.7. All four plots indicate that the MLE method is better than the moment method, as the absolute errors from the oracle method of the Bayes estimators and the PESLs by the MLE method are much smaller than those of the moment method.

The averages and proportions of the absolute errors from the oracle method by the moment method, and the MLE method for the Bayes estimators and the PESLs are summarized in table 2.7. See Subsection 1.8.3 for details. From the table, we observe that the averages of the absolute errors from the oracle method by the MLE method are much smaller than those by the moment method. Moreover, the proportions of the absolute errors from the oracle method by the MLE method are much larger than those by the moment method. In summary, the table illustrates that the MLE method is better than the moment method in terms of the averages and proportions of the absolute errors from the oracle method.

FIG. 2.6 — IG-IG: Comparisons of the Bayes estimators and the PESLs of the three methods for sample size and number of simulations . (a) . (b) . (c) . (d) .

FIG. 2.7 — IG-IG: The boxplots of the absolute errors from the oracle method by the moment method, and the MLE method for sample size and number of simulations . (a) Absolute errors for . (b) Absolute errors for . (c) Absolute errors for . (d) Absolute errors for .

TAB. 2.7 — IIG-IG: The averages and proportions of the absolute errors from the oracle method by the moment method and the MLE method for the Bayes estimators and the PESLs.

	Averages		Proportions
	Moment	MLE	Moment	MLE
	0.0491	0.0210	0.20	0.80
	0.0626	0.0275	0.21	0.79
	0.0056	0.0022	0.29	0.71
	0.0075	0.0030	0.29	0.71

2.3.5　Marginal Densities for Various Hyperparameters

In this subsection, we will plot the marginal densities of the hierarchical inverse gamma and inverse gamma model (2.1) for various hyperparameters , , and . The motivation of this subsection is that we want to know what kind of data can be modeled by the hierarchical inverse gamma and inverse gamma model (2.1). Note that the marginal density of is given by (2.3) specified by three hyperparameters , , and . We will explore how the marginal densities change around the marginal density with hyperparameters specified by , , and . Other numerical values of the hyperparameters can also be specified.

Figure 2.8 plots the marginal densities for varied , holding and fixed. From the figure, we see that as increases, the peak value of the curve decreases. In other words, the variance of the marginal density increases as

(2.16)

is an increasing function of . Moreover, all the marginal densities are right-skewed.

FIG. 2.8 — IG-IG: The marginal densities for varied , holding and fixed.

Figure 2.9 plots the marginal densities for varied , holding and fixed. From the figure, we also see that as increases, the peak value of the curve decreases. In other words, the variance of the marginal density increases, as (2.16) is an increasing function of . Moreover, all the marginal densities are also right-skewed.

Figure 2.10 plots the marginal densities for varied , holding and fixed. From the figure, we see that as increases, the peak value of the curve increases. In other words, the variance of the marginal density decreases, as (2.16) is a decreasing function of . Moreover, all the marginal densities are also right-skewed.

FIG. 2.9 — IG-IG: The marginal densities for varied , holding and fixed.

2.4　Conclusions and Discussions

For the hierarchical inverse gamma and inverse gamma model (2.1), we calculate the posterior density and the marginal density in Theorem 2.1. Since is a rate parameter in (2.1), the Bayes estimator of under the squared error loss function is not appropriate. In contrast, we should choose Stein’s loss function because it penalizes gross overestimation and gross underestimation equally. After that, we calculate the Bayes estimators of , and , and the PESLs of , and .

FIG. 2.10 — IG-IG: The marginal densities for varied , holding and fixed.

In order to calculate the empirical Bayes estimator of the rate parameter , we must calculate the estimators of the hyperparameters of model (2.1). The estimators of the hyperparameters of model (2.1) by the moment method and their consistencies are summarized in Theorem 2.2. Moreover, the estimators of the hyperparameters of model (2.1) by the MLE method and their consistencies are summarized in Theorem 2.3. Finally, the empirical Bayes estimators of the rate parameter of the model (2.1) under Stein’s loss function by the moment method and the MLE method are summarized in Theorem 2.4.

Note that in Theorem 2.3, we only stated that the estimators of the hyperparameters of model (2.1) by the MLE method , , and are the solutions to the equations (2.13)–(2.15). We can exploit Newton’s method to solve the equations (2.13)–(2.15) and to numerically obtain the MLEs of , , and . However, we can not prove the existence and uniqueness of the solutions to our system. The interested readers who have such kind of knowledge and skills are encouraged to solve this issue.

Numerical simulations illustrate that the moment estimators (, , and ) and the MLEs (, , and ) are consistent estimators of the hyperparameters (, , and ), as reported in table 2.5. Moreover, table 2.6 indicates that the hierarchical inverse gamma and inverse gamma model (2.1) fits the simulated data well in terms of the KS test goodness-of-fit by the moment method and the MLE method.

Chapter 3　The Empirical Bayes Estimators of the Rate Parameter of the Gamma Distribution with a Conjugate Gamma Prior under Stein’s Loss Function

3.1　Introduction

The hierarchical gamma and gamma model (3.1) has been considered in table 3.3.1 (p. 121) and table 4.2.1 (p. 176) of Robert (2007). However, he only calculated the Bayes estimator of under the squared error loss function. Since the rate parameter is a positive parameter, the Bayes estimator of under the squared error loss function is not appropriate. In contrast, we should choose Stein’s loss function (James and Stein (1961); see also Brown (1990)) because it penalizes gross overestimation and gross underestimation equally. To determine the unknown hyperparameters, we adopt the empirical Bayes method. In this chapter, we calculate the Bayes estimator of under Stein’s loss function and the corresponding PESL. We also obtain the Bayes estimator of under the squared error loss function and the corresponding PESL. Moreover, we obtain the empirical Bayes estimators of the rate parameter of the gamma distribution with a conjugate gamma prior by the moment method and the MLE method under Stein’s loss function.

3.2　Theoretical Results

Suppose that we observe from the hierarchical gamma and gamma model:

(3.1)

where , , and are hyperparameters to be estimated, is the unknown parameter of interest, is the gamma distribution with shape parameter and rate parameter , and is the gamma distribution with shape parameter and rate parameter . As described in Deely and Lindley (1981), the statistician observes data and wishes to make an inference about . Therefore, provides direct information about the parameter , while supplementary information is also available. The connection between the prime data and the supplementary information is provided by the common distributions and . The pdfs of and can be found in section 1.2.

Now we give the justifications of why is the only parameter of interest. The pdf of is

for and . We can only handle the case of an unknown rate parameter , letting the shape parameter be a hyperparameter to be determined. If the shape parameter is also an unknown parameter of interest, then we have to deal with the part in the posterior distribution, which is very complicated and has no analytical solutions. It seems that the Bayesian community avoids dealing with such a situation by letting be a known constant or assuming be a hyperparameter to be determined.

3.2.1　The Bayes Estimators and the PESLs

For the hierarchical gamma and gamma model (3.1), the posterior density of and the marginal density of are given by the following theorem, whose proof can be found in appendix A.5.

Theorem 3.1. For the hierarchical gamma and gamma model (3.1), the posterior density of is

where

(3.2)

Moreover, the marginal density of is

(3.3)

for and .

From Theorem 3.1, we have

Since is a rate parameter of the gamma distribution, the Bayes estimator of under the squared error loss function is not appropriate. In contrast, we should choose Stein’s loss function because it penalizes gross overestimation and gross underestimation equally. From (1.12), the Bayes estimator of under Stein’s loss function is given by

(3.4)

for , where and are given by (3.2). Moreover, from (1.13), the Bayes estimator of under the usual squared error loss function is given by

(3.5)

It is easy to see that

(3.6)

which exemplifies the theoretical study of (1.16). Furthermore, from (1.14) and (1.15), the PESLs at and are respectively given by

(3.7)

and

(3.8)

where

is the digamma function. It is worth noting that the two PESLs and ) depend on , which is given by

The calculations of and the two PESLs can be found in appendix A.6. It is easy to show that

(3.9)

which exemplifies the theoretical study of (1.17). The numerical simulations will exemplify (3.6) and (3.9).

It is worth noting that the Bayes estimators and the PESLs in this subsection assume that the hyperparameters , , and are known. In other words, the Bayes estimators and the PESLs in this subsection are calculated by the oracle method, which will be further discussed in subsection 3.2.3.

3.2.2　The Empirical Bayes Estimators of θ_n+1

To obtain the empirical Bayes estimators of , we need to estimate the hyperparameters from the supplementary information . There are two common methods to estimate the hyperparameters: the moment method and the MLE method.

The estimators of the hyperparameters of the model (3.1) by the moment method , , and and their consistencies are summarized in the following theorem, whose proof can be found in appendix A.7.

Theorem 3.2. The estimators of the hyperparameters of the model (3.1) by the moment method are

(3.10)

(3.11)

(3.12)

where , , is the sample th moment of . Moreover, the moment estimators are consistent estimators of the hyperparameters.

The estimators of the hyperparameters of the model (3.1) by the MLE method , , and and their consistencies are summarized in the following theorem whose proof can be found in appendix A.8.

Theorem 3.3. The estimators of the hyperparameters of the model (3.1) by the MLE method , , and are the solutions to the following equations:

(3.13)

(3.14)

(3.15)

The analytical calculations of the MLEs of , , and by solving the equations (3.13)–(3.15) are impossible, and thus we have to resort to numerical solutions. We can exploit Newton’s method to solve the equations (3.13)–(3.15) and to obtain the MLEs of , , and . Note that the MLEs of , , and are very sensitive to the initial estimators, and the moment estimators are usually proven to be good initial estimators.

Theorem 3.4. The empirical Bayes estimator of the parameter of the model (3.1) under Stein’s loss function by the moment method is given by (3.4) with the hyperparameters estimated by in Theorem 3.2. Alternatively, the empirical Bayes estimator of the parameter of the model (3.1) under Stein’s loss function by the MLE method is given by (3.4) with the hyperparameters estimated by numerically determined in Theorem 3.3

3.2.3　Theoretical Comparisons of the Bayes Estimators and the PESLs of Three Methods

Note that the subscripts 0, 1, and 2 below are for the oracle method, the moment method, and the MLE method, respectively. The PESL functions of the three methods are respectively given by

where

, , and are unknown hyperparameters, , , and given in Theorem 3.2 are the moment estimators of the hyperparameters, and , , and numerically determined in Theorem 3.3 are the MLEs of the hyperparameters.

The Bayes estimators of under Stein’s loss function are given by

The Bayes estimators of under the squared error loss function are given by

The PESLs evaluated at the Bayes estimators are given by

The PESLs evaluated at the Bayes estimators are given by

3.3　Simulations

The simulated data are generated according to the model (3.1) with the hyperparameters specified by , , and . The reason why we choose these values is that , , and are required in the moment estimations of the hyperparameters. Other numerical values of the hyperparameters can also be specified.

3.3.1　Two Inequalities of the Bayes Estimators and the PESLs

First, we fix , , and . Then we set a seed number 1 in R software and draw from . After that, we draw from . Figure 3.1 shows the histogram of and the density estimation curve of . It is that we find to minimize the PESL. Numerical results show that

and

which exemplify the theoretical studies of (3.6) and (3.9).

FIG. 3.1 — G-G: The histogram of and the density estimation curve of .

In figure 3.2, we fix , , and , but allow to change from 1 to 10. From the figure, we see that the Bayes estimators and the PESLs are functions of . The numerical values of the Bayes estimators and the PESLs in figure 3.2 are displayed in table 3.1. We see from plot (a) or the first two lines of table 3.1 that the Bayes estimators are decreasing functions of , and are unanimously smaller than , and thus (3.6) is exemplified. Plot (b) or the last two lines of table 3.1 exhibit that the PESLs do not depend on , and are unanimously smaller than , and thus (3.9) is exemplified.

FIG. 3.2 — G-G: The Bayes estimators and the PESLs as functions of . (a) Bayes estimators. (b) PESLs.

TAB. 3.1 — G-G: The numerical values of the Bayes estimators and the PESLs in figure 3.2: changes.

	1	2	3	4	5	6	7	8	9	10
	1.6667	1.2500	1.0000	0.8333	0.7143	0.6250	0.5556	0.5000	0.4545	0.4167
	2.0000	1.5000	1.2000	1.0000	0.8571	0.7500	0.6667	0.6000	0.5455	0.5000
	0.0967	0.0967	0.0967	0.0967	0.0967	0.0967	0.0967	0.0967	0.0967	0.0967
	0.1144	0.1144	0.1144	0.1144	0.1144	0.1144	0.1144	0.1144	0.1144	0.1144

Now we allow one of the three parameters , , and to change, holding other parameters fixed. Moreover, we also assume that the datum is fixed, as is the case for the real data. Figure 3.3 shows the Bayes estimators and the PESLs as functions of , , and . We see from the left plots of the figure that the Bayes estimators depend on , , and , and (3.6) is exemplified. Moreover, the Bayes estimators are increasing functions of and , and they are decreasing functions of . The right plots of the figure exhibit that the PESLs depend on and , but not on , and (3.9) is exemplified. In addition, the PESLs are decreasing functions of and . Furthermore, tables 3.2–3.4 display the numerical values of the Bayes estimators and the PESLs in figure 3.3. In summary, the results of figure 3.3 and tables 3.2–3.4 exemplify the theoretical studies of (3.6) and (3.9).

TAB. 3.2 — G-G: The numerical values of the Bayes estimators and the PESLs in figure 3.3: changes.

	4	5	6	7	8	9	10	11	12	13
	2.0700	2.4840	2.8980	3.3120	3.7260	4.1400	4.5540	4.9681	5.3821	5.7961
	2.4840	2.8980	3.3120	3.7260	4.1400	4.5540	4.9681	5.3821	5.7961	6.2101
	0.0967	0.0810	0.0697	0.0612	0.0545	0.0492	0.0448	0.0411	0.0380	0.0353
	0.1144	0.0935	0.0791	0.0684	0.0603	0.0539	0.0487	0.0444	0.0408	0.0377

Since the Bayes estimators and and the PESLs and depend on and , where and , we can plot the surfaces of the Bayes estimators and the PESLs on the domain via the R function persp3d() in the R package rgl (see Sun et al. (2021); Zhang et al. (2017, 2019b); Adler et al. (2017)). We remark that the R function persp() in the R package graphics can not add another surface to the existing surface, but persp3d() can. Moreover, persp3d() allows one to rotate the perspective plots of the surface according to one’s wishes. Figure 3.4 plots the surfaces of the Bayes estimators and the PESLs, and the surfaces of the differences of the Bayes estimators and the PESLs. The domain for is for all the plots. a is for and b is for in the axes of all the plots. The red surface is for and the blue surface is for in the upper two plots. From the left two plots of the figure, we see that for all on . From the right two plots of the figure, we see that for all on . The results of the figure exemplify the theoretical studies of (3.6) and (3.9).

FIG. 3.3 — G-G: The Bayes estimators and the PESLs as functions of , , and . (a) Bayes estimators vs. . (b) PESLs vs. . (c) Bayes estimators vs. . (d) PESLs vs. . (e) Bayes estimators vs. . (f) PESLs vs. .

TAB. 3.3 — G-G: The numerical values of the Bayes estimators and the PESLs in figure 3.3: changes.

	1	2	3	4	5	6	7	8	9	10
	3.5325	2.0700	1.4639	1.1324	0.9233	0.7794	0.6743	0.5941	0.5310	0.4801
	4.2390	2.4840	1.7567	1.3589	1.1079	0.9352	0.8091	0.7130	0.6373	0.5761
	0.0967	0.0967	0.0967	0.0967	0.0967	0.0967	0.0967	0.0967	0.0967	0.0967
	0.1144	0.1144	0.1144	0.1144	0.1144	0.1144	0.1144	0.1144	0.1144	0.1144

TAB. 3.4 — G-G: The numerical values of the Bayes estimators and the PESLs in figure 3.3: changes.

	1	2	3	4	5	6	7	8	9	10
	1.4624	1.8279	2.1935	2.5591	2.9247	3.2903	3.6559	4.0215	4.3871	4.7527
	1.8279	2.1935	2.5591	2.9247	3.2903	3.6559	4.0215	4.3871	4.7527	5.1182
	0.1198	0.0967	0.0810	0.0697	0.0612	0.0545	0.0492	0.0448	0.0411	0.0380
	0.1467	0.1144	0.0935	0.0791	0.0684	0.0603	0.0539	0.0487	0.0444	0.0408

3.3.2　Consistencies of the Moment Estimators and the MLEs

In this subsection, similar to subsection 1.8.1, we will numerically exemplify that the moment estimators and the MLEs are consistent estimators of the hyperparameters , , and of the hierarchical gamma and gamma model (3.1). The motivation of this subsection is that in Theorems 3.2 and 3.3, we have theoretically shown that the moment estimators and the MLEs are consistent estimators of the hyperparameters. Note that only are used in this subsection.

The frequencies of the moment estimators (, , and ) and the MLEs (, , and ) of the hyperparameters (, , and ) as varies for and , 0.5, and 0.1 are reported in table 3.5. From the table, we observe the following facts.

1. Given , 0.5, or 0.1, the frequencies of the estimators (, , , or , ) tend to 0 as increases to infinity, which means that the moment estimators and the MLEs are consistent estimators of the hyperparameters. For , the frequencies of the estimators are still very large. However, we observe the tendencies of declining to 0 as increases to infinity.

FIG. 3.4 — G-G: (a) The Bayes estimators as functions of and . (b) The PESLs as functions of and . (c) The surface of which is positive for all on . (d) The surface of which is also positive for all on .

2. Comparing the frequencies corresponding to , 0.5, and 0.1, we observe that as gets smaller, the frequencies tend to be larger, since the constraints

are easier to meet.

3. Comparing the moment estimators and the MLEs of the hyperparameters , , and , we see that the frequencies of the MLEs are smaller than those of the moment estimators, which means that the MLEs are better than the moment estimators when estimating the hyperparameters in terms of consistency.

3.3.3　Goodness-of-Fit of the Model: KS Test

In this subsection, similar to subsection 1.8.2, we will calculate the goodness-of-fit of the hierarchical gamma and gamma model (3.1) to the simulated data. The motivation of this subsection is that we want to check whether the marginal distribution of the hierarchical gamma and gamma model (3.1) fits the simulated data well. Note that only are used in this subsection.

In our problem, the null hypothesis specifies that

where is the marginal distribution of the hierarchical gamma and gamma model (3.1). The marginal density of the distribution is given by (3.3), which is obviously one-dimensional continuous. Hence, the KS test can be utilized as a measure of the goodness-of-fit.

The results of the KS test goodness-of-fit of the model (3.1) to the simulated data are reported in table 3.6. Note that the data are simulated according to the hierarchical gamma and gamma model (3.1) with , , and . In the table, the hyperparameters , , and are estimated by three methods. The first method is the oracle method, in that we know the hyperparameters , , and . The second method is the moment method, in that the hyperparameters , , and are estimated by their moment estimators (see Theorem 3.2). The third method is the MLE method, in that the hyperparameters , , and are estimated by their MLEs (see Theorem 3.3). In the table, the sample size is , and the number of simulations is .

TAB. 3.5 — G-G: The frequencies of the moment estimators and the MLEs of the hyperparameters as varies for and , 0.5, and 0.1.

		Moment estimators			MLEs
	1e4	0.27	0.56	0.01	0	0.01	0.01
	2e4	0.21	0.46	0	0	0	0
	4e4	0.11	0.32	0	0	0	0
	8e4	0.04	0.21	0	0	0	0
	1e4	0.70	0.78	0.15	0.04	0.04	0.02
	2e4	0.63	0.75	0.04	0	0.01	0
	4e4	0.50	0.69	0.03	0.01	0.03	0.01
	8e4	0.44	0.60	0.03	0.01	0.01	0
	1e4	0.95	0.99	0.92	0.60	0.55	0.07
	2e4	0.92	0.93	0.91	0.44	0.36	0.01
	4e4	0.93	0.96	0.93	0.27	0.19	0.05
	8e4	0.91	0.93	0.90	0.13	0.08	0.01

From table 3.6, we observe the following facts.

1. The values for the three methods are respectively given by 0.2674, 0.0292, and 0.0088, which means that the MLE method is the best method, the moment method is the second-best method, and the oracle method is the worst method. A possible explanation for this phenomenon is that in (1.23), the empirical cdf is based on data, and the population cdfs for the MLE method and the moment method are also based on data, while the population cdf for the oracle method is not based on data.

2. The values for the three methods are respectively given by 0.0161, 0.0679, and 0.7956, which also means that the MLE method ranks first, the moment method ranks second, and the oracle method ranks third. The possible explanation for this phenomenon is described in the previous paragraph.

3. The values for the three methods are respectively given by 0.01, 0.06, and 0.93. The value for the MLE method accounts for nearly all of the simulations. The order of preference for the three methods is the MLE method, the moment method, and the oracle method.

4. The values for the three methods are respectively given by 0.02, 0.07, and 0.91. The value for the MLE method accounts for nearly all of the simulations. The order of preference for the three methods is the MLE method, the moment method, and the oracle method.

5. The values for the three methods are respectively given by 0.04, 0.20, and 0.95. Once again, the order of preference for the three methods is the MLE method, the moment method, and the oracle method. The value for the MLE method is over , which means that the MLE method has good performance in terms of goodness-of-fit.

6. In summary, for the five indices (, , , , ), the order of preference for the three methods is the MLE method, the moment method, and the oracle method. Comparing the moment method and the MLE method, we find that the MLE method has a better performance than the moment method in terms of all five indices.

TAB. 3.6 — G-G: The results of the KS test goodness-of-fit of the model (3.1) to the simulated data.

	0.2674	0.0292	0.0088
	0.0161	0.0679	0.7956
	0.01	0.06	0.93
%	0.02	0.07	0.91
%	0.04	0.20	0.95

The boxplots of the values and the p-values for the three methods are displayed in figure 3.5. From the figure, we observe the following facts.

1. The values of the oracle method are significantly larger than those of the other two methods. Since for the value, the smaller the better, the order of preference for the three methods is the MLE method, the moment method, and the oracle method.

2. The p-values of the MLE method are significantly larger than those of the other two methods. Since for the p-value, the larger the better, the order of preference for the three methods is the MLE method, the moment method, and the oracle method.

3. Small values correspond to large p-values, and large values correspond to small p-values.

4. The MLE method has a better performance than the moment method in terms of the values and the p-values.

FIG. 3.5 — G-G: The boxplots of the values and the p-values for the three methods. (a) values. (b) p-values.

3.3.4　Numerical Comparisons of the Bayes Estimators and the PESLs of Three Methods

In this subsection, we will numerically compare the Bayes estimators and the PESLs of the three methods (the oracle method, the moment method, and the MLE method). The motivation of this subsection is that the theoretical comparisons of the Bayes estimators and the PESLs of the three methods can be found in subsection 3.2.3. Note that the full data are used in this subsection.

Note that the data are simulated according to the hierarchical gamma and gamma model (3.1) with hyperparameters , , and . Moreover, the oracle method knows the hyperparameters , , and in simulations.

Comparisons of the Bayes estimators and the PESLs of the three methods for sample size and number of simulations are displayed in figure 3.6. From the figure, we observe the following facts.

1. Plot (a): For the Bayes estimators of under Stein’s loss function , the MLE method is slightly closer to the oracle method than the moment method.

2. Plot (b): For the Bayes estimators of under the squared error loss function , the MLE method is also slightly closer to the oracle method than the moment method.

3. Plot (c): For the PESLs , the MLE method is much closer to the oracle method than the moment method.

4. Plot (d): For the PESLs , the MLE method is also much closer to the oracle method than the moment method.

5. All four plots indicate that the MLE method is better than the moment method, as the Bayes estimators and the PESLs of the MLE method are closer to those of the oracle method than those of the moment method.

The boxplots of the absolute errors from the oracle method by the moment method and the MLE method for sample size and number of simulations are displayed in figure 3.7. All four plots indicate that the MLE method is better than the moment method, as the absolute errors from the oracle method of the Bayes estimators and the PESLs by the MLE method are much smaller than those by the moment method.

FIG. 3.6 — G-G: Comparisons of the Bayes estimators and the PESLs of the three methods for sample size and number of simulations . (a) . (b) . (c) . (d) .

FIG. 3.7 — G-G: The boxplots of the absolute errors from the oracle method by the moment method and the MLE method for sample size and number of simulations . (a) Absolute errors for . (b) Absolute errors for . (c) Absolute errors for . (d) Absolute errors for .

The averages and proportions of the absolute errors from the oracle method by the moment method and the MLE method for the Bayes estimators and the PESLs are summarized in table 3.7. See subsection 1.8.3 for details. From the table, we observe that the averages of the absolute errors from the oracle method by the MLE method are much smaller than those by the moment method. Moreover, the proportions of the absolute errors from the oracle method by the MLE method are much larger than those by the moment method. In summary, the table illustrates that the MLE method is better than the moment method in terms of the averages and proportions of the absolute errors from the oracle method.

TAB. 3.7 — G-G: The averages and proportions of the absolute errors from the oracle method by the moment method and the MLE method for the Bayes estimators and the PESLs.

	Averages		Proportions
	Moment	MLE	Moment	MLE
	0.3332	0.0485	0.06	0.94
	0.4136	0.0621	0.06	0.94
	0.0079	0.0025	0.15	0.85
	0.0107	0.0034	0.15	0.85

3.3.5　Marginal Densities for Various Hyperparameters

In this subsection, we will plot the marginal densities of the hierarchical gamma and gamma model (3.1) for various hyperparameters , , and . The motivation of this subsection is that we want to know what kind of data can be modeled by the hierarchical gamma and gamma model (3.1). Note that the marginal density of is given by (3.3) specified by three hyperparameters , , and . We will explore how the marginal densities change around the marginal density with hyperparameters specified by , , and . Other numerical values of the hyperparameters can also be specified.

Figure 3.8 plots the marginal densities for varied , holding and fixed. From the figure, we see that as increases, the peak value of the curve increases. In other words, the variance of the marginal density decreases, as

(3.16)

is a decreasing function of . Moreover, all the marginal densities are right-skewed.

FIG. 3.8 — G-G: The marginal densities for varied , holding and fixed.

Figure 3.9 plots the marginal densities for varied , holding and fixed. From the figure, we see that as increases, the peak value of the curve decreases. In other words, the variance of the marginal density increases, as (3.16) is an increasing function of . Moreover, all the marginal densities are also right-skewed.

FIG. 3.9 — G-G: The marginal densities for varied , holding and fixed.

Figure 3.10 plots the marginal densities for varied , holding and fixed. From the figure, we see that as increases, the peak value of the curve decreases. In other words, the variance of the marginal density increases, as (3.16) is an increasing function of . Moreover, all the marginal densities are also right-skewed.

3.4　Conclusions and Discussions

For the hierarchical gamma and gamma model (3.1), we calculate the posterior density and the marginal density in Theorem 3.1. Since is a rate parameter in (3.1), the Bayes estimator of under the squared error loss function is not appropriate. In contrast, we should choose Stein’s loss function because it penalizes gross overestimation and gross underestimation equally. After that, we calculate the Bayes estimators of , and , and the PESLs of , and .

FIG. 3.10 — G-G: The marginal densities for varied , holding and fixed.

In order to calculate the empirical Bayes estimator of the rate parameter , we must calculate the estimators of the hyperparameters of model (3.1). The estimators of the hyperparameters of model (3.1) by the moment method and their consistencies are summarized in Theorem 3.2. Moreover, the estimators of the hyperparameters of model (3.1) by the MLE method and their consistencies are summarized in Theorem 3.3. Finally, the empirical Bayes estimators of the rate parameter of the model (3.1) under Stein’s loss function by the moment method and the MLE method are summarized in Theorem 3.4.

Note that in Theorem 3.3, we only stated that the estimators of the hyperparameters of the model (3.1) by the MLE method , , and are the solutions to the equations (3.13)–(3.15). We can exploit Newton’s method to solve the equations (3.13)–(3.15), and to numerically obtain the MLEs of , , and . However, we can not prove the existence and uniqueness of the solutions to our system. The interested readers who have such knowledge and skills are encouraged to solve this issue.

Chapter 4　The Empirical Bayes Estimators of the Mean Parameter of the Exponential Distribution with a Conjugate Inverse Gamma Prior under Stein’s Loss Function

4.1　Introduction

In the hierarchical exponential and inverse gamma model (4.1), our parameter of interest is the mean which is a positive parameter. Therefore, we will choose Stein’s loss function because it penalizes gross overestimation and gross underestimation equally.

4.2　Theoretical Results

Suppose that we observe from the hierarchical exponential and inverse gamma model:

(4.1)

where and are hyperparameters to be estimated, is the unknown parameter of interest, is the exponential distribution with mean parameter , and is the inverse gamma distribution with shape parameter and scale parameter β. As described in Deely and Lindley (1981), the statistician observes data and wishes to make an inference about . Therefore, provides direct information about the parameter , while supplementary information is also available. The connection between the prime data and the supplementary information is provided by the common distributions and . The pdfs of and can be found in section 1.2.

There are two pdf forms for the exponential distribution. One form uses the mean (or scale) as the parameter (see Casella and Berger (2002)), and with pdf , for and . Another form utilizes the rate as the parameter (see Gelman et al. (2013); Mao and Tang (2012)), and with pdf , for and . The two pdfs are the same with a relationship of the parameters .

The exponential-gamma model is assumed to generate observations :

(4.2)

Firstly, the exponential-inverse gamma model (4.1) and the exponential-gamma model (4.2) are equivalent in the sense of their marginal pdfs. For convenience, now let , , , , , and be random variables having the corresponding distributions. For example, is a random variable having the distribution. It is easy to see that

Therefore, the two hierarchical models (4.1) and (4.2) are equivalent. It is straightforward to derive that the two marginal pdfs of the two hierarchical models are the same, and they are equal to

Since the two marginal pdfs are the same, the moment estimators (displayed in Theorem 4.2) and the MLEs (see Theorem 4.3) of the hyperparameters and for the two hierarchical models (4.1) and (4.2) are the same.

Another reason to use (4.1) is that it motivates us to consider 16 hierarchical models of the gamma and inverse gamma distributions (see Zhang and Zhang (2022)). It is easy to see that is a conjugate prior for the distribution. Writing in the form of likelihood-prior, it is . Similarly, the is a conjugate prior for the distribution. Writing in the form of likelihood-prior, it is . The expression motivates us to consider , , , and as the likelihood, and , , , and as the prior, leading to 16 combinations of the likelihood-prior.

4.2.1　The Bayes Estimators and the PESLs

The posterior distribution of and the marginal density of of the hierarchical exponential and inverse gamma model (4.1) are summarized in the following theorem whose proof can be found in appendix A.9.

Theorem 4.1. For the hierarchical exponential and inverse gamma model (4.1), the posterior distribution of is

(4.3)

where

(4.4)

Moreover, the marginal density of is

(4.5)

for and .

Now, let us analytically calculate the Bayes estimators and , and the PESLs and under the hierarchical exponential and inverse gamma model (4.1).

From (4.3), we have

From (1.12), the Bayes estimator of under Stein’s loss function, it is given by

(4.6)

where and are given by (4.4). Moreover, from (1.13), the Bayes posterior estimator of under the usual squared error loss function, it is given by

(4.7)

for . It is easy to see that

(4.8)

which exemplifies the theoretical study of (1.16).

To analytically calculate the PESLs and , we need to analytically calculate . For the sake of simplicity, the *’s are dropped from and . We have

where

is the digamma function.

Therefore, from (1.14) and (1.15), after some algebraic operations, the PESLs at and are respectively given by

(4.9)

and

(4.10)

for . It is easy to show that

(4.11)

which exemplifies the theoretical study of (1.17). It is worth noting that the PESLs and depend only on , but not on . Therefore, the PESLs depend only on , but not on and .

In the simulations section and the real data section, we will exemplify the two inequalities (4.8) and (4.11). Moreover, we will exemplify that the PESLs depend only on , but not on and .

It is worth noting that the Bayes estimators and the PESLs in this subsection assume that the hyperparameters and are known. In other words, the Bayes estimators and the PESLs in this subsection are calculated by the oracle method, which will be further discussed in subsection 4.2.3.

4.2.2　The Empirical Bayes Estimators of θ_n+1

To obtain the empirical Bayes estimators of , we need to estimate the hyperparameters from the supplementary information . There are two common methods to estimate the hyperparameters, that is, the moment method and the MLE method.

The estimators of the hyperparameters of the model (4.1) by the moment method and , and their consistencies are summarized in the following theorem whose proof can be found in appendix A.10.

Theorem 4.2. The estimators of the hyperparameters of the model (4.1) by the moment method are

(4.12)

(4.13)

where , , is the sample th moment of . Moreover, the moment estimators are consistent estimators of the hyperparameters.

The estimators of the hyperparameters of the model (4.1) by the MLE method and , and their consistencies are summarized in the following theorem whose proof can be found in appendix A.11.

Theorem 4.3. The estimators of the hyperparameters of the model (4.1) by the MLE method and are the solutions to the following equations:

(4.14)

(4.15)

The analytical calculations of the MLEs of and by solving the equations (4.14) and (4.15) are impossible, and thus we have to resort to numerical solutions. We can exploit Newton’s method to solve the equations (4.14) and (4.15) and to obtain the MLEs of and . Note that the MLEs of and are very sensitive to the initial estimators, and the moment estimators are usually proved to be good initial estimators.

Theorem 4.4. The empirical Bayes estimator of the parameter of the model (4.1) under Stein’s loss function by the moment method is given by (4.6) with the hyperparameters estimated by in Theorem 4.2. Alternatively, the empirical Bayes estimator of the parameter of the model (4.1) under Stein’s loss function by the MLE method is given by (4.6) with the hyperparameters estimated by numerically determined in Theorem 4.3.

4.2.3　Theoretical Comparisons of the Bayes Estimators and the PESLs of Three Methods

Note that the subscripts 0, 1, and 2 below are for the oracle method, the moment method, and the MLE method, respectively. The PESL functions of the three methods are respectively given by ()

where

and are unknown hyperparameters, and given in Theorem 4.2 are the moment estimators of the hyperparameters, and and numerically determined in Theorem 4.3 are the MLEs of the hyperparameters.

The Bayes estimators of under Stein’s loss function are given by

The Bayes estimators of under the squared error loss function are given by

for .

The PESLs evaluated at the Bayes estimators are given by

The PESLs evaluated at the Bayes estimators are given by

4.3　Simulations

The simulated data are generated according to the model (4.1) with the hyperparameters specified by and . The reason why we choose these values is that and are required in the moment estimations of the hyperparameters. Other numerical values of the hyperparameters can also be specified.

4.3.1　Two Inequalities of the Bayes Estimators and the PESLs

First, we fix , , and . Figure 4.1 shows the histogram of and the density estimation curve of . It is that we find to minimize the PESL. Numerical results show that

and

which exemplify the theoretical studies of (4.8) and (4.11).

FIG. 4.1 — Exp-IG: The histogram of and the density estimation curve of .