Objective and subjective prior distributions for the Gompertz distribution

Abstract

This paper takes into account the estimation for the unknown parameters of the Gompertz distribution from the frequentist and Bayesian view points by using both objective and subjective prior distributions. We first derive non-informative priors using formal rules, such as Jefreys prior and maximal data information prior (MDIP), based on Fisher information and entropy, respectively. We also propose a prior distribution that incorporate the expert’s knowledge about the issue under study. In this regard, we assume two independent gamma distributions for the parameters of the Gompertz distribution and it is employed for an elicitation process based on the predictive prior distribution by using Laplace approximation for integrals. We suppose that an expert can summarize his/her knowledge about the reliability of an item through statements of percentiles. We also present a set of priors proposed by Singpurwala assuming a truncated normal prior distribution for the median of distribution and a gamma prior for the scale parameter. Next, we investigate the effects of these priors in the posterior estimates of the parameters of the Gompertz distribution. The Bayes estimates are computed using Markov Chain Monte Carlo (MCMC) algorithm. An extensive numerical simulation is carried out to evaluate the performance of the maximum likelihood estimates and Bayes estimates based on bias, mean-squared error and coverage probabilities. Finally, a real data set have been analyzed for illustrative purposes.

Key words
Gompertz distribution; objective prior; Jeffreys prior; subjective prior; maximal data information prior; elicitation

INTRODUCTION

Gompertz distribution was introduced in connection with human mortality and actuarial sciences by Benzamin Gompertz (¹⁸²⁵). Right from the time of its introduction, this distribution has been receiving great attention from demographers and actuarist. This distribution is a generalization of the exponential distribution and is applied in various fields especially in reliability and life testing studies, actuarial science, epidemiological and biomedical studies. Gompertz distribution has some interesting relations with some of the well-known distributions such as exponential, double exponential, Weibull, extreme value (Gumbel Distribution) or generalized logistic distribution (^{Willekens 2002}). An important characteristic of the Gompertz distribution is that it has an exponentially increasing failure rate for the life of the systems and is often used to model highly negatively skewed data in survival analysis (^{Elandt-Johnson and Johnson 1979}). In recent past, many authors have contributed to the studies of statistical methodology and characterization of this distribution; for example, Garg et al. (¹⁹⁷⁰), Read (¹⁹⁸³), Makany (¹⁹⁹¹), Rao and Damaraju (¹⁹⁹²), Franses (¹⁹⁹⁴), Chen (¹⁹⁹⁷) and Wu and Lee (¹⁹⁹⁹). Jaheen (^2003a, ^b) studied this distribution based on progressive type-II censoring and record values using Bayesian approach. Wu et al. (²⁰⁰³) derived the point and interval estimators for the parameters of the Gompertz distribution based on progressive type II censored samples. Wu et al. (²⁰⁰⁴) used least squared method to estimate the parameters of the Gompertz distribution. Wu et al. (²⁰⁰⁶) also studied this distribution under progressive censoring with binomial removals. Ismail (²⁰¹⁰) obtained Bayes estimators under partially accelerated life tests with type-I censoring. Ismail (²⁰¹¹) also discussed the point and interval estimations of a two-parameter Gompertz distribution under partially accelerated life tests with Type-II censoring. Asgharzadeh and Abdi (²⁰¹¹) studied different types of exact confidence intervals and exact joint confidence regions for the parameters of the two-parameter Gompertz distribution based on record values. Kiani et al. (²⁰¹²) studied the performance of the Gompertz model with time-dependent covariate in the presence of right censored data. Moreover, they compared the performance of the model under different censoring proportions (CP) and sample sizes. Shanubhogue and Jain (²⁰¹³) studied uniformly minimum variance unbiased estimation for the parameter of the Gompertz distribution based on progressively Type II censored data with binomial removals. Lenart (²⁰¹⁴) obtained moments of the Gompertz distribution and maximum likelihood estimators of its parameters. Lenart and Missov (²⁰¹⁶) studied Goodness-of-fit tests for the Gompertz distribution. Recently, Singh et al. (²⁰¹⁶) studied different methods of estimation for the parameters of Gompertz distribution when the available data are in the form of fuzzy numbers. They also obtained Bayes estimators of the parameters under different symmetric and asymmetric loss functions.

In this paper, we present a Bayesian analysis when there is a limited prior knowledge about the parameter of interest. In this regard, it is important to use noninformative priors, however, it can be difficult to choose a prior distribution that represents this situation, because there is hardly any precise definition of the concept of noninformative prior. Nevertheless, we have many noninformative priors, for instance, Jeffreys prior (^{Jeffreys 1967}), MDIP prior (^{Zellner 1977}, ¹⁹⁸⁴, ¹⁹⁹⁰), Tibshirani prior (^{Tibshirani 1989}), reference prior (^{Bernardo 1979}) and many others which seemingly appropriate for a number of inference problems. It is to be noted that lack of enough information on the part of analysts often forces them to choose noninformative priors and this consideration ensures that the inferences are mostly data driven. In Bayesian analysis, many authors consider independent gamma priors for the estimation of parameters of the model, representing weak information as the use of a priori independence assumption simplifies the computations. Our main interest in the Bayesian analysis is to select a prior distribution that represents better dependence structure of the parameters in which the information regarding the parameters is not considered substantial as compared with information from the data. The focus is on the comparison of independent gamma prior, Jeffreys prior, maximal data information prior (MDIP), Singpurwalla’s prior and elicited prior. Jeffreys (1967) proposed a noninformative prior resulting from an argument based on the Fisher Information Measure and Zellner (1977, 1984) proposed an alternative prior, named maximal data information prior (MDIP) based on the Entropy Measure. The prior proposed by Singpurwalla (1988) for estimation of the parameters of Weibull distribution is also considered in this paper to estimate the parameters of Gompertz distribution.

There are many methods for eliciting parameters of prior distributions. In this paper, we also consider an elicitation method to specify the values of hyperparameters of the two gamma priors assigned to the parameters of the Gompertz distribution. The method requires the derivation of predictive prior distribution and it is assumed that the expert is able to provide some percentiles values. Thus, the main aim of this paper is to propose noninformative and informative prior distributions for the parameters $c$ and $\lambda$ of the Gompertz distribution and to study the effects of these different priors in the resulting posterior distributions, especially in situations of small sample sizes, a common situation in applications.

The paper is organized as follows. Some probability properties of the Gompertz distribution such as quantiles, moments, moment generating function are reviewed in Section 2. Section 3 describes the maximum likelihood estimation method. The Bayesian approach with proposed informative and noninformative priors is presented in section 4. In Section 5, simulation study is carried out to evaluate the performance of several estimation procedures along with coverage percentages is provided. The methodology developed in this paper and the usefulness of the Gompertz distribution is illustrated by using a real data example in Section 6. Finally, concluding remarks are provided in Section 7.

MODEL AND ITS BASIC PROPERTIES

A random variable X has the Gompertz distribution with parameters $c$ and $\lambda$, say $GM(c,\lambda )$, if its density function is

$f\left(x\right)=\lambda e^{cx}{e}^{-\frac{\lambda }{c}(e^{cx}-1)};x>0,c,\lambda >0,$ (1)

and the corresponding c.d.f is given by

$F\left(x\right)=1-e^{-\frac{\lambda }{c}(e^{cx}-1)};x>0,c,\lambda >0.$ (2)

The basic tools for studying the ageing and reliability characteristics of the system are the hazard rate $h\left(x\right)$. The hazard function gives the rate of failure of the system immediately after time x. Thus the hazard rate function of the Gompertz distribution is given by

$h\left(x\right)=\frac{f\left(x\right)}{1-F\left(x\right)}=\lambda e^{cx}.$ (3)

Note that the hazard rate function is increasing function if $c>0$ or constant if $c=0$. Figure 1b shows the shapes of the hazard function for different selected values of the parameters $c$ and $\lambda$. From the plot, it is quite evident that the Gompertz distribution has increasing hazard rate function.

