Bayesian Estimation of a New Pareto-Type Distribution Based on Mixed Gibbs Sampling Algorithm

Abstract

In this paper, based on the mixed Gibbs sampling algorithm, a Bayesian estimation procedure is proposed for a new Pareto-type distribution in the case of complete and type II censored samples. Simulation studies show that the proposed method is consistently superior to the maximize likelihood estimation in the context of small samples. Also, an analysis of some real data is provided to test the Bayesian estimation.

1. Introduction

As an important statistical inference method, the Bayesian method has been widely applied and developed in the fields of statistics and artificial intelligence. Especially with the rise of big data and machine learning in the early 21st century, the Bayesian method has received more and more attention. Researchers have improved the efficiency and accuracy of the Bayesian method by developing new algorithms and tools and applying them to various fields, such as computer vision, natural language processing, bioinformatics, and so on. Recently, many scholars [1,2,3] have studied it and obtained some valuable results.

The Pareto distribution [4] is a failure function with a decreasing value and has been widely applied in reliability, personal income, stock price fluctuation, and some other models. The origin and some applications of the Pareto distribution can be seen in Johnson et al. [5]. Bourguignon et al. [6] proposed a new Pareto-type (NP) distribution, which is a generalization of the Pareto distribution and often applied for income and reliability data analysis. Raof et al. [7] used the NP distribution model to simulate the income of upper class groups in Malaysia. Karakaya et al. [8] proposed a new generalization of the Pareto distribution, incorporating a truncation parameter. Sarabia et al. [9] acquired further significant properties of this distribution and derived simpler expressions for relevant inequality indices and risk measures. Nik et al. [10] estimated the parameters of the NP distribution from two different perspectives, namely, classical and Bayesian statistics, in which the Bayesian estimation used importance sampling method. Nik et al. [11] studied the parameter estimation as well as predicting the failure times of the removed units in multiple stages of the progressively censored sample coming from the NP distribution, in which a Bayesian method (i.e., Markov chain Monte Carlo, MCMC) was applied to estimate the unknown parameters involved in the model. Soliman [12] studied the estimation of the reliability function in a generalized life model. The Bayesian estimations of the symmetric loss function (i.e., secondary loss) and asymmetric loss function (i.e., LINEX loss and GE loss) were obtained. These estimations were compared with the maximum likelihood estimation (MLE) with the Burr-XII model using the Bayesian approximation based on Lindley. Furthermore, Soliman [12] studied the Bayesian estimation of the parameter of interest on the premise that some other parameters are known. However, all of the parameters are often unknown in practice. Soliman [13] gave an approximate calculation formula for Bayesian estimation; however, the derivation of the formula is complicated and difficult to understand. Soliman’s work shows that the Bayesian method may encounter high dimensions, complex models, and big data in practical use. Thus, scholars can greatly improve the efficiency and accuracy of the Bayesian method by improving the MCMC algorithm.

The Gibbs sampling algorithm was originally proposed by Geman [14]. In recent years, with the development of technologies such as big data analysis, artificial intelligence, and machine learning, the applications of the Gibbs sampling algorithm to solve Bayesian estimation problems has become increasingly common. Ahmed [15] proposed two prediction methods for predicting the censored units with progressive type II censored samples. The lifetimes under consideration are set to follow a new two-parameter Pareto distribution. Moreover, point and interval estimation of the unknown parameters of the NP distribution can be obtained using the maximum likelihood and Bayesian estimation methods. Since Bayesian estimators cannot be expressed explicitly, the Gibbs sampling and MCMC techniques are utilized for Bayesian calculation. Ibrahim et al. [16] proposed a new discrete analogue of the Weibull class, and gave a Bayesian estimation procedure under the square error loss function. Furthermore, they compared the non-Bayesian estimation with the Bayesian estimation using Gibbs sampling and the Metropolis–Hastings algorithm.

For the NP distribution, Nik et al. [10] proposed some classical estimation methods and a Bayesian estimation based on the importance sampling method. However, the sampling efficiency of the importance sampling method is often low in high-dimensional situations. Thus, some other efficient methods are preferred. In our work, we use the mixed Gibbs sampling algorithm to estimate the parameters of the NP distribution with complete samples and type II censored samples.

The Gibbs sampling algorithm can solve some complex Bayes analysis problems directly, and has become a useful tool in Bayesian statistics. By establishing a Markov chain and iterating to the equilibrium state repeatedly, the Gibbs sampling algorithm can obtain a posterior distribution sample. However, the Gibbs sampling algorithm also has some shortcomings. For instance, the fully conditional posterior distribution of a parameter may be difficult to be directly sampled for some special prior distribution. Thus, the Gibbs sampling algorithm can be combined with other sampling methods to achieve the stability of the efficiency. Nowadays, there are some sampling methods, such as the Metropolis algorithm, importance sampling, accept–reject sampling, and so on. Among them, the Metropolis algorithm is a relatively easy and flexible sampling method. The Metropolis algorithm was improved by Hastings [17], and the result is called the Metropolis–Hastings algorithm. However, the Metropolis–Hastings algorithm is time consuming and inefficient in high-dimensional cases. Sampling may even be rejected after time-consuming calculations. A mixture of the Gibbs sampling algorithm and Metropolis–Hastings algorithm, i.e., the mixed Gibbs sampling algorithm, can overcome these shortcomings of the Metropolis–Hastings algorithm. Furthermore, the mixed Gibbs sampling algorithm can make the Gibbs sampling algorithm and Metropolis–Hastings algorithm perform fully to their own advantages. Thus, in our work, we use the mixed Gibbs sampling algorithm to propose a Bayesian estimation procedure for the NP distribution with complete and type II censored samples.

The rest of this paper is organized as follows. We state the MLE of the NP distribution in Section 2. In Section 3, we propose a Bayesian estimation procedure for the NP distribution with complete and type II censored samples using the mixed Gibbs sampling algorithm. Simulation studies and real data analysis are presented in Section 4. In Section 5, we present a brief discussion and conclusion of the results and methods.

