The Empirical Bayes Estimators of the Variance Parameter of the Normal Distribution with a Normal-Inverse-Gamma Prior under Stein's Loss Function

1. Introduction

The empirical Bayes method depends on a conjugate prior modeling, where the hyperparameters are estimated from the observations and the “estimated prior" is then used as a regular prior in the subsequent inference. See [1,2,3] and the references therein. The empirical Bayes method is introduced in [4,5,6]. From Bayesian perspective, it means that the sampling distribution is known, but the prior distribution is not. The marginal distribution is then used to restore the prior distribution from the observations. More literature on empirical Bayes method can be found, for instance, in [7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24].

The motivations of this paper are summarized as follows. Example 1.5.1 (p.20) of [25], Part I (pp.69-70) of [26], and [19] have considered the following hierarchical normal and normal-inverse-gamma model:

$$\left\{\begin{array}{c}{X}_i|\left(\mu ,\theta \right)\stackrel{\mathrm{iid}}{\sim }N\left(\mu ,\theta \right),i=1,2,\dots ,n,\hfill \\ \mu |\theta \sim N({\mu }_0,\theta /{\kappa }_0),\theta \sim \widetilde{IG}\left({\nu }_0/2,{\nu }_0{\sigma }_0^2/2\right),\hfill \end{array}\right.$$

(1)

where $-\infty <{\mu }_0<\infty$, ${\kappa }_0>0$, ${\nu }_0>0$, and ${\sigma }_0^2>0$ are known hyperparameters, and iid means independent and identically distributed. The $N-\widetilde{IG}({\mu }_0,{\kappa }_0,{\nu }_0,{\sigma }_0^2)$ distribution is a joint conjugate prior for $\left(\mu ,\theta \right)$ of the normal distribution $N\left(\mu ,\theta \right)$, so that the posterior distribution of $\left(\mu ,\theta \right)$ is an $N-\widetilde{IG}({\mu }_n,{\kappa }_n,{\nu }_n,{\sigma }_n^2)$ distribution with updated hyperparameters. However, in reality, the hyperparameters are unknown. [19] have estimated the hyperparameters of the model (1) by the moment method and the Maximum Likelihood Estimation (MLE) method. Moreover, they obtained the Bayes estimators of the mean and variance parameters of the model (1) under the squared error loss function. Finally, they obtained the empirical Bayes estimators of the mean and variance parameters of the model (1) under the squared error loss function by the moment method and the MLE method. However, in their empirical Bayes estimators, the sample ${\mathit{x}}^{\left(n\right)}=\left(x_1,x_2,\dots ,x_n\right)$ have been used twice. First, the sample ${\mathit{x}}^{\left(n\right)}$ are utilized to estimate the hyperparameters ${\mu }_0$, ${\kappa }_0$, ${\nu }_0$, and ${\sigma }_0^2$. Second, the sample ${\mathit{x}}^{\left(n\right)}$ are used to obtain the Bayes estimators. To avoid using the sample ${\mathit{x}}^{\left(n\right)}$ twice, and to be compatible with the usual empirical Bayes analysis, we will use the following hierarchical normal and normal-inverse-gamma model in this paper:

$$\left\{\begin{array}{c}{X}_i|\left({\mu }_i,{\theta }_i\right)\stackrel{\mathrm{independent}}{\sim }N\left({\mu }_i,{\theta }_i\right),i=1,2,\dots ,n,n+1,\hfill \\ {\mu }_i|{\theta }_i\stackrel{\mathrm{independent}}{\sim }N({\mu }_0,{\theta }_i/{\kappa }_0),{\theta }_i\stackrel{\mathrm{iid}}{\sim }\widetilde{IG}\left({\nu }_0/2,{\nu }_0{\sigma }_0^2/2\right),\hfill \end{array}\right.$$

(2)

where $-\infty <{\mu }_0<\infty$, ${\kappa }_0>0$, ${\nu }_0>0$, and ${\sigma }_0^2>0$ are hyperparameters to be determined, $-\infty <{\mu }_{n+1}<\infty$ and ${\theta }_{n+1}>0$ are the unknown parameters of interest, $N\left({\mu }_i,{\theta }_i\right)$ is a normal distribution with an unknown mean ${\mu }_i$ and an unknown variance ${\theta }_i$, the conditional conjugate prior distribution of ${\mu }_i$ given ${\theta }_i$ is $N({\mu }_0,{\theta }_i/{\kappa }_0)$ which is a normal distribution with mean ${\mu }_0$ and an unknown variance ${\theta }_i/{\kappa }_0$, the marginal conjugate prior distribution of ${\theta }_i$ is $\widetilde{IG}\left({\nu }_0/2,{\nu }_0{\sigma }_0^2/2\right)$ which is an inverse gamma distribution with shape parameter ${\nu }_0/2$ and scale parameter ${\nu }_0{\sigma }_0^2/2$. Note that the joint conjugate prior $\pi (\mu ,\theta )\sim N-\widetilde{IG}({\mu }_0,{\kappa }_0,{\nu }_0,{\sigma }_0^2)$ is a normal-inverse-gamma distribution. As described in [7], the statistician observes data $\left(x_1,x_2,\dots ,x_{n+1}\right)={\mathit{x}}^{\left(n+1\right)}$ and wishes to make an inference about ${\mu }_{n+1}$ and ${\theta }_{n+1}$. Therefore, $x_{n+1}$ provides direct information about the parameters ${\mu }_{n+1}$ and ${\theta }_{n+1}$, while supplementary information $\left(x_1,x_2,\dots ,x_n\right)={\mathit{x}}^{\left(n\right)}$ is also available. The connection between the prime data $x_{n+1}$ and the supplementary information ${\mathit{x}}^{\left(n\right)}$ is provided by the common distributions $N\left({\mu }_i,{\theta }_i\right)$, $N({\mu }_0,{\theta }_i/{\kappa }_0)$, and $\widetilde{IG}\left({\nu }_0/2,{\nu }_0{\sigma }_0^2/2\right)$. Moreover, since the variance parameter of the normal distribution is a positive restricted parameter, the squared error loss function is not appropriate. In contrast, we will choose Stein’s loss function because it penalizes gross overestimation and gross underestimation equally, that is, an action a will incur an infinite loss when it tends to 0 or ∞. Note that the squared error loss function does not have this property. For more literature on Stein’s loss function, we refer readers to [20,27,28,29,30].

The key contributions of this work are outlined below. For the hierarchical normal and normal-inverse-gamma (N-NIG) framework (2), we begin by deriving the posterior distributions and marginal density, as stated in Theorem 1. Subsequently, we derive Bayes estimators for the variance parameter of the normal distribution under two loss criteria: Stein’s loss and the squared error loss. We further compute the corresponding Posterior Expected Stein’s Losses (PESLs) associated with these two estimators. It is shown that the two Bayes estimators and their respective PESLs satisfy specific inequalities, namely (15) and (16). Additionally, we derive moment-based estimators for the hyperparameters of the N-NIG model (2) and establish their consistency in Theorem 2. Similarly, Maximum Likelihood Estimators (MLEs) for the same hyperparameters are derived, and their consistency is proven in Theorem 3. Building on these results, we construct empirical Bayes estimators for the variance parameter under Stein’s loss using both the method of moments and MLE, as presented in Theorem 4. The theoretical findings are supported by numerical experiments and an application to real-world data.

Comparing models (1) and (2) carefully, we find that the samples ${\mathit{x}}_1^{\left(n\right)}=\left(x_{11},x_{12},\dots ,x_{1n}\right)$ and ${\mathit{x}}_2^{\left(n\right)}=\left(x_{21},x_{22},\dots ,x_{2n}\right)$ generated from the two models are iid from different distributions. On the one hand, the sample ${\mathit{x}}_1^{\left(n\right)}$ generated from (1) are iid from $N\left(\mu ,\theta \right)$. Although the marginal densities of $x_{11},x_{12},\dots ,x_{1n}$ are $t_{{\nu }_0}\left({\mu }_0,{\sigma }_0^2\left({\kappa }_0+1\right)/{\kappa }_0\right)$, and thus $x_{11},x_{12},\dots ,x_{1n}$ can be thought to be from the t distribution. However, $x_{11},x_{12},\dots ,x_{1n}$ are not iid from the t distribution. In other words, $x_{11},x_{12},\dots ,x_{1n}$ are dependent from the t distribution. On the other hand, the sample ${\mathit{x}}_2^{\left(n\right)}$ generated from (2) are iid from the t distribution. That is, $x_{21},x_{22},\dots ,x_{2n}$ are independent and identically distributed from the t distribution. The sample ${\mathit{x}}_1^{\left(n\right)}$ can be used to estimate the parameters $\mu$ and $\theta$ from $N\left(\mu ,\theta \right)$, while the sample ${\mathit{x}}_2^{\left(n\right)}$ can be used to estimate the hyperparameters ${\mu }_0$, ${\nu }_0$, and $u_0$ from $t_{{\nu }_0}\left({\mu }_0,u_0\right)$, where $u_0={\sigma }_0^2\left({\kappa }_0+1\right)/{\kappa }_0$.

The structure of this paper is outlined as follows. Section 2 begins with the derivation of posterior distributions and the marginal likelihood for the hierarchical normal and normal-inverse-gamma model. Subsequently, we derive the corresponding Bayes estimators and PESLs, which are shown to fulfill two key inequalities: (15) and (16). Additionally, Theorem 4 presents the empirical Bayes estimators for the model’s variance component under Stein’s loss function, obtained via both the method of moments and MLE. In Section 3, we conduct numerical experiments from five different perspectives to evaluate performance. Section 4 illustrates the methodology using real-world data, specifically poverty rate data reflecting the proportion of individuals living below the poverty line. Finally, Section 5 offers concluding remarks and a discussion of the findings.

2. Theoretical Results

Suppose that $X\sim \tilde{G}\left(\alpha ,\beta \right)$ and $Y=1/X\sim \widetilde{IG}\left(\alpha ,\beta \right)$. The probability density functions (pdfs) of X and Y are respectively given by

$$\begin{array}{cc}\hfill f_X\left(x|\alpha ,\beta \right)& ={\displaystyle \frac{{\beta }^{\alpha }}{\Gamma \left(\alpha \right)}}{x}^{\alpha -1}exp\left(-x\beta \right),x>0,\alpha ,\beta >0,\hfill \\ \hfill f_Y\left(y|\alpha ,\beta \right)& ={\displaystyle \frac{{\beta }^{\alpha }}{\Gamma \left(\alpha \right)}}{\left({\displaystyle \frac{1}{y}}\right)}^{\alpha +1}exp\left(-{\displaystyle \frac{\beta }{y}}\right),y>0,\alpha ,\beta >0.\hfill \end{array}$$

It is easy to calculate $\mathrm{E}X=\alpha /\beta$ and $\mathrm{E}Y=\beta /\left(\alpha -1\right)$, for $\alpha >1$ and $\beta >0$.

2.1. The Bayes Estimators and the PESLs

For the hierarchical normal and normal-inverse-gamma model (2), we have the following theorem which calculates the posterior densities $\pi \left({\mu }_{n+1},{\theta }_{n+1}|x_{n+1}\right)$, $\pi \left({\theta }_{n+1}|x_{n+1}\right)$, $\pi \left({\mu }_{n+1}|x_{n+1}\right)$, $\pi \left({\mu }_{n+1}|{\theta }_{n+1},x_{n+1}\right)$, and the marginal density $m\left(x_{n+1}|{\mu }_0,{\kappa }_0,{\nu }_0,{\sigma }_0^2\right)$. The proof of the theorem can be found in the Appendix A.

Theorem 1. 

For the hierarchical normal and normal-inverse-gamma model (2), the joint posterior density of $\left({\mu }_{n+1},{\theta }_{n+1}\right)$ is

$$\pi \left({\mu }_{n+1},{\theta }_{n+1}|x_{n+1}\right)\sim N-\widetilde{IG}({\mu }_1,{\kappa }_1,{\nu }_1,{\sigma }_1^2),$$

the marginal posterior density of ${\theta }_{n+1}$ is

$$\pi \left({\theta }_{n+1}|x_{n+1}\right)\sim \widetilde{IG}\left({\alpha }^{*},{\beta }^{*}\right),$$

the marginal posterior density of ${\mu }_{n+1}$ is

$$\pi \left({\mu }_{n+1}|x_{n+1}\right)\sim t_{{\nu }_1}\left({\mu }_1,{\displaystyle \frac{{\sigma }_1^2}{{\kappa }_1}}\right),$$

the conditional posterior density of ${\mu }_{n+1}$ is

$$\pi \left({\mu }_{n+1}|{\theta }_{n+1},x_{n+1}\right)\sim N({\mu }_1,{\theta }_{n+1}/{\kappa }_1),$$

where

$$\left\{\begin{array}{c}{\mu }_1=\left(x_{n+1}+{\kappa }_0{\mu }_0\right)/\left({\kappa }_0+1\right),\hfill \\ {\kappa }_1={\kappa }_0+1,\hfill \\ {\nu }_1={\nu }_0+1,\hfill \\ {\nu }_1{\sigma }_1^2={\nu }_0{\sigma }_0^2+{\displaystyle \frac{{\kappa }_0}{{\kappa }_0+1}}{\left(x_{n+1}-{\mu }_0\right)}^2,\hfill \end{array}\right.$$

(3)

$${\alpha }^{*}={\displaystyle \frac{{\nu }_1}{2}}={\displaystyle \frac{{\nu }_0+1}{2}},$$

(4)

and

$${\beta }^{*}={\displaystyle \frac{{\nu }_1{\sigma }_1^2}{2}}={\displaystyle \frac{1}{2}}\left[{\nu }_0{\sigma }_0^2+{\displaystyle \frac{{\kappa }_0}{{\kappa }_0+1}}{\left(x_{n+1}-{\mu }_0\right)}^2\right].$$

(5)

Moreover, the marginal density of $x_{n+1}$ is given by

$$m\left(x_{n+1}|{\mu }_0,{\kappa }_0,{\nu }_0,{\sigma }_0^2\right)\sim t_{{\nu }_0}\left({\mu }_0,{\displaystyle \frac{{\sigma }_0^2\left({\kappa }_0+1\right)}{{\kappa }_0}}\right),$$

with probability density function (pdf) given by

$$\begin{array}{cc}\hfill & m\left(x_{n+1}|{\mu }_0,{\kappa }_0,{\nu }_0,{\sigma }_0^2\right)\hfill \\ \hfill & ={\displaystyle \frac{\Gamma \left(\frac{{\nu }_0+1}{2}\right)}{\Gamma \left(\frac{{\nu }_0}{2}\right)\sqrt{\pi {\nu }_0{\sigma }_0^2\left({\kappa }_0+1\right)/{\kappa }_0}}}{\left[1+{\displaystyle \frac{1}{{\nu }_0}}{\displaystyle \frac{{\left(x_{n+1}-{\mu }_0\right)}^2}{{\sigma }_0^2\left({\kappa }_0+1\right)/{\kappa }_0}}\right]}^{-{\displaystyle \frac{{\nu }_0+1}{2}}},\hfill \end{array}$$

(6)

for $-\infty <x_{n+1}<\infty$, $-\infty <{\mu }_0<\infty$, ${\kappa }_0>0$, ${\nu }_0>0$, and ${\sigma }_0^2>0$.

In the following, we will calculate the Bayes estimator of ${\theta }_{n+1}$ under Stein’s loss function ${\delta }_s^{\pi ,{\theta }_{n+1}}\left(x_{n+1}\right)$, the Bayes estimator of ${\theta }_{n+1}$ under the usual squared error loss function ${\delta }_2^{\pi ,{\theta }_{n+1}}\left(x_{n+1}\right)$, and the PESLs at ${\delta }_s^{\pi ,{\theta }_{n+1}}\left(x_{n+1}\right)$ and ${\delta }_2^{\pi ,{\theta }_{n+1}}\left(x_{n+1}\right)$ ($PES{L}_s^{\pi ,{\theta }_{n+1}}\left(x_{n+1}\right)$ and $PES{L}_2^{\pi ,{\theta }_{n+1}}\left(x_{n+1}\right)$).

Similar to [28], the two Bayes estimators and the two PESLs are respectively given by

$${\delta }_s^{\pi ,{\theta }_{n+1}}\left(x_{n+1}\right)={\displaystyle \frac{1}{\mathrm{E}\left(\frac{1}{{\theta }_{n+1}}|x_{n+1}\right)}},$$

(7)

$${\delta }_2^{\pi ,{\theta }_{n+1}}\left(x_{n+1}\right)=\mathrm{E}\left({\theta }_{n+1}|x_{n+1}\right),$$

(8)

$$PES{L}_s^{\pi ,{\theta }_{n+1}}\left(x_{n+1}\right)=log\mathrm{E}\left({\displaystyle \frac{1}{{\theta }_{n+1}}}|x_{n+1}\right)+\mathrm{E}\left(log{\theta }_{n+1}|x_{n+1}\right),$$

(9)

$$\begin{array}{cc}\hfill PES{L}_2^{\pi ,{\theta }_{n+1}}\left(x_{n+1}\right)& =\mathrm{E}\left({\theta }_{n+1}|x_{n+1}\right)\mathrm{E}\left({\displaystyle \frac{1}{{\theta }_{n+1}}}|x_{n+1}\right)-1\hfill \\ \hfill & -log\mathrm{E}\left({\theta }_{n+1}|x_{n+1}\right)+\mathrm{E}\left(log{\theta }_{n+1}|x_{n+1}\right).\hfill \end{array}$$

(10)

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

$$E_1=\mathrm{E}\left({\theta }_{n+1}|x_{n+1}\right),E_2=\mathrm{E}\left({\displaystyle \frac{1}{{\theta }_{n+1}}}|x_{n+1}\right),\mathrm{and}{E}_3=\mathrm{E}\left(log{\theta }_{n+1}|x_{n+1}\right).$$

Similar to [28], it can be shown that

$${\delta }_s^{\pi ,{\theta }_{n+1}}\left(x_{n+1}\right)\le {\delta }_2^{\pi ,{\theta }_{n+1}}\left(x_{n+1}\right)$$

(11)

by exploiting Jensen’s inequality. Moreover,

$$PES{L}_s^{\pi ,{\theta }_{n+1}}\left(x_{n+1}\right)\le PES{L}_2^{\pi ,{\theta }_{n+1}}\left(x_{n+1}\right),$$

(12)

which is a direct consequence of the general methodology for finding a Bayes estimator. According to construction, ${\delta }_s^{\pi ,{\theta }_{n+1}}\left(x_{n+1}\right)$ minimizes the Posterior Expected Stein’s Loss (PESL).

Now let us analytically calculate the Bayes estimators ${\delta }_s^{\pi ,{\theta }_{n+1}}\left(x_{n+1}\right)$ and ${\delta }_2^{\pi ,{\theta }_{n+1}}\left(x_{n+1}\right)$, and the PESLs $PES{L}_s^{\pi ,{\theta }_{n+1}}\left(x_{n+1}\right)$ and $PES{L}_2^{\pi ,{\theta }_{n+1}}\left(x_{n+1}\right)$ under the hierarchical normal and normal-inverse-gamma model (2) from (7), (8), (9), and (10). The three expectations are calculated as

$$\begin{array}{cc}\hfill \mathrm{E}\left({\theta }_{n+1}|x_{n+1}\right)& ={\displaystyle \frac{{\beta }^{*}}{{\alpha }^{*}-1}},\hfill \\ \hfill \mathrm{E}\left({\displaystyle \frac{1}{{\theta }_{n+1}}}|x_{n+1}\right)& ={\displaystyle \frac{{\alpha }^{*}}{{\beta }^{*}}},\hfill \\ \hfill \mathrm{E}\left(log{\theta }_{n+1}|x_{n+1}\right)& =log{\beta }^{*}-\psi \left({\alpha }^{*}\right),\hfill \end{array}$$

for ${\alpha }^{*}>1$, where ${\alpha }^{*}$ and ${\beta }^{*}$ are given by (4) and (5). From (7), the Bayes estimator of ${\theta }_{n+1}$ under Stein’s loss function is given by

$${\delta }_s^{\pi ,{\theta }_{n+1}}\left(x_{n+1}\right)={\displaystyle \frac{1}{\mathrm{E}\left(\frac{1}{{\theta }_{n+1}}|x_{n+1}\right)}}={\displaystyle \frac{{\beta }^{*}}{{\alpha }^{*}}}.$$

(13)

From (8), the Bayes estimator of ${\theta }_{n+1}$ under the usual squared error loss function is given by

$${\delta }_2^{\pi ,{\theta }_{n+1}}\left(x_{n+1}\right)=\mathrm{E}\left({\theta }_{n+1}|x_{n+1}\right)={\displaystyle \frac{{\beta }^{*}}{{\alpha }^{*}-1}},$$

(14)

for ${\alpha }^{*}>1$. It is easy to show that

$${\delta }_s^{\pi ,{\theta }_{n+1}}\left(x_{n+1}\right)={\displaystyle \frac{{\beta }^{*}}{{\alpha }^{*}}}<{\displaystyle \frac{{\beta }^{*}}{{\alpha }^{*}-1}}={\delta }_2^{\pi ,{\theta }_{n+1}}\left(x_{n+1}\right),$$

(15)

which exemplifies the theoretical study of (11). Furthermore, from (9) and (10), the PESL at ${\delta }_s^{\pi ,{\theta }_{n+1}}\left(x_{n+1}\right)$ and ${\delta }_2^{\pi ,{\theta }_{n+1}}\left(x_{n+1}\right)$ are respectively given by

$$PES{L}_s^{\pi ,{\theta }_{n+1}}\left(x_{n+1}\right)=log{\alpha }^{*}-\psi \left({\alpha }^{*}\right)$$

and

$$PES{L}_2^{\pi ,{\theta }_{n+1}}\left(x_{n+1}\right)={\displaystyle \frac{1}{{\alpha }^{*}-1}}+log\left({\alpha }^{*}-1\right)-\psi \left({\alpha }^{*}\right),$$

where

$$\psi \left(x\right)={\displaystyle \frac{{\Gamma }^{\prime }\left(x\right)}{\Gamma \left(x\right)}}={\displaystyle \frac{d}{dx}}log\Gamma \left(x\right)=\mathrm{digamma}\left(x\right)$$

is the digamma function. It can be shown that

$$\begin{array}{cc}\hfill & PES{L}_s^{\pi ,{\theta }_{n+1}}\left(x_{n+1}\right)\hfill \\ \hfill & =log{\alpha }^{*}-\psi \left({\alpha }^{*}\right)\hfill \\ \hfill & \le {\displaystyle \frac{1}{{\alpha }^{*}-1}}+log\left({\alpha }^{*}-1\right)-\psi \left({\alpha }^{*}\right)\hfill \\ \hfill & =PES{L}_2^{\pi ,{\theta }_{n+1}}\left(x_{n+1}\right),\hfill \end{array}$$

(16)

which exemplifies the theoretical study of (12). It is worth noting that the PESLs $PES{L}_s^{\pi ,{\theta }_{n+1}}\left(x_{n+1}\right)$ and $PES{L}_2^{\pi ,{\theta }_{n+1}}\left(x_{n+1}\right)$ depend only on ${\alpha }^{*}$, but not on ${\beta }^{*}$. Therefore, the PESLs depend only on ${\nu }_0$, but not on ${\mu }_0$, ${\kappa }_0$, ${\sigma }_0^2$, and $x_{n+1}$.

We will exemplify the two inequalities (15) and (16) in the simulations section and the real data section. Moreover, we will exemplify that the PESLs depend only on ${\nu }_0$, but not on ${\mu }_0$, ${\kappa }_0$, ${\sigma }_0^2$, and $x_{n+1}$.

2.2. The Empirical Bayes Estimators of ${\theta }_{n+1}$

The hyperparameters of model (2) are $-\infty <{\mu }_0<\infty$, ${\kappa }_0>0$, ${\nu }_0>0$, and ${\sigma }_0^2>0$. However, we can not directly obtain the estimators of the four hyperparameters of model (2) by the moment method. Let

$$u_0={\displaystyle \frac{{\sigma }_0^2\left({\kappa }_0+1\right)}{{\kappa }_0}}.$$

(17)

Since ${\kappa }_0$ and ${\sigma }_0$ appear together in $u_0$, we can not directly obtain the estimators of ${\kappa }_0$ and ${\sigma }_0$ by the moment method. In other words, ${\kappa }_0$ and ${\sigma }_0$ are unidentifiable. In the empirical Bayesian statistical literature, common approaches to addressing the issue of unidentifiability of hyperparameters include the following two. One is to estimate the prior distribution through non-parametric or semi-parametric methods, avoiding strong assumptions about the functional form of the prior distribution, thereby circumventing the problem of unidentifiability of hyperparameters ([31,32]). The other is to use auxiliary data, model structure constraints, or specific assumptions (such as sparsity, spatial correlation) to provide additional information, making the unidentifiable hyperparameters identifiable ([33,34,35]). We adopt the second approach to make hyperparameters (${\kappa }_0$ and ${\sigma }_0^2$) identifiable. More specifically, when ${\kappa }_0$ is fixed to be a known constant (our recommendation is ${\kappa }_0=1$, which will be made clear later in this subsection), then ${\kappa }_0$ and ${\sigma }_0^2$ are identifiable. Otherwise, ${\kappa }_0$ and ${\sigma }_0^2$ are unidentifiable.

However, we can obtain the estimator of $u_0$ by the moment method. In the following, we are interested in the hyperparameters $-\infty <{\mu }_0<\infty$, ${\nu }_0>0$, and $u_0>0$. Using hyperparameters ${\mu }_0$, ${\nu }_0$, and $u_0$, the marginal density (6) changes to

$$\begin{array}{cc}\hfill & m\left(x_{n+1}|{\mu }_0,{\nu }_0,u_0\right)\hfill \\ \hfill & ={\displaystyle \frac{\Gamma \left(\frac{{\nu }_0+1}{2}\right)}{\Gamma \left(\frac{{\nu }_0}{2}\right)\sqrt{\pi {\nu }_0{u}_0}}}{\left[1+{\displaystyle \frac{1}{{\nu }_0}}{\displaystyle \frac{{\left(x_{n+1}-{\mu }_0\right)}^2}{u_0}}\right]}^{-{\displaystyle \frac{{\nu }_0+1}{2}}}\hfill \\ \hfill & \sim t_{{\nu }_0}\left({\mu }_0,u_0\right).\hfill \end{array}$$

(18)

The moment-based estimators for the hyperparameters ${\mu }_0$, ${\nu }_0$, and $u_0$ in model (2), denoted as ${\mu }_{01}\left(n\right)$, ${\nu }_{01}\left(n\right)$, and $u_{01}\left(n\right)$, along with their consistency properties, are presented in the subsequent theorem. The proof of this theorem is provided in the Appendix A.

Theorem 2. 

The estimators of the hyperparameters ${\mu }_0$, ${\nu }_0$, and $u_0$ of model (2) by the moment method are

$$\begin{array}{cc}\hfill {\mu }_{01}\left(n\right)& =A_1,\hfill \\ \hfill {\nu }_{01}\left(n\right)& ={\displaystyle \frac{-\frac{14}{3}{A}_1^4+4A_1^2{A}_2+2A_2^2-\frac{4}{3}{A}_4}{-\frac{2}{3}{A}_1^4+A_2^2-\frac{1}{3}{A}_4}},\hfill \\ \hfill u_{01}\left(n\right)& ={\displaystyle \frac{-\frac{1}{3}\left(A_2-A_1^2\right)\left(5A_1^4-6A_1^2{A}_2+A_4\right)}{-\frac{7}{3}{A}_1^4+2A_1^2{A}_2+A_2^2-\frac{2}{3}{A}_4}},\hfill \end{array}$$

where

$$A_k={\displaystyle \frac{1}{n}}\sum _{i=1}^n{X}_i^k,k=1,2,4,$$

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

We remark that the moment estimators ${\mu }_{01}\left(n\right)$, ${\nu }_{01}\left(n\right)$, and $u_{01}\left(n\right)$ in Theorem 2 are the same as those in Theorem 2 in [19]. The reason for the same moment estimators is that for the two hierarchical normal and normal-inverse-gamma models (1) and (2), the marginal distributions are the same, and the population moments of X are the same. Moreover, in Theorem 2 of this paper, we have shown that the moment estimators are consistent estimators of the hyperparameters, and this result has not been derived in [19].

The MLE method is used to derive estimators for the hyperparameters ${\mu }_0$, ${\nu }_0$, and $u_0$ of the model (2), denoted as ${\mu }_{02}\left(n\right)$, ${\nu }_{02}\left(n\right)$, and $u_{02}\left(n\right)$, respectively. The consistency properties of these estimators are presented in the subsequent theorem, with the detailed proof provided in the Appendix A.

Theorem 3. 

The estimators of the hyperparameters ${\mu }_0$, ${\nu }_0$, and $u_0$ of the model (2) by the MLE method ${\mu }_{02}\left(n\right)$, ${\nu }_{02}\left(n\right)$, and $u_{02}\left(n\right)$ are the solutions to the following equations:

$$f_1\left({\mu }_0,{\nu }_0,u_0\right)={\displaystyle \frac{1}{n}}\sum _{i=1}^n{\displaystyle \frac{\left({\nu }_0+1\right)\left(x_i-{\mu }_0\right)}{{\nu }_0{u}_0+{\left(x_i-{\mu }_0\right)}^2}}=0,$$

(19)

$$\begin{array}{cc}\hfill f_2\left({\mu }_0,{\nu }_0,u_0\right)& ={\displaystyle \frac{1}{2}}\psi \left({\displaystyle \frac{{\nu }_0+1}{2}}\right)-{\displaystyle \frac{1}{2}}\psi \left({\displaystyle \frac{{\nu }_0}{2}}\right)-{\displaystyle \frac{1}{2{\nu }_0}}\hfill \\ \hfill & -{\displaystyle \frac{1}{2}}{\displaystyle \frac{1}{n}}\sum _{i=1}^n\left\{log\left[1+{\displaystyle \frac{{\left(x_i-{\mu }_0\right)}^2}{{\nu }_0{u}_0}}\right]-{\displaystyle \frac{\left({\nu }_0+1\right){\left(x_i-{\mu }_0\right)}^2}{\left[{\nu }_0{u}_0+{\left(x_i-{\mu }_0\right)}^2\right]{\nu }_0}}\right\}=0,\hfill \end{array}$$

(20)

$$f_3\left({\mu }_0,{\nu }_0,u_0\right)=-{\displaystyle \frac{1}{2u_0}}+{\displaystyle \frac{1}{2}}{\displaystyle \frac{1}{n}}\sum _{i=1}^n{\displaystyle \frac{\left({\nu }_0+1\right){\left(x_i-{\mu }_0\right)}^2}{\left[{\nu }_0{u}_0+{\left(x_i-{\mu }_0\right)}^2\right]u_0}}=0.$$

(21)

Moreover, the MLEs are consistent estimators of the hyperparameters.

The MLEs of the hyperparameters ${\mu }_0$, ${\nu }_0$, and $u_0$ cannot be derived analytically by solving equations (19), (20), and (21). As a result, numerical methods must be employed. Newton’s method can be applied to solve these equations and thereby obtain the MLEs of the hyperparameters. It should be emphasized that the accuracy of the resulting MLEs is highly sensitive to the choice of initial values, and moment-based estimators have generally been found to serve as effective starting points.

Finally, the empirical Bayes estimates of the model (2)’s variance parameter under Stein’s loss function, derived using the method of moments and the MLE, are presented in the subsequent theorem.

Theorem 4. 

The empirical Bayes estimator of the variance parameter of the model (2) under Stein’s loss function by the moment method is given by (13) with the hyperparameters estimated by $\left({\mu }_{01}\left(n\right),{\nu }_{01}\left(n\right),u_{01}\left(n\right)\right)$ in Theorem 2. Alternatively, the empirical Bayes estimator of the variance parameter of the model (2) under Stein’s loss function by the MLE method is given by (13) with the hyperparameters estimated by $\left({\mu }_{02}\left(n\right),{\nu }_{02}\left(n\right),u_{02}\left(n\right)\right)$ numerically determined in Theorem 3.

Now let us discuss the selection of ${\kappa }_0$. We recommend choosing ${\kappa }_0=1$, and the reason is given as follows. From (13), (4), and (5), we have

$$\begin{array}{cc}\hfill {\delta }_s^{\pi ,{\theta }_{n+1}}\left(x_{n+1}\right)& ={\displaystyle \frac{{\beta }^{*}}{{\alpha }^{*}}}\hfill \\ \hfill & ={\displaystyle \frac{\frac{1}{2}\left[{\nu }_0{\sigma }_0^2+\frac{{\kappa }_0}{{\kappa }_0+1}{\left(x_{n+1}-{\mu }_0\right)}^2\right]}{\frac{{\nu }_0+1}{2}}}\hfill \\ \hfill & ={\displaystyle \frac{1}{{\nu }_0+1}}\left[{\nu }_0{\sigma }_0^2+{\displaystyle \frac{{\kappa }_0}{{\kappa }_0+1}}{\left(x_{n+1}-{\mu }_0\right)}^2\right]\hfill \\ \hfill & ={\displaystyle \frac{{\kappa }_0}{{\kappa }_0+1}}{\displaystyle \frac{1}{{\nu }_0+1}}\left[{\nu }_0{\displaystyle \frac{{\sigma }_0^2\left({\kappa }_0+1\right)}{{\kappa }_0}}+{\left(x_{n+1}-{\mu }_0\right)}^2\right]\hfill \\ \hfill & ={\displaystyle \frac{{\kappa }_0}{{\kappa }_0+1}}{\displaystyle \frac{1}{{\nu }_0+1}}\left[{\nu }_0{u}_0+{\left(x_{n+1}-{\mu }_0\right)}^2\right].\hfill \end{array}$$

(22)

It is easy to show that the factor

$${\displaystyle \frac{{\kappa }_0}{{\kappa }_0+1}}\in \left(0,1\right),$$

for ${\kappa }_0>0$. Because we have little information about ${\kappa }_0$, we will choose

$${\displaystyle \frac{{\kappa }_0}{{\kappa }_0+1}}={\displaystyle \frac{1}{2}},$$

which is in the middle of the above range. Hence,

$${\kappa }_0=1.$$

(23)

Therefore, (22) reduces to

$${\delta }_s^{\pi ,{\theta }_{n+1}}\left(x_{n+1}\right)={\displaystyle \frac{1}{2}}{\displaystyle \frac{1}{{\nu }_0+1}}\left[{\nu }_0{u}_0+{\left(x_{n+1}-{\mu }_0\right)}^2\right],$$

which can then be estimated once the hyperparameters $\left({\mu }_0,{\nu }_0,u_0\right)$ are estimated. From (17) and (23), we have

$${\sigma }_0^2={\displaystyle \frac{{\kappa }_0}{{\kappa }_0+1}}{u}_0={\displaystyle \frac{1}{2}}{u}_0,$$

which can then be estimated once the hyperparameter $u_0$ is estimated.

Another reason to choose ${\kappa }_0=1$ is given below. Since $-\infty <{\mu }_{n+1}<\infty$, the squared error loss function is appropriate. The Bayes estimator of ${\mu }_{n+1}$ under the squared error loss function is given by

$${\delta }_2^{\pi ,{\mu }_{n+1}}\left(x_{n+1}\right)=\mathrm{E}\left({\mu }_{n+1}|x_{n+1}\right)={\mu }_1={\displaystyle \frac{x_{n+1}+{\kappa }_0{\mu }_0}{{\kappa }_0+1}}.$$

(24)

In (24), ${\kappa }_0$ represents a strength of belief on ${\mu }_0$. If one harbors no belief on ${\mu }_0$, then ${\kappa }_0=0$, and thus ${\delta }_2^{\pi ,{\mu }_{n+1}}\left(x_{n+1}\right)=x_{n+1}$, which depends only on the datum $x_{n+1}$. In contrast, if one harbors complete belief on ${\mu }_0$, then ${\kappa }_0=\infty$, and thus ${\delta }_2^{\pi ,{\mu }_{n+1}}\left(x_{n+1}\right)={\mu }_0$, which depends only on ${\mu }_0$. However, if one believes that $x_{n+1}$ and ${\mu }_0$ are equally important, then ${\kappa }_0=1$ is a reasonable choice, and thus ${\delta }_2^{\pi ,{\mu }_{n+1}}\left(x_{n+1}\right)=\left(x_{n+1}+{\mu }_0\right)/2$, which is a balanced combination of $x_{n+1}$ and ${\mu }_0$.

We remark that ${\kappa }_0$ and ${\sigma }_0^2$ affect ${\beta }^{*}$, ${\delta }_s^{\pi ,{\theta }_{n+1}}\left(x_{n+1}\right)$, ${\delta }_2^{\pi ,{\theta }_{n+1}}\left(x_{n+1}\right)$, and ${\delta }_2^{\pi ,{\mu }_{n+1}}\left(x_{n+1}\right)$, and they do not affect ${\alpha }^{*}$, $PES{L}_s^{\pi ,{\theta }_{n+1}}\left(x_{n+1}\right)$, and $PES{L}_2^{\pi ,{\theta }_{n+1}}\left(x_{n+1}\right)$.

3. Simulations

This section presents numerical simulations for the hierarchical normal and normal-inverse-gamma model (2). We examine five key aspects: the consistency of moment estimators and MLEs, the influence of ${\kappa }_0$ on quantities of interest as functions of $x_{n+1}$, the validity of two inequalities concerning Bayes estimators and PESLs, the model’s goodness-of-fit to the generated data, and the marginal density distributions of the model parameters.

The simulated data are generated according to the hierarchical normal and normal-inverse-gamma model (2) with the hyperparameters specified by ${\mu }_0=3$, ${\kappa }_0=1$, ${\nu }_0=6$, and ${\sigma }_0^2=1$. The reason why we choose these values is that $-\infty <{\mu }_0<\infty$, ${\kappa }_0>0$, ${\nu }_0>0$, and ${\sigma }_0^2>0$. Moreover, ${\nu }_0>4$ is required in moment estimations of the hyperparameters. Additional numerical values for the hyperparameters can likewise be defined.

3.1. Consistencies of the Moment Estimators and the MLEs

Currently, both the moment estimators $\left({\mu }_{01}(n),{\nu }_{01}(n),u_{01}(n)\right)$ and the MLEs $\left({\mu }_{02}(n),{\nu }_{02}(n),u_{02}(n)\right)$ serve as consistent estimates for the hyperparameters $\left({\mu }_0,{\nu }_0,u_0\right)$ of the hierarchical normal and normal-inverse-gamma model (2). The purpose of this subsection is to emphasize that, as established in Theorems 2 and 3, we have rigorously proven the consistency of both the moment-based and MLE-based estimators with respect to the true hyperparameters. It should be noted that all analyses here are based solely on the observed data ${\mathit{x}}^{(n)}=(x_1,x_2,\dots ,x_n)$.

First, we will numerically exemplify that the sample ${\mathit{x}}_1^{\left(n\right)}=\left(x_{11},x_{12},\dots ,x_{1n}\right)$ generated from the model (1) can not be used to estimate the hyperparameters $\left({\mu }_0,{\nu }_0,u_0\right)$, while the sample ${\mathit{x}}_2^{\left(n\right)}=\left(x_{21},x_{22},\dots ,x_{2n}\right)$ generated from the model (2) can be used to estimate the hyperparameters $\left({\mu }_0,{\nu }_0,u_0\right)$, where $u_0={\sigma }_0^2\left({\kappa }_0+1\right)/{\kappa }_0$. Moreover, we will exemplify that the moment estimators and the MLEs of $u_0$ can correctly estimate the true hyperparameter $u_0$ regardless of the ${\kappa }_0$ and ${\sigma }_0^2$ values.

The histograms of the samples and their density curves are plotted in Figure 1. From the figure, we observe the following facts.

• Plot (a): the sample ${\mathit{x}}_1^{\left(n\right)}$ generated from (1) are iid from $N\left(\mu ,\theta \right)$ with $\mu =4.817$ and $\theta =1.866$.
• Plots (b) and (c): the sample ${\mathit{x}}_2^{\left(n\right)}$ are generated from (2) with $\left({\mu }_0,{\kappa }_0,{\nu }_0,{\sigma }_0^2\right)=\left(3,1,6,1\right)$ and $u_0=2$. The sample ${\mathit{x}}_3^{\left(n\right)}$ are generated from (2) with $\left({\mu }_0,{\kappa }_0,{\nu }_0,{\sigma }_0^2\right)=\left(3,2,6,4/3\right)$ and $u_0=2$. Although $\left({\kappa }_0,{\sigma }_0^2\right)=\left(1,1\right)$ and $\left(2,4/3\right)$ are different in the two plots, $u_0=2$ is the same. Therefore, the two samples ${\mathit{x}}_2^{\left(n\right)}$ and ${\mathit{x}}_3^{\left(n\right)}$ are from the same marginal distribution $t_{{\nu }_0}\left({\mu }_0,u_0\right)$ with $\left({\mu }_0,{\nu }_0,u_0\right)=\left(3,6,2\right)$.

The moment estimators and the MLEs of the hyperparameters $\left({\mu }_0,{\nu }_0,u_0\right)=\left(3,6,2\right)$ for the samples ${\mathit{x}}_1^{\left(n\right)}$, ${\mathit{x}}_2^{\left(n\right)}$, and ${\mathit{x}}_3^{\left(n\right)}$ are summarized in Table 1. From the table, we observe the following facts.

• The moment estimators of the hyperparameters $\left({\mu }_0,{\nu }_0,u_0\right)$ for sample ${\mathit{x}}_1^{\left(n\right)}$ are far away from the true hyperparameters $\left(3,6,2\right)$, and thus the samples generated from the model (1) can not be used to estimate the hyperparameters $\left({\mu }_0,{\nu }_0,u_0\right)$.
• For ${\mathit{x}}_1^{\left(n\right)}$, since the moment estimator of ${\nu }_0$ is $-159.725$ which is negative, the MLE method fails to iterate, and thus the MLEs of the hyperparameters $\left({\mu }_0,{\nu }_0,u_0\right)$ are equal to the moment estimators.
• For ${\mathit{x}}_2^{\left(n\right)}$ and ${\mathit{x}}_3^{\left(n\right)}$, the moment estimators and the MLEs of the hyperparameters $\left({\mu }_0,{\nu }_0,u_0\right)$ are close to the true hyperparameters $\left(3,6,2\right)$, and thus the samples generated from the model (2) can be used to estimate the hyperparameters $\left({\mu }_0,{\nu }_0,u_0\right)$.
• The sample ${\mathit{x}}_2^{\left(n\right)}$ is generated from the model (2) with $\left({\kappa }_0,{\sigma }_0^2\right)=\left(1,1\right)$, while the sample ${\mathit{x}}_3^{\left(n\right)}$ is generated from the model (2) with $\left({\kappa }_0,{\sigma }_0^2\right)=\left(2,4/3\right)$. Although $\left({\kappa }_0,{\sigma }_0^2\right)=\left(1,1\right)$ and $\left(2,4/3\right)$ are different for the two samples ${\mathit{x}}_2^{\left(n\right)}$ and ${\mathit{x}}_3^{\left(n\right)}$, $u_0=2$ is the same. We find that both the moment method and the MLE method correctly estimate the true hyperparameter $u_0=2$ for the two samples ${\mathit{x}}_2^{\left(n\right)}$ and ${\mathit{x}}_3^{\left(n\right)}$.
• For ${\mathit{x}}_2^{\left(n\right)}$ and ${\mathit{x}}_3^{\left(n\right)}$, the MLEs are closer to the true hyperparameters $\left({\mu }_0,{\nu }_0,u_0\right)=\left(3,6,2\right)$ than the moment estimators for this simulation.

Following the approach in [21], let $\xi$ represent one of the hyperparameters ${\mu }_0$, ${\nu }_0$, or $u_0$. Consistency is defined as the property that

$${\xi }_i(n)\stackrel{p}{⟶}\xi \mathrm{as}n\to \infty ,$$

for $i=1,2$, where $i=1$ corresponds to the moment estimator and $i=2$ to the MLE, with $\stackrel{p}{⟶}$ denoting convergence in probability and n being the sample size. Equivalently, consistency can be expressed as

$$\underset{n\to \infty }{lim}P\left(\left|{\xi }_i\left(n\right)-\xi \right|\ge \epsilon \right)=0,$$

for any $\epsilon >0$ and for all $\xi \in \Xi$. These probabilities are estimated using empirical frequencies from simulations:

$$P\left(\left|{\xi }_i\left(n\right)-\xi \right|\ge \epsilon \right)\approx {\displaystyle \frac{1}{M}}\sum _{m=1}^{M}I\left(\left|{\xi }_i^{\left(m\right)}\left(n\right)-\xi \right|\ge \epsilon \right)\triangleq F_{{\xi }_i}^{\epsilon ,M}\left(n\right),$$

where $I(A)$ is the indicator function taking value 1 if event A occurs and 0 otherwise, and M denotes the number of simulation replications. Hence, a decreasing trend in the simulated frequencies $F_{{\xi }_i}^{\epsilon ,M}(n)$ toward zero as n increases provides empirical evidence of estimator consistency, for both estimation methods ($i=1,2$) and across all hyperparameters $\xi \in \{{\mu }_0,{\nu }_0,u_0\}$.

We now demonstrate that the moment estimators $\left({\mu }_{01}(n),{\nu }_{01}(n),u_{01}(n)\right)$ and the MLEs $\left({\mu }_{02}(n),{\nu }_{02}(n),u_{02}(n)\right)$ are consistent estimators of the hyperparameters $\left({\mu }_0,{\nu }_0,u_0\right)$. Table 2 presents the frequency distributions of these estimators for both moment-based and MLE for varying sample sizes n, under settings where $\epsilon =1$, $0.5$, and $0.1$. The simulated data are generated from the hierarchical normal and normal-inverse-gamma model (2), with hyperparameter values set to $\left({\mu }_0,{\kappa }_0,{\nu }_0,{\sigma }_0^2\right)=\left(3,1,6,1\right)$ and $u_0=2$. While other configurations of these hyperparameters could also be considered, the current choice serves to illustrate the behavior effectively. Analysis of the table reveals several key observations.

• With $\epsilon =1$, $0.5$, or $0.1$, the frequencies of the estimators approach zero as n tends to infinity, indicating that both the moment estimators and the MLEs are consistent for estimating the hyperparameters. In the case of $\epsilon =0.1$, the frequencies of ${\nu }_{01}(n)$ and ${\nu }_{02}(n)$ remain relatively high (at least $0.2$ across all scenarios). Nevertheless, a decreasing trend toward zero is evident as n grows large, suggesting eventual convergence.
• By comparing the frequencies corresponding to $\epsilon =1$, $0.5$, and $0.1$, we find that as $\epsilon$ decreases, the frequencies generally increase. This occurs because the constraints

  $$\left|{\xi }_i^{\left(m\right)}\left(n\right)-\xi \right|\ge \epsilon ,\mathrm{for}\xi ={\mu }_0,{\nu }_0,\mathrm{and}{u}_0,i=1,2,$$

  become easier to satisfy when $\epsilon$ is smaller.
• When comparing the moment estimators with the MLEs of the hyperparameters ${\mu }_0$, ${\nu }_0$, and $u_0$, it is observed that for large sample sizes n, the MLEs exhibit smaller values than the moment estimators. This indicates that, in terms of consistency, the MLEs perform more reliably in estimating the hyperparameters.

3.2. The Effect of ${\kappa }_0$ for Quantities of Interest as Functions of $x_{n+1}$

In this subsection, we will investigate the effect of ${\kappa }_0$ for quantities of interest (${\beta }^{*}$, ${\delta }_s^{\pi ,{\theta }_{n+1}}\left(x_{n+1}\right)$, ${\delta }_2^{\pi ,{\theta }_{n+1}}\left(x_{n+1}\right)$, ${\delta }_2^{\pi ,{\mu }_{n+1}}\left(x_{n+1}\right)$, ${\alpha }^{*}$, $PES{L}_s^{\pi ,{\theta }_{n+1}}\left(x_{n+1}\right)$, and $PES{L}_2^{\pi ,{\theta }_{n+1}}\left(x_{n+1}\right)$) as functions of $x_{n+1}$. The motivation of this subsection is that ${\kappa }_0$ is unidentifiable, and we want to know whether using ${\kappa }_0=1$ (our recommendation) for the computation will overestimate or underestimate quantities of interest as functions of $x_{n+1}$ when the true ${\kappa }_0>0$ (e.g., 2 or $0.5$) is different from ${\kappa }_0=1$.

Quantities of interest as functions of $x_{n+1}$ for ${\kappa }_0=2,1,0.5$ in the simulations are plotted in Figure 2. In all the plots, $u_0={\sigma }_0^2\left({\kappa }_0+1\right)/{\kappa }_0=2$ is fixed. From the figure, we observe the following facts.

• From the left column of the figure, we see that ${\beta }^{*}$, ${\delta }_s^{\pi ,{\theta }_{n+1}}\left(x_{n+1}\right)$, ${\delta }_2^{\pi ,{\theta }_{n+1}}\left(x_{n+1}\right)$, and ${\delta }_2^{\pi ,{\mu }_{n+1}}\left(x_{n+1}\right)$ are affected by ${\kappa }_0$. Moreover, ${\beta }^{*}$, ${\delta }_s^{\pi ,{\theta }_{n+1}}\left(x_{n+1}\right)$, and ${\delta }_2^{\pi ,{\theta }_{n+1}}\left(x_{n+1}\right)$ are quadratic functions of $x_{n+1}$ (note that their y ranges are different), while ${\delta }_2^{\pi ,{\mu }_{n+1}}\left(x_{n+1}\right)$ is a linear function of $x_{n+1}$.
• For plot (a), we observe that ${\beta }^{*}$ when using ${\kappa }_0=1$ will underestimate ${\beta }^{*}$ when using ${\kappa }_0=2$, and it will overestimate ${\beta }^{*}$ when using ${\kappa }_0=0.5$. In other words, ${\beta }^{*}$ is an increasing function of ${\kappa }_0$ for a fixed $x_{n+1}$. Similar phenomena are observed for ${\delta }_s^{\pi ,{\theta }_{n+1}}\left(x_{n+1}\right)$ (plot (c)) and ${\delta }_2^{\pi ,{\theta }_{n+1}}\left(x_{n+1}\right)$ (plot (e)).
• For ${\delta }_2^{\pi ,{\mu }_{n+1}}\left(x_{n+1}\right)$ (plot (g)), we observe an interesting phenomenon that when $x_{n+1}<{\mu }_0$ ($=3$ in the simulation), ${\delta }_2^{\pi ,{\mu }_{n+1}}\left(x_{n+1}\right)$ when using ${\kappa }_0=1$ will underestimate that when using ${\kappa }_0=2$, and it will overestimate that when using ${\kappa }_0=0.5$. In other words, ${\delta }_2^{\pi ,{\mu }_{n+1}}\left(x_{n+1}\right)$ is an increasing function of ${\kappa }_0$ when $x_{n+1}<{\mu }_0$. However, when $x_{n+1}>{\mu }_0$, we observe a reversed phenomenon that ${\delta }_2^{\pi ,{\mu }_{n+1}}\left(x_{n+1}\right)$ when using ${\kappa }_0=1$ will overestimate that when using ${\kappa }_0=2$, and it will underestimate that when using ${\kappa }_0=0.5$. In other words, ${\delta }_2^{\pi ,{\mu }_{n+1}}\left(x_{n+1}\right)$ is a decreasing function of ${\kappa }_0$ when $x_{n+1}>{\mu }_0$. The phenomena can be explained by the theoretical analysis of

  $${\displaystyle \frac{\partial {\delta }_2^{\pi ,{\mu }_{n+1}}\left(x_{n+1}\right)}{\partial {\kappa }_0}}={\displaystyle \frac{{\mu }_0-x_{n+1}}{{\left({\kappa }_0+1\right)}^2}}\left\{\begin{array}{cc}>0,& \mathrm{if}{x}_{n+1}<{\mu }_0,\\ =0,& \mathrm{if}{x}_{n+1}={\mu }_0,\\ <0,& \mathrm{if}{x}_{n+1}>{\mu }_0.\end{array}\right.$$
• From the right column of the figure, we see that ${\alpha }^{*}$, $PES{L}_s^{\pi ,{\theta }_{n+1}}\left(x_{n+1}\right)$, and $PES{L}_2^{\pi ,{\theta }_{n+1}}\left(x_{n+1}\right)$ are not affected by ${\kappa }_0$. Moreover, they also do not depend on $x_{n+1}$.

3.3. Two Inequalities of the Bayes Estimators and the PESLs

This subsection aims to provide numerical illustrations of the two inequalities related to the Bayes estimators and the PESLs—specifically, (15) and (16)—under the oracle scenario where the hyperparameters $\left({\mu }_0,{\kappa }_0,{\nu }_0,{\sigma }_0^2\right)$ are known. The rationale for this section lies in the theoretical foundation established by these two inequalities, which we seek to validate through numerical examples.

First, we fix ${\mu }_0=3$, ${\kappa }_0=1$, ${\nu }_0=6$, and ${\sigma }_0^2=1$. Then we set a seed number 1 in R software and draw ${\theta }_{n+1}=1.866135$ from $\widetilde{IG}\left({\nu }_0/2,{\nu }_0{\sigma }_0^2/2\right)$. Next, we draw ${\mu }_{n+1}=4.816593$ from $N({\mu }_0,{\theta }_{n+1}/{\kappa }_0)$. After that, we draw $x_{n+1}=6.554815$ from $N\left({\mu }_{n+1},{\theta }_{n+1}\right)$. Figure 3 shows the histogram of ${\theta }_{n+1}|x_{n+1}$ and the density estimation curve of $\pi ({\theta }_{n+1}|x_{n+1})$. It is $\pi ({\theta }_{n+1}|x_{n+1})$ that we find ${\delta }_s^{\pi ,{\theta }_{n+1}}\left(x_{n+1}\right)$ to minimize the PESL. Numerical results show that

$${\delta }_s^{\pi ,{\theta }_{n+1}}\left(x_{n+1}\right)=1.759765<2.463671={\delta }_2^{\pi ,{\theta }_{n+1}}\left(x_{n+1}\right)$$

and

$$PES{L}_s^{\pi ,{\theta }_{n+1}}\left(x_{n+1}\right)=0.1496063<0.2131341=PES{L}_2^{\pi ,{\theta }_{n+1}}\left(x_{n+1}\right),$$

which exemplify the theoretical studies of (15) and (16).

We now consider the scenario where one of the five parameters—${\mu }_0$, ${\kappa }_0$, ${\nu }_0$, ${\sigma }_0^2$, and $x_{n+1}$—is allowed to vary while the others remain constant. Specifically, our focus lies on conducting a sensitivity analysis of the Bayes estimators and the PESLs with respect to these five parameters.

Figure 4 illustrates the Bayes estimators and the PESLs as varying with ${\nu }_0$. As shown in the left panel, the Bayes estimators exhibit dependence on ${\nu }_0$, demonstrating the validity of inequality (15). Specifically, ${\delta }_s^{\pi ,{\theta }_{n+1}}\left(x_{n+1}\right)$ increases as ${\nu }_0$ increases, whereas ${\delta }_2^{\pi ,{\theta }_{n+1}}\left(x_{n+1}\right)$ decreases with rising ${\nu }_0$. The right panel reveals that the PESLs are also influenced by ${\nu }_0$, supporting the inequality (16), and both PESLs decrease monotonically as ${\nu }_0$ grows. Additionally, numerical values corresponding to these quantities are provided in Table 3, which aligns with the graphical representation in Figure 4. Overall, both the figure and the table confirm the two theoretical inequalities: (15) and (16).

Figure 5 shows the Bayes estimators and the PESLs as functions of ${\mu }_0$, ${\kappa }_0$, ${\sigma }_0^2$, and $x_{n+1}$. We see from the left plots of the figure that the Bayes estimators depend on ${\mu }_0$, ${\kappa }_0$, ${\sigma }_0^2$, and $x_{n+1}$, and (15) is exemplified. Moreover, the Bayes estimators are first decreasing and then increasing functions of ${\mu }_0$ and $x_{n+1}$, and they are increasing functions of ${\kappa }_0$ and ${\sigma }_0^2$. The right plots of the figure exhibit that the PESLs do not depend on ${\mu }_0$, ${\kappa }_0$, ${\sigma }_0^2$, and $x_{n+1}$, and (16) is exemplified. Furthermore, Table 4, Table 5, Table 6 and Table 7 display the numerical values of the Bayes estimators and the PESLs in Figure 5. In summary, the results of Figure 5 and Table 4, Table 5, Table 6 and Table 7 exemplify the two inequalities (15) and (16).

In brief, the results of Figure 4 and Figure 5 exemplify that the PESLs depend only on ${\nu }_0$, but not on ${\mu }_0$, ${\kappa }_0$, ${\sigma }_0^2$, and $x_{n+1}$.

Since the estimators ${\delta }_2^{\pi ,{\theta }_{n+1}}\left(x_{n+1}\right)$ and ${\delta }_s^{\pi ,{\theta }_{n+1}}\left(x_{n+1}\right)$, along with their corresponding PESLs—$PES{L}_2^{\pi ,{\theta }_{n+1}}\left(x_{n+1}\right)$ and $PES{L}_s^{\pi ,{\theta }_{n+1}}\left(x_{n+1}\right)$—are functions of the parameters ${\alpha }^{*}$ and ${\beta }^{*}$, where ${\alpha }^{*}>1$ and ${\beta }^{*}>0$, we visualize their behavior across the domain $D=(1,10]\times (0,10]$ using the persp3d() function from the R package rgl (refer to [20,21,36,37]). A key advantage of persp3d() over the standard persp() function in the graphics package is its ability to overlay multiple surfaces within the same plot and enable interactive rotation of 3D visualizations, thereby offering greater flexibility in graphical analysis. The resulting surface plots for the estimators, the PESLs, and their pairwise differences are displayed in Figure 6. As shown in the two leftmost panels, ${\delta }_2^{\pi ,{\theta }_{n+1}}\left(x_{n+1}\right)$ consistently exceeds ${\delta }_s^{\pi ,{\theta }_{n+1}}\left(x_{n+1}\right)$ across all values of $\left({\alpha }^{*},{\beta }^{*}\right)\in D$, providing empirical support for the inequality stated in (15). Similarly, the rightmost panels indicate that $PES{L}_2^{\pi ,{\theta }_{n+1}}\left(x_{n+1}\right)$ remains greater than $PES{L}_s^{\pi ,{\theta }_{n+1}}\left(x_{n+1}\right)$ throughout D, which visually confirms the theoretical result in (16). Overall, the numerical illustrations in the figure corroborate the analytical findings presented in (15) and (16).

3.4. Goodness-of-Fit of the Model: KS Test

This section focuses on evaluating the goodness-of-fit of the hierarchical normal and normal-inverse-gamma model (2) to the generated synthetic data. The primary objective is to assess how well the marginal distribution of the N-NIG model captures the characteristics of the simulated dataset. It should be noted that the analysis here is based solely on the observed data points ${\mathit{x}}^{\left(n\right)}=\left(x_1,x_2,\dots ,x_n\right)$.

The Kolmogorov-Smirnov (KS) statistic measures the maximum discrepancy, denoted as D, between the empirical cumulative distribution function (cdf), $F_n\left(x\right)$, and the theoretical population cdf, $F_0\left(x\right)$. Mathematically, this is expressed as

$$D=\underset{-\infty <x<\infty }{sup}\left|F_n\left(x\right)-F_0\left(x\right)\right|.$$

(25)

In R, the KS test can be carried out using the built-in function ks.test() (refer to [38,39,40]). While the standard KS test is primarily designed for one-dimensional continuous distributions, various adaptations known as KS-type tests have also been proposed in the literature for handling discrete data ([41,42,43]).

The R function ks.test() performs the KS test, and its return value is a list containing the components statistic (the value of the test statistic, or the D value) and p.value (the p-value of the test). In the literature, it is well known that a smaller D value or a larger p-value indicates a better fit of the model to the simulated data. Based on and inspired by the D value and p-value, [21] propose five indices ($\overline{D}$, $\overline{\mathrm{p}-\mathrm{value}}$, $\left(minD\right)\%$, $\left(max\mathrm{p}-\mathrm{value}\right)\%$, $\left(\mathrm{Accept}{H}_0\right)\%$) to compare the three methods, that is, the oracle method, the moment method, and the MLE method in simulations. $\overline{D}$ is the average D values (25) of M simulations (the smaller the better). $\overline{\mathrm{p}-\mathrm{value}}$ is the average p-values of M simulations (the larger the better). $\left(minD\right)\%$ is the percentage of M simulations which attain the minimum D value in the three methods (the larger the better). The three $\left(minD\right)\%$ values should sum to 1. $\left(max\mathrm{p}-\mathrm{value}\right)\%$ is the percentage of M simulations which attain the maximum p-value in the three methods (the larger the better). The three $\left(max\mathrm{p}-\mathrm{value}\right)\%$ values should sum to 1. $\left(\mathrm{Accept}{H}_0\right)\%$ is the percentage of accepting $H_0$ (defined as p-value $>0.05$) in M simulations for each method (the larger the better). Each $\left(\mathrm{Accept}{H}_0\right)\%$ value should between $0\%-100\%$.

In our problem, the null hypothesis specifies that

$$H_0:X\sim N-N-\widetilde{IG}({\mu }_0,{\nu }_0,u_0),$$

where $N-N-\widetilde{IG}({\mu }_0,{\nu }_0,u_0)$ is the marginal distribution of the hierarchical normal and normal-inverse-gamma model (2). The marginal density of the $N-N-\widetilde{IG}({\mu }_0,{\nu }_0,u_0)$ distribution is given by (18) which is obviously one-dimensional and continuous. Therefore, the KS test can serve as an indicator of how well the model fits the data.

It should be noted that the data presented in this section are generated based on the hierarchical normal and normal-inverse-gamma model (2), where the hyperparameters are set to $\left({\mu }_0,{\kappa }_0,{\nu }_0,{\sigma }_0^2\right)=\left(3,1,6,1\right)$ and $u_0=2$. Alternative values for these hyperparameters may also be chosen depending on the context.

Table 8 presents the goodness-of-fit results from the KS test for model (2) applied to simulated data. The simulations are based on a sample size of $n=1,000$ and repeated $M=100$ times. Hyperparameter estimation for ${\mu }_0$, ${\nu }_0$, and $u_0$ is conducted using three distinct approaches. The first is the oracle method, which assumes prior knowledge of the true values: ${\mu }_0=3$, ${\nu }_0=6$, and $u_0=2$. The second approach employs the method of moments, utilizing moment estimators derived in Theorem 2. The third approach relies on MLE, where the hyperparameters are numerically computed using the MLEs provided in Theorem 3.

As shown in Table 8, the following observations can be made.

• The $\overline{D}$ values corresponding to the three approaches are $0.0268$, $0.0239$, and $0.0204$, respectively. This indicates that the MLE method performs best, followed by the moment method, while the oracle method yields the least satisfactory results. A plausible explanation for this outcome is that, as shown in (25), the empirical cdf $F_n\left(x\right)$ relies on observed data, and similarly, the theoretical cdf $F_0\left(x\right)$ used in both the moment and MLE methods are derived from data. In contrast, the $F_0\left(x\right)$ associated with the oracle method does not depend on actual data.
• The $\overline{\mathrm{p}-\mathrm{value}}$ corresponding to the three methods are $0.5230$, $0.6245$, and $0.7674$, respectively, indicating that the MLE method performs best, followed by the moment method, with the oracle method ranking last. A plausible explanation for this outcome is provided in the preceding paragraph.
• The $\left(minD\right)\%$ values for the three approaches are $0.20$, $0.16$, and $0.64$, respectively. For the MLE method, the $\left(minD\right)\%$ value exceeds half of the M simulations. In terms of performance ranking, the methods are ordered as follows: MLE method, oracle method, and moment method.
• The three methods yield $\left(max\mathrm{p}-\mathrm{value}\right)\%$ values of $0.20$, $0.16$, and $0.64$, respectively. Since a lower D value is associated with a higher p-value, the method producing the minimum D value will also produce the maximum p-value. As a result, the $\left(minD\right)\%$ and $\left(max\mathrm{p}-\mathrm{value}\right)\%$ values are identical across the three approaches. Based on performance, the methods can be ranked in descending order of preference as follows: MLE method, oracle method, and moment method.
• The $\left(\mathrm{Accept}{H}_0\right)\%$ values for the three methods are $98\%$, $99\%$, and $99\%$, respectively, indicating that these values are close to $100\%$. This suggests that all three methods exhibit strong performance in terms of model fit.
• To summarize, the MLE method consistently achieves the top ranking across all five evaluation metrics—$\overline{D}$, $\overline{\mathrm{p}-\mathrm{value}}$, $\left(minD\right)\%$, $\left(max\mathrm{p}-\mathrm{value}\right)\%$, and $\left(\mathrm{Accept}{H}_0\right)\%$. When comparing the moment method with the MLE method, the results indicate that the MLE method outperforms the moment method on each of these five criteria.

Figure 7 presents the boxplots of the D values and p-values obtained from the three methods. The following observations can be made based on this figure.

• The D values obtained by the oracle method are higher than those of the other two approaches. Since lower D values indicate better performance, the preferred ranking of the three methods is as follows: MLE method, moment method, and then the oracle method.
• The oracle method yields the lowest p-values compared to the other two approaches. Since higher p-values are preferable, the methods can be ranked in descending order of preference as follows: MLE method, moment method, and finally the oracle method.
• Large D values are associated with small p-values, whereas small D values are linked to large p-values.
• Based on the D values and p-values, the MLE method outperforms the method of moments.

3.5. Marginal Densities for Various Hyperparameters

This subsection presents the marginal density plots of the hierarchical normal and normal-inverse-gamma model (2) under different settings of the hyperparameters ${\mu }_0$, ${\kappa }_0$, ${\nu }_0$, and ${\sigma }_0^2$. The primary objective is to investigate the types of data distributions that can be captured by this hierarchical model. Recall that the marginal density of X is defined by equation (6), which depends on four hyperparameters: ${\mu }_0$ (ranging over the real line), ${\kappa }_0>0$, ${\nu }_0>0$, and ${\sigma }_0^2>0$. Our analysis focuses on how variations in these hyperparameters affect the shape of the marginal density, using the baseline configuration ${\mu }_0=3$, ${\kappa }_0=1$, ${\nu }_0=6$, and ${\sigma }_0^2=1$ as a reference point. Alternative values for the hyperparameters may also be considered to further explore the model’s flexibility.

Figure 8 displays the marginal density functions under different values of ${\mu }_0=-3,0,3,6,9$, with ${\kappa }_0=1$, ${\nu }_0=6$, and ${\sigma }_0^2=1$ held constant. It can be observed that as ${\mu }_0$ increases, the entire distribution shifts rightward without altering its shape, indicating that ${\mu }_0$ acts as a location parameter. Furthermore, each marginal density is symmetric around its corresponding mean value ${\mu }_0$.

Figure 9 illustrates the marginal density functions under different values of ${\kappa }_0=0.25,0.5,1,2,4$, with ${\mu }_0=3$, ${\nu }_0=6$, and ${\sigma }_0^2=1$ held constant. It can be observed that as ${\kappa }_0$ grows larger, the peak of the marginal density becomes more pronounced, indicating a reduction in its spread. This is consistent with the theoretical expression for the variance of X, given by

$$\mathrm{Var}\left(X\right)={\displaystyle \frac{{\nu }_0}{{\nu }_0-2}}{\displaystyle \frac{{\sigma }_0^2\left({\kappa }_0+1\right)}{{\kappa }_0}},$$

(26)

which decreases as ${\kappa }_0$ increases. Additionally, each marginal density remains symmetric around the common mean ${\mu }_0$.

Figure 10 displays the marginal density functions under different values of ${\nu }_0=2,4,6,12,24$, with ${\mu }_0=3$, ${\kappa }_0=1$, and ${\sigma }_0^2=1$ held constant. It can be observed that as ${\nu }_0$ grows larger, the height of the density peak rises, indicating a reduction in the spread of the distribution. This is consistent with (26), which shows that the variance of the marginal density decreases as ${\nu }_0$ increases. Additionally, each marginal density exhibits symmetry around the common mean ${\mu }_0$.

Figure 11 displays the marginal density functions under different values of ${\sigma }_0^2=0.25,0.5,1,2,4$, with ${\mu }_0=3$, ${\kappa }_0=1$, and ${\nu }_0=6$ held constant. As observed in the plot, an increase in ${\sigma }_0^2$ leads to a reduction in the height of the density peak, indicating a broader distribution. This reflects a larger variance in the marginal density, which aligns with (26) showing that variance increases with ${\sigma }_0^2$. Additionally, each marginal density remains symmetric around the mean value ${\mu }_0$.

4. A Real Data Example

In this section, we exploit the poverty level data. The data represent percentages of all persons below the poverty level. The sample is from a random collection of $80\left(=n+1\right)$ cities in the Western U.S. Source: County and City Data Book 12th edition, U.S. Department of Commerce. The data can be downloaded from

https://college.cengage.com/mathematics/brase/understandable_statistics/7e/students/datasets/svls/frames/frame.html

under section title “Poverty Level". Note that there is an “_" between understandable and statistics.

The histogram of the sample ${\mathit{x}}^{\left(n+1\right)}$ and its density estimation curveis depicted in Figure 12. From the figure, we see that the data are roughly symmetric about $0.15$.

The estimators of the hyperparameters ${\mu }_0$, ${\nu }_0$, and $u_0$, the goodness-of-fit of the model, the empirical Bayes estimators of the variance parameter of the normal distribution with a conjugate normal-inverse-gamma prior and the PESLs, and the mean and variance of the real data by the moment method and the MLE method are summarized in Table 9. From the table, we observe the following facts.

• The moment estimator of the hyperparameter ${\mu }_0$ is equal to the sample mean $\overline{x}=0.1669620$ of the first n observations. It is worth noting that the MLE of the hyperparameter ${\mu }_0$ is equal to $0.1619865$, which is very close to the moment estimator of the hyperparameter ${\mu }_0$. Furthermore, the moment estimator and the MLE of the hyperparameter $u_0$ are also very similar, and they are close to $0.006$. However, the moment estimator and the MLE of the hyperparameter ${\nu }_0$ are quite different. This does not mean that the moment estimator and the MLE are not consistent estimators of the hyperparameter ${\nu }_0$, nor mean that the hierarchical normal and normal-inverse-gamma model (2) does not fit the real data. For the big difference between the two estimators, the reason is that the sample size $n=79$ is small. Certainly, the MLE of the hyperparameter ${\nu }_0$ is more reliable, as guaranteed from the previous figures and tables in the simulations section.
• The KS test is used as a measure of the goodness-of-fit. The p-value of the moment method is $0.5313>0.05$, which implies that the $t_{{\nu }_0}\left({\mu }_0,u_0\right)$ distribution with ${\mu }_0$, ${\nu }_0$, and $u_0$ estimated by their moment estimators fits the sample ${\mathit{x}}^{\left(n\right)}=\left(x_1,x_2,\dots ,x_n\right)$ well. Moreover, the p-value of the MLE method is $0.8233>0.05$, which implies that the $t_{{\nu }_0}\left({\mu }_0,u_0\right)$ distribution with ${\mu }_0$, ${\nu }_0$, and $u_0$ estimated by their MLEs fits the sample ${\mathit{x}}^{\left(n\right)}$ even better. Comparing the moment method and the MLE method, we observe that the D value of the MLEs is smaller, and the p-value of the MLEs is larger, which means that the $t_{{\nu }_0}\left({\mu }_0,u_0\right)$ distribution with the hyperparameters ${\mu }_0$, ${\nu }_0$, and $u_0$ estimated by the MLEs has a better fit to the sample ${\mathit{x}}^{\left(n\right)}$ than that estimated by the moment estimators.
• When the hyperparameters are estimated by the MLE method, we see that

  $${\delta }_s^{\pi ,{\theta }_{n+1}}\left({\mathit{x}}^{\left(n+1\right)}\right)=0.0032576<0.0038348={\delta }_2^{\pi ,{\theta }_{n+1}}\left({\mathit{x}}^{\left(n+1\right)}\right)$$

  and

  $$PES{L}_s^{\pi ,{\theta }_{n+1}}\left({\mathit{x}}^{\left(n+1\right)}\right)=0.0771456<0.0912061=PES{L}_2^{\pi ,{\theta }_{n+1}}\left({\mathit{x}}^{\left(n+1\right)}\right).$$
  We observe a similar phenomenon for the Bayes estimators and the PESLs when the hyperparameters are estimated by the moment method. Therefore, the two inequalities (15) and (16) are exemplified.
• The mean of X (the poverty level data) is estimated by $\mathrm{E}X\approx {\widehat{\mu }}_0$. The variance of X is estimated by $\mathrm{Var}\left(X\right)\approx {\displaystyle \frac{{\widehat{\nu }}_0}{{\widehat{\nu }}_0-2}}{\widehat{u}}_0$. It is useful to note that the mean and variance of X by the two methods are very similar, although the estimators of the hyperparameters are very different. Moreover, it is worth mentioning that

  $$\begin{array}{cc}\hfill \mathrm{E}X& \approx 0.1619865=16.19865\%,\hfill \\ \hfill \mathrm{Var}\left(X\right)& \approx 0.0080536=80.536{\%}^2\hfill \end{array}$$

  for the MLE method. For the moment method, the mean and variance of X are similar. Hence, the variance of X is very small, not large!

5. Conclusions and Discussions

For the hierarchical normal and normal-inverse-gamma (2) framework, we begin by deriving the posterior distributions $\pi \left({\mu }_{n+1},{\theta }_{n+1}|x_{n+1}\right)$, $\pi \left({\theta }_{n+1}|x_{n+1}\right)$, $\pi \left({\mu }_{n+1}|x_{n+1}\right)$, and $\pi \left({\mu }_{n+1}|{\theta }_{n+1},x_{n+1}\right)$, along with the marginal distribution $m\left(x_{n+1}|{\mu }_0,{\kappa }_0,{\nu }_0,{\sigma }_0^2\right)$, as presented in Theorem 1. Subsequently, we obtain the Bayes estimators ${\delta }_s^{\pi ,{\theta }_{n+1}}\left(x_{n+1}\right)$ and ${\delta }_2^{\pi ,{\theta }_{n+1}}\left(x_{n+1}\right)$, as well as the corresponding PESLs, denoted by $PES{L}_s^{\pi ,{\theta }_{n+1}}\left(x_{n+1}\right)$ and $PES{L}_2^{\pi ,{\theta }_{n+1}}\left(x_{n+1}\right)$. These quantities fulfill two key inequalities, given in (15) and (16). Additionally, Theorem 2 provides the moment-based estimates for the model’s hyperparameters—specifically ${\mu }_{01}\left(n\right)$, ${\nu }_{01}\left(n\right)$, and $u_{01}\left(n\right)$. In parallel, Theorem 3 summarizes the MLEs of the same hyperparameters, namely ${\mu }_{02}\left(n\right)$, ${\nu }_{02}\left(n\right)$, and $u_{02}\left(n\right)$. Lastly, Theorem 4 compiles the empirical Bayes estimators for the model’s variance component under Stein’s loss function, derived using both estimation approaches.

In the simulation study, we conduct numerical experiments on the hierarchical normal and normal-inverse-gamma model (2) from five perspectives: the consistency of moment estimators and MLEs, the influence of ${\kappa }_0$ on quantities of interest as functions of $x_{n+1}$, the validation of two inequalities related to Bayes estimators and PESLs, the model’s goodness-of-fit to the generated data, and the examination of the model’s marginal density distributions.

In the real data analysis, we utilize poverty level data reflecting the proportion of individuals living below the poverty threshold. Table 9 presents a summary of several key results: estimates of the hyperparameters ${\mu }_0$, ${\nu }_0$, and $u_0$; model goodness-of-fit measures; empirical Bayes estimates for the variance parameter under a normal distribution with a conjugate normal-inverse-gamma prior; PESL values; as well as the mean and variance estimates of the poverty data computed via both the method of moments and MLE. The MLE-based estimates demonstrate a superior fit to the observed sample ${\mathit{x}}^{\left(n\right)}$ compared to moment-based estimates, as evidenced by a smaller D statistic and a higher corresponding p-value, indicating that the $t_{{\nu }_0}\left({\mu }_0,u_0\right)$ distribution using MLE-derived hyperparameters better captures the characteristics of the data.

In empirical Bayes approaches, hyperparameters are treated as unknown quantities, and their estimation relies on the marginal distribution derived from observed data. Two widely used techniques for this purpose are the method of moments and MLE. This study applies both methods to estimate the hyperparameters within a hierarchical normal and normal-inverse-gamma model (2).

The simulation results presented in the preceding figures and tables demonstrate that MLEs outperform moment-based estimators in estimating the hyperparameters ${\mu }_0$, ${\nu }_0$, and $u_0$, particularly with respect to consistency and model fit. Nevertheless, these advantages come at a cost. Unlike moment estimators, MLEs are prone to numerical instability, impose higher computational demands, necessitate the constraints ${\nu }_0>0$ and $u_0>0$ throughout all iterations, and lack closed-form solutions. Importantly, the performance of MLEs is highly sensitive to the choice of initial values, and moment estimators have been empirically shown to serve as reliable starting points. Furthermore, in scenarios where MLEs fail to converge or do not exist, moment estimators provide a justifiable and robust alternative.

Finally, we outline potential directions for future research. One possible extension involves adapting the hierarchical normal and normal-inverse-gamma model (2) to accommodate various non-conjugate prior distributions for the parameters of the normal distribution (refer to [1,44,45] and related works). In these cases, closed-form solutions may not be feasible, necessitating the development of numerical methods to compute the corresponding estimators.

Appendix A.1. The Proof of Theorem 1

Let ${\mathit{\eta }}_0={\left({\mu }_0,{\kappa }_0,{\nu }_0,{\sigma }_0^2\right)}^{\prime }$ be the hyperparameters. The marginal distribution of $x_{n+1}$ in model (2) is

$$m\left(x_{n+1}|{\mathit{\eta }}_0\right)={\displaystyle \frac{f\left(x_{n+1}|{\mu }_{n+1},{\theta }_{n+1}\right)\pi \left({\mu }_{n+1},{\theta }_{n+1}|{\mathit{\eta }}_0\right)}{\pi \left({\mu }_{n+1},{\theta }_{n+1}|x_{n+1},{\mathit{\eta }}_0\right)}}.$$

To lighten notations, the ${\mathit{\eta }}_0$ will be dropped in the densities. Some of the following derivations are quoted from Example 1.5.1 (p.20) of [25] and [19].

For the random variables, parameters, and hyperparameters, their domains are respectively given by

$$-\infty <x_{n+1},{\mu }_{n+1},{\mu }_0<\infty$$

and

$${\theta }_{n+1},{\kappa }_0,{\nu }_0,{\sigma }_0^2>0.$$

By the Bayes Theorem, the joint posterior distribution of ${\mu }_{n+1}$ and ${\theta }_{n+1}$ is

$$\pi \left({\mu }_{n+1},{\theta }_{n+1}|x_{n+1}\right)\propto f\left(x_{n+1}|{\mu }_{n+1},{\theta }_{n+1}\right)\pi \left({\mu }_{n+1},{\theta }_{n+1}\right).$$

The joint conjugate prior distribution of ${\mu }_{n+1}$ and ${\theta }_{n+1}$ is decomposed as

$$\pi \left({\mu }_{n+1},{\theta }_{n+1}\right)=\pi \left({\mu }_{n+1}|{\theta }_{n+1}\right)\pi \left({\theta }_{n+1}\right),$$

which is a normal-inverse-gamma distribution. Hence,

$$\pi \left({\mu }_{n+1},{\theta }_{n+1}|x_{n+1}\right)\propto f\left(x_{n+1}|{\mu }_{n+1},{\theta }_{n+1}\right)\pi \left({\mu }_{n+1}|{\theta }_{n+1}\right)\pi \left({\theta }_{n+1}\right).$$

It is easy to see that

$$\begin{array}{cc}\hfill f\left(x_{n+1}|{\mu }_{n+1},{\theta }_{n+1}\right)& ={\displaystyle \frac{1}{\sqrt{2\pi }\sqrt{{\theta }_{n+1}}}}exp\left[-{\displaystyle \frac{1}{2{\theta }_{n+1}}}{\left(x_{n+1}-{\mu }_{n+1}\right)}^2\right]\hfill \\ \hfill & \propto {\theta }_{n+1}^{-{\displaystyle \frac{1}{2}}}exp\left[-{\displaystyle \frac{1}{2{\theta }_{n+1}}}{\left(x_{n+1}-{\mu }_{n+1}\right)}^2\right],\hfill \end{array}$$

$$\begin{array}{cc}\hfill \pi \left({\mu }_{n+1}|{\theta }_{n+1}\right)& ={\displaystyle \frac{1}{\sqrt{2\pi }\sqrt{{\theta }_{n+1}/{\kappa }_0}}}exp\left[-{\displaystyle \frac{{\kappa }_0}{2{\theta }_{n+1}}}{\left({\mu }_{n+1}-{\mu }_0\right)}^2\right]\hfill \\ \hfill & \propto {\theta }_{n+1}^{-{\displaystyle \frac{1}{2}}}exp\left[-{\displaystyle \frac{{\kappa }_0}{2{\theta }_{n+1}}}{\left({\mu }_{n+1}-{\mu }_0\right)}^2\right],\hfill \end{array}$$

and

$$\begin{array}{cc}\hfill \pi \left({\theta }_{n+1}\right)& ={\displaystyle \frac{{\left(\frac{{\nu }_0{\sigma }_0^2}{2}\right)}^{\frac{{\nu }_0}{2}}}{\Gamma \left(\frac{{\nu }_0}{2}\right)}}{\left({\displaystyle \frac{1}{{\theta }_{n+1}}}\right)}^{{\displaystyle \frac{{\nu }_0}{2}}+1}exp\left(-{\displaystyle \frac{{\nu }_0{\sigma }_0^2}{2}}{\displaystyle \frac{1}{{\theta }_{n+1}}}\right)\hfill \\ \hfill & \propto {\left({\displaystyle \frac{1}{{\theta }_{n+1}}}\right)}^{{\displaystyle \frac{{\nu }_0}{2}}+1}exp\left(-{\displaystyle \frac{{\nu }_0{\sigma }_0^2}{2}}{\displaystyle \frac{1}{{\theta }_{n+1}}}\right).\hfill \end{array}$$

Therefore,

$$\begin{array}{cc}\hfill & \pi \left({\mu }_{n+1},{\theta }_{n+1}|x_{n+1}\right)\hfill \\ \hfill & \propto {\theta }_{n+1}^{-{\displaystyle \frac{1}{2}}}exp\left[-{\displaystyle \frac{1}{2{\theta }_{n+1}}}{\left(x_{n+1}-{\mu }_{n+1}\right)}^2\right]{\theta }_{n+1}^{-{\displaystyle \frac{1}{2}}}exp\left[-{\displaystyle \frac{{\kappa }_0}{2{\theta }_{n+1}}}{\left({\mu }_{n+1}-{\mu }_0\right)}^2\right]\hfill \\ \hfill & \times {\theta }_{n+1}^{-\left({\displaystyle \frac{{\nu }_0}{2}}+1\right)}exp\left(-{\displaystyle \frac{{\nu }_0{\sigma }_0^2}{2{\theta }_{n+1}}}\right)\hfill \\ \hfill & ={\theta }_{n+1}^{-{\displaystyle \frac{1}{2}}}{\theta }_{n+1}^{-\left({\displaystyle \frac{{\nu }_0+1}{2}}+1\right)}exp\left\{-{\displaystyle \frac{1}{2{\theta }_{n+1}}}\left[{\left(x_{n+1}-{\mu }_{n+1}\right)}^2+{\kappa }_0{\left({\mu }_{n+1}-{\mu }_0\right)}^2+{\nu }_0{\sigma }_0^2\right]\right\}\hfill \end{array}$$

(A1)

The expression in the square brackets of (A1) changes to

$$\begin{array}{cc}\hfill & {\left(x_{n+1}-{\mu }_{n+1}\right)}^2+{\kappa }_0{\left({\mu }_{n+1}-{\mu }_0\right)}^2+{\nu }_0{\sigma }_0^2\hfill \\ \hfill & ={\mu }_{n+1}^2-2x_{n+1}{\mu }_{n+1}+x_{n+1}^2+{\kappa }_0{\mu }_{n+1}^2-2{\kappa }_0{\mu }_0{\mu }_{n+1}+{\kappa }_0{\mu }_0^2+{\nu }_0{\sigma }_0^2\hfill \\ \hfill & =\left({\kappa }_0+1\right){\mu }_{n+1}^2-2\left(x_{n+1}+{\kappa }_0{\mu }_0\right){\mu }_{n+1}+x_{n+1}^2+{\kappa }_0{\mu }_0^2+{\nu }_0{\sigma }_0^2\hfill \\ \hfill & =\left({\kappa }_0+1\right)\left[{\mu }_{n+1}^2-2{\displaystyle \frac{x_{n+1}+{\kappa }_0{\mu }_0}{{\kappa }_0+1}}{\mu }_{n+1}\right]+x_{n+1}^2+{\kappa }_0{\mu }_0^2+{\nu }_0{\sigma }_0^2\hfill \\ \hfill & =\left({\kappa }_0+1\right){\left({\mu }_{n+1}-{\displaystyle \frac{x_{n+1}+{\kappa }_0{\mu }_0}{{\kappa }_0+1}}\right)}^2-{\displaystyle \frac{{\left(x_{n+1}+{\kappa }_0{\mu }_0\right)}^2}{{\kappa }_0+1}}+x_{n+1}^2+{\kappa }_0{\mu }_0^2+{\nu }_0{\sigma }_0^2\hfill \\ \hfill & =\left({\kappa }_0+1\right){\left({\mu }_{n+1}-{\displaystyle \frac{x_{n+1}+{\kappa }_0{\mu }_0}{{\kappa }_0+1}}\right)}^2+{\displaystyle \frac{{\kappa }_0{\left(x_{n+1}-{\mu }_0\right)}^2}{{\kappa }_0+1}}+{\nu }_0{\sigma }_0^2,\hfill \end{array}$$

where

$$\begin{array}{cc}\hfill & -{\displaystyle \frac{{\left(x_{n+1}+{\kappa }_0{\mu }_0\right)}^2}{{\kappa }_0+1}}+x_{n+1}^2+{\kappa }_0{\mu }_0^2\hfill \\ \hfill & ={\displaystyle \frac{\left({\kappa }_0+1\right)x_{n+1}^2+\left({\kappa }_0+1\right){\kappa }_0{\mu }_0^2-{\left(x_{n+1}+{\kappa }_0{\mu }_0\right)}^2}{{\kappa }_0+1}}\hfill \\ \hfill & ={\displaystyle \frac{\left({\kappa }_0+1\right)x_{n+1}^2+\left({\kappa }_0+1\right){\kappa }_0{\mu }_0^2-x_{n+1}^2-{\kappa }_0^2{\mu }_0^2-2{\kappa }_0{\mu }_0{x}_{n+1}}{{\kappa }_0+1}}\hfill \\ \hfill & ={\displaystyle \frac{{\kappa }_0{x}_{n+1}^2+{\kappa }_0{\mu }_0^2-2{\kappa }_0{\mu }_0{x}_{n+1}}{{\kappa }_0+1}}\hfill \\ \hfill & ={\displaystyle \frac{{\kappa }_0\left(x_{n+1}^2+{\mu }_0^2-2x_{n+1}{\mu }_0\right)}{{\kappa }_0+1}}\hfill \\ \hfill & ={\displaystyle \frac{{\kappa }_0{\left(x_{n+1}-{\mu }_0\right)}^2}{{\kappa }_0+1}}.\hfill \end{array}$$

Let

$$\left\{\begin{array}{c}{\mu }_1=\left(x_{n+1}+{\kappa }_0{\mu }_0\right)/\left({\kappa }_0+1\right),\hfill \\ {\kappa }_1={\kappa }_0+1,\hfill \\ {\nu }_1={\nu }_0+1,\hfill \\ {\nu }_1{\sigma }_1^2={\nu }_0{\sigma }_0^2+{\displaystyle \frac{{\kappa }_0}{{\kappa }_0+1}}{\left(x_{n+1}-{\mu }_0\right)}^2.\hfill \end{array}\right.$$

(A2)

Then (A1) reduces to

$$\begin{array}{cc}\hfill & \pi \left({\mu }_{n+1},{\theta }_{n+1}|x_{n+1}\right)\hfill \\ \hfill & \propto {\theta }_{n+1}^{-{\displaystyle \frac{1}{2}}}{\theta }_{n+1}^{-\left({\displaystyle \frac{{\nu }_1}{2}}+1\right)}exp\left\{-{\displaystyle \frac{1}{2{\theta }_{n+1}}}\left[{\kappa }_1{\left({\mu }_{n+1}-{\mu }_1\right)}^2+{\nu }_1{\sigma }_1^2\right]\right\}.\hfill \end{array}$$

It is shown that $\pi \left({\mu }_{n+1},{\theta }_{n+1}|x_{n+1}\right)$ is a normal-inverse-gamma distribution as follows. The joint posterior distribution $\pi \left({\mu }_{n+1},{\theta }_{n+1}|x_{n+1}\right)$ can be written as $\pi \left({\mu }_{n+1}|{\theta }_{n+1},x_{n+1}\right)\pi \left({\theta }_{n+1}|x_{n+1}\right)$, where

$$\left\{\begin{array}{c}\pi \left({\mu }_{n+1}|{\theta }_{n+1},x_{n+1}\right){\left({\displaystyle \frac{{\theta }_{n+1}}{{\kappa }_1}}\right)}^{-1/2}exp\left[-{\displaystyle \frac{{\kappa }_1}{2{\theta }_{n+1}}}{\left({\mu }_{n+1}-{\mu }_1\right)}^2\right],\hfill \\ \pi \left({\theta }_{n+1}|x_{n+1}\right){\theta }_{n+1}^{-\left({\displaystyle \frac{{\nu }_1}{2}}+1\right)}exp\left(-{\displaystyle \frac{1}{{\theta }_{n+1}}}{\displaystyle \frac{{\nu }_1{\sigma }_1^2}{2}}\right).\hfill \end{array}\right.$$

That is,

$$\left\{\begin{array}{c}{\mu }_{n+1}|{\theta }_{n+1},x_{n+1}\sim N\left({\mu }_1,{\theta }_{n+1}/{\kappa }_1\right),\hfill \\ {\theta }_{n+1}|x_{n+1}\sim IG\left({\alpha }^{*}={\nu }_1/2,{\beta }^{*}={\nu }_1{\sigma }_1^2/2\right).\hfill \end{array}\right.$$

Therefore, the joint posterior distribution $\pi \left({\mu }_{n+1},{\theta }_{n+1}|x_{n+1}\right)$ is

$$\begin{array}{cc}\hfill & \pi \left({\mu }_{n+1},{\theta }_{n+1}|x_{n+1}\right)\hfill \\ \hfill & =\pi \left({\mu }_{n+1}|{\theta }_{n+1},x_{n+1}\right)\pi \left({\theta }_{n+1}|x_{n+1}\right)\hfill \\ \hfill & ={\displaystyle \frac{1}{\sqrt{2\pi }\sqrt{{\theta }_{n+1}/{\kappa }_1}}}exp\left[-{\displaystyle \frac{{\kappa }_1}{2{\theta }_{n+1}}}{\left({\mu }_{n+1}-{\mu }_1\right)}^2\right]\hfill \\ \hfill & \times {\displaystyle \frac{{\left(\frac{{\nu }_1{\sigma }_1^2}{2}\right)}^{\frac{{\nu }_1}{2}}}{\Gamma \left(\frac{{\nu }_1}{2}\right)}}{\theta }_{n+1}^{-\left({\displaystyle \frac{{\nu }_1}{2}}+1\right)}exp\left(-{\displaystyle \frac{{\nu }_1{\sigma }_1^2}{2}}{\displaystyle \frac{1}{{\theta }_{n+1}}}\right)\hfill \\ \hfill & ={\displaystyle \frac{\sqrt{{\kappa }_1}}{\sqrt{2\pi }}}{\displaystyle \frac{{\left(\frac{{\nu }_1{\sigma }_1^2}{2}\right)}^{\frac{{\nu }_1}{2}}}{\Gamma \left(\frac{{\nu }_1}{2}\right)}}{\theta }_{n+1}^{-{\displaystyle \frac{1}{2}}}{\theta }_{n+1}^{-\left({\displaystyle \frac{{\nu }_1}{2}}+1\right)}exp\left\{-{\displaystyle \frac{1}{2{\theta }_{n+1}}}\left[{\kappa }_1{\left({\mu }_{n+1}-{\mu }_1\right)}^2+{\nu }_1{\sigma }_1^2\right]\right\}\hfill \\ \hfill & \sim N-\widetilde{IG}({\mu }_1,{\kappa }_1,{\nu }_1,{\sigma }_1^2).\hfill \end{array}$$

(A3)

The joint prior distribution $\pi \left({\mu }_{n+1},{\theta }_{n+1}\right)$ is

$$\begin{array}{cc}\hfill & \pi \left({\mu }_{n+1},{\theta }_{n+1}\right)\hfill \\ \hfill & =\pi \left({\mu }_{n+1}|{\theta }_{n+1}\right)\pi \left({\theta }_{n+1}\right)\hfill \\ \hfill & ={\displaystyle \frac{1}{\sqrt{2\pi }\sqrt{{\theta }_{n+1}/{\kappa }_0}}}exp\left[-{\displaystyle \frac{{\kappa }_0}{2{\theta }_{n+1}}}{\left({\mu }_{n+1}-{\mu }_0\right)}^2\right]\hfill \\ \hfill & \times {\displaystyle \frac{{\left(\frac{{\nu }_0{\sigma }_0^2}{2}\right)}^{\frac{{\nu }_0}{2}}}{\Gamma \left(\frac{{\nu }_0}{2}\right)}}{\left({\displaystyle \frac{1}{{\theta }_{n+1}}}\right)}^{{\displaystyle \frac{{\nu }_0}{2}}+1}exp\left(-{\displaystyle \frac{{\nu }_0{\sigma }_0^2}{2}}{\displaystyle \frac{1}{{\theta }_{n+1}}}\right)\hfill \\ \hfill & ={\displaystyle \frac{\sqrt{{\kappa }_0}}{\sqrt{2\pi }}}{\displaystyle \frac{{\left(\frac{{\nu }_0{\sigma }_0^2}{2}\right)}^{\frac{{\nu }_0}{2}}}{\Gamma \left(\frac{{\nu }_0}{2}\right)}}{\theta }_{n+1}^{-{\displaystyle \frac{1}{2}}}{\theta }_{n+1}^{-\left({\displaystyle \frac{{\nu }_0}{2}}+1\right)}exp\left\{-{\displaystyle \frac{1}{2{\theta }_{n+1}}}\left[{\kappa }_0{\left({\mu }_{n+1}-{\mu }_0\right)}^2+{\nu }_0{\sigma }_0^2\right]\right\}\hfill \\ \hfill & \sim N-\widetilde{IG}({\mu }_0,{\kappa }_0,{\nu }_0,{\sigma }_0^2).\hfill \end{array}$$

(A4)

Comparing (A4) and (A3), we find that the normal-inverse-gamma distribution is a conjugate prior for $\left({\mu }_{n+1},{\theta }_{n+1}\right)$ of $N\left({\mu }_{n+1},{\theta }_{n+1}\right)$.

Now let us derive the marginal posterior density of ${\mu }_{n+1}$. We have

$$\begin{array}{cc}\hfill & \pi \left({\mu }_{n+1}|x_{n+1}\right)\hfill \\ \hfill & ={\int }_0^{\infty }\pi \left({\mu }_{n+1},{\theta }_{n+1}|x_{n+1}\right)d{\theta }_{n+1}\hfill \\ \hfill & ={\int }_0^{\infty }{\displaystyle \frac{\sqrt{{\kappa }_1}}{\sqrt{2\pi }}}{\displaystyle \frac{{\left(\frac{{\nu }_1{\sigma }_1^2}{2}\right)}^{\frac{{\nu }_1}{2}}}{\Gamma \left(\frac{{\nu }_1}{2}\right)}}{\theta }_{n+1}^{-\left({\displaystyle \frac{{\nu }_1+1}{2}}+1\right)}exp\left\{-{\displaystyle \frac{1}{2{\theta }_{n+1}}}\left[{\kappa }_1{\left({\mu }_{n+1}-{\mu }_1\right)}^2+{\nu }_1{\sigma }_1^2\right]\right\}d{\theta }_{n+1}\hfill \\ \hfill & ={\displaystyle \frac{\sqrt{{\kappa }_1}}{\sqrt{2\pi }}}{\displaystyle \frac{{\left(\frac{{\nu }_1{\sigma }_1^2}{2}\right)}^{\frac{{\nu }_1}{2}}}{\Gamma \left(\frac{{\nu }_1}{2}\right)}}{\int }_0^{\infty }{\theta }_{n+1}^{-\left({\displaystyle \frac{{\nu }_1+1}{2}}+1\right)}exp\left\{-{\displaystyle \frac{1}{2{\theta }_{n+1}}}\left[{\kappa }_1{\left({\mu }_{n+1}-{\mu }_1\right)}^2+{\nu }_1{\sigma }_1^2\right]\right\}d{\theta }_{n+1}\hfill \\ \hfill & ={\displaystyle \frac{\sqrt{{\kappa }_1}}{\sqrt{2\pi }}}{\displaystyle \frac{{\left(\frac{{\nu }_1{\sigma }_1^2}{2}\right)}^{\frac{{\nu }_1}{2}}}{\Gamma \left(\frac{{\nu }_1}{2}\right)}}{\displaystyle \frac{\Gamma \left(\frac{{\nu }_1+1}{2}\right)}{{\left[\frac{{\kappa }_1{\left({\mu }_{n+1}-{\mu }_1\right)}^2+{\nu }_1{\sigma }_1^2}{2}\right]}^{\frac{{\nu }_1+1}{2}}}},\hfill \end{array}$$

by noting that the integrand of the above integral is the kernel of an inverse gamma distribution $\widetilde{IG}\left({\alpha }_1,{\beta }_1\right)$ with

$$\begin{array}{cc}\hfill {\alpha }_1& ={\displaystyle \frac{{\nu }_1+1}{2}},\hfill \\ \hfill {\beta }_1& ={\displaystyle \frac{{\kappa }_1{\left({\mu }_{n+1}-{\mu }_1\right)}^2+{\nu }_1{\sigma }_1^2}{2}}.\hfill \end{array}$$

Therefore,

$$\begin{array}{cc}\hfill & \pi \left({\mu }_{n+1}|x_{n+1}\right)\hfill \\ \hfill & ={\displaystyle \frac{\Gamma \left(\frac{{\nu }_1+1}{2}\right)}{\Gamma \left(\frac{{\nu }_1}{2}\right)}}{\displaystyle \frac{\sqrt{{\kappa }_1}}{\sqrt{2\pi }}}{\left({\displaystyle \frac{{\nu }_1{\sigma }_1^2}{2}}\right)}^{{\displaystyle \frac{{\nu }_1}{2}}}{\displaystyle \frac{1}{{\left(\frac{{\nu }_1{\sigma }_1^2}{2}\right)}^{\frac{{\nu }_1+1}{2}}{\left[1+\frac{1}{{\nu }_1}\frac{{\left({\mu }_{n+1}-{\mu }_1\right)}^2}{{\sigma }_1^2/{\kappa }_1}\right]}^{\frac{{\nu }_1+1}{2}}}}\hfill \\ \hfill & ={\displaystyle \frac{\Gamma \left(\frac{{\nu }_1+1}{2}\right)}{\Gamma \left(\frac{{\nu }_1}{2}\right)\sqrt{\pi {\nu }_1{\sigma }_1^2/{\kappa }_1}}}{\left[1+{\displaystyle \frac{1}{{\nu }_1}}{\displaystyle \frac{{\left({\mu }_{n+1}-{\mu }_1\right)}^2}{{\sigma }_1^2/{\kappa }_1}}\right]}^{-{\displaystyle \frac{{\nu }_1+1}{2}}}\hfill \\ \hfill & \sim t_{{\nu }_1}\left({\mu }_1,{\displaystyle \frac{{\sigma }_1^2}{{\kappa }_1}}\right).\hfill \end{array}$$

Now let us calculate

$$\begin{array}{cc}\hfill & f\left(x_{n+1}|{\mu }_{n+1},{\theta }_{n+1}\right)\pi \left({\mu }_{n+1},{\theta }_{n+1}\right)\hfill \\ \hfill & ={\displaystyle \frac{1}{\sqrt{2\pi }\sqrt{{\theta }_{n+1}}}}exp\left[-{\displaystyle \frac{1}{2{\theta }_{n+1}}}{\left(x_{n+1}-{\mu }_{n+1}\right)}^2\right]\hfill \\ \hfill & \times {\displaystyle \frac{\sqrt{{\kappa }_0}}{\sqrt{2\pi }}}{\displaystyle \frac{{\left(\frac{{\nu }_0{\sigma }_0^2}{2}\right)}^{\frac{{\nu }_0}{2}}}{\Gamma \left(\frac{{\nu }_0}{2}\right)}}{\theta }_{n+1}^{-\left({\displaystyle \frac{{\nu }_0+1}{2}}+1\right)}exp\left\{-{\displaystyle \frac{1}{2{\theta }_{n+1}}}\left[{\kappa }_0{\left({\mu }_{n+1}-{\mu }_0\right)}^2+{\nu }_0{\sigma }_0^2\right]\right\}\hfill \\ \hfill & ={\left(2\pi \right)}^{-1}\sqrt{{\kappa }_0}{\displaystyle \frac{{\left(\frac{{\nu }_0{\sigma }_0^2}{2}\right)}^{\frac{{\nu }_0}{2}}}{\Gamma \left(\frac{{\nu }_0}{2}\right)}}{\theta }_{n+1}^{-{\displaystyle \frac{1}{2}}}{\theta }_{n+1}^{-\left({\displaystyle \frac{{\nu }_0+1}{2}}+1\right)}\hfill \\ \hfill & \times exp\left\{-{\displaystyle \frac{1}{2{\theta }_{n+1}}}\left[{\left(x_{n+1}-{\mu }_{n+1}\right)}^2+{\kappa }_0{\left({\mu }_{n+1}-{\mu }_0\right)}^2+{\nu }_0{\sigma }_0^2\right]\right\}\hfill \\ \hfill & ={\left(2\pi \right)}^{-1}\sqrt{{\kappa }_0}{\displaystyle \frac{{\left(\frac{{\nu }_0{\sigma }_0^2}{2}\right)}^{\frac{{\nu }_0}{2}}}{\Gamma \left(\frac{{\nu }_0}{2}\right)}}{\theta }_{n+1}^{-\left({\displaystyle \frac{{\nu }_1+1}{2}}+1\right)}exp\left\{-{\displaystyle \frac{1}{2{\theta }_{n+1}}}\left[{\kappa }_1{\left({\mu }_{n+1}-{\mu }_1\right)}^2+{\nu }_1{\sigma }_1^2\right]\right\}.\hfill \end{array}$$

(A5)

Finally, let us derive the marginal density of $x_{n+1}$. Combining (A5) and (A3), we have

$$\begin{array}{cc}\hfill & m\left(x_{n+1}|{\mu }_0,{\kappa }_0,{\nu }_0,{\sigma }_0^2\right)\hfill \\ \hfill & ={\displaystyle \frac{f\left(x_{n+1}|{\mu }_{n+1},{\theta }_{n+1}\right)\pi \left({\mu }_{n+1},{\theta }_{n+1}\right)}{\pi \left({\mu }_{n+1},{\theta }_{n+1}|x_{n+1}\right)}}\hfill \\ \hfill & ={\displaystyle \frac{{\left(2\pi \right)}^{-1}\sqrt{{\kappa }_0}\frac{{\left(\frac{{\nu }_0{\sigma }_0^2}{2}\right)}^{\frac{{\nu }_0}{2}}}{\Gamma \left(\frac{{\nu }_0}{2}\right)}{\theta }_{n+1}^{-\left(\frac{{\nu }_1+1}{2}+1\right)}exp\left\{-\frac{1}{2{\theta }_{n+1}}\left[{\kappa }_1{\left({\mu }_{n+1}-{\mu }_1\right)}^2+{\nu }_1{\sigma }_1^2\right]\right\}}{{\left(2\pi \right)}^{-\frac{1}{2}}\sqrt{{\kappa }_1}\frac{{\left(\frac{{\nu }_1{\sigma }_1^2}{2}\right)}^{\frac{{\nu }_1}{2}}}{\Gamma \left(\frac{{\nu }_1}{2}\right)}{\theta }_{n+1}^{-\left(\frac{{\nu }_1+1}{2}+1\right)}exp\left\{-\frac{1}{2{\theta }_{n+1}}\left[{\kappa }_1{\left({\mu }_{n+1}-{\mu }_1\right)}^2+{\nu }_1{\sigma }_1^2\right]\right\}}}\hfill \end{array}$$

$$\begin{array}{cc}\hfill & ={\displaystyle \frac{{\left(2\pi \right)}^{-\frac{1}{2}}\sqrt{{\kappa }_0}{\left(\frac{{\nu }_0{\sigma }_0^2}{2}\right)}^{\frac{{\nu }_0}{2}}\Gamma \left(\frac{{\nu }_1}{2}\right)}{\sqrt{{\kappa }_1}{\left(\frac{{\nu }_1{\sigma }_1^2}{2}\right)}^{\frac{{\nu }_1}{2}}\Gamma \left(\frac{{\nu }_0}{2}\right)}}\hfill \\ \hfill & ={\displaystyle \frac{1}{\sqrt{2\pi }}}{\displaystyle \frac{\sqrt{{\kappa }_0}}{\sqrt{{\kappa }_0+1}}}{\displaystyle \frac{{\left(\frac{{\nu }_0{\sigma }_0^2}{2}\right)}^{\frac{{\nu }_0}{2}}\Gamma \left(\frac{{\nu }_0+1}{2}\right)}{\Gamma \left(\frac{{\nu }_0}{2}\right){\left(\frac{{\nu }_0{\sigma }_0^2+\frac{{\kappa }_0}{{\kappa }_0+1}{\left(x_{n+1}-{\mu }_0\right)}^2}{2}\right)}^{\frac{{\nu }_0+1}{2}}}}\hfill \\ \hfill & ={\displaystyle \frac{\Gamma \left(\frac{{\nu }_0+1}{2}\right)}{\Gamma \left(\frac{{\nu }_0}{2}\right)\sqrt{2\pi \left({\kappa }_0+1\right)/{\kappa }_0}}}{\displaystyle \frac{{\left(\frac{{\nu }_0{\sigma }_0^2}{2}\right)}^{\frac{{\nu }_0}{2}}}{{\left(\frac{{\nu }_0{\sigma }_0^2}{2}\right)}^{\frac{{\nu }_0+1}{2}}}}{\left[1+{\displaystyle \frac{1}{{\nu }_0}}{\displaystyle \frac{{\left(x_{n+1}-{\mu }_0\right)}^2}{{\sigma }_0^2\left({\kappa }_0+1\right)/{\kappa }_0}}\right]}^{-{\displaystyle \frac{{\nu }_0+1}{2}}}\hfill \\ \hfill & ={\displaystyle \frac{\Gamma \left(\frac{{\nu }_0+1}{2}\right)}{\Gamma \left(\frac{{\nu }_0}{2}\right)\sqrt{\pi {\nu }_0{\sigma }_0^2\left({\kappa }_0+1\right)/{\kappa }_0}}}{\left[1+{\displaystyle \frac{1}{{\nu }_0}}{\displaystyle \frac{{\left(x_{n+1}-{\mu }_0\right)}^2}{{\sigma }_0^2\left({\kappa }_0+1\right)/{\kappa }_0}}\right]}^{-{\displaystyle \frac{{\nu }_0+1}{2}}}\hfill \\ \hfill & \sim t_{{\nu }_0}\left({\mu }_0,{\displaystyle \frac{{\sigma }_0^2\left({\kappa }_0+1\right)}{{\kappa }_0}}\right).\hfill \end{array}$$

The proof of the theorem is complete. □

Appendix A.2. The Proof of Theorem 2

The hyperparameters of model (2) are $-\infty <{\mu }_0<\infty$, ${\kappa }_0>0$, ${\nu }_0>0$, and ${\sigma }_0>0$. By Theorem 1, we know that the marginal distribution of X in model (2) is a non-standardized student-t distribution, that is,

$$X\sim t_{{\nu }_0}\left({\mu }_0,{\displaystyle \frac{{\sigma }_0^2\left({\kappa }_0+1\right)}{{\kappa }_0}}\right).$$

Since there are four hyperparameters, if we want to obtain the estimators of the hyperparameters of model (2) by the moment method, we need to calculate the first four moments of X at least. By Lemma 4 in [19], we obtain the first six population moments of X as follows. Furthermore, letting the population moments be equal to the sample moments, we obtain

$$\begin{array}{cc}\hfill \mathrm{E}X& ={\mu }_0=A_1,\hfill \\ \hfill \mathrm{E}{X}^2& ={\mu }_0^2+{\displaystyle \frac{{\nu }_0}{{\nu }_0-2}}{\displaystyle \frac{{\sigma }_0^2\left({\kappa }_0+1\right)}{{\kappa }_0}}=A_2,\mathrm{if}{\nu }_0>2,\hfill \\ \hfill \mathrm{E}{X}^3& ={\mu }_0^3+3{\mu }_0{\displaystyle \frac{{\nu }_0}{{\nu }_0-2}}{\displaystyle \frac{{\sigma }_0^2\left({\kappa }_0+1\right)}{{\kappa }_0}}=A_3,\mathrm{if}{\nu }_0>2,\hfill \\ \hfill \mathrm{E}{X}^4& ={\mu }_0^4+6{\mu }_0^2{\displaystyle \frac{{\nu }_0}{{\nu }_0-2}}{\displaystyle \frac{{\sigma }_0^2\left({\kappa }_0+1\right)}{{\kappa }_0}}+3{\displaystyle \frac{{\nu }_0^2}{\left({\nu }_0-2\right)\left({\nu }_0-4\right)}}{\left[{\displaystyle \frac{{\sigma }_0^2\left({\kappa }_0+1\right)}{{\kappa }_0}}\right]}^2=A_4,\mathrm{if}{\nu }_0>4,\hfill \\ \hfill \mathrm{E}{X}^5& ={\mu }_0^5+10{\mu }_0^3{\displaystyle \frac{{\nu }_0}{{\nu }_0-2}}{\displaystyle \frac{{\sigma }_0^2\left({\kappa }_0+1\right)}{{\kappa }_0}}+15{\mu }_0{\displaystyle \frac{{\nu }_0^2}{\left({\nu }_0-2\right)\left({\nu }_0-4\right)}}{\left[{\displaystyle \frac{{\sigma }_0^2\left({\kappa }_0+1\right)}{{\kappa }_0}}\right]}^2=A_5,\mathrm{if}{\nu }_0>4,\hfill \\ \hfill \mathrm{E}{X}^6& =\left\{\begin{array}{c}{\mu }_0^6+15{\mu }_0^4{\displaystyle \frac{{\nu }_0}{{\nu }_0-2}}{\displaystyle \frac{{\sigma }_0^2\left({\kappa }_0+1\right)}{{\kappa }_0}}+45{\mu }_0^2{\displaystyle \frac{{\nu }_0^2}{\left({\nu }_0-2\right)\left({\nu }_0-4\right)}}{\left[{\displaystyle \frac{{\sigma }_0^2\left({\kappa }_0+1\right)}{{\kappa }_0}}\right]}^2\\ +15{\displaystyle \frac{{\nu }_0^3}{\left({\nu }_0-2\right)\left({\nu }_0-4\right)\left({\nu }_0-6\right)}}{\left[{\displaystyle \frac{{\sigma }_0^2\left({\kappa }_0+1\right)}{{\kappa }_0}}\right]}^3\end{array}\right\}=A_6,\mathrm{if}{\nu }_0>6.\hfill \end{array}$$

From the first moment of X, we obtain the moment estimator of ${\mu }_0$ as

$${\mu }_{01}\left(n\right)=A_1.$$

Let

$$\begin{array}{cc}\hfill w_1& ={\displaystyle \frac{{\nu }_0}{{\nu }_0-2}}{\displaystyle \frac{{\sigma }_0^2\left({\kappa }_0+1\right)}{{\kappa }_0}},\hfill \\ \hfill w_2& ={\displaystyle \frac{{\nu }_0^2}{\left({\nu }_0-2\right)\left({\nu }_0-4\right)}}{\left[{\displaystyle \frac{{\sigma }_0^2\left({\kappa }_0+1\right)}{{\kappa }_0}}\right]}^2,\hfill \\ \hfill w_3& ={\displaystyle \frac{{\nu }_0^3}{\left({\nu }_0-2\right)\left({\nu }_0-4\right)\left({\nu }_0-6\right)}}{\left[{\displaystyle \frac{{\sigma }_0^2\left({\kappa }_0+1\right)}{{\kappa }_0}}\right]}^3.\hfill \end{array}$$

From the second moment of X, we obtain the moment estimator of $w_1$ as

$$w_1\left(n\right)=A_2-A_1^2.$$

From the third moment of X, we obtain the moment estimator of $w_1$ as

$$w_1\left(n\right)={\displaystyle \frac{A_3-A_1^3}{3A_1}}.$$

Obviously, the two moment estimators of $w_1$ are not equal. Therefore, we choose one of them as the moment estimator of $w_1$. For simplicity, we use the moment estimator of $w_1$ calculated from $\mathrm{E}{X}^2$, and ignore the third equation involving $\mathrm{E}{X}^3$. Similarly, the equations involving $\mathrm{E}{X}^4$ and $\mathrm{E}{X}^5$ both have $w_1$ and $w_2$. For simplicity, we use the equation involving $\mathrm{E}{X}^4$, and ignore the equation involving $\mathrm{E}{X}^5$. To have four equations, we will use the equation involving $\mathrm{E}{X}^6$. Therefore, the moment equations become

$$\begin{array}{cc}\hfill \mathrm{E}X& ={\mu }_0=A_1,\hfill \\ \hfill \mathrm{E}{X}^2& ={\mu }_0^2+w_1=A_2,\hfill \\ \hfill \mathrm{E}{X}^4& ={\mu }_0^4+6{\mu }_0^2{w}_1+3w_2=A_4,\hfill \\ \hfill \mathrm{E}{X}^6& ={\mu }_0^6+15{\mu }_0^4{w}_1+45{\mu }_0^2{w}_2+15w_3=A_6.\hfill \end{array}$$

Solving the above moment equations, we obtain

$$\begin{array}{cc}\hfill {\mu }_{01}\left(n\right)& =A_1,\hfill \\ \hfill w_1\left(n\right)& ={\displaystyle \frac{{\nu }_0}{{\nu }_0-2}}{\displaystyle \frac{{\sigma }_0^2\left({\kappa }_0+1\right)}{{\kappa }_0}}\hfill \\ \hfill & =A_2-A_1^2\triangleq B_1,\hfill \\ \hfill w_2\left(n\right)& ={\displaystyle \frac{{\nu }_0^2}{\left({\nu }_0-2\right)\left({\nu }_0-4\right)}}{\left[{\displaystyle \frac{{\sigma }_0^2\left({\kappa }_0+1\right)}{{\kappa }_0}}\right]}^2\hfill \\ \hfill & ={\displaystyle \frac{1}{3}}\left(A_4-6A_1^2{A}_2+5A_1^4\right)\triangleq B_2,\hfill \\ \hfill w_3\left(n\right)& ={\displaystyle \frac{{\nu }_0^3}{\left({\nu }_0-2\right)\left({\nu }_0-4\right)\left({\nu }_0-6\right)}}{\left[{\displaystyle \frac{{\sigma }_0^2\left({\kappa }_0+1\right)}{{\kappa }_0}}\right]}^3\hfill \\ \hfill & ={\displaystyle \frac{1}{15}}\left(-61A_1^6+75A_2{A}_1^4-15A_4{A}_1^2+A_6\right)\triangleq B_3.\hfill \end{array}$$

(A6)

Let

$$u_0={\displaystyle \frac{{\sigma }_0^2\left({\kappa }_0+1\right)}{{\kappa }_0}}.$$

Since ${\kappa }_0$ and ${\sigma }_0$ appear together in $u_0$, we can not directly obtain the estimators of ${\kappa }_0$ and ${\sigma }_0$ by the moment method. But we can obtain the estimator of $u_0$ by the moment method. In the following, we are interested in obtaining the moment estimators of $-\infty <{\mu }_0<\infty$, ${\nu }_0>0$, and $u_0>0$. The moment equations involving $w_1\left(n\right)$, $w_2\left(n\right)$, and $w_3\left(n\right)$ become

$$\begin{array}{cc}\hfill {\displaystyle \frac{{\nu }_0}{{\nu }_0-2}}{u}_0& =B_1,\hfill \\ \hfill {\displaystyle \frac{{\nu }_0^2}{\left({\nu }_0-2\right)\left({\nu }_0-4\right)}}{u}_0^2& =B_2,\hfill \\ \hfill {\displaystyle \frac{{\nu }_0^3}{\left({\nu }_0-2\right)\left({\nu }_0-4\right)\left({\nu }_0-6\right)}}{u}_0^3& =B_3.\hfill \end{array}$$

Since there are three equations and only two parameters ${\nu }_0$ and $u_0$, for simplicity, we will only use the first two equations, and ignore the third equation. Solving the above first two equations for ${\nu }_0$ and $u_0$, we obtain the moment estimators of ${\nu }_0$ and $u_0$ as

$$\begin{array}{cc}\hfill {\nu }_{01}\left(n\right)& ={\displaystyle \frac{2B_1^2-4B_2}{B_1^2-B_2}},\hfill \\ \hfill u_{01}\left(n\right)& ={\displaystyle \frac{-B_1{B}_2}{B_1^2-2B_2}}.\hfill \end{array}$$

Substituting the expressions of $B_1$ and $B_2$, and after some algebra, we obtain the expressions of ${\nu }_{01}\left(n\right)$ and $u_{01}\left(n\right)$ in terms of $A_i\left(i=1,2,4\right)$ as

$$\begin{array}{cc}\hfill {\nu }_{01}\left(n\right)& ={\displaystyle \frac{-\frac{14}{3}{A}_1^4+4A_1^2{A}_2+2A_2^2-\frac{4}{3}{A}_4}{-\frac{2}{3}{A}_1^4+A_2^2-\frac{1}{3}{A}_4}},\hfill \end{array}$$

(A7)

$$\begin{array}{cc}\hfill u_{01}\left(n\right)& ={\displaystyle \frac{-\frac{1}{3}\left(A_2-A_1^2\right)\left(5A_1^4-6A_1^2{A}_2+A_4\right)}{-\frac{7}{3}{A}_1^4+2A_1^2{A}_2+A_2^2-\frac{2}{3}{A}_4}}.\hfill \end{array}$$

(A8)

Therefore, the moment estimators of ${\mu }_0$, ${\nu }_0$, and $u_0$ are given by (A6), (A7), and (A8).

Now let us show that the moment estimators are consistent estimators of the hyperparameters. It is easy to show that

$$A_k={\displaystyle \frac{1}{n}}\sum _{i=1}^n{X}_i^k\stackrel{p}{\to }\mathrm{E}{X}^k$$

for $k=1,2,4$, where $\stackrel{p}{\to }$ means convergence in probability. Hence,

$$\begin{array}{cc}\hfill A_1\stackrel{p}{\to }\mathrm{E}X& ={\mu }_0,\hfill \\ \hfill A_2\stackrel{p}{\to }\mathrm{E}{X}^2& ={\mu }_0^2+{\displaystyle \frac{{\nu }_0}{{\nu }_0-2}}{u}_0,\hfill \\ \hfill A_4\stackrel{p}{\to }\mathrm{E}{X}^4& ={\mu }_0^4+6{\mu }_0^2{\displaystyle \frac{{\nu }_0}{{\nu }_0-2}}{u}_0+3{\displaystyle \frac{{\nu }_0^2}{\left({\nu }_0-2\right)\left({\nu }_0-4\right)}}{u}_0^2.\hfill \end{array}$$

Therefore,

$${\mu }_{01}\left(n\right)=A_1\triangleq f\left(A_1\right)\stackrel{p}{\to }f\left(\mathrm{E}X\right)={\mu }_0,$$

$$\begin{array}{cc}\hfill & {\nu }_{01}\left(n\right)\hfill \\ \hfill & ={\displaystyle \frac{-\frac{14}{3}{A}_1^4+4A_1^2{A}_2+2A_2^2-\frac{4}{3}{A}_4}{-\frac{2}{3}{A}_1^4+A_2^2-\frac{1}{3}{A}_4}}\hfill \\ \hfill & \triangleq g\left(A_1,A_2,A_4\right)\hfill \\ \hfill & \stackrel{p}{\to }g\left(\mathrm{E}X,\mathrm{E}{X}^2,\mathrm{E}{X}^4\right)\hfill \\ \hfill & ={\displaystyle \frac{\begin{array}{c}-\frac{14}{3}{\mu }_0^4+4{\mu }_0^2\left({\mu }_0^2+\frac{{\nu }_0}{{\nu }_0-2}{u}_0\right)+2{\left({\mu }_0^2+\frac{{\nu }_0}{{\nu }_0-2}{u}_0\right)}^2\\ -\frac{4}{3}\left({\mu }_0^4+6{\mu }_0^2\frac{{\nu }_0}{{\nu }_0-2}{u}_0+3\frac{{\nu }_0^2}{\left({\nu }_0-2\right)\left({\nu }_0-4\right)}{u}_0^2\right)\end{array}}{-\frac{2}{3}{\mu }_0^4+{\left({\mu }_0^2+\frac{{\nu }_0}{{\nu }_0-2}{u}_0\right)}^2-\frac{1}{3}\left({\mu }_0^4+6{\mu }_0^2\frac{{\nu }_0}{{\nu }_0-2}{u}_0+3\frac{{\nu }_0^2}{\left({\nu }_0-2\right)\left({\nu }_0-4\right)}{u}_0^2\right)}}\hfill \\ \hfill & ={\nu }_0,\hfill \end{array}$$

and

$$\begin{array}{cc}\hfill & u_{01}\left(n\right)\hfill \\ \hfill & ={\displaystyle \frac{-\frac{1}{3}\left(A_2-A_1^2\right)\left(5A_1^4-6A_1^2{A}_2+A_4\right)}{-\frac{7}{3}{A}_1^4+2A_1^2{A}_2+A_2^2-\frac{2}{3}{A}_4}}\hfill \\ \hfill & \triangleq h\left(A_1,A_2,A_4\right)\hfill \\ \hfill & \stackrel{p}{\to }h\left(\mathrm{E}X,\mathrm{E}{X}^2,\mathrm{E}{X}^4\right)\hfill \\ \hfill & ={\displaystyle \frac{-\frac{1}{3}\left(\left({\mu }_0^2+\frac{{\nu }_0}{{\nu }_0-2}{u}_0\right)-{\mu }_0^2\right)\left(\begin{array}{c}5{\mu }_0^4-6{\mu }_0^2\left({\mu }_0^2+\frac{{\nu }_0}{{\nu }_0-2}{u}_0\right)\\ +\left({\mu }_0^4+6{\mu }_0^2\frac{{\nu }_0}{{\nu }_0-2}{u}_0+3\frac{{\nu }_0^2}{\left({\nu }_0-2\right)\left({\nu }_0-4\right)}{u}_0^2\right)\end{array}\right)}{\begin{array}{c}-\frac{7}{3}{\mu }_0^4+2{\mu }_0^2\left({\mu }_0^2+\frac{{\nu }_0}{{\nu }_0-2}{u}_0\right)+{\left({\mu }_0^2+\frac{{\nu }_0}{{\nu }_0-2}{u}_0\right)}^2\\ -\frac{2}{3}\left({\mu }_0^4+6{\mu }_0^2\frac{{\nu }_0}{{\nu }_0-2}{u}_0+3\frac{{\nu }_0^2}{\left({\nu }_0-2\right)\left({\nu }_0-4\right)}{u}_0^2\right)\end{array}}}\hfill \\ \hfill & =u_0.\hfill \end{array}$$

The proof of the theorem is complete. □

Appendix A.3. The Proof of Theorem 3

Now we derive the Maximum Likelihood Estimators (MLEs) of ${\mu }_0$, ${\nu }_0$, and $u_0$. By Theorem 1, we know that the marginal distribution of X of the model (2) is

$$m\left(x|{\mu }_0,{\nu }_0,u_0\right)={\displaystyle \frac{\Gamma \left(\frac{{\nu }_0+1}{2}\right)}{\Gamma \left(\frac{{\nu }_0}{2}\right)\sqrt{\pi {\nu }_0{u}_0}}}{\left[1+{\displaystyle \frac{1}{{\nu }_0}}{\displaystyle \frac{{\left(x-{\mu }_0\right)}^2}{u_0}}\right]}^{-{\displaystyle \frac{{\nu }_0+1}{2}}},$$

for $-\infty <x<\infty$,$-\infty <{\mu }_0<\infty$, ${\nu }_0>0$, and $u_0>0$. Then the likelihood function of ${\mu }_0$, ${\nu }_0$, and $u_0$ is

$$\begin{array}{cc}\hfill & L\left({\mu }_0,{\nu }_0,u_0|\mathit{x}\right)\hfill \\ \hfill & =m\left(\mathit{x}|{\mu }_0,{\nu }_0,u_0\right)\hfill \\ \hfill & =\prod _{i=1}^{n}m\left(x_i|{\mu }_0,{\nu }_0,u_0\right)\hfill \\ \hfill & =\prod _{i=1}^n{\displaystyle \frac{\Gamma \left(\frac{{\nu }_0+1}{2}\right)}{\Gamma \left(\frac{{\nu }_0}{2}\right)\sqrt{\pi {\nu }_0{u}_0}}}{\left[1+{\displaystyle \frac{1}{{\nu }_0}}{\displaystyle \frac{{\left(x_i-{\mu }_0\right)}^2}{u_0}}\right]}^{-{\displaystyle \frac{{\nu }_0+1}{2}}}\hfill \\ \hfill & ={\displaystyle \frac{{\Gamma }^n\left(\frac{{\nu }_0+1}{2}\right)}{{\Gamma }^n\left(\frac{{\nu }_0}{2}\right)}}{\left(\pi {\nu }_0{u}_0\right)}^{-{\displaystyle \frac{n}{2}}}\prod _{i=1}^n{\left[1+{\displaystyle \frac{{\left(x_i-{\mu }_0\right)}^2}{{\nu }_0{u}_0}}\right]}^{-{\displaystyle \frac{{\nu }_0+1}{2}}}.\hfill \end{array}$$

Consequently, the log-likelihood function of ${\mu }_0$, ${\nu }_0$, and $u_0$ is

$$\begin{array}{cc}\hfill & logL\left({\mu }_0,{\nu }_0,u_0|\mathit{x}\right)\hfill \\ \hfill & =nlog\Gamma \left({\displaystyle \frac{{\nu }_0+1}{2}}\right)-nlog\Gamma \left({\displaystyle \frac{{\nu }_0}{2}}\right)-{\displaystyle \frac{n}{2}}log\left(\pi {\nu }_0{u}_0\right)\hfill \\ \hfill & +\sum _{i=1}^n\left(-{\displaystyle \frac{{\nu }_0+1}{2}}\right)log\left[1+{\displaystyle \frac{{\left(x_i-{\mu }_0\right)}^2}{{\nu }_0{u}_0}}\right]\hfill \\ \hfill & =nlog\Gamma \left({\displaystyle \frac{{\nu }_0+1}{2}}\right)-nlog\Gamma \left({\displaystyle \frac{{\nu }_0}{2}}\right)-{\displaystyle \frac{n}{2}}\left(log\pi +log{\nu }_0+log{u}_0\right)\hfill \\ \hfill & -\sum _{i=1}^n{\displaystyle \frac{{\nu }_0+1}{2}}log\left[1+{\displaystyle \frac{{\left(x_i-{\mu }_0\right)}^2}{{\nu }_0{u}_0}}\right].\hfill \end{array}$$

Taking partial derivatives with respect to ${\mu }_0$, ${\nu }_0$, and $u_0$ and setting them to zeros, we obtain

$$\begin{array}{cc}\hfill & {\displaystyle \frac{\partial }{\partial {\mu }_0}}logL\hfill \\ \hfill & =-\sum _{i=1}^n{\displaystyle \frac{{\nu }_0+1}{2}}{\displaystyle \frac{1}{1+\frac{{\left(x_i-{\mu }_0\right)}^2}{{\nu }_0{u}_0}}}{\displaystyle \frac{1}{{\nu }_0{u}_0}}2\left(x_i-{\mu }_0\right)\left(-1\right)\hfill \\ \hfill & =\sum _{i=1}^n{\displaystyle \frac{\left({\nu }_0+1\right)\left(x_i-{\mu }_0\right)}{{\nu }_0{u}_0+{\left(x_i-{\mu }_0\right)}^2}}\hfill \\ \hfill & =0,\hfill \end{array}$$

$$\begin{array}{cc}\hfill & {\displaystyle \frac{\partial }{\partial {\nu }_0}}logL\hfill \\ \hfill & =n\psi \left({\displaystyle \frac{{\nu }_0+1}{2}}\right){\displaystyle \frac{1}{2}}-n\psi \left({\displaystyle \frac{{\nu }_0}{2}}\right){\displaystyle \frac{1}{2}}-{\displaystyle \frac{n}{2}}{\displaystyle \frac{1}{{\nu }_0}}\hfill \\ \hfill & -{\displaystyle \frac{1}{2}}\sum _{i=1}^n\left\{1·log\left[1+{\displaystyle \frac{{\left(x_i-{\mu }_0\right)}^2}{{\nu }_0{u}_0}}\right]+\left({\nu }_0+1\right){\displaystyle \frac{1}{1+\frac{{\left(x_i-{\mu }_0\right)}^2}{{\nu }_0{u}_0}}}{\displaystyle \frac{{\left(x_i-{\mu }_0\right)}^2}{u_0}}\left(-{\displaystyle \frac{1}{{\nu }_0^2}}\right)\right\}\hfill \\ \hfill & ={\displaystyle \frac{n}{2}}\psi \left({\displaystyle \frac{{\nu }_0+1}{2}}\right)-{\displaystyle \frac{n}{2}}\psi \left({\displaystyle \frac{{\nu }_0}{2}}\right)-{\displaystyle \frac{n}{2}}{\displaystyle \frac{1}{{\nu }_0}}\hfill \\ \hfill & -{\displaystyle \frac{1}{2}}\sum _{i=1}^n\left\{log\left[1+{\displaystyle \frac{{\left(x_i-{\mu }_0\right)}^2}{{\nu }_0{u}_0}}\right]-{\displaystyle \frac{\left({\nu }_0+1\right){\left(x_i-{\mu }_0\right)}^2}{\left[{\nu }_0{u}_0+{\left(x_i-{\mu }_0\right)}^2\right]{\nu }_0}}\right\}\hfill \\ \hfill & =0,\hfill \end{array}$$

and

$$\begin{array}{cc}\hfill & {\displaystyle \frac{\partial }{\partial u_0}}logL\hfill \\ \hfill & =-{\displaystyle \frac{n}{2}}{\displaystyle \frac{1}{u_0}}-\sum _{i=1}^n{\displaystyle \frac{{\nu }_0+1}{2}}{\displaystyle \frac{1}{1+\frac{{\left(x_i-{\mu }_0\right)}^2}{{\nu }_0{u}_0}}}{\displaystyle \frac{{\left(x_i-{\mu }_0\right)}^2}{{\nu }_0}}\left(-{\displaystyle \frac{1}{u_0^2}}\right)\hfill \\ \hfill & =-{\displaystyle \frac{n}{2}}{\displaystyle \frac{1}{u_0}}+{\displaystyle \frac{1}{2}}\sum _{i=1}^n{\displaystyle \frac{\left({\nu }_0+1\right){\left(x_i-{\mu }_0\right)}^2}{\left[{\nu }_0{u}_0+{\left(x_i-{\mu }_0\right)}^2\right]u_0}}\hfill \\ \hfill & =0,\hfill \end{array}$$

where

$$\psi \left(x\right)={\displaystyle \frac{{\Gamma }^{\prime }\left(x\right)}{\Gamma \left(x\right)}}={\displaystyle \frac{d}{dx}}log\Gamma \left(x\right)=\mathrm{digamma}\left(x\right)$$

which can be directly calculated in R software by digamma(x) ([46]). Let

$$g_i={\nu }_0{u}_0+{\left(x_i-{\mu }_0\right)}^2.$$

Thus,

$${\left(x_i-{\mu }_0\right)}^2=g_i-{\nu }_0{u}_0.$$

After some algebra, the above equations reduce to

$$\begin{array}{cc}\hfill f_1\left({\mu }_0,{\nu }_0,u_0\right)& ={\displaystyle \frac{1}{n}}\sum _{i=1}^n{\displaystyle \frac{\left({\nu }_0+1\right)\left(x_i-{\mu }_0\right)}{{\nu }_0{u}_0+{\left(x_i-{\mu }_0\right)}^2}}\hfill \\ \hfill & ={\displaystyle \frac{1}{n}}\sum _{i=1}^n{\displaystyle \frac{\left({\nu }_0+1\right)\left(x_i-{\mu }_0\right)}{g_i}}\hfill \\ \hfill & =0,\hfill \end{array}$$

(A9)

$$\begin{array}{cc}\hfill & f_2\left({\mu }_0,{\nu }_0,u_0\right)\hfill \\ \hfill & ={\displaystyle \frac{1}{2}}\psi \left({\displaystyle \frac{{\nu }_0+1}{2}}\right)-{\displaystyle \frac{1}{2}}\psi \left({\displaystyle \frac{{\nu }_0}{2}}\right)-{\displaystyle \frac{1}{2{\nu }_0}}\hfill \\ \hfill & -{\displaystyle \frac{1}{2}}{\displaystyle \frac{1}{n}}\sum _{i=1}^n\left\{log\left[1+{\displaystyle \frac{{\left(x_i-{\mu }_0\right)}^2}{{\nu }_0{u}_0}}\right]-{\displaystyle \frac{\left({\nu }_0+1\right){\left(x_i-{\mu }_0\right)}^2}{\left[{\nu }_0{u}_0+{\left(x_i-{\mu }_0\right)}^2\right]{\nu }_0}}\right\}\hfill \\ \hfill & ={\displaystyle \frac{1}{2}}\psi \left({\displaystyle \frac{{\nu }_0+1}{2}}\right)-{\displaystyle \frac{1}{2}}\psi \left({\displaystyle \frac{{\nu }_0}{2}}\right)-{\displaystyle \frac{1}{2{\nu }_0}}\hfill \\ \hfill & -{\displaystyle \frac{1}{2}}{\displaystyle \frac{1}{n}}\sum _{i=1}^n\left\{log{g}_i-log\left({\nu }_0{u}_0\right)-{\displaystyle \frac{\left({\nu }_0+1\right)\left(g_i-{\nu }_0{u}_0\right)}{g_i{\nu }_0}}\right\}\hfill \\ \hfill & =0,\hfill \end{array}$$

(A10)

$$\begin{array}{cc}\hfill f_3\left({\mu }_0,{\nu }_0,u_0\right)& =-{\displaystyle \frac{1}{2u_0}}+{\displaystyle \frac{1}{2}}{\displaystyle \frac{1}{n}}\sum _{i=1}^n{\displaystyle \frac{\left({\nu }_0+1\right){\left(x_i-{\mu }_0\right)}^2}{\left[{\nu }_0{u}_0+{\left(x_i-{\mu }_0\right)}^2\right]u_0}}\hfill \\ \hfill & =-{\displaystyle \frac{1}{2u_0}}+{\displaystyle \frac{1}{2}}{\displaystyle \frac{1}{n}}\sum _{i=1}^n{\displaystyle \frac{\left({\nu }_0+1\right)\left(g_i-{\nu }_0{u}_0\right)}{g_i{u}_0}}\hfill \\ \hfill & =0.\hfill \end{array}$$

(A11)

In (A9), (A10), and (A11), the expressions involving $g_i$ are used for simplifying the R coding.

The MLEs of the hyperparameters ${\mu }_0$, ${\nu }_0$, and $u_0$ are the solutions to the equations (A9), (A10), and (A11). The analytical calculations of the MLEs of ${\mu }_0$, ${\nu }_0$, and $u_0$ by solving the above equations are impossible, and thus we have to resort to numerical solutions. We can exploit Newton’s method to solve the above equations and to obtain the MLEs of ${\mu }_0$, ${\nu }_0$, and $u_0$. The iterative scheme of Newton’s method is

$${\mathit{p}}^{\left(k+1\right)}={\mathit{p}}^{\left(k\right)}-{\left[\mathit{J}\left({\mathit{p}}^{\left(k\right)}\right)\right]}^{-1}\mathit{f}\left({\mathit{p}}^{\left(k\right)}\right),k=0,1,\dots ,$$

where $\mathit{J}\left(\mathit{p}\right)$ is the Jacobian matrix of $\mathit{f}\left(\mathit{p}\right)={\left(f_1\left(\mathit{p}\right),f_2\left(\mathit{p}\right),f_3\left(\mathit{p}\right)\right)}^{\prime }$ and $\mathit{p}={\left({\mu }_0,{\nu }_0,u_0\right)}^{\prime }$. Note that the MLEs of ${\mu }_0$, ${\nu }_0$, and $u_0$ are very sensitive to the initial estimators, and the moment estimators are usually proved to be good initial estimators. The Jacobian matrix of ${\mu }_0$, ${\nu }_0$, and $u_0$ is given by

$$J=\left(\begin{array}{ccc}{J}_{11}& J_{12}& J_{13}\\ J_{21}& J_{22}& J_{23}\\ J_{31}& J_{32}& J_{33}\end{array}\right)=\left(\begin{array}{ccc}{\displaystyle \frac{\partial f_1}{\partial {\mu }_0}}& {\displaystyle \frac{\partial f_1}{\partial {\nu }_0}}& {\displaystyle \frac{\partial f_1}{\partial u_0}}\\ {\displaystyle \frac{\partial f_2}{\partial {\mu }_0}}& {\displaystyle \frac{\partial f_2}{\partial {\nu }_0}}& {\displaystyle \frac{\partial f_2}{\partial u_0}}\\ {\displaystyle \frac{\partial f_3}{\partial {\mu }_0}}& {\displaystyle \frac{\partial f_3}{\partial {\nu }_0}}& {\displaystyle \frac{\partial f_3}{\partial u_0}}\end{array}\right),$$

where

$$\begin{array}{cc}\hfill J_{11}& ={\displaystyle \frac{\partial f_1}{\partial {\mu }_0}}\hfill \\ \hfill & ={\displaystyle \frac{1}{n}}\sum _{i=1}^n\left({\nu }_0+1\right){\displaystyle \frac{\left(-1\right)\left[{\nu }_0{u}_0+{\left(x_i-{\mu }_0\right)}^2\right]-\left(x_i-{\mu }_0\right)2\left(x_i-{\mu }_0\right)\left(-1\right)}{{\left[{\nu }_0{u}_0+{\left(x_i-{\mu }_0\right)}^2\right]}^2}}\hfill \\ \hfill & ={\displaystyle \frac{1}{n}}\sum _{i=1}^n\left({\nu }_0+1\right){\displaystyle \frac{2{\left(x_i-{\mu }_0\right)}^2-{\nu }_0{u}_0-{\left(x_i-{\mu }_0\right)}^2}{{\left[{\nu }_0{u}_0+{\left(x_i-{\mu }_0\right)}^2\right]}^2}}\hfill \\ \hfill & ={\displaystyle \frac{1}{n}}\sum _{i=1}^n\left({\nu }_0+1\right){\displaystyle \frac{{\left(x_i-{\mu }_0\right)}^2-{\nu }_0{u}_0}{{\left[{\nu }_0{u}_0+{\left(x_i-{\mu }_0\right)}^2\right]}^2}}\hfill \\ \hfill & ={\displaystyle \frac{1}{n}}\sum _{i=1}^n\left({\nu }_0+1\right){\displaystyle \frac{g_i-2{\nu }_0{u}_0}{g_i^2}},\hfill \end{array}$$

$$\begin{array}{cc}\hfill J_{12}& ={\displaystyle \frac{\partial f_1}{\partial {\nu }_0}}\hfill \\ \hfill & ={\displaystyle \frac{1}{n}}\sum _{i=1}^n\left(x_i-{\mu }_0\right){\displaystyle \frac{1·\left[{\nu }_0{u}_0+{\left(x_i-{\mu }_0\right)}^2\right]-\left({\nu }_0+1\right)u_0}{{\left[{\nu }_0{u}_0+{\left(x_i-{\mu }_0\right)}^2\right]}^2}}\hfill \\ \hfill & ={\displaystyle \frac{1}{n}}\sum _{i=1}^n{\displaystyle \frac{\left(x_i-{\mu }_0\right)\left[{\left(x_i-{\mu }_0\right)}^2-u_0\right]}{{\left[{\nu }_0{u}_0+{\left(x_i-{\mu }_0\right)}^2\right]}^2}}\hfill \\ \hfill & ={\displaystyle \frac{1}{n}}\sum _{i=1}^n{\displaystyle \frac{\left(x_i-{\mu }_0\right)\left(g_i-{\nu }_0{u}_0-u_0\right)}{g_i^2}},\hfill \end{array}$$

$$\begin{array}{cc}\hfill J_{13}& ={\displaystyle \frac{\partial f_1}{\partial u_0}}\hfill \\ \hfill & ={\displaystyle \frac{1}{n}}\sum _{i=1}^n\left({\nu }_0+1\right)\left(x_i-{\mu }_0\right)\left(-1\right){\displaystyle \frac{1}{{\left[{\nu }_0{u}_0+{\left(x_i-{\mu }_0\right)}^2\right]}^2}}{\nu }_0\hfill \\ \hfill & =-{\displaystyle \frac{1}{n}}\sum _{i=1}^n{\displaystyle \frac{{\nu }_0\left({\nu }_0+1\right)\left(x_i-{\mu }_0\right)}{{\left[{\nu }_0{u}_0+{\left(x_i-{\mu }_0\right)}^2\right]}^2}}\hfill \\ \hfill & =-{\displaystyle \frac{1}{n}}\sum _{i=1}^n{\displaystyle \frac{{\nu }_0\left({\nu }_0+1\right)\left(x_i-{\mu }_0\right)}{g_i^2}},\hfill \end{array}$$

$$\begin{array}{cc}\hfill J_{21}& ={\displaystyle \frac{\partial f_2}{\partial {\mu }_0}}\hfill \\ \hfill & =-{\displaystyle \frac{1}{2}}{\displaystyle \frac{1}{n}}\sum _{i=1}^n\left\{\begin{array}{c}{\displaystyle \frac{1}{1+\frac{{\left(x_i-{\mu }_0\right)}^2}{{\nu }_0{u}_0}}}{\displaystyle \frac{1}{{\nu }_0{u}_0}}2\left(x_i-{\mu }_0\right)\left(-1\right)\\ -{\displaystyle \frac{{\nu }_0+1}{{\nu }_0}}{\displaystyle \frac{2\left(x_i-{\mu }_0\right)\left(-1\right)g_i-{\left(x_i-{\mu }_0\right)}^{2}2\left(x_i-{\mu }_0\right)\left(-1\right)}{g_i^2}}\end{array}\right\}\hfill \\ \hfill & =-{\displaystyle \frac{1}{2}}{\displaystyle \frac{1}{n}}\sum _{i=1}^n\left\{{\displaystyle \frac{-2\left(x_i-{\mu }_0\right)}{g_i}}+{\displaystyle \frac{{\nu }_0+1}{{\nu }_0}}{\displaystyle \frac{2\left(x_i-{\mu }_0\right)g_i-2{\left(x_i-{\mu }_0\right)}^3}{g_i^2}}\right\}\hfill \\ \hfill & =-{\displaystyle \frac{1}{2}}{\displaystyle \frac{1}{n}}\sum _{i=1}^n\left\{{\displaystyle \frac{-2\left(x_i-{\mu }_0\right)}{g_i}}+{\displaystyle \frac{{\nu }_0+1}{{\nu }_0}}{\displaystyle \frac{2\left(x_i-{\mu }_0\right){\nu }_0{u}_0}{g_i^2}}\right\}\hfill \end{array}$$

$$\begin{array}{cc}\hfill & ={\displaystyle \frac{1}{2}}{\displaystyle \frac{1}{n}}\sum _{i=1}^n\left\{{\displaystyle \frac{2\left(x_i-{\mu }_0\right)}{g_i}}-{\displaystyle \frac{2\left({\nu }_0+1\right)u_0\left(x_i-{\mu }_0\right)}{g_i^2}}\right\}\hfill \\ \hfill & ={\displaystyle \frac{1}{2}}{\displaystyle \frac{1}{n}}\sum _{i=1}^n{\displaystyle \frac{2\left(x_i-{\mu }_0\right)g_i-2\left(x_i-{\mu }_0\right)\left({\nu }_0{u}_0+u_0\right)}{g_i^2}}\hfill \\ \hfill & ={\displaystyle \frac{1}{2}}{\displaystyle \frac{1}{n}}\sum _{i=1}^n{\displaystyle \frac{2\left(x_i-{\mu }_0\right)\left[g_i-{\nu }_0{u}_0-u_0\right]}{g_i^2}}\hfill \\ \hfill & ={\displaystyle \frac{1}{n}}\sum _{i=1}^n{\displaystyle \frac{\left(x_i-{\mu }_0\right)\left(g_i-{\nu }_0{u}_0-u_0\right)}{g_i^2}}\hfill \\ \hfill & =J_{12},\hfill \end{array}$$

$$\begin{array}{cc}\hfill J_{22}& ={\displaystyle \frac{\partial f_2}{\partial {\nu }_0}}\hfill \\ \hfill & ={\displaystyle \frac{1}{2}}{\psi }^{\prime }\left({\displaystyle \frac{{\nu }_0+1}{2}}\right){\displaystyle \frac{1}{2}}-{\displaystyle \frac{1}{2}}{\psi }^{\prime }\left({\displaystyle \frac{{\nu }_0}{2}}\right){\displaystyle \frac{1}{2}}-{\displaystyle \frac{1}{2}}\left(-{\displaystyle \frac{1}{{\nu }_0^2}}\right)\hfill \\ \hfill & -{\displaystyle \frac{1}{2}}{\displaystyle \frac{1}{n}}\sum _{i=1}^n\left\{{\displaystyle \frac{1}{1+\frac{{\left(x_i-{\mu }_0\right)}^2}{{\nu }_0{u}_0}}}{\displaystyle \frac{{\left(x_i-{\mu }_0\right)}^2}{u_0}}\left(-{\displaystyle \frac{1}{{\nu }_0^2}}\right)-{\left(x_i-{\mu }_0\right)}^2{\displaystyle \frac{-\frac{1}{{\nu }_0^2}{g}_i-\left(1+\frac{1}{{\nu }_0}\right)u_0}{g_i^2}}\right\}\hfill \\ \hfill & ={\displaystyle \frac{1}{4}}{\psi }^{\prime }\left({\displaystyle \frac{{\nu }_0+1}{2}}\right)-{\displaystyle \frac{1}{4}}{\psi }^{\prime }\left({\displaystyle \frac{{\nu }_0}{2}}\right)+{\displaystyle \frac{1}{2{\nu }_0^2}}\hfill \\ \hfill & +{\displaystyle \frac{1}{2}}{\displaystyle \frac{1}{n}}\sum _{i=1}^n\left\{{\displaystyle \frac{{\left(x_i-{\mu }_0\right)}^2}{g_i{\nu }_0}}-{\left(x_i-{\mu }_0\right)}^2{\displaystyle \frac{\frac{1}{{\nu }_0^2}{g}_i+\left(1+\frac{1}{{\nu }_0}\right)u_0}{g_i^2}}\right\}\hfill \\ \hfill & ={\displaystyle \frac{1}{4}}{\psi }^{\prime }\left({\displaystyle \frac{{\nu }_0+1}{2}}\right)-{\displaystyle \frac{1}{4}}{\psi }^{\prime }\left({\displaystyle \frac{{\nu }_0}{2}}\right)+{\displaystyle \frac{1}{2{\nu }_0^2}}\hfill \\ \hfill & +{\displaystyle \frac{1}{2}}{\displaystyle \frac{1}{n}}\sum _{i=1}^n{\displaystyle \frac{\left(g_i-{\nu }_0{u}_0\right)\left[g_i\left({\nu }_0-1\right)-u_0{\nu }_0\left({\nu }_0+1\right)\right]}{g_i^2{\nu }_0^2}},\hfill \end{array}$$

where

$$\begin{array}{cc}\hfill & {\displaystyle \frac{{\left(x_i-{\mu }_0\right)}^2}{g_i{\nu }_0}}-{\left(x_i-{\mu }_0\right)}^2{\displaystyle \frac{\frac{1}{{\nu }_0^2}{g}_i+\left(1+\frac{1}{{\nu }_0}\right)u_0}{g_i^2}}\hfill \\ \hfill & ={\displaystyle \frac{g_i-{\nu }_0{u}_0}{g_i{\nu }_0}}-{\displaystyle \frac{\left(g_i-{\nu }_0{u}_0\right)\left[\frac{1}{{\nu }_0^2}{g}_i+\left(1+\frac{1}{{\nu }_0}\right)u_0\right]}{g_i^2}}\hfill \\ \hfill & ={\displaystyle \frac{g_i-{\nu }_0{u}_0}{g_i{\nu }_0}}-{\displaystyle \frac{\left(g_i-{\nu }_0{u}_0\right)\left[g_i+u_0{\nu }_0\left({\nu }_0+1\right)\right]}{g_i^2{\nu }_0^2}}\hfill \\ \hfill & ={\displaystyle \frac{\left(g_i-{\nu }_0{u}_0\right)g_i{\nu }_0-\left(g_i-{\nu }_0{u}_0\right)\left[g_i+u_0{\nu }_0\left({\nu }_0+1\right)\right]}{g_i^2{\nu }_0^2}}\hfill \\ \hfill & ={\displaystyle \frac{\left(g_i-{\nu }_0{u}_0\right)\left[g_i{\nu }_0-g_i-u_0{\nu }_0\left({\nu }_0+1\right)\right]}{g_i^2{\nu }_0^2}}\hfill \\ \hfill & ={\displaystyle \frac{\left(g_i-{\nu }_0{u}_0\right)\left[g_i\left({\nu }_0-1\right)-u_0{\nu }_0\left({\nu }_0+1\right)\right]}{g_i^2{\nu }_0^2}},\hfill \end{array}$$

$$\begin{array}{cc}\hfill J_{23}& ={\displaystyle \frac{\partial f_2}{\partial u_0}}\hfill \\ \hfill & =-{\displaystyle \frac{1}{2}}{\displaystyle \frac{1}{n}}\sum _{i=1}^n\left\{\begin{array}{c}{\displaystyle \frac{1}{1+\frac{{\left(x_i-{\mu }_0\right)}^2}{{\nu }_0{u}_0}}}{\displaystyle \frac{{\left(x_i-{\mu }_0\right)}^2}{{\nu }_0}}\left(-{\displaystyle \frac{1}{u_0^2}}\right)\\ -{\displaystyle \frac{\left({\nu }_0+1\right){\left(x_i-{\mu }_0\right)}^2}{{\nu }_0}}\left(-{\displaystyle \frac{1}{g_i^2}}\right){\nu }_0\end{array}\right\}\hfill \\ \hfill & ={\displaystyle \frac{1}{2}}{\displaystyle \frac{1}{n}}\sum _{i=1}^n\left\{{\displaystyle \frac{{\left(x_i-{\mu }_0\right)}^2}{g_i{u}_0}}-{\displaystyle \frac{\left({\nu }_0+1\right){\left(x_i-{\mu }_0\right)}^2}{g_i^2}}\right\}\hfill \\ \hfill & ={\displaystyle \frac{1}{2}}{\displaystyle \frac{1}{n}}\sum _{i=1}^n{\displaystyle \frac{{\left(x_i-{\mu }_0\right)}^2\left[{\nu }_0{u}_0+{\left(x_i-{\mu }_0\right)}^2\right]-u_0\left({\nu }_0+1\right){\left(x_i-{\mu }_0\right)}^2}{g_i^2{u}_0}}\hfill \end{array}$$

$$\begin{array}{cc}\hfill & ={\displaystyle \frac{1}{2}}{\displaystyle \frac{1}{n}}\sum _{i=1}^n{\displaystyle \frac{{\nu }_0{u}_0{\left(x_i-{\mu }_0\right)}^2+{\left(x_i-{\mu }_0\right)}^4-{\nu }_0{u}_0{\left(x_i-{\mu }_0\right)}^2-u_0{\left(x_i-{\mu }_0\right)}^2}{g_i^2{u}_0}}\hfill \\ \hfill & ={\displaystyle \frac{1}{2}}{\displaystyle \frac{1}{n}}\sum _{i=1}^n{\displaystyle \frac{{\left(x_i-{\mu }_0\right)}^4-u_0{\left(x_i-{\mu }_0\right)}^2}{g_i^2{u}_0}}\hfill \\ \hfill & ={\displaystyle \frac{1}{2}}{\displaystyle \frac{1}{n}}\sum _{i=1}^n{\displaystyle \frac{{\left(x_i-{\mu }_0\right)}^2\left[{\left(x_i-{\mu }_0\right)}^2-u_0\right]}{g_i^2{u}_0}}\hfill \\ \hfill & ={\displaystyle \frac{1}{2}}{\displaystyle \frac{1}{n}}\sum _{i=1}^n{\displaystyle \frac{\left(g_i-{\nu }_0{u}_0\right)\left(g_i-{\nu }_0{u}_0-u_0\right)}{g_i^2{u}_0}},\hfill \end{array}$$

$$\begin{array}{cc}\hfill J_{31}& ={\displaystyle \frac{\partial f_3}{\partial {\mu }_0}}\hfill \\ \hfill & ={\displaystyle \frac{1}{2}}{\displaystyle \frac{1}{n}}\sum _{i=1}^n{\displaystyle \frac{{\nu }_0+1}{u_0}}{\displaystyle \frac{2\left(x_i-{\mu }_0\right)\left(-1\right)g_i-{\left(x_i-{\mu }_0\right)}^{2}2\left(x_i-{\mu }_0\right)\left(-1\right)}{g_i^2}}\hfill \\ \hfill & ={\displaystyle \frac{1}{2}}{\displaystyle \frac{1}{n}}\sum _{i=1}^n{\displaystyle \frac{{\nu }_0+1}{u_0}}{\displaystyle \frac{2\left(x_i-{\mu }_0\right)\left[{\left(x_i-{\mu }_0\right)}^2-g_i\right]}{g_i^2}}\hfill \\ \hfill & ={\displaystyle \frac{1}{2}}{\displaystyle \frac{1}{n}}\sum _{i=1}^n{\displaystyle \frac{{\nu }_0+1}{u_0}}{\displaystyle \frac{2\left(x_i-{\mu }_0\right)\left(-{\nu }_0{u}_0\right)}{g_i^2}}\hfill \\ \hfill & =-{\displaystyle \frac{1}{n}}\sum _{i=1}^n{\displaystyle \frac{{\nu }_0\left({\nu }_0+1\right)\left(x_i-{\mu }_0\right)}{g_i^2}}\hfill \\ \hfill & =J_{13},\hfill \end{array}$$

$$\begin{array}{cc}\hfill J_{32}& ={\displaystyle \frac{\partial f_3}{\partial {\nu }_0}}\hfill \\ \hfill & ={\displaystyle \frac{1}{2}}{\displaystyle \frac{1}{n}}\sum _{i=1}^n{\displaystyle \frac{{\left(x_i-{\mu }_0\right)}^2}{u_0}}{\displaystyle \frac{1·g_i-\left({\nu }_0+1\right)u_0}{g_i^2}}\hfill \\ \hfill & ={\displaystyle \frac{1}{2}}{\displaystyle \frac{1}{n}}\sum _{i=1}^n{\displaystyle \frac{\left(g_i-{\nu }_0{u}_0\right)\left(g_i-{\nu }_0{u}_0-u_0\right)}{g_i^2{u}_0}}\hfill \\ \hfill & =J_{23},\hfill \end{array}$$

$$\begin{array}{cc}\hfill J_{33}& ={\displaystyle \frac{\partial f_3}{\partial u_0}}\hfill \\ \hfill & =-{\displaystyle \frac{1}{2}}\left(-{\displaystyle \frac{1}{u_0^2}}\right)+{\displaystyle \frac{1}{2}}{\displaystyle \frac{1}{n}}\sum _{i=1}^n\left({\nu }_0+1\right){\left(x_i-{\mu }_0\right)}^2\left(-1\right){\displaystyle \frac{{\nu }_{0}2u_0+{\left(x_i-{\mu }_0\right)}^2}{g_i^2{u}_0^2}}\hfill \\ \hfill & ={\displaystyle \frac{1}{2u_0^2}}-{\displaystyle \frac{1}{2}}{\displaystyle \frac{1}{n}}\sum _{i=1}^n\left({\nu }_0+1\right){\left(x_i-{\mu }_0\right)}^2{\displaystyle \frac{{\nu }_0{u}_0+g_i}{g_i^2{u}_0^2}}\hfill \\ \hfill & ={\displaystyle \frac{1}{2u_0^2}}-{\displaystyle \frac{1}{2}}{\displaystyle \frac{1}{n}}\sum _{i=1}^n\left({\nu }_0+1\right){\displaystyle \frac{\left(g_i-{\nu }_0{u}_0\right)\left(g_i+{\nu }_0{u}_0\right)}{g_i^2{u}_0^2}}\hfill \\ \hfill & ={\displaystyle \frac{1}{2u_0^2}}-{\displaystyle \frac{1}{2}}{\displaystyle \frac{1}{n}}\sum _{i=1}^n{\displaystyle \frac{\left({\nu }_0+1\right)\left(g_i^2-{\nu }_0^2{u}_0^2\right)}{g_i^2{u}_0^2}}.\hfill \end{array}$$

Note that

$${\psi }^{\prime }\left(x\right)={\displaystyle \frac{d^2}{d{x}^2}}log\Gamma \left(x\right)=\mathrm{trigamma}\left(x\right)$$

which can be directly calculated in R software by trigamma(x) ([46]). In $J_{ij}\left(i=1,2,3,j=1,2,3\right)$, the expressions involving $g_i$ are used for simplifying the R coding.

Now let us show that the MLEs are consistent estimators of the hyperparameters. From Theorem 10.1.6 in [47], we know that the MLEs are consistent estimators of the hyperparameters under some regularity conditions in Miscellanea 10.6.2 in [47]. The regularity conditions are listed below:

(C1). 

We observe $X_1,X_2,\dots ,X_n$, where $X_i\sim m\left(x|{\mu }_0,{\nu }_0,u_0\right)$ are iid.

(C2). 

The parameter is identifiable; that is, if $\left({\mu }_0,{\nu }_0,u_0\right)\ne \left({\mu }_0^{\prime },{\nu }_0^{\prime },u_0^{\prime }\right)$, then $m\left(x|{\mu }_0,{\nu }_0,u_0\right)\ne m\left(x|{\mu }_0^{\prime },{\nu }_0^{\prime },u_0^{\prime }\right)$.

(C3). 

The densities $m\left(x|{\mu }_0,{\nu }_0,u_0\right)$ have common support, and $m\left(x|{\mu }_0,{\nu }_0,u_0\right)$ is differentiable in ${\mu }_0$, ${\nu }_0$, and $u_0$.

(C4). 

The parameter space $\Omega$ contains an open set $\omega$ of which the true parameter value $\left({\mu }_{00},{\nu }_{00},u_{00}\right)$ is an interior point.

It remains to show that the marginal distribution $m\left(x|{\mu }_0,{\nu }_0,u_0\right)$ satisfies all the regularity conditions. Let $\mathcal{X}=\left(-\infty ,\infty \right)$ be the support set of $m\left(x|{\mu }_0,{\nu }_0,u_0\right)$.

First, (C1) is satisfied, as $X_1,X_2,\dots ,X_n$ is a random sample from $m\left(x|{\mu }_0,{\nu }_0,u_0\right)$.

Second, let us show that (C2) is satisfied. The parameter is identifiable

$$\iff \left({\mu }_0,{\nu }_0,u_0\right)\ne \left({\mu }_0^{\prime },{\nu }_0^{\prime },u_0^{\prime }\right)\Rightarrow m\left(x|{\mu }_0,{\nu }_0,u_0\right)\ne m\left(x|{\mu }_0^{\prime },{\nu }_0^{\prime },u_0^{\prime }\right)\mathrm{for}\mathrm{some}x\in \mathcal{X}$$

$$\iff m\left(x|{\mu }_0,{\nu }_0,u_0\right)=m\left(x|{\mu }_0^{\prime },{\nu }_0^{\prime },u_0^{\prime }\right)\mathrm{for}\mathrm{all}x\in \mathcal{X}\Rightarrow \left({\mu }_0,{\nu }_0,u_0\right)=\left({\mu }_0^{\prime },{\nu }_0^{\prime },u_0^{\prime }\right).$$

(A12)

We have

$$m\left(x|{\mu }_0,{\nu }_0,u_0\right)=m\left(x|{\mu }_0^{\prime },{\nu }_0^{\prime },u_0^{\prime }\right)\mathrm{for}\mathrm{all}x\in \mathcal{X}$$

$$\iff {\displaystyle \frac{m\left(x|{\mu }_0,{\nu }_0,u_0\right)}{m\left(x|{\mu }_0^{\prime },{\nu }_0^{\prime },u_0^{\prime }\right)}}=1\mathrm{for}\mathrm{all}x\in \mathcal{X}.$$

(A13)

Note that

$$m\left(x|{\mu }_0,{\nu }_0,u_0\right)={\displaystyle \frac{\Gamma \left(\frac{{\nu }_0+1}{2}\right)}{\Gamma \left(\frac{{\nu }_0}{2}\right)\sqrt{\pi {\nu }_0{u}_0}}}{\left[1+{\displaystyle \frac{1}{{\nu }_0}}{\displaystyle \frac{{\left(x-{\mu }_0\right)}^2}{u_0}}\right]}^{-{\displaystyle \frac{{\nu }_0+1}{2}}},$$

and thus

$$m\left(x|{\mu }_0^{\prime },{\nu }_0^{\prime },u_0^{\prime }\right)={\displaystyle \frac{\Gamma \left(\frac{{\nu }_0^{\prime }+1}{2}\right)}{\Gamma \left(\frac{{\nu }_0^{\prime }}{2}\right)\sqrt{\pi {\nu }_0^{\prime }{u}_0^{\prime }}}}{\left[1+{\displaystyle \frac{1}{{\nu }_0^{\prime }}}{\displaystyle \frac{{\left(x-{\mu }_0^{\prime }\right)}^2}{u_0^{\prime }}}\right]}^{-{\displaystyle \frac{{\nu }_0^{\prime }+1}{2}}}.$$

Therefore, (A13) is equivalent to

$${\displaystyle \frac{\Gamma \left(\frac{{\nu }_0+1}{2}\right)}{\Gamma \left(\frac{{\nu }_0}{2}\right)\sqrt{\pi {\nu }_0{u}_0}}}{\displaystyle \frac{\Gamma \left(\frac{{\nu }_0^{\prime }}{2}\right)\sqrt{\pi {\nu }_0^{\prime }{u}_0^{\prime }}}{\Gamma \left(\frac{{\nu }_0^{\prime }+1}{2}\right)}}{\displaystyle \frac{{\left[1+\frac{1}{{\nu }_0}\frac{{\left(x-{\mu }_0\right)}^2}{u_0}\right]}^{-\frac{{\nu }_0+1}{2}}}{{\left[1+\frac{1}{{\nu }_0^{\prime }}\frac{{\left(x-{\mu }_0^{\prime }\right)}^2}{u_0^{\prime }}\right]}^{-\frac{{\nu }_0^{\prime }+1}{2}}}}=1\mathrm{for}\mathrm{all}x\in \mathcal{X}$$

$$\iff {\displaystyle \frac{\Gamma \left(\frac{{\nu }_0+1}{2}\right)\Gamma \left(\frac{{\nu }_0^{\prime }}{2}\right)\sqrt{\pi {\nu }_0^{\prime }{u}_0^{\prime }}}{\Gamma \left(\frac{{\nu }_0^{\prime }+1}{2}\right)\Gamma \left(\frac{{\nu }_0}{2}\right)\sqrt{\pi {\nu }_0{u}_0}}}{\displaystyle \frac{{\left[1+\frac{1}{{\nu }_0^{\prime }}\frac{{\left(x-{\mu }_0^{\prime }\right)}^2}{u_0^{\prime }}\right]}^{\frac{{\nu }_0^{\prime }+1}{2}}}{{\left[1+\frac{1}{{\nu }_0}\frac{{\left(x-{\mu }_0\right)}^2}{u_0}\right]}^{\frac{{\nu }_0+1}{2}}}}=1\mathrm{for}\mathrm{all}x\in \mathcal{X},$$

which implies

$$\begin{array}{cc}\hfill {\displaystyle \frac{\Gamma \left(\frac{{\nu }_0+1}{2}\right)\Gamma \left(\frac{{\nu }_0^{\prime }}{2}\right)\sqrt{\pi {\nu }_0^{\prime }{u}_0^{\prime }}}{\Gamma \left(\frac{{\nu }_0^{\prime }+1}{2}\right)\Gamma \left(\frac{{\nu }_0}{2}\right)\sqrt{\pi {\nu }_0{u}_0}}}& ={\displaystyle \frac{1}{c_1}},\hfill \\ \hfill {\displaystyle \frac{{\left[1+\frac{1}{{\nu }_0^{\prime }}\frac{{\left(x-{\mu }_0^{\prime }\right)}^2}{u_0^{\prime }}\right]}^{\frac{{\nu }_0^{\prime }+1}{2}}}{{\left[1+\frac{1}{{\nu }_0}\frac{{\left(x-{\mu }_0\right)}^2}{u_0}\right]}^{\frac{{\nu }_0+1}{2}}}}& =c_1\mathrm{for}\mathrm{all}x\in \mathcal{X},\hfill \end{array}$$

(A14)

where $c_1$ is a constant which does not depend on x but may depend on ${\mu }_0$, ${\nu }_0$, $u_0$, ${\mu }_0^{\prime }$, ${\nu }_0^{\prime }$, and $u_0^{\prime }$. From (A14), we obtain

$${\mu }_0={\mu }_0^{\prime }$$

(A15)

and

$${\nu }_0{u}_0={\nu }_0^{\prime }{u}_0^{\prime }.$$

(A16)

Substituting (A15) and (A16) into (A14), we obtain

$${\left[1+{\displaystyle \frac{{\left(x-{\mu }_0\right)}^2}{{\nu }_0{u}_0}}\right]}^{{\displaystyle \frac{{\nu }_0^{\prime }-{\nu }_0}{2}}}=c_1\mathrm{for}\mathrm{all}x\in \mathcal{X},$$

which implies

$${\nu }_0={\nu }_0^{\prime }.$$

(A17)

Substituting (A17) into (A16), we obtain

$$u_0=u_0^{\prime }.$$

(A18)

Therefore,

$$\left({\mu }_0,{\nu }_0,u_0\right)=\left({\mu }_0^{\prime },{\nu }_0^{\prime },u_0^{\prime }\right).$$

Consequently, (A12) is correct, and (C2) is satisfied.

Third, (C3) is satisfied, as the densities $m\left(x|{\mu }_0,{\nu }_0,u_0\right)$ have common support $\mathcal{X}=\left(-\infty ,\infty \right)$, and $m\left(x|{\mu }_0,{\nu }_0,u_0\right)$ is differentiable in ${\mu }_0$, ${\nu }_0$, and $u_0$.

Finally, (C4) is satisfied, as the true parameter value

$$\left({\mu }_{00},{\nu }_{00},u_{00}\right)\in \omega =\Omega =\left(-\infty ,\infty \right)\times \left(0,\infty \right)\times \left(0,\infty \right),$$

which is an open set.

Therefore, the marginal densities $m\left(x|{\mu }_0,{\nu }_0,u_0\right)$ satisfy all the regularity conditions, and the MLEs are consistent estimators of the hyperparameters.

The proof of the theorem is complete. □