Figure 1a shows the shapes of the pdf of the Gompertz distribution for different values of the parameters $c$ and $\lambda$ and from the plot, it is quite evident that the Gompertz distribution is positively skewed distribution.

Figure 1
Pdf function (a) and hazard function (b).

The quantile function $x_p=Q\left(p\right)=F^{-1}\left(p\right)$, for $0<p<1$, of the Gompertz distribution is obtained from ([eq.2]), thus the quantile function $x_p$ is

$x_p=\frac{1}{c}ln(1-\frac{c}{\lambda }ln(1-p\left)\right).$ (4)

In particular, the median of the Gompertz distribution can be written as

$Md\left(X\right)=M_d=\frac{1}{c}ln(1-\frac{c}{\lambda }ln(1-0.5\left)\right).$ (5)

If the random variable $X$ is distributed $GM(c,\lambda )$, then its $n$th moment around zero can be expressed as

$E\left(X^n\right)=\frac{\lambda e^{\frac{\lambda }{c}}}{c}{\int }_1^{\infty }\frac{1}{c^n}{e}^{-\frac{\lambda }{c}x}[ln(x){]}^{n}dx.$ (6)

On simplification, we get

$E\left(X^n\right)=\frac{n!}{c^n}{e}^{\frac{\lambda }{c}}{E}_1^{n-1}\left(\frac{\lambda }{c}\right),$ (7)

where

$E_s^n\left(z\right)=\frac{1}{n!}{\int }_1^{\infty }(ln(x){)}^n{x}^{-s}{e}^{-zx}dx,$

$E_n(x)={\int }_1^{\mathrm{\infty }}\frac{e^{-xt}}{t^n}dt$

and

$E_s^0\left(z\right)=E_s\left(z\right),$

is the generalized integro-exponential function (^{Milgram 1985}).

The mean and variance of the random variable X of the Gompertz distribution are respectively, given by

$E\left(X\right)=\frac{1}{c}{e}^{\frac{\lambda }{c}}{E}_1\left(\frac{\lambda }{c}\right)$ (8)

and

$var(X)=\frac{2}{c^2}{e}^{\frac{\lambda }{c}}{E}_1^1(\frac{\lambda }{c})-(\frac{1}{c}{e}^{\frac{\lambda }{c}}{E}_1(\frac{\lambda }{c}){)}^2$ (9)

Many of the interesting characteristics and features of a distribution can be obtained via its moment generating function and moments. Let X denote a random variable with the probability density function (1). By definition of moment generating function of X and using (1), we have

$M_x\left(t\right)=E\left(e^{tx}\right)={\int }_0^{\infty }{e}^{tx}f\left(x\right)dx=\frac{\lambda }{c}\sum _{p=0}^{\infty }(-1{)}^p\left(\begin{array}{c}\hfill \frac{t}{\lambda }\hfill \\ \hfill p\hfill \end{array}\right)\frac{\Gamma (p+1)}{{\left(\frac{\lambda }{c}\right)}^{p+1}}.$ (10)

MAXIMUM LIKELIHOOD ESTIMATION

The method of maximum likelihood is the most frequently used method of parameter estimation (^{Casella and Berger 2001}). The success of the method stems no doubt from its many desirable properties including consistency, asymptotic efficiency, invariance property as well as its intuitive appeal. Let $x_1,\cdots ,x_n$ be a random sample of size $n$ from (1), then the log-likelihood function of (1) without constant terms is given by

$\begin{array}{ccc}\hfill \ell (c,\lambda ;x)& =\hfill & \hfill logL(c,\lambda ;x)\\ \hfill & =\hfill & \hfill nlog\lambda +c\sum _{i=1}^n{x}_i-\frac{\lambda }{c}\sum _{i=1}^n(e^{c{x}_i}-1).\end{array}$

For ease of notation, we denote the first partial derivatives of any function $f(x,y)$ by $f_x$ and $f_y$. Now setting

${\ell }_c=0\mathrm{\text{and}}{\ell }_{\lambda }=0,$

we have

${\ell }_{\lambda }=\frac{n}{\lambda }-\frac{1}{c}\sum _{i=1}^n(e^{c{x}_i}-1)=0,$ (11)

and

${\ell }_c=\sum _{i=1}^n{x}_i-\frac{\lambda }{c}\sum _{i=1}^n{x}_i{e}^{c{x}_i}+\frac{\lambda }{c^2}\sum _{i=1}^n{e}^{c{x}_i}-\frac{n\lambda }{c^2}=0.$ (12)

From (11) and (12), we find the MLE for $\lambda$ given by,

$\hat{\lambda }=\frac{n\hat{c}}{\sum _{i=1}^n(e^{\hat{c}{x}_i}-1)}$

The MLE for "$c$" is obtained by solving the non-linear equation,

$\sum _{i=1}^n{x}_i-\frac{n\sum _{i=1}^n{x}_i{e}^{c{x}_i}}{\sum _{i=1}^n(e^{c{x}_i}-1)}+\frac{n}{c}=0$

The asymptotic distribution of the MLE $\hat{\theta }$ is

$(\hat{\theta }-\theta )\to N_2(0,I^{-1}(\theta \left)\right)$

(^{Lawless 2003}), where $I^{-1}\left(\theta \right)$ is the inverse of the observed information matrix of the unknown parameters $\theta =(c,\lambda )$.

$I^{-1}(\theta )={\left(\begin{array}{c}-\frac{{\mathrm{\partial }}^2\log L}{\mathrm{\partial }{c}^2}-\frac{{\mathrm{\partial }}^2\log L}{\mathrm{\partial }c\mathrm{\partial }\lambda }\\ -\frac{{\mathrm{\partial }}^2\log L}{\mathrm{\partial }\lambda \mathrm{\partial }c}-\frac{{\mathrm{\partial }}^2\log L}{\mathrm{\partial }{\lambda }^2}\end{array}\right)|}_{(c,\lambda )=(\hat{c},\hat{\lambda })}^{-1}$ (13)

$=\left(\begin{array}{c}{var(\hat{c}}_{MLE}){cov(\hat{c}}_{MLE}{,\hat{\lambda }}_{MLE})\\ {cov(\hat{\lambda }}_{MLE}{,\hat{c}}_{MLE}){var(\hat{\lambda }}_{MLE})\end{array}\right)=\left(\begin{array}{c}{\sigma }_{cc}{\sigma }_{c\lambda }\\ {\sigma }_{\lambda c}{\sigma }_{\lambda \lambda }\end{array}\right)$

The derivatives in $I\left(\theta \right)$ are given as follows

${\frac{{\mathrm{\partial }}^2\log L}{\mathrm{\partial }{c}^2}|}_{c={\hat{c}}_{MLE}}=\frac{\lambda }{c^2}[-c\underset{i=1}{\overset{n}{\sum }}{x}_i^2{e}^{c{x}_i}+2\underset{i=1}{\overset{n}{\sum }}{x}_i{e}^{c{x}_i}-\frac{2}{c}\underset{i=1}{\overset{n}{\sum }}{e}^{c{x}_i}+\frac{2n}{c}]$ (14)

${\frac{{\partial }^{2}logL}{\partial {\lambda }^2}|}_{\lambda ={\hat{\lambda }}_{MLE}}=-\frac{n}{{\lambda }^2},$ (15)

${\frac{{\partial }^{2}logL}{\partial c\partial \lambda }|}_{c={\hat{c}}_{MLE},\lambda ={\hat{\lambda }}_{MLE}}=\left[\frac{1}{c^2}\sum _{i=1}^n\right(e^{c{x}_i}-1)-\frac{1}{c}\sum _{i=1}^n{x}_i{e}^{c{x}_i}].$ (16)

Therefore, the above approach is used to derive the approximate $100(1-\tau )\%$ confidence intervals of the parameters $\theta =(c,\lambda )$ as in the following forms

${\hat{c}}_{MLE}\pm z_{\frac{\tau }{2}}\sqrt{Var\left({\hat{c}}_{MLE}\right)},{\hat{\lambda }}_{MLE}\pm z_{\frac{\tau }{2}}\sqrt{Var\left({\hat{\lambda }}_{MLE}\right)}.$