2. Maximum Likelihood Estimation of NP Distribution

Suppose X is a random variable that follows the NP distribution with a shape parameter $\alpha$ and a scale parameter $\beta$, that is, $X\sim \mathrm{NP}(\alpha ,\beta )$. The cumulative distribution function (CDF) $F\left(x\right|\alpha ,\beta )$ of X is given as

$$F\left(x\right|\alpha ,\beta )=\frac{1-{(\beta /x)}^{\alpha }}{1+{(\beta /x)}^{\alpha }}=\frac{x^{\alpha }-{\beta }^{\alpha }}{x^{\alpha }+{\beta }^{\alpha }}=1-\frac{2{(\beta /x)}^{\alpha }}{1+{(\beta /x)}^{\alpha }}=1-\frac{2{\beta }^{\alpha }}{x^{\alpha }+{\beta }^{\alpha }},x\ge \beta .$$

(1)

The probability density function (PDF) $f\left(x\right|\alpha ,\beta )$ of X is

$$f\left(x\right|\alpha ,\beta )=\frac{2\alpha {(\beta /x)}^{\alpha }}{x{[1+{(\beta /x)}^{\alpha }]}^2}=\frac{2\alpha {(\beta /x)}^{\alpha +1}}{\beta {[1+{(\beta /x)}^{\alpha }]}^2}=\frac{2\alpha {\beta }^{\alpha }{x}^{\alpha -1}}{{(x^{\alpha }+{\beta }^{\alpha })}^2},x\ge \beta .$$

(2)

The PDF and CDF diagrams are shown in Figure 1.

It should be highlighted that compared with the Pareto distribution, the NP distribution has the following advantages: (1) It can exhibit either an upside-down bathtub or a decreasing hazard rate function, depending on the values of its parameters. (2) It offers mathematical simplicity, as the probability and distribution functions of the new distribution have a simple form, in contrast to some generalizations of the Pareto distribution which involve special functions such as log, beta, and gamma functions. (3) The proposed distribution has only two parameters, unlike some generalizations of the Pareto distribution which have three or four parameters. Thus, in our work, it is assumed that the life of the product X follows the NP distribution. We present an MLE of the NP distribution below for comparison with the Bayesian estimation.

2.1. Complete Samples Case

In the case of complete samples, assuming that there are n products for a life test, and the order failure data $x_1\le x_2\le \cdots \le x_n$ can be obtained, the log-likelihood function $\ell (\alpha ,\beta )$ can be written as

$$\ell (\alpha ,\beta )=nlog\left(2\alpha \right)-nlog\beta +(\alpha +1)\sum _{i=1}^{n}log\left(\frac{\beta }{x_i}\right)-2\sum _{i=1}^{n}log\left(1+{\left(\frac{\beta }{x_i}\right)}^{\alpha }\right).$$

It can be seen that $\ell (\alpha ,\beta )$ increases monotonically with $\beta$, that is to say, the greater the value of $\beta$, the greater the value of $\ell (\alpha ,\beta )$. Since $x\ge \beta$, we obtain the MLE $\widehat{\beta }=x_1$. Thus, the MLE $\widehat{\alpha }$ of $\alpha$ can be obtained from the solution of Equation (3):

$$\begin{array}{c}\hfill {\displaystyle \frac{\partial \ell (\alpha ,x_1)}{\partial \alpha }}=\frac{n}{\alpha }+\sum _{i=2}^{n}log\left(\frac{x_1}{x_i}\right)-2\sum _{i=2}^n\frac{{\left(\frac{x_1}{x_i}\right)}^{\alpha }log\left(\frac{x_1}{x_i}\right)}{1+{\left(\frac{x_1}{x_i}\right)}^{\alpha }}=0.\end{array}$$

(3)

2.2. Type II Censored Samples Case

In the case of type II censored samples, assume that there are n products for a life test, then the test stops when there are r product failures, and the order failure data $x_1\le x_2\le \cdots \le x_r$ can be obtained. MLE is used to estimate the unknown parameters $\alpha ,\beta$ below. In order to obtain the likelihood function, it is necessary to know the probability of occurrence as follows.

(1) The probability of a product failing at $(x_i,x_i+d{x}_i]$ is approximately

$$\begin{array}{c}\hfill f\left(x_i\right)d{x}_i={\displaystyle \frac{2\alpha {\beta }^{\alpha }{x}_i^{\alpha -1}}{{(x_i^{\alpha }+{\beta }^{\alpha })}^2}}d{x}_i,i=1,2,\cdots ,r.\end{array}$$

(2) The probability that the lifetime of the remaining $n-r$ products exceed $x_r$ is

$${\left({\int }_{x_r}^{\infty }{\displaystyle \frac{2\alpha {\beta }^{\alpha }{x}_i^{\alpha -1}}{{(x_i^{\alpha }+{\beta }^{\alpha })}^2}}dx\right)}^{n-r}={\left(\left(1-\frac{2{\beta }^{\alpha }}{x_i^{\alpha }+{\beta }^{\alpha }}\right){|}_{x_r}^{\infty }\right)}^{n-r}={\left(\frac{2{\beta }^{\alpha }}{x_r^{\alpha }+{\beta }^{\alpha }}\right)}^{n-r}.$$

Thus, the probability of the above observation $\{x_i,i=1,2,\cdots ,r\}$ is approximately