Here, $Z_{\frac{\tau }{2}}$ is the upper ($\frac{\tau }{2}$)th percentile of the standard normal distribution.

BAYESIAN ANALYSIS

In this section, we consider Bayesian inference of the unknown parameters of the $GM(c,\lambda )$. First, we assume that $c$ and $\lambda$ has the independent gamma prior distributions with probability density functions

$\pi \left(c\right)\propto c^{a_1-1}{e}^{-b_{1}c}c>0$ (17)

and

$\pi \left(\lambda \right)\propto {\lambda }^{a_2-1}{e}^{-b_2\lambda }\lambda >0.$ (18)

The hyperparameters $a_1$, $a_2$, $b_1$ and $b_2$ are known and positives. If both parameters $c$ and $\lambda$ are unknown, we cannot see an easy way to work with conjugation, since the expression of the likelihood function does not suggest any known form for the joint density of $(c,\lambda )$. It is not unreasonable to assume independent gamma priors on the shape and scale parameters for a two-parameter $GM(c,\lambda )$, because gamma distributions are very flexible. The joint prior distribution for both parameters in this case is given by

$\pi (c,\lambda )\propto c^{a_1-1}exp(-b_{1}c){\lambda }^{a_2-1}exp(-b_2\lambda ).$ (19)

Thus, the joint posterior distribution is given by

$p(c,\lambda |𝐱)\propto {\lambda }^{n+a_2-1}{c}^{a_1-1}{e}^{c(\sum _{i=1}^n{x}_i-b_1)}{e}^{-\lambda \left[\frac{1}{c}\sum _{i=1}^n\right(e^{c{x}_i}-1)+b_2]}.$ (20)

The conditional distribution of $c$ given $\lambda$ and data is given by

$p\left(c\right|\lambda ,𝐱)\propto c^{a_1-1}{e}^{c(\sum _{i=1}^n{x}_i-b_1)}{e}^{-\frac{\lambda }{c}\sum _{i=1}^n(e^{c{x}_i}-1)}.$ (21)

Similarly, the conditional distribution of $\lambda$ given $c$ and data is given by

$p\left(\lambda \right|c,𝐱)\propto {\lambda }^{n+a_2-1}{e}^{-\lambda \left[\frac{1}{c}\sum _{i=1}^n\right(e^{c{x}_i}-1)+b_2]}.$ (22)

Note that the although the conditional $p\left(\lambda \right|c,𝐱)$ is a gamma distribution, the conditional distributions $p\left(c\right|\lambda ,𝐱)$ is not identified as known distributions that are easy to simulate. In this way, Bayesian inference for the parameters $c$and $\lambda$ can be performed by Metropolis-Hastings (MH) algorithm considering the gamma distribution as the target density for $\lambda$, and $c$ can be generated from the conditional $p\left(c\right|\lambda ,𝐱)$ by using the rejection method.

JEFFREYS PRIOR

A well known non-informative prior, which represents a situation with little a priori information on the parameters was introduced by Jeffreys (1967), also known as the Jeffreys rule. The Jeffreys prior has been widely used due to the invariance property for one to one transformations of the parameters. Since then Jeffreys prior has played an important role in Bayesian inference. This prior is derived from Fisher Information matrix $I$($c$, $\lambda$) as

$\pi (c,\lambda )\propto \sqrt{det\left(I\right(c,\lambda \left)\right)}.$ (23)

However, $I$($c$, $\lambda$) can not be analytically obtained for the parameters of Gompertz distribution. A possible simplification is to consider a noninformative prior given by $\pi (c,\lambda )=\pi \left(\lambda \right|c\left)\pi \right(c)$. Using the Jeffreys’rule, we have,

$\pi (c,\lambda )=[E(-\frac{{\mathrm{\partial }}^2}{\mathrm{\partial }{\lambda }^2}logL(c,\lambda )){]}^{\frac{1}{2}}\pi (c)$ (24)

where $E(-\frac{{\mathrm{\partial }}^2}{\mathrm{\partial }{\lambda }^2}logL(c,\lambda ))$ is given by (15) and $\pi \left(c\right)$ is a noninformative prior, for instance, a gamma distribution with hyper-parameters equal to 0.01.

In this way, from (15) and (24) the non-informative prior for ($c$, $\lambda$) parameters is given by:

$\pi (c,\lambda )=\frac{1}{\lambda }\pi \left(c\right).$ (25)

Let us denote the prior (25) as "Jeffreys prior".

Thus, the corresponding posterior distribution is given by

$p(\lambda ,c|\mathit{x})={\lambda }^{n-1}exp(c\underset{i=1}{\overset{n}{\sum }}{x}_i-\frac{\lambda }{c}\underset{i=1}{\overset{n}{\sum }}{e}^{c{x}_i}+\frac{n\lambda }{c})\pi (c)$ (26)

Proposition 1: For the parameters of the Gompertz distribution, the posterior distribution given in (26) under Jeffreys prior $\pi$($\lambda$, $c$) given in (25) is proper.

We need to prove that ${\int }_0^{\infty }{\int }_0^{\infty }p(\lambda$, $c|𝐱$)$d\lambda dc$ is finite.

Indeed,

${\int }_0^{\mathrm{\infty }}{\int }_0^{\mathrm{\infty }}p(\lambda ,c|\mathit{x})d\lambda dc=\mathrm{\Gamma }(n){\int }_0^{\mathrm{\infty }}\frac{c^{n}exp(c\underset{i=1}{\overset{n}{\sum }}{x}_i)}{(\underset{i=1}{\overset{n}{\sum }}(e^{c{x}_i}-1){)}^n}\pi (c)dc$ (27)

The function $h(c)=\frac{c^{n}exp(c\underset{i=1}{\overset{n}{\sum }}{x}_i)}{(\underset{i=1}{\overset{n}{\sum }}(e^{c{x}_i}-1){)}^n}$ is unimodal with maximum point $\hat{c}$ as the solution of the nonlinear equation given by

$\frac{c\sum _{i=1}^n{x}_i{e}^{c{x}_i}-\sum _{i=1}^n{e}^{c{x}_i}+n}{c\sum _{i=1}^n(e^{c{x}_i}-1)}=\overline{X}$

Therefore, from (27) we have

${\int }_0^{\infty }{\int }_0^{\infty }p(\lambda ,c|𝐱)d\lambda dc\le \Gamma (n\left)h\right(\hat{c}\left){\int }_0^{\infty }\pi \right(c)dc<\infty ,$

where ${\int }_0^{\infty }\pi$($c$)$dc=1$. This completes the proof.

MAXIMAL DATA INFORMATION PRIOR (MDIP)

It is interesting to note that the data gives more information about the parameter than the information from the prior density, otherwise, there would not be justification for the realization of the experiment. Let $X$ be a random variable with density $f\left(x\right|\varphi )$, $x\in R_X$ ($R_X\subseteq \Re$), parameter $\varphi \in$a, b. Thus, we wish a prior distribution $\pi \left(\varphi \right)$ that provides the gain in the information supplied by the data as much as possible relative to the prior information of the parameter, that is, maximizes the information on the data. With this idea, Zellner (1977, 1984, 1990) and Zellner and Min (¹⁹⁹³) derived a prior distribution which maximize the information from the data in relation to the prior information on the parameters. Let

$H\left(\varphi \right)={\int }_{R_X}f\left(x\right|\varphi \left)lnf\right(x\left|\varphi \right)dx$ (28)

be a negative entropy of $f$($x|\varphi$), the measure of the information in $f$($x|\varphi$). Thus, the following functional form is employed in the MDIP approach:

$G[\pi (\varphi )]={\int }_{\text{a}}^{\text{b}}H(\varphi )\pi (\varphi )d\varphi -{\int }_{\text{a}}^{\text{b}}\pi (\varphi )ln\pi (\varphi )d\varphi$ (29)

which is the prior average information in the data density minus the information in the prior density. $G[\pi$($\varphi$)$]$ is maximized by selection of $\pi$($\varphi$) subject to ${\int }_{\mathrm{\text{a}}}^{\mathrm{\text{b}}}\pi$($\varphi$)$d\varphi =1$.

The following theorem proposed by Zellner provides the formula for the MDIP prior.

Theorem: The MDIP prior is given by:

$\pi (\varphi )=kexp(H(\varphi ))a\le \varphi \le b$ (30)

where $k^{-1}={\int }_{\text{a}}^{\text{b}}exp(H(\varphi ))d\varphi$ is the normalizing constant.

Proof. We have to maximize the function $U=G[\pi$($\varphi$)]$-\lambda ({\int }_{\text{a}}^{\text{b}}\pi$($\varphi$)$d\varphi -1)$ where $\lambda$ is the Lagrange multiplier.$$

Thus, $\frac{\partial U}{\partial \pi }=0$ $\iff$ $H\left(X\right)=-ln(\pi$($\varphi$)$)-1-\lambda =0$ and the solution is given by $\pi$($\varphi$)$=k$ $exp(H(X))$.

Therefore, the MDIP is a prior that leads to an emphasis on the information in the data density or likelihood function, that is, its information is weak in comparison with data information.

Zellner (1984) shows several interesting properties of MDIP and additional conditions that can also be imposed to the approach refleting given initial information.

Suppose that we do not have much prior information available about $c$ and $\lambda$. Under this condition, the prior distribution MDIP for the parameters ($c$, $\lambda$) of Gompertz distribution (1) is obtained as follows. Firstly, we have to evaluate the measure of information $H$($c,$ $\lambda$)$=E(logf(x))$, that is

$E(logf(x))=log(\lambda )+cE(X)-\frac{\lambda }{c}(E(e^{cX})-1)$ (31)

where $E\left(X\right)$is obtained from (8) given by $E(X)=\frac{1}{c}{e}^{\frac{\lambda }{c}}{E}_1(\frac{\lambda }{c})$ and $E_1\left(x\right)={\int }_1^{\infty }\frac{e^{-xu}}{u}du$ is the Exponential Integral. After some algebra, we also obtain $E\left(e^{cX}\right)=\frac{c}{\lambda }+1$. Therefore,

$H(c,\lambda )=log(\lambda )+e^{\frac{\lambda }{c}}{E}_1(\frac{\lambda }{c})-\frac{\lambda }{c}(\frac{c}{\lambda }+1)+\frac{\lambda }{c}$ (32)

Hence the MDIP prior is given by

${\pi }_Z(c,\lambda )\propto \lambda exp(e^{\frac{\lambda }{c}}{E}_1(\frac{\lambda }{c}))$ (33)

Now combining the likelihood function given by

$L(\lambda ,c|\mathit{x})={\lambda }^{n}exp(c\underset{i=1}{\overset{n}{\sum }}{x}_i-\frac{\lambda }{c}\underset{i=1}{\overset{n}{\sum }}{e}^{c{x}_i}+\frac{n\lambda }{c})$ (34)

and the MDIP prior in (33), the posterior densitiy for the parameters $c$ and $\lambda$ is given by

$p(\lambda ,c|\mathit{x})={\lambda }^{n+1}exp(c\underset{i=1}{\overset{n}{\sum }}{x}_i+e^{\frac{\lambda }{c}}{E}_1(\frac{\lambda }{c})-\frac{\lambda }{c}\underset{i=1}{\overset{n}{\sum }}(e^{c{x}_i}-1))$ (35)

Proposition 2: For the parameters of the Gompertz distribution, the posterior distribution given in (35) under the corresponding MDIP prior $\pi$($\lambda$, $c$) given in (33) is proper.

Proof. Indeed,

${\int }_0^{\mathrm{\infty }}{\int }_0^{\mathrm{\infty }}p(\lambda ,c|\mathit{x})d\lambda dc={\int }_0^{\mathrm{\infty }}{\int }_0^{\mathrm{\infty }}{\lambda }^{n+1}exp(c\underset{i=1}{\overset{n}{\sum }}{x}_i+e^{\frac{\lambda }{c}}{E}_1(\frac{\lambda }{c})-\frac{\lambda }{c}\underset{i=1}{\overset{n}{\sum }}(e^{c{x}_i}-1))d\lambda dc$

Now, we consider a substituition of variables in the integral above as