$$\begin{array}{cc}\hfill & \left(\genfrac{}{}{0pt}{}{n}{r}\right)\left({\displaystyle \frac{2\alpha {\beta }^{\alpha }{x}_1^{\alpha -1}}{{(x_1^{\alpha }+{\beta }^{\alpha })}^2}}d{x}_1\right)\left({\displaystyle \frac{2\alpha {\beta }^{\alpha }{x}_2^{\alpha -1}}{{(x_2^{\alpha }+{\beta }^{\alpha })}^2}}d{x}_2\right)\cdots \left({\displaystyle \frac{2\alpha {\beta }^{\alpha }{x}_r^{\alpha -1}}{{(x_r^{\alpha }+{\beta }^{\alpha })}^2}}d{x}_r\right){\left(\frac{2{\beta }^{\alpha }}{x_r^{\alpha }+{\beta }^{\alpha }}\right)}^{n-r}\hfill \\ \hfill & ={\alpha }^r{\left(2{\beta }^{\alpha }\right)}^n\left(\genfrac{}{}{0pt}{}{n}{r}\right)\frac{{\prod }_{i=1}^r{x}_i^{\alpha -1}}{{\prod }_{i=1}^r{(x_i^{\alpha }+{\beta }^{\alpha })}^2}{\left(\frac{1}{x_r^{\alpha }+{\beta }^{\alpha }}\right)}^{n-r}d{x}_{1}d{x}_2\cdots d{x}_r,\hfill \end{array}$$

where $d{x}_{1}d{x}_2\cdots d{x}_r$ is a constant. Ignoring a constant factor does not affect the MLE of $\alpha$ and $\beta$, and the likelihood function $L(\alpha ,\beta )$ can be given as

$$\begin{array}{c}\hfill L(\alpha ,\beta )={\alpha }^r{\left(2{\beta }^{\alpha }\right)}^n\frac{{\prod }_{i=1}^r{x}_i^{\alpha -1}}{{\prod }_{i=1}^r{(x_i^{\alpha }+{\beta }^{\alpha })}^2}{\left(\frac{1}{x_r^{\alpha }+{\beta }^{\alpha }}\right)}^{n-r}.\end{array}$$

Thus, the log-likelihood function $\ell (\alpha ,\beta )$ is as

$$\begin{array}{cc}\hfill \ell (\alpha ,\beta )=& rlog\alpha +nlog2+n\alpha log\beta +(\alpha -1)\sum _{i=1}^{r}log{x}_i-2\sum _{i=1}^{r}log(x_i^{\alpha }+{\beta }^{\alpha })\hfill \\ \hfill & -(n-r)log(x_r^{\alpha }+{\beta }^{\alpha }).\hfill \end{array}$$

Similarly, we can obtain the MLE $\widehat{\beta }$ of $\beta$ as

$$\begin{array}{c}\hfill \widehat{\beta }=x_1.\end{array}$$

(4)

Therefore, the MLE $\widehat{\alpha }$ of $\alpha$ can be obtained from the solution of Equation (5)

$$\begin{array}{cc}\hfill \frac{\partial \ell (\alpha ,x_1)}{\partial \alpha }=& \frac{r}{\alpha }+nlog{x}_1+\sum _{i=1}^{r}log{x}_i-2\sum _{i=1}^r\frac{x_i^{\alpha }log{x}_i+x_1^{\alpha }log{x}_1}{x_i^{\alpha }+x_1^{\alpha }}\hfill \\ \hfill & -(n-r)\frac{x_r^{\alpha }log{x}_r+x_1^{\alpha }log{x}_1}{x_r^{\alpha }+x_1^{\alpha }}=0.\hfill \end{array}$$

(5)

3. Bayesian Estimation of NP Distribution

The Bayesian estimations of the NP distribution under complete samples and type II censored samples are given below.

3.1. Complete Samples Case

In the case of complete samples, assume that there are n products for a life test, and the order failure data $x_1\le x_2\le \cdots \le x_n$ can be obtained. When the parameters $\alpha$ and $\beta$ are unknown, we assume that the joint prior distribution $\pi (\alpha ,\beta )$ of $\alpha$ and $\beta$ is

$$\pi (\alpha ,\beta )\propto 1/\left(\alpha \beta \right).$$

(6)

Then, the posterior distribution $\pi (\alpha ,\beta |x_1,x_2,\cdots ,x_n)$ is given as (7):

$$\begin{array}{cc}\hfill \pi (\alpha ,\beta |x_1,x_2,\cdots ,x_n)& \propto {\displaystyle \frac{{\left(2\alpha \right)}^n{\beta }^{n\alpha }{\prod }_{i=1}^n{x}_i^{\alpha -1}}{{\prod }_{i=1}^n{(x_i^{\alpha }+{\beta }^{\alpha })}^2}}·\frac{1}{\alpha \beta }\hfill \\ \hfill & =\frac{2^n{\alpha }^{n-1}{\beta }^{n\alpha -1}{\prod }_{i=1}^n{x}_i^{\alpha -1}}{{\prod }_{i=1}^n{(x_i^{\alpha }+{\beta }^{\alpha })}^2}.\hfill \end{array}$$

(7)

Thus, the fully conditional posterior distribution $\pi \left(\alpha \right|x_1,x_2,\cdots ,x_n,\beta )$ of $\alpha$ is

$$\pi \left(\alpha \right|x_1,x_2,\cdots ,x_n,\beta )\propto \frac{{\alpha }^{n-1}{\beta }^{n\alpha -1}{\prod }_{i=1}^n{x}_i^{\alpha -1}}{{\prod }_{i=1}^n{(x_i^{\alpha }+{\beta }^{\alpha })}^2},$$

(8)

and the fully conditional posterior distribution $\pi \left(\beta \right|x_1,x_2,\cdots ,x_n,\alpha )$ of $\beta$ is

$$\pi \left(\beta \right|x_1,x_2,\cdots ,x_n,\alpha )\propto \frac{{\beta }^{n\alpha -1}}{{\prod }_{i=1}^n{(x_i^{\alpha }+{\beta }^{\alpha })}^2}.$$

(9)

According to (8) and (9), both $\pi \left(\alpha \right|x_1,x_2,\cdots ,x_r,\beta )$ and $\pi \left(\beta \right|x_1,x_2,\cdots ,x_r,\alpha )$ do not have explicit expressions in the case of complete samples. Therefore, it is difficult to conduct sampling directly in practice. However, the mixed Gibbs sampling algorithm does not require an explicit expression of the fully conditional posterior distribution, and can be used to sampling different prior distributions. Thus, in our work, we consider the case of the prior distribution denoted as (6) using the mixed Gibbs sampling algorithm in Section 3.1. For Bayesian estimations for parameters of the NP distribution in the case of complete samples, we state the iterative steps of the mixed Gibbs sampling algorithm as follows.