$\{\begin{array}{cc}w=\frac{\lambda }{c}\hfill & \hfill \\ u=c\hfill & \hfill \end{array}\Rightarrow J=\mid \begin{array}{cc}\hfill w\hfill & \hfill u\hfill \\ \hfill 1\hfill & \hfill 0\hfill \end{array}\mid =-u$

resulting in

${\int }_0^{\mathrm{\infty }}{\int }_0^{\mathrm{\infty }}p(\lambda ,c|\mathit{x})d\lambda dc={\int }_0^{\mathrm{\infty }}[{\int }_0^{\mathrm{\infty }}{u}^{n+2}exp(u\underset{i=1}{\overset{n}{\sum }}{x}_i-w\underset{i=1}{\overset{n}{\sum }}{e}^{u{x}_i})du]h(w)dw$

where $h\left(w\right)=w^{n+1}exp(nw+e^w{E}_1(w\left)\right)$.

Let us denote $x_{\left(n\right)}=max(x_1,x_2,...,x_n)$then ${\sum }_{i=1}^n{e}^{u{x}_i}<n{e}^{u{x}_{\left(n\right)}}$. Hence,

${\int }_0^{\mathrm{\infty }}{\int }_0^{\mathrm{\infty }}p(\lambda ,c|\mathit{x})d\lambda dc<{\int }_0^{\mathrm{\infty }}[{\int }_0^{\mathrm{\infty }}{u}^{n+2}exp(n\bar{x}u-nw{e}^{u{x}_{(n)}})du]h(w)dw$

where $\overline{x}=\frac{{\sum }_{i=1}^n{x}_i}{n}$. As $\overline{x}<x_{\left(n\right)}$we have

${\int }_0^{\mathrm{\infty }}{u}^{n+2}exp(n\bar{x}u-nw{e}^{u{x}_{(n)}})du<{\int }_0^{\mathrm{\infty }}{u}^{n+2}exp(n{x}_{(n)}u-nw{e}^{u{x}_{(n)}})du<$

$<{\int }_0^{\mathrm{\infty }}{u}^{n+2}exp(n{x}_{(n)}u-w{e}^{u{x}_{(n)}})du$

Now consider $e^{u{x}_{\left(n\right)}}=z$ then by substitution process the integral above becomes

${\displaystyle {\int }_0^{\mathrm{\infty }}{u}^{n+2}exp(n{x}_{(n)}u-w{e}^{u{x}_{(n)}})du=\frac{1}{x_{(n)}^{n+3}}{\int }_1^{\mathrm{\infty }}(logz{)}^{n+2}{z}^{n-1}exp(-wz)dz{\displaystyle =\frac{1}{x_{(n)}^{n+3}}\Gamma (n+2)G_{n\text{+}3,n\text{+4}}^{0,0}(\begin{array}{c}\text{-(}n-1\text{), . . ., -(}n-1\text{)}\\ \text{-n, . . ., -n, 0}\end{array}|w)}}$ (36)

where $G_{p,q}^{m,n}(\begin{array}{c}{\text{a}}_1,...,{\text{a}}_p\\ b_1,...,b_q\end{array}|w)$ is the Meijer G-function introduced by Meijer (¹⁹³⁶) and given by

$G_{p,q}^{m,n}(\begin{array}{c}{\text{a}}_1,...,{\text{a}}_p\\ b_1,...,b_q\end{array}|w)=\frac{1}{2\pi i}{\int }_L\frac{\underset{j=1}{\overset{m}{\prod }}\mathrm{\Gamma }(b_j-s)\underset{j=1}{\overset{n}{\prod }}\mathrm{\Gamma }(1-{\text{a}}_j+s)}{\underset{j=m+1}{\overset{q}{\prod }}\mathrm{\Gamma }(1-b_j+s)\underset{j=n+1}{\overset{p}{\prod }}\mathrm{\Gamma }({\text{a}}_j-s)}{z}^{s}ds$

for 0$\le m\le q$ and 0$\le n\le p$, where $m$, $n$, $p$ and $q$ are integer numbers.

From (36) we have

${\int }_0^{\mathrm{\infty }}{\int }_0^{\mathrm{\infty }}p(\lambda ,c|\mathit{x})d\lambda dc<\frac{1}{x_{(n)}^{n+3}}\Gamma (n+2){\int }_0^{\mathrm{\infty }}{G}_{n\text{+}3,n\text{+4}}^{0,0}(\begin{array}{c}\text{-(}n-1\text{), . . ., -(}n-1\text{)}\\ \text{-n, . . ., -n, 0}\end{array}|w)h(w)dw$

which is not possible to obtain an analytical expression for this integral. However, the software Mathematica gives a convergence result of the integral.

PRIORS PROPOSED BY SINGPURWALLA

Singpurwalla (¹⁹⁸⁸) presented a procedure for the construction of the prior distribution with the use of expert opinion in order to estimate the parameters $\alpha$ and $\beta$ of Weibull distribution. Expert’s opinion about measures of central tendency such as the median can be easily found, since most people are accustomed to this term. Singpurwalla introduces the median life $M$ and elicit expert opinion on $M$ and $\beta$ through the priors $\pi$($M$) and $\pi$($\beta$). He focuses attention on $\beta$ and on $M=\alpha exp(\frac{k}{\beta })$ where $k=ln\left(ln2\right)$. Since $M$ is restricted to being nonnegative, it is assumed a Gaussian distribution truncated at 0 with parameters $\mu$ and $\sigma$. A gamma prior distribution with parameters $a$ and $b$ is chosen to model the uncertainty about $\beta$. After this reparametrization, the prior $\pi (\alpha ,\beta )$ is derived by transformation of variables.

Our aim is to derive the prior $\pi$($c$, $\lambda$) applying a similar procedure proposed by Singpurwalla in order to estimate the parameters $c$ and $\lambda$ of Gompertz distribution. Differently of Singpurwalla’s priors who considered elicitation from expert for the parameters, we assume absence of information, hence the hyperparameters of the priors are chosen to provide noninformative prior and we use the information from the data to the parameter $\mu$ through the median of the data.

Consider the median of $X$ is given by $M=\frac{1}{c}log(1+0.6931\frac{c}{\lambda })$. Thus, it is assumed a Gaussian distribution truncated at 0 for the parameter $M$, making the density as

$\pi (M|\mu ,\sigma )\propto exp(-\frac{1}{2}(\frac{M-\mu }{\sigma }{)}^2)$

where $0\le M<\infty$ with parameter $\mu$ equal to median of the data and standard deviation $\sigma =100$.

A gamma prior distribution is chosen to model the uncertainty about $\lambda$ with density

$\pi \left(\lambda \right|a,b)\propto {\lambda }^{\mathrm{\text{a}}-1}exp(-b\lambda ),$ (38)

with the parameters $a$ and $b$ specified as 0.01 representing a noninformative prior for $\lambda$.

Thus, we can determine the conditional prior distribution $\pi$($c$ $|\lambda$, $\mu$, $\sigma$) for the parameter $c$ given $\lambda$ through the reparametrization $M=\frac{1}{c}log(1+0.6931\frac{c}{\lambda })$ with the Jacobian given by $|\frac{dM}{dc}|=|\frac{0.6931}{c\lambda +0.6931c^2}-\frac{1}{c^2}log(1+0.6931\frac{c}{\lambda })|$. Therefore, the conditional prior $\pi$($c$ $|\lambda$, $\mu$, $\sigma$) is given by

$\pi (c|\lambda ,\mu ,\sigma )\propto exp(-\frac{1}{2}(\frac{\frac{1}{c}log(1+0.6931\frac{c}{\lambda })-\mu }{\sigma }{)}^2)|\frac{0.6931}{c\lambda +0.6931c^2}-\frac{1}{c^2}log(1+0.6931\frac{c}{\lambda })|$ (39)

Finally, the joint prior for $c$ and $\lambda$, obtained as $\pi$($c$, $\lambda$)$=\pi$($c|\lambda$)$\pi$($\lambda$), that is, by the product of (38) and (39), is given by

$\pi (c,\lambda |\mathrm{\Theta })\propto {\lambda }^{\text{a}-1}exp(-\frac{1}{2}(\frac{\frac{1}{c}log(1+0.6931\frac{c}{\lambda })-\mu }{\sigma }{)}^2-b\lambda )|\frac{0.6931}{c\lambda +0.6931c^2}-\frac{1}{c^2}log(1+0.6931\frac{c}{\lambda })|$ (40)

where the vector of parameters $\Theta =($a$,b,\mu ,\sigma )$ is known.

ELICITED PRIOR

In this Section, we provide a methodology that permits the experts to use their knowledges about the reliability of an item through statements of percentiles. This method requires the derivation of prior predictive distribution for elicitation. Suppose that joint prior $\pi (c,\lambda )$ is given, then the reliability based on the predictive prior distribution is given by:

$R\left(t_p\right)=P\{T\ge t_p\}={\int }_0^{\infty }{\int }_0^{\infty }{\int }_{t_p}^{\infty }f\left(t\right|c,\lambda \left)\pi \right(c,\lambda )dtdcd\lambda ,$ (41)

for a fixed mission time $t_p$.

In order to elicit the four hyperparameters $(a_1,a_2,b_1,b_2)$ of the prior $\pi (c,\lambda |a_1,a_2,b_1,b_2)$ the integral have been considered as follows

$p={\int }_0^{\infty }{\int }_0^{\infty }R\left(t_p\right|c,\lambda \left)\pi \right(c,\lambda |a_1,a_2,b_1,b_2)dcd\lambda ,$ (42)

for a given p-th percentile elicited from the expert where $p=R\left(t_p\right)$.

By considering Gompertz distribution, the reliability function $R\left(t_p\right|c$, $\lambda )$ is given by

$R(t_p|c,\lambda )=exp(-\frac{\lambda }{c}(e^{c{t}_p}-1))$ (43)

and assuming a joint prior $\pi (c$, $\lambda |$ a${}_1,$a${}_2,b_1,b_2)$ given by the product of gamma priors, we have

$\pi (c,\lambda |a_1,a_2,b_1,b_2)=k{c}^{a_1-1}{\lambda }^{a_2-1}exp(-(b_{1}c+b_2\lambda ))$ (44)

where$k=\frac{b_1^{a_1}{b}_2^{a_2}}{\Gamma \left(a_1\right)\Gamma \left(a_2\right)}$.

Using (42), (43) and (44), the probability in (42) becomes

$p=k{\int }_0^{\mathrm{\infty }}{\int }_0^{\mathrm{\infty }}{c}^{a_1-1}{\lambda }^{a_2-1}exp(-\frac{\lambda }{c}(e^{c{t}_p}-1)-(b_{1}c+\lambda b_2))dcd\lambda$ (45)

Let $d=\frac{1}{c}(e^{c{t}_p}-1)$, then the integral for $\lambda$ in equation (45) takes the gamma shape resulting in

$p=k{\int }_0^{\mathrm{\infty }}[{\int }_0^{\mathrm{\infty }}{\lambda }^{a_2-1}exp(-\lambda (b_2+d))d\lambda ]c^{a_1-1}{e}^{b_{1}c}dc$ (46)

that is,

$p={\displaystyle \frac{b_1^{a_1}{b}_2^{a_2}}{\mathrm{\Gamma }(a_1)}{\int }_0^{\mathrm{\infty }}\frac{c^{a_1-1}{e}^{b_{1}c}}{(b_2+\frac{1}{c}(e^{c{t}_p}-1){)}^{a_2}}dc}$ (47)

Since it is not possible to obtain a closed form for the integral (47), one possibility to work around this problem is to use the Laplace approximation.

Assuming $h$ is a smooth function of an one-dimensional parameter $\varphi$ with $-h$ having a maximum at $\hat{\varphi }$, the Laplace approach asymptotically approximates an integral of the form,

$I={\int }_{-\infty }^{+\infty }exp(-nh(\varphi \left)\right)d\varphi$ (48)

by expanding $h$ in a Taylor series about $\hat{\varphi }$ (Tierney and Kadane 1986). The Laplace’s method gives the approximation

${\int }_{-\infty }^{+\infty }exp(-nh(\varphi \left)\right)d\varphi \approx \sqrt{\frac{2\pi {\sigma }^2}{n}}exp(-nh(\hat{\varphi }\left)\right)(1+O(n^{-1}\left)\right)$ (49)

where $\hat{\varphi }$ is the root of equation $h^{\left(1\right)}\left(\varphi \right)=0$ and $\hat{\sigma }=\frac{1}{\sqrt{h^{\left(2\right)}\mathrm{\text{(}}\hat{\varphi }\mathrm{\text{)}}}}$.

We can write (47) as

${\int }_0^{\mathrm{\infty }}\frac{c^{a_1-1}{e}^{b_1-1}}{(b_2+\frac{1}{c}(e^{c{t}_p}-1){)}^{a_2}}dc={\int }_0^{\mathrm{\infty }}exp(-a_{2}log(b_2+\frac{1}{c}(e^{c{t}_p}-1))+(a_1-1)log(c)-b_{1}c)dc$ (50)

Thus, the function $h\left(c\right)$ in (50) is given by

$h(c)=-a_{2}log(b_2+\frac{1}{c}(e^{c{t}_p}-1))+(a_1-1)log(c)-b_{1}c$ (51)

By applying Laplace approximation to the integral in (47) we have

$p=\frac{b_1^{a_1}{b}_2^{a_2}}{\Gamma \left(a_1\right)}\hat{\sigma }\sqrt{2\pi }exp(-h(\hat{c}\left)\right)$ (52)

where $\hat{c}$ is the root of the equation $h^{\left(1\right)}\left(c\right)=0$ and $\hat{\sigma }=\frac{1}{h^{\mathrm{\text{(}}2\mathrm{\text{)}}}\mathrm{\text{(}}\hat{c}\mathrm{\text{)}}}$ .

We suppose that an expert can summarize his/her knowledge about the reliability of an item through statements of percentiles. Thus, we ask for

expert’s information in the form of four distinct percentiles $t_p$ for given $p$ be provided to generate four equations in (52). In particular, the expert needs to specify $t_p$ for $p=0.25,0.50,0.75,0.90$.

The nonlinear system composed by the equation (52) under the four pair of values $(t_p,p)$ is solved numerically to obtain the required values of the hyperparameter a${}_1,$a${}_2,b_1$ and $b_2$ of the joint prior $\pi (c$, $\lambda |$ a${}_1,$a${}_2,b_1,b_2)$. A program has been developed in R package to solve the system.

SIMULATIONS

In this section, we perform a simulation study to examine the behavior of the proposed methods under different conditions. We considered three different sample sizes; $n=10$, 50, 100, and used several values of $(c,\lambda$). We computed MLEs of the unknown parameters of the Gompertz distribution along with the confidence intervals using the method described in Section 3. All results of the simulation study are based on 1,000 samples. The performance of the estimates is compared with respect to the average biases and the mean squared errors (MSE). To obtain the Bayes estimates and credible intervals, we need to appeal MCMC algorithm in order to obtain a sample of values of $c$ and $\lambda$ from the joint posterior distribution. To conduct the MCMC procedure, Markov chains of size 20,000 are generated from both conditional distributions $p\left(c\right|\lambda ,𝐱)$ and $p\left(c\right|\lambda ,𝐱)$ corresponding to the joint posteriors obtained under each proposed prior distribution in this paper, using MH algorithm and the first 5,000 of the observations are removed to eliminate the efect of the starting distribution. Then, in order to reduce the dependence among the generated samples, we take every 5^th sampled value which result in final chains of size 10,000, and subsequently obtained Bayes estimates based on mean of the chain, and credible intervals. The rejection rate for is around 43 and 41 over 5,000 iterations. This ensures that the choice of proposal distribution works reasonably well in sampling posterior. subsequently obtained Bayes estimates, and credible intervals.

To investigate the convergence of the MCMC sampling via MH algorithm, we have used the Gelman-Rubin multiple sequence diagnostics. For computation, we have used R package coda. For each case of “$c$” and $\lambda$, we ran two different chains with two distinct starting values of the Monte Carlo samples. Then we get two potential scale reduction factor (psrf ) values for “c” and $\lambda$. If psrf values are close to 1, we say that samples converge to the stationary distribution. For both the cases, using 10,000 samples, we get psrf value equal to 1, which suffices the convergence of the MCMC sampling procedure.

For the informative gamma priors, the elicited percentiles provided by the expert and the corresponding elicited values of the hyperparameters have been found to be: for Table I, we have $t_{25}=1.2625,t_{50}=1.5512,t_{75}=1.7806$ and $t_{90}=1.9491$ resulting $(a_1,a_2,b_1,b_2)=(38.5645,1.615,11.6416,78.2737)$, for Table II, $t_{25}=0.1083,t_{50}=0.2010,t_{75}=0.2992$ and $t_{90}=0.3820$ resulting $(a_1,a_2,b_1,b_2)=(23.3800,4.2410,23.3795,11.9681)$ and for Table III with $t_{25}=0.2272,t_{50}=0.4348,t_{75}=0.6638$ and $t_{90}=0.8618$ provides $(a_1,a_2,b_1,b_2)=(38.7505,17.2753,13.2237,13.5219)$. Frequentist property of coverage probabilities for the parameters $c$ and $\lambda$ have also been obtained to compare the Bayes estimators with different priors and MLE. Tables IVIV, VV and VI summarize the simulated coverage probabilities of 95% confidence/ credible intervals.

From the simulation results, we reach to the following conclusions: 1. With increase in sample size, biases and MSEs of the estimators decrease for given values of $n$, $c$ and $\lambda$. 2. The performance of the MLEs are quite satisfactory. The Bayes’ estimates using noninformative prior works quite well, and in most of the cases it performs better, in terms of MSE than the MLE when the value of $\lambda$ is very small. 3. Bayes estimates based on elicited prior produces much smaller bias and MSE than using the other assumed priors. 4. For the three sample sizes considered here, the elicited prior produces an over-coverage probability for small sample sizes while MLE and independent gamma priors seem to have an under-coverage for some cases. Coverage probabilities are very close to the nominal value when n increases.

TABLE I
Average bias of the estimates of $c$ and $\lambda$ and their associated MSEs (in parenthesis) for the different methods with $𝐜=3$ and $𝛌=0.02$.

Method	$c=3$			$\lambda =0.02$
	$n=10$	$n=50$	$n=100$	$n=10$	$n=50$	$n=100$
MLE	0.8476	0.3080	0.2077	0.0196	0.0088	0.0060
	(1.4217)	(0.1556)	(0.0685)	(0.0009)	(0.0001)	(6.0e-05)
Gamma prior	0.6607	0.2935	0.1997	0.0363	0.0097	0.0062
	(0.6696)	(0.1378)	(0.0649)	(0.0042)	(0.0002)	(6.9e-05)
Jeffrey’s prior	0.5204	0.2907	0.1978	0.0360	0.0096	0.0061
	(0.4204)	(0.1341)	(0.0613)	(0.0037)	(0.0001)	(6.8e-05)
MDIP	0.5639	0.2749	0.1913	0.0573	0.0114	0.0068
	(0.4732)	(0.1178)	(0.0570)	(0.0063)	(0.0002)	(8.4e-05)
Singpurwalla’s prior	0.4687	0.2898	0.1780	0.0248	0.0091	0.0058
	(0.3387)	(0.1327)	(0.0494)	(0.0015)	(0.0001)	(6.2e-05)
Elicited prior	0.2146	0.1822	0.1528	0.0040	0.0045	0.0040
	(0.0680)	(0.0538)	(0.0367)	(2.5e-05)	(3.2e-05)	(2.5e-05)

TABLE II
Average bias of the estimates of $c$ and $\lambda$ and their associated MSEs (in parenthesis) for the different methods with $𝐜=5$ and $𝛌=2$.

Method	$c=5$			$\lambda =2$
	$n=10$	$n=50$	$n=100$	$n=10$	$n=50$	$n=100$
MLE	2.8437	0.9839	0.6354	1.0200	0.4439	0.3116
	(15.2918)	(1.5513)	(0.6406)	(1.8551)	(0.3017)	(0.1536)
Gamma prior	2.5298	1.0271	0.6519	1.2956	0.5039	0.3304
	(10.2579)	(1.6401)	(0.6689)	(2.8684)	(0.4233)	(0.1806)
Jeffrey’s prior	2.2830	1.0219	0.6508	1.1710	0.5007	0.3308
	(8.5106)	(1.6199)	(0.6672)	(2.4183)	(0.4134)	(0.1805)
MDIP	1.7570	0.8304	0.6015	0.6181	0.4049	0.3135
	(6.6596)	(1.1011)	(0.5643)	(0.6339)	(0.2464)	(0.1575)
Singpurwalla’s prior	2.5279	0.9786	0.6323	0.7457	0.4248	0.3042
	(12.1530)	(1.5391)	(0.6358)	(0.8266)	(0.2741)	(0.1458)
Elicited prior	0.5518	0.4103	0.3436	0.1301	0.1467	0.1445
	(0.4509)	(0.2713)	(0.1920)	(0.0255)	(0.0326)	(0.0324)

TABLE III
Average bias of the estimates of $c$ and $\lambda$ and their associated MSEs (in parenthesis) for the different methods with $𝐜=2$ and $𝛌=1$.

Method	$c=2$			$\lambda =1$
	$n=10$	$n=50$	$n=100$	$n=10$	$n=50$	$n=100$
MLE	1.2171	0.4331	0.2967	0.4908	0.2139	0.1575
	(2.9378)	(0.3166)	(0.1414)	(0.4166)	(0.0720)	(0.0393)
Gamma prior	1.0368	0.4557	0.3030	0.5919	0.24807	0.1659
	(1.8688)	(0.3422)	(0.1449)	(0.5873)	(0.1055)	(0.0460)
Jeffrey’s prior	0.9578	0.4536	0.3020	0.5465	0.2466	0.1661
	(1.6174)	(0.3380)	(0.1443)	(0.5124)	(0.1039)	(0.0461)
MDIP	0.7924	0.3410	0.2636	0.2828	0.1757	0.1454
	(1.3917)	(0.2110)	(0.1102)	(0.1328)	(0.0482)	(0.0331)
Singpurwalla’s prior	1.0689	0.4324	0.2965	0.3593	0.2054	0.1543
	(2.2494)	(0.3157)	(0.1415)	(0.1916)	(0.0653)	(0.0374)
Elicited prior	0.2306	0.1922	0.1674	0.1071	0.0932	0.0858
	(0.0664)	(0.0539)	(0.0432)	(0.0171)	(0.0130)	(0.0112)

TABLE IV
Coverage probabilities for the parameters $𝐜$ and $𝛌$.

Method	$c=3$			$\lambda =0.02$
	$n=10$	$n=50$	$n=100$	$n=10$	$n=50$	$n=100$
MLE	0.95	0.95	0.95	0.73	0.87	0.91
Gamma prior	0.91	0.92	0.92	0.91	0.92	0.92
Jeffrey’s prior	0.97	0.96	0.96	0.97	0.96	0.96
MDIP	0.94	0.95	0.96	0.93	0.95	0.96
Singpurwalla’s prior	0.97	0.96	0.98	0.98	0.96	0.98
Elicited prior	1.00	0.98	0.98	1.00	0.99	0.98

TABLE V
Coverage probabilities for the parameters $𝐜$ and $𝛌$.

Method	$c=5$			$\lambda =2$
	$n=10$	$n=50$	$n=100$	$n=10$	$n=50$	$n=100$
MLE	0.90	0.97	0.96	0.82	0.92	0.95
Gamma prior	0.89	0.96	0.95	0.89	0.96	0.94
Jeffrey’s prior	0.95	0.95	0.95	0.95	0.95	0.95
MDIP	0.97	0.97	0.96	0.98	0.97	0.95
Singpurwalla’s prior	0.94	0.95	0.96	0.92	0.93	0.95
Elicited prior	1.00	0.99	0.99	1.00	1.00	0.99

TABLE VI
Coverage probabilities for the parameters $𝐜$ and $𝛌$.

Method	$c=2$			$\lambda =1$
	$n=10$	$n=50$	$n=100$	$n=10$	$n=50$	$n=100$
MLE	0.92	0.96	0.95	0.84	0.93	0.94
Gamma prior	0.90	0.95	0.94	0.92	0.93	0.94
Jeffrey’s prior	0.96	0.94	0.94	0.95	0.94	0.94
MDIP	0.97	0.96	0.95	0.98	0.96	0.95
Singpurwalla’s prior	0.94	0.94	0.94	0.92	0.94	0.94
Elicited prior	1.00	0.99	0.99	1.00	0.99	0.98

AN EXAMPLE WITH LITERATURE DATA

In this section, we use a real data set to illustrate the proposed estimation methods discussed in the previous sections.

Let us consider the following data set provided in King et al. (¹⁹⁷⁹):

112, 68, 84, 109, 153, 143, 60, 70, 98, 164, 63, 63, 77, 91, 91, 66, 70, 77, 63, 66, 66, 94, 101, 105, 108, 112, 115, 126, 161, 178.

These data represent the numbers of tumor-days of 30 rats fed with unsaturated diet. Chen (1997) and Asgharzadeh and Abdi (2011) used the Gompertz distribution for these data set in order to obtain exact confidence intervals and joint confidence regions for the parameters based on two different statistical analysis. Let us also assume the Gompertz distribution with density (1) fitted to the data and to compare the performance of the methods discussed in this paper.

For a Bayesian analysis, we assume independent Gamma prior distributions for the parameters $c$ and $\lambda$, with the hyper parameter values $a=b=\alpha =\beta =0.01$. The Bayes estimates cannot be obtained in closed form therefore we use MCMC procedure to compute the Bayes estimates and also to construct credible intervals. Using the software R, we simulated 50,000 MCMC samples (5,000 “burn-in-samples”) for the joint posterior distribution. The convergence of the chains was monitored from trace plots of the simulated samples. The estimates and 95% confidence intervals under classical method and Bayesian estimates with 95% credible intervals for the parameters $c$ and $\lambda$ of the Gompertz distribution are given in Table VII. The results show that among the Bayes estimators, Bayes estimate based on Singpurwalla’s prior performs the best in terms of credible intervals for both the parameters.

The marginal posterior distributions for the parameters $c$ and $\lambda$ considering the proposed prior distributions are shown in Figures 2 and 3, respectively.

TABLE VII
Estimators and 95% confidence/credible intervals of $𝐜$ and $𝛌$ of the Gompertz distribution for different estimation methods.

Method	$\hat{c}$	95% CI	$\hat{\lambda }$	95% CI
MLE	0.0241	(0.0160, 0.0322)	0.0016	(0.0002, 0.0031)
Gamma Prior	0.0232	(0.0150, 0.0312)	0.0019	(0.0007, 0.0038)
Jeffrey’s prior	0.0234	(0.0152, 0.0318)	0.0018	(0.0007, 0.0038)
MDIP	0.0226	(0.0147, 0.0304)	0.0020	(0.0008, 0.0045)
Singpurwalla’s Prior	0.0242	(0.0165, 0.0317)	0.0016	(0.0006, 0.0033)

Figure 2
The posterior densities for the parameter c of the Gompertz distribution fitted by the data.

Figure 3
The posterior densities for the parameter λ of the Gompertz distribution fitted by the data.

CONCLUSIONS

In this paper, we have considered estimation of the parameters of the Gompertz distribution using frequentist and Bayesian methods. In Bayesian methods, we have consider objective priors (Jeffreys and MDIP), gamma prior, Singpurwalla’s prior and Elicited prior. We have performed an extensive simulation study to compare these methods. From the simulation study regarding the bias, MSE and CP we observe that in general the MDIP provides best results for both parameters and in some cases, with MDIP and Jeffreys priors the results are quite similar. The real data application shows the same situation. It is worth remembering that both forms result from formal procedures for representing absence of information, that is, they are noninformative. The commonly assumption used of independent gamma priors and the priors proposed by Singpurwalla do not present as good results as the objective priors. The independent gamma priors are generally used in situations where no objective priors are possible to obtain or provide improper posterior distributions and mainly due to computational ease. Elicited prior produces much smaller bias and MSE than using the other assumed priors and also provides an over-coverage probability than their counterparts. Hence, we can conclude that, in the situation of the absence of information, the MDIP prior is more indicate for a Bayesian estimation of the two-parameter Gompertz distribution. On the other hand, in the situation where we have available expert’s information, the Elicited prior will perform the best estimators.

REFERENCES

• ¹
  ASGHARZADEH A AND ABDI M. 2011. Exact Confidence Intervals and Joint Confidence Regions for the Parameters of the Gompertz Distribution based on Records. Pak J Statist 271: 55-64.
• ²
  BERNARDO JM. 1979. Reference Posterior Distributions for Bayesian Inference. J R Stat Soc Series B Stat Methodol 41: 113-147.
• ³
  CASELLA G AND BERGER RL. 2001. Statistical Inference. Cengage Learning, 660 p.
• ⁴
  CHEN Z. 1997. Parameter estimation of the Gompertz population. Biom J 39: 117-124.
• ⁵
  ELANDT-JOHNSON RC AND JOHNSON NL. 1979. Survival Models and Data Analysis. J Wiley & Sons: NY, 457 p.
• ⁶
  FRANSES PH. 1994. Fitting a Gompertz curve. J Oper Res Soc 45: 109-113.
• ⁷
  GARG M, RAO B AND REDMOND C. 1970. Maximum-likelihood estimation of the parameters of the Gompertz survival function. J R Stat Soc Ser C Appl Stat 19: 152-159.
• ⁸
  GOMPERTZ B. 1825. On the nature of the function expressive of the law of human mortality and on a new mode of determining the value of life contingencies. Philos Trans R Soc Lond 115: 513-583.
• ⁹
  ISMAIL AA. 2010. Bayes estimation of Gompertz distribution parameters and acceleration factor under partially accelerated life tests with type-I censoring. J Stat Comput Simul 80: 1253-1264.
• ¹⁰
  ISMAIL AA. 2011. Planning step-stress life tests with type-II censored Data. Sci Res Essays 6: 4021-4028.
• ¹¹
  JAHEEN ZF. 2003a. Prediction of Progressive Censored Data from the Gompertz Model. Commun Stat Simul Comput 32: 663-676.
• ¹²
  JAHEEN ZF. 2003b. A Bayesian analysis of record statistics from the Gompertz model. Appl Math Comput 145: 307-320.
• ¹³
  JEFFREYS SIR HAROLD. 1967. Theory of probability. Oxford U Press: London, 470 p.
• ¹⁴
  KIANI K, ARASAN J AND MIDI H. 2012. Interval estimations for parameters of Gompertz model with time-dependent covariate and right censored data. Sains Malays 414: 471-480.
• ¹⁵
  KING M, BAILEY DM, GIBSON DG, PITHA JV AND MCCAY PB. 1979. Incidence and growth of mammary tumors induced by 7, l2-Dimethylbenz antheacene as related to the dietary content of fat and antioxidant. J Natl Cancer Inst 63: 656-664.
• ¹⁶
  LAWLESS JF. 2003. Statistical Models and Methods for Lifetime Data. J Wiley & Sons: New Jersey, 664 p.
• ¹⁷
  LENART A. 2014. The moments of the Gompertz distribution and maximum likelihood estimation of its parameters. Scand Actuar J 3: 255-277.
• ¹⁸
  LENART A AND MISSOV TI. 2016. Goodness-of-fit tests for the Gompertz distribution. Commun Stat Theory Methods 45: 2920-2937.
• ¹⁹
  MAKANY RA. 1991. Theoretical basis of Gompertz’s curve. Biom J 33: 121-128.
• ²⁰
  MEIJER CS. 1936. Über Whittakersche bzw. Besselsche Funktionen und deren Produkte. Nieuw Archief voor Wiskunde 18: 10-39.
• ²¹
  MILGRAM M. 1985. The generalized integro-exponential function. Math Comp 44: 443-458.
• ²²
  RAO BR AND DAMARAJU CV. 1992. New better than used and other concepts for a class of life distribution. Biom J 34: 919-935.
• ²³
  READ CB. 1983. Gompertz Distribution. Encyclopedia of Statistical Sciences. J Wiley & Sons, NY.
• ²⁴
  SHANUBHOGUE A AND JAIN NR. 2013. Minimum Variance Unbiased Estimation in the Gompertz Distribution under Progressive Type II Censored Data with Binomial Removals. Int Sch Res Notices 2013: 1-7.
• ²⁵
  SINGH N, YADAV KK AND RAJASEKHARAN R. 2016. ZAP1-mediated modulation of triacylglycerol levels in yeast by transcriptional control of mitochondrial fatty acid biosynthesis. Mol Microbiol 1001: 55-75.
• ²⁶
  SINGPURWALLA ND. 1988. An interactive PC-based procedure for reliability assessment incorporating expert opinion and survival data. J Am Stat Assoc 83: 43-51.
• ²⁷
  TIBSHIRANI R. 1989. Noninformative Priors for One Parameters of Many. Biometrika 76: 604-608.
• ²⁸
  WILLEKENS F. 2002. Gompertz in context: the Gompertz and related distributions. In: Tabeau E, Van den Berg JA and Heathcote C (Eds), Forecasting mortality in developed countries - insights from a statistical demographic and epidemiological perspective European studies of population. Dordrecht: Kluwer Academic Publishers 9: 105-126.
• ²⁹
  WU JW, HUNG WL AND TSAI CH. 2004. Estimation of parameters of the Gompertz distribution using the least squares method. Appl Math Comput 158: 133-147.
• ³⁰
  WU JW AND LEE WC. 1999. Characterization of the mixtures of Gompertz distributions by conditional expectation of order statistics. Biom J 41: 371-381.
• ³¹
  WU JW AND TSENG HC. 2006. Statistical inference about the shape parameter of the Weibull distribution by upper record values. Stat Pap 48: 95-129.
• ³²
  WU SJ, CHANG CT AND TSAI TR. 2003. Point and interval estimations for the Gompertz distribution under progressive type-II censoring. Metron LXI: 403-418.
• ³³
  ZELLNER A. 1977. Maximal Data Information Prior Distributions. In: Aykac A and Brumat C (Eds), New Methods in the applications of Bayesian Methods. North-Holland Amsterdam, p. 211-232.
• ³⁴
  ZELLNER A. 1984. Maximal Data Information Prior Distributions. Basic Issues in Econometrics, U. of Chicago Press.
• ³⁵
  ZELLNER A. 1990. Bayesian Methods and Entropy in Economics and Econometrics. In: Grandy Junior WT and Schick LH (Eds), Maximum Entropy and Bayesian Methods Dordrecht Netherlands: Kluwer Academic Publishers, p. 17-31.
• ³⁶
  ZELLNER A AND MIN C. 1993. Bayesian analysis model selection and prediction. In: Grandy Junior WT and Milonni PW (Eds), Physics and Probability: Essays in honour of Edwin T Jaynes, Cambridge University Press, Cambridge, UK.

Publication Dates

• Publication in this collection
  Sept 2018

History

• Received
  29 Dec 2017
• Accepted
  26 Mar 2018