Step 0 

Choose an initial value $({\alpha }^{\left(0\right)},{\beta }^{\left(0\right)})$ of $(\alpha ,\beta )$. Denote the value of the ith iteration as $({\alpha }^{\left(i\right)},{\beta }^{\left(i\right)})$. Then, the value of the $(i+1)$th iteration $({\alpha }^{(i+1)},{\beta }^{(i+1)})$ can be obtained using the following Step 1 and Step 2.

Step 1 

Obtain ${\alpha }^{(i+1)}$ from $\pi \left(\alpha \right|x_1,x_2,\cdots ,x_n,{\beta }^{\left(i\right)})$ using the Metropolis–Hastings algorithm in this step, which consists of the following four steps (1)∼(4).

(1)

Set $\alpha \sim \mathrm{N}({\alpha }_k,{\sigma }_{\alpha }^2)$, where ${\alpha }_k$ is the current state and ${\sigma }_{\alpha }$ is the standard deviation. Obtain a sample ${\alpha }^{\prime }$ of $\alpha$; if ${\alpha }^{\prime }\le 0$, then a new sampling is required.

(2)

Calculate the acceptance probability ${\mathrm{P}}_{\mathrm{accept}}({\alpha }_k,{\alpha }^{\prime })$ as (10) below:

$$\begin{array}{cc}\hfill {\mathrm{P}}_{\mathrm{accept}}& ({\alpha }_k,{\alpha }^{\prime })\hfill \\ & =min\left\{1,\frac{\pi \left({\alpha }^{\prime }\right|x_1,x_2,\cdots ,x_n,{\beta }^{\left(i\right)})}{\pi \left({\alpha }_k\right|x_1,x_2,\cdots ,x_n,{\beta }^{\left(i\right)})}\right\}\hfill \\ & =min\left\{1,\frac{{{\alpha }^{\prime }}^{n-1}{\beta }^{n{\alpha }^{\prime }-1}{\prod }_{i=1}^n{x}_i^{{\alpha }^{\prime }-1}}{{\alpha }_k^{n-1}{\beta }^{n{\alpha }_k-1}{\prod }_{i=1}^n{x}_i^{{\alpha }_k-1}}·\frac{{\prod }_{i=1}^n{(x_i^{{\alpha }_k}+{\beta }^{{\alpha }_k})}^2}{{\prod }_{i=1}^n{(x_i^{{\alpha }^{\prime }}+{\beta }^{{\alpha }^{\prime }})}^2}\right\}\hfill \\ & =min\left\{1,{\left(\frac{{\alpha }^{\prime }}{{\alpha }_k}\right)}^{n-1}·{\beta }^{n·\Delta {\alpha }_k}·{\left(\prod _{i=1}^n{x}_i\right)}^{\Delta {\alpha }_k}·\frac{{\prod }_{i=1}^n{(x_i^{{\alpha }_k}+{\beta }^{{\alpha }_k})}^2}{{\prod }_{i=1}^n{(x_i^{{\alpha }^{\prime }}+{\beta }^{{\alpha }^{\prime }})}^2}\right\},\hfill \end{array}$$

(10)

where $\Delta {\alpha }_k={\alpha }^{\prime }-{\alpha }_k$.

(3)

Generate a random number u from the uniform distribution $U(0,1)$, and then obtain ${\alpha }_{k+1}$ according to the following conditions (11):

$${\alpha }_{k+1}=\left\{\begin{array}{c}{\alpha }^{\prime },u\le {\mathrm{P}}_{\mathrm{accept}}({\alpha }_k,{\alpha }^{\prime })\hfill \\ {\alpha }_k,u>{\mathrm{P}}_{\mathrm{accept}}({\alpha }_k,{\alpha }^{\prime })\hfill \end{array}.\right.$$

(11)

(4)

Set $k=k+1$, then return to Step 1 (1).

The resulting $\{{\alpha }^{\left(i\right)},i=1,2,\cdots ,k\}$ is a Markov chain. Step 1 ends when the Markov chain reaches equilibrium, then we can obtain ${\alpha }^{(i+1)}$.

Step 2 

Obtain ${\beta }^{(i+1)}$ from $\pi \left(\beta \right|x_1,x_2,\cdots ,x_n,{\alpha }^{(i+1)})$ using the Metropolis–Hastings algorithm in this step.

(1)

Set $\beta \sim \mathrm{N}({\beta }_k,{\sigma }_{\beta }^2)$, where ${\beta }_k$ is the current state and ${\sigma }_{\beta }$ is the standard deviation. Obtain a sample ${\beta }^{\prime }$ of $\beta$; if ${\beta }^{\prime }\le 0$ or ${\beta }^{\prime }>x_1$, then a new sampling is required.

(2)

Calculate the acceptance probability ${\mathrm{P}}_{\mathrm{accept}}({\beta }_k,{\beta }^{\prime })$ as (12):

$$\begin{array}{cc}\hfill {\mathrm{P}}_{\mathrm{accept}}({\beta }_k,{\beta }^{\prime })& =min\left\{1,\frac{\pi \left({\beta }^{\prime }\right|x_1,x_2,\cdots ,x_n,{\alpha }^{(i+1)})}{\pi \left({\beta }_k\right|x_1,x_2,\cdots ,x_n,{\alpha }^{i+1})}\right\}\hfill \\ \hfill & =min\left\{1,{\left(\frac{{\beta }^{\prime }}{{\beta }_k}\right)}^{n{\alpha }^{(i+1)}-1}·\frac{{\prod }_{i=1}^n{(x_i^{{\alpha }^{(i+1)}}+{\beta }_k^{{\alpha }^{(i+1)}})}^2}{{\prod }_{i=1}^n{(x_i^{{\alpha }^{(i+1)}}+{{\beta }^{\prime }}^{{\alpha }^{(i+1)}})}^2}\right\}.\hfill \end{array}$$

(12)

(3)

Generate a random number u from the uniform distribution of $U(0,1)$, and then obtain ${\beta }_{k+1}$ according to the following conditions (13):

$${\beta }_{k+1}=\left\{\begin{array}{c}{\beta }^{\prime },u\le {\mathrm{P}}_{\mathrm{accept}}({\beta }_k,{\beta }^{\prime })\hfill \\ {\beta }_k,u>{\mathrm{P}}_{\mathrm{accept}}({\beta }_k,{\beta }^{\prime })\hfill \end{array}.\right.$$

(13)

(4)

Set $k=k+1$, then return to Step 2 (1).

Similarly, the resulting $\{{\beta }^{\left(i\right)},i=1,2,\cdots ,k\}$ is a Markov chain. We can obtain ${\beta }^{(i+1)}$ when the Markov chain reaches equilibrium.

Step 3 

Let $i=i+1$, and repeat Step 1 and Step 2 in sequence until the Markov chains reach equilibrium.

It is worth noting that the sample trajectory of the Markov chain can be used as samples of the posterior distributions of $\alpha ,\beta$ when the Markov chain reaches equilibrium, further leading to the Bayesian estimation of $\alpha ,\beta$. In order to ensure the approximate independence of the posterior distribution samples, subsequent determinations will be made based on the sample autocorrelation coefficients.

3.2. Type II Censored Samples Case

In the case of type II censored samples, assuming that there are n products for a life test, the test stops when there are r products failure, and the order failure data $x_1\le x_2\le \cdots \le x_r$ can be obtained. We still take the joint prior distribution of $\alpha ,\beta$ as $\pi (\alpha ,\beta )\propto 1/\left(\alpha \beta \right)$, then the posterior distribution $\pi (\alpha ,\beta |x_1,x_2,\cdots ,x_r)$ is the following (14):

$$\begin{array}{cc}\hfill \pi (\alpha ,\beta |x_1,x_2,\cdots ,x_r)& \propto \prod _{i=1}^{r}f\left(x_i\right|\alpha ,\beta )·\prod _{i=r+1}^n(1-F\left(x_r\right|\alpha ,\beta ))·\frac{1}{\alpha \beta }\hfill \\ & =\frac{{\left(2\alpha \right)}^r{\beta }^{r\alpha }{\prod }_{i=1}^r{x}_i^{\alpha -1}}{{\prod }_{i=1}^r{(x_i^{\alpha }+{\beta }^{\alpha })}^2}·\frac{{\left(2{\beta }^{\alpha }\right)}^{n-r}}{{(x_r^{\alpha }+{\beta }^{\alpha })}^{n-r}}·\frac{1}{\alpha \beta }\hfill \\ & =\frac{2^n{\alpha }^{r-1}{\beta }^{n\alpha -1}{\prod }_{i=1}^r{x}_i^{\alpha -1}}{{\prod }_{i=1}^r{(x_i^{\alpha }+{\beta }^{\alpha })}^2{(x_r^{\alpha }+{\beta }^{\alpha })}^{n-r}}.\hfill \end{array}$$

(14)

Thus, the fully conditional posterior distribution $\pi \left(\alpha \right|x_1,x_2,\cdots ,x_r,\beta )$ of $\alpha$, denoted as (15), is

$$\pi \left(\alpha \right|x_1,x_2,\cdots ,x_r,\beta )\propto \frac{{\alpha }^{r-1}{\beta }^{n\alpha -1}{\prod }_{i=1}^r{x}_i^{\alpha -1}}{{\prod }_{i=1}^r{(x_i^{\alpha }+{\beta }^{\alpha })}^2{(x_r^{\alpha }+{\beta }^{\alpha })}^{n-r}},$$

(15)

and the fully conditional posterior distribution $\pi \left(\beta \right|x_1,x_2,\cdots ,x_r,\alpha )$ of $\beta$, denoted as (16), is

$$\pi \left(\beta \right|x_1,x_2,\cdots ,x_r,\alpha )\propto \frac{{\beta }^{n\alpha -1}}{{\prod }_{i=1}^r{(x_i^{\alpha }+{\beta }^{\alpha })}^2{(x_r^{\alpha }+{\beta }^{\alpha })}^{n-r}}.$$

(16)

As can be seen from the above, in the case of type II censored samples, neither (15) nor (16) has an explicit expression, and so they are difficult to sample directly. Similarly, as previously discussed, the mixed Gibbs sampling algorithm does not require a display expression for the complete conditional posterior distribution, so we continue to use the mixed Gibbs algorithm. Similarly, we still only consider the same prior distribution as in Section 3.1. The steps of the mixed Gibbs sampling algorithm for Bayesian estimation of the NP distribution parameters in the case of type II censored samples are stated as follows.

Step 0 

Give an initial value $({\alpha }^{\left(0\right)},{\beta }^{\left(0\right)})$ using the Gibbs sampling. For simplicity, the value of the ith iteration is still denoted as $({\alpha }^{\left(i\right)},{\beta }^{\left(i\right)})$. Then, we can use the following Step 1 and Step 2 to obtain $({\alpha }^{(i+1)},{\beta }^{(i+1)})$.

Step 1 

Obtain a sample ${\alpha }^{(i+1)}$ from $\pi \left(\alpha \right|x_1,x_2,\cdots ,x_r,{\beta }^{\left(i\right)})$ using the Metropolis–Hastings algorithm in this step.

(1)

Set $\alpha \sim \mathrm{N}({\alpha }_k,{\sigma }_{\alpha }^2)$, where ${\alpha }_k$ is the current state and ${\sigma }_{\alpha }$ is the standard deviation. Obtain a sample ${\alpha }^{\prime }$ of $\alpha$; if ${\alpha }^{\prime }\le 0$, then a new sampling is required.

(2)

Calculate the acceptance probability ${\mathrm{P}}_{\mathrm{accept}}({\alpha }_k,{\alpha }^{\prime })$, denoted as (17):

$$\begin{array}{cc}& {\mathrm{P}}_{\mathrm{accept}}({\alpha }_k,{\alpha }^{\prime })\hfill \\ & =min\left\{1,\frac{\pi \left({\alpha }^{\prime }\right|x_1,x_2,\cdots ,x_r,{\beta }^{\left(i\right)})}{\pi \left({\alpha }_k\right|x_1,x_2,\cdots ,x_r,{\beta }^{\left(i\right)})}\right\}\hfill \\ & =min\left\{1,\frac{{{\alpha }^{\prime }}^{r-1}{\beta }^{n{\alpha }^{\prime }-1}{\prod }_{i=1}^r{x}_i^{{\alpha }^{\prime }-1}}{{\alpha }_k^{r-1}{\beta }^{n{\alpha }_k-1}{\prod }_{i=1}^r{x}_i^{{\alpha }_k-1}}·\frac{{\prod }_{i=1}^r{(x_i^{{\alpha }_k}+{\beta }^{{\alpha }_k})}^2·{(x_r^{{\alpha }_k}+{\beta }^{{\alpha }_k})}^{n-r}}{{\prod }_{i=1}^r{(x_i^{{\alpha }^{\prime }}+{\beta }^{{\alpha }^{\prime }})}^2·{(x_r^{{\alpha }^{\prime }}+{\beta }^{{\alpha }^{\prime }})}^{n-r}}\right\}\hfill \\ & =min\left\{1,{\left(\frac{{\alpha }^{\prime }}{{\alpha }_k}\right)}^{r-1}{\beta }^{n·\Delta {\alpha }_k}{\left(\prod _{i=1}^r{x}_i\right)}^{\Delta {\alpha }_k}\frac{{\prod }_{i=1}^r{(x_i^{{\alpha }_k}+{\beta }^{{\alpha }_k})}^2}{{\prod }_{i=1}^r{(x_i^{{\alpha }^{\prime }}+{\beta }^{{\alpha }^{\prime }})}^2}{\left(\frac{x_r^{{\alpha }_k}+{\beta }^{{\alpha }_k}}{x_r^{{\alpha }^{\prime }}+{\beta }^{{\alpha }^{\prime }}}\right)}^{n-r}\right\},\hfill \end{array}$$

(17)

where $\Delta {\alpha }_k={\alpha }^{\prime }-{\alpha }_k$.

(3)

Generate a random number u from the uniform distribution of $U(0,1)$, and then obtain ${\alpha }_{k+1}$ according to the following (18):

$${\alpha }_{k+1}=\left\{\begin{array}{c}{\alpha }^{\prime },u\le {\mathrm{P}}_{\mathrm{accept}}({\alpha }_k,{\alpha }^{\prime })\hfill \\ {\alpha }_k,u>{\mathrm{P}}_{\mathrm{accept}}({\alpha }_k,{\alpha }^{\prime })\hfill \end{array}.\right.$$

(18)

(4)

Set $k=k+1$, then return to Step 1(1).

The resulting $\{{\alpha }^{\left(i\right)},i=1,2,\cdots ,k\}$ is a Markov chain, and we can obtain ${\alpha }^{(i+1)}$ when the Markov chain reaches equilibrium.

Step 2 

Obtain a sample ${\beta }^{(i+1)}$ from $\pi \left(\beta \right|x_1,x_2,\cdots ,x_r,{\alpha }^{(i+1)})$ using the Metropolis–Hastings algorithm in this step.

(1)

Set $\beta \sim \mathrm{N}({\beta }_k,{\sigma }_{\beta }^2)$, where ${\beta }_k$ is the current state and ${\sigma }_{\beta }$ is the standard deviation. Obtain a sample ${\beta }^{\prime }$ of $\beta$; if ${\beta }^{\prime }\le 0$ or ${\beta }^{\prime }>x_1$, then a new sampling is required.

(2)

Calculate the acceptance probability ${\mathrm{P}}_{\mathrm{accept}}({\beta }_k,{\beta }^{\prime })$ as (19):

$$\begin{array}{cc}\hfill {\mathrm{P}}_{\mathrm{accept}}({\beta }_k,{\beta }^{\prime })& =min\left\{1,\frac{\pi \left({\beta }^{\prime }\right|x_1,x_2,\cdots ,x_r,{\alpha }^{(i+1)})}{\pi \left({\beta }_k\right|x_1,x_2,\cdots ,x_r,{\alpha }^{(i+1)})}\right\}\hfill \\ & =min\left\{1,{\left(\frac{{\beta }^{\prime }}{{\beta }_k}\right)}^{n{\alpha }^{(i+1)}-1}·\frac{{\prod }_{i=1}^r{(x_i^{{\alpha }^{(i+1)}}+{\beta }_k^{{\alpha }^{(i+1)}})}^2}{{\prod }_{i=1}^r{(x_i^{{\alpha }^{(i+1)}}+{{\beta }^{\prime }}^{{\alpha }^{(i+1)}})}^2}\right\}.\hfill \end{array}$$

(19)

(3)

Generate a random number u from the uniform distribution of $U(0,1)$, and then obtain ${\beta }_{K+1}$ according to the following conditions (20):

$${\beta }_{k+1}=\left\{\begin{array}{c}{\beta }^{\prime },u\le {\mathrm{P}}_{\mathrm{accept}}({\beta }_k,{\beta }^{\prime })\hfill \\ {\beta }_k,u>{\mathrm{P}}_{\mathrm{accept}}({\beta }_k,{\beta }^{\prime })\hfill \end{array}.\right.$$

(20)

(4)

Set $k=k+1$, then return to Step 2(1).

Similarly, the resulting $\{{\beta }^{\left(i\right)},i=1,2,\cdots ,k\}$ is a Markov chain, and we can obtain ${\beta }^{(i+1)}$ when the Markov chain reaches equilibrium.

Step 3 

Let $i=i+1$, and repeat Step 1 and Step 2 in sequence until the Markov chains reach equilibrium.

The specific approach to obtaining the Bayesian estimation of the target parameters is similar to that in the case of complete samples and will not be further elaborated here.

4. Numerical Studies

In Section 4, we apply the R4.1.3 software in our work.

4.1. Simulation Studies

Case 1. Denote X as the life of a product, and set $X\sim \mathrm{NP}(\alpha ,\beta )$ with three groups of $\alpha$ and $\beta$. A Monte Carlo simulation was carried out in the complete samples case and 20% type II censored samples case. Generate 100 NP random variables using the inverse sampling method in each simulation.

To illustrate the complete implementation process of the mixed Gibbs sampling algorithm, we first present the Bayesian estimation and 95% credible interval of three groups of parameter simulations in both the complete samples case and the 20% type II censored samples case (refer to

Table 1. Bayesian estimation and 95% credible interval for parameters in complete samples case under different loss functions.

Loss Function	$(\mathit{\alpha },\mathit{\beta })$	$\widehat{\mathit{\alpha }}$	$\widehat{\mathit{\beta }}$	95% Credible Interval for $\mathit{\alpha }$	95% Credible Interval for $\mathit{\beta }$
LINEX loss	(2, 1)	1.8271	0.9938	[1.5392, 2.1614]	[0.9631, 1.0052]
	(3, 1.5)	2.9647	1.4964	[2.5138, 3.4898]	[1.4691, 1.5067]
	(4, 1)	3.9284	0.9955	[3.3583, 4.6654]	[0.9819, 1.0006]
Square loss	(2, 1)	2.0986	0.9977	[1.7768, 2.4561]	[0.9718, 1.0067]
	(3, 1.5)	2.8759	1.4952	[2.4206, 3.3609]	[1.4685, 1.5056]
	(4, 1)	3.9840	0.9977	[3.3579, 4.6719]	[0.9845, 1.0027]

and

Table 2. Bayesian estimation and 95% credible interval for parameters in 20% type II censored samples case under different loss functions.

Loss Function	$(\mathit{\alpha },\mathit{\beta })$	$\widehat{\mathit{\alpha }}$	$\widehat{\mathit{\beta }}$	95% Credible Interval for $\mathit{\alpha }$	95% Credible Interval for $\mathit{\beta }$
LINEX loss	(2, 1)	1.9896	0.9908	[1.2957, 3.0046]	[0.9630, 1.0008]
	(3, 1.5)	2.9361	1.4974	[1.9785, 4.5394]	[1.4681, 1.5077]
	(4, 1)	4.1126	0.9988	[2.8326, 6.5555]	[0.9855, 1.0037]
Square loss	(2, 1)	1.9576	0.9946	[1.2310, 2.8416]	[0.9637, 1.0051]
	(3, 1.5)	3.1365	1.4977	[1.9785, 4.5394]	[1.4681, 1.5077]
	(4, 1)	4.3141	0.9955	[3.5397, 6.2427]	[0.9816, 1.0002]

). The parameters are estimated using Bayesian estimation based on the sample separately under the square loss function and the LINEX loss function (with a loss parameter of 1). Furthermore, when the parameters $(\alpha ,\beta )=(4,1)$, we demonstrate the process of mixed Gibbs sampling in both the complete samples case and the 20% type II censored samples case. This includes the trace plot, histogram, and autocorrelation coefficient diagram of corresponding parameters (refer to Figure 2 and Figure 3).

As stated in Table 1 and Table 2, the Bayesian estimation of parameters using the mixed Gibbs sampling algorithm has produced relatively satisfactory results for the three groups of selected parameters. The results show that the estimated value is close to the true value, which means that the proposed method has high accuracy and reliability concerning parameter estimation. In addition, the credible interval of the parameters also covers the true value, and the interval length is relatively short, indicating high precision in the estimation of the parameters.

From Figure 2 and Figure 3, the trace plots show that the values of $\alpha$ and $\beta$ are randomly scattered around the average. Furthermore, the traversal mean of the parameters $(\alpha ,\beta )$ tends to stabilize after 200 iterations. It can be considered that the Gibbs sampling process reaches the equilibrium state at this point, and the remaining samples can be regarded as the observed value of the target parameter. As a precaution, the first 500 iterations are discarded, and then 4500 iterations are performed to obtain 4500 samples of each parameter. Additionally, the autocorrelation function indicates that the autocorrelation coefficient of all samples approaches 0 at a lag of 1 order. Therefore, the Bayesian estimation of the parameters is ultimately obtained through the 4500 samples that have been obtained.

Case 2. It is assumed that the life X of a product follows the NP distribution of two parameters $(\alpha ,\beta )=(4,1)$. To investigate the strengths and weaknesses of the two estimation methods with small sample sizes, 20 random variables from the NP distribution were generated in each simulation in the complete samples case. In the type II censored samples case, 100 random variables of the NP distribution were generated for each simulation. Furthermore, 20% of the order failure data for the complete random variables are selected in each simulation. The mean value and mean square error (MSE) of MLE and Bayesian estimation (BE) are given in

Table 3. Mean and mean square error of Bayesian estimation and maximum likelihood estimation.

Sample Case	n	Estimation Method	(Mean of $\mathit{\alpha }$, ${\mathbf{MSE}}_{\mathit{\alpha }}$)	(Mean of $\mathit{\beta }$, ${\mathbf{MSE}}_{\mathit{\beta }}$)
Complete samples	20	BE	(4.0457, 0.0389)	(0.9896, 0.0002)
		MLE	(4.4374, 0.4789)	(1.0287, 0.0014)
Type II censored samples (20%)	100	BE	(4.1136, 0.1671)	(0.9989, 1.6214 × 10${}^{-5}$)
		MLE	(4.0985, 0.4564)	(1.0056, 5.3148 × 10${}^{-5}$)

. The simulation shows that there is no significant difference between maximum likelihood estimation and Bayesian estimation in the large sample case, so we do not report the results in our article.

As shown in Table 3, the MSEs for Bayesian estimation are consistently smaller than the MLEs. Consequently, we can infer that the BE method is evidently superior to the MLE in the context of small samples. In summary, the findings highlight the superior accuracy and robustness of Bayesian estimation utilizing the mixed Gibbs sampling algorithm compared to MLE.

4.2. Real Data Analysis

Take two reliability datasets from Bourguignon [6]. The first dataset (DataSet I), given by Linhart and Zucchini [18], represents the number of failures of an aircraft’s air conditioning system: 23, 261, 87, 7, 120, 14, 62, 47, 225, 71, 246, 21,42, 20, 5, 12, 120, 11, 3, 14, 71, 11, 14, 11, 11, 90, 1, 16, 52, 95. The second dataset (DataSet II) comes from Murthy et al. [19], which is composed of the number of failures of 20 mechanical components with the following data: 0.067, 0.068, 0.076, 0.081, 0.084, 0.085, 0.085, 0.086, 0.089, 0.098, 0.098, 0.114, 0.115, 0.121, 0.125, 0.131, 0.149, 0.160, 0.485, 0.485. It can be observed that both datasets are positive. Due to the origin of the NP distribution, positive data are ideally modeled from this distribution. Therefore, it is reasonable to use the NP distribution to fit these two datasets. In this section, the Bayesian estimation of two datasets by the NP distribution is presented using the mixed Gibbs sampling algorithm.

DataSet I. For DataSet I, Bourguignon [6] provided the MLEs (i.e., ${\widehat{\alpha }}_{MLE}=0.433$, ${\widehat{\beta }}_{MLE}=1$) of parameters $\alpha$ and $\beta$. Subsequently, Bayesian analysis was conducted for both type II censored data ($r=24$) and the complete data ($r=n=30$). In our work, 2000 sampling iterations were performed for each case. Testing the sampling mean revealed that the traversal mean of the parameters stabilized within 100 iterations, and the mixed Gibbs sampling algorithm converged in both cases. To mitigate the influence of the transition process, the sampled data from 201 to 2000 iterations were utilized as the mixed Gibbs sampling values for Bayesian estimation (BE), and the posterior estimation of the distribution parameters was performed using these sampling values (refer to

Table 4. Numerical results of the real data analyses.

r	Parameters	BE	2.5% Percentile	50% Percentile	97.5% Percentile
$r=24$	$\alpha$	0.3329	0.2152	0.3304	0.4645
	$\beta$	0.7609	0.3180	0.7963	0.9923
$r=n=30$	$\alpha$	0.4279	0.3152	0.4264	0.5590
	$\beta$	0.8723	0.5896	0.9016	0.9956

).

In Table 4, the posterior expectation and median represent the Bayesian point estimates of the parameters under square loss and absolute loss. The 95% credible intervals of the parameters are determined by the 2.5th and 97.5th percentiles. The results in Table 4 demonstrate the high accuracy of Bayesian estimation using the mixed Gibbs sampling algorithm for both complete data and type II censored data. Compared to frequency estimation methods, Bayesian estimation proves to be more effective in handling small sample sizes and provides probability distribution estimates of parameters, rather than just point estimates, for reliability data with limited observations. Consequently, Bayesian estimation offers a more precise reflection of parameter uncertainties, particularly when dealing with small datasets.

DataSet II. In the case of complete data, we conducted Bayesian analysis for DataSet II with 2000 sampling iterations. After testing the sampling mean, we found that the traversal mean of the parameters tended to be stable within 100 iterations, and the mixed Gibbs sampling algorithm converged. To eliminate the influence of the transition process, the sampled data of 201 to 2000 times were used as the mixed Gibbs sampling values for Bayesian estimation, and the posterior estimation of the distribution parameters was performed using these sampling values. Then, the Bayesian estimation (BE) results under the square loss were compared with the MLE results of Bourguignon [6] (refer to

Table 5. Comparison of two estimation methods.

Estimation Method	Sample Case	Parameter	Estimation Results
MLE	Complete samples	$\alpha$	2.871
		$\beta$	0.067
BE	Complete samples	$\alpha$	2.823
		$\beta$	0.065

).

It can be seen from the above Table 5 that the Bayesian estimation is very close to the MLE result, which further verifies the stability and reliability of the mixed Gibbs sampling algorithm.

5. Discussion and Conclusions

Some scholars have conducted valuable research for the NP distribution in the literature, such as [10,11,15]. However, Bayesian methods may encounter challenges with high-dimensional, complex models, and big data in practice. Contrary to the low efficiency of importance sampling [10] and the complexity of the Lindley approximation [12] methods, a Bayesian estimation procedure is proposed for the NP distribution using the mixed Gibbs sampling algorithm. When the fully conditional distributions in Gibbs sampling steps lack an explicit form, direct sampling from these distributions becomes challenging. In such instances, alternative sampling methods must be proposed. As each step within a loop of Gibbs sampling is itself a Metropolis–Hastings iteration, the algorithm can be used to combine Gibbs sampling with the Metropolis–Hastings algorithm if necessary.

In this paper, we proposed a Bayesian estimation procedure for the NP distribution with complete and type II censored samples using the mixed Gibbs sampling algorithm. For instance, the mean and MSEs of the resulting parameters estimations are given, and the stability and effectiveness of the mixed Gibbs sampling algorithm are demonstrated with simulation studies and real data analysis. It can be seen that the mixed Gibbs sampling algorithm is stable, feasible, and high precision.

In reliability-related product tests, it is often challenging to obtain a sufficient number of samples to assess the quality of the products. Bayesian estimation, however, can assist in obtaining reliable results even in cases of limited sample sizes. By incorporating prior information, Bayesian estimation can improve the accuracy of parameter estimation and provide valuable uncertainty information about the parameters, which is particularly crucial when sample data is scarce. Therefore, in the context of product reliability tests, Bayesian estimation can serve as an effective tool for reliably evaluating and predicting product quality in situations with small sample sizes. Thus, in view of the merits of the Bayesian method with the mixed Gibbs sampling algorithm in this paper, it will become an important Bayesian inference method.