Type-I Heavy-Tailed Burr XII Distribution with Applications to Quality Control, Skewed Reliability Engineering Systems and Lifetime Data

1  Introduction

In the analysis of real-world phenomena, rare events and complex hazard behaviors are common, necessitating the use of flexible statistical distributions. Classical models such as the normal, exponential, Weibull, and Rayleigh often prove inadequate in modeling data characterized by significant skewness, heavy tails, or non-monotonic hazard functions. Consequently, the development of more versatile distributions, particularly those with heavy-tailed behavior, has been a central focus of modern applied statistics.

Heavy-tailed distributions, which allocate a higher probability mass to extreme outcomes than Gaussian models, have been observed across a wide range of disciplines. In finance, for example, asset prices and stock returns frequently exhibit fat-tailed tendencies, with the probability of large gains or losses exceeding what would be predicted by a normal distribution [1]. This has led to the adoption of stable laws to better explain high market volatility. Similarly, in actuarial science, catastrophic claims in insurance are a classic example of heavy-tailed behavior, particularly when modeling loss severity [2]. Environmental and natural systems also display these characteristics; hydrological records, rain intensities, and seismic events often follow distributions with slowly decaying tails [3,4]. In medicine and biology, skewed and extreme-value data, such as tumor sizes and survival times, similarly require heavy-tailed models for accurate analysis [5–7].

Formally, a distribution with distribution function F(x) is defined as heavy-tailed if its survival function F¯(x)=1−F(x) converges to zero more slowly than any exponential rate, that is,

limx→∞F¯(x)e−λx=∞for all λ>0.

This implies that heavy-tailed distributions place significantly more probability mass in the tails than exponential-type distributions. Classic examples include the Cauchy distribution by Cauchy [8], Student’s t distribution by Student [9], Fréchet distribution by Fréchet [10], and the Pareto family of distributions by Pareto [11], and Arnold [12], which includes the Pareto II (Lomax) distribution by Lomax [13]. A specific and widely-studied type of heavy-tailed distribution is the power-law distribution, characterized by a survival function that decays as a power of x for large x:

F¯(x)∼x−αas x→∞,

where α>0 is the tail index. This slow, polynomial decay of the tail is a defining feature that distinguishes power-law distributions from those with exponential tails. Given that empirical data frequently violate the assumptions of standard models like the exponential, gamma, and Weibull distributions, researchers have increasingly turned to advanced generalization strategies such as transformation techniques, compounding, and mixture models to construct distributions with superior tail behavior and greater structural flexibility, see [14–16] as well as [17–20].

Among these flexible models, the Burr Type XII (BXII) distribution has been particularly influential. Introduced by Burr [21] as part of a system of distributions designed to fit a wide range of empirical data, it is a continuous probability distribution for a non-negative random variable x≥0. Singh [22] demonstrated the applicability of Burr-type distributions in different data scenarios. The conventional Burr Type XII distribution is characterized by two positive shape parameters, commonly represented as c and k, and its high degree of flexibility is well-known. By altering the shape parameters, the BXII distribution can take on many forms and approximate numerous other common distributions, such as the log-logistic, Lomax (Pareto Type II), and Weibull, and can even approximate the normal distribution under certain conditions. It offers a parsimonious yet powerful framework for analyzing data across diverse fields, including reliability, survival analysis, and econometrics.

The ongoing demand for models capable of capturing complex data characteristics has led to a proliferation of BXII generalizations. These extensions often enhance the original BXII’s flexibility by adding new parameters or by embedding it within a larger family of distributions. For example, Gad et al. [23] introduced the Burr XII-Burr XII (BXI-BXII) distribution, a compound model that demonstrated a superior fit to certain datasets compared to the baseline BXII. The literature is rich with such constructions, which can be broadly categorized by their methodological approach.

One common strategy is the use of general families of distributions. The T-X family by Alzaatreh et al. [24], for instance, has been used to derive the Burr XII-moment exponential (BXII-ME) distribution, which offers greater flexibility for modeling lifetime data [25]. Similarly, the Kumaraswamy-G family has led to the Kumaraswamy Burr XII (KBXII) distribution, a four-parameter model with a remarkable ability to accommodate diverse functional shapes and hazard rate behaviors [26]. Another important approach involves the Marshall-Olkin extended family, which produced the Marshall-Olkin extended Burr Type XII (MOEBXII) distribution to provide enhanced flexibility for reliability and survival analysis [27]. For more details see [28–30] as well as [31–33].

Other specialized extensions of the BXII distribution have been developed to address specific challenges. Ocloo [34] introduced a novel extension of the Burr XII distribution, demonstrating its enhanced capability to model data with complex characteristics such as bimodality and heavy tails, and providing a better fit than many existing BXII extensions. Recent innovations also include the Sine-Burr XII distribution, which utilizes a trigonometric transformation to add flexibility without introducing new parameters [35], and the Odd Burr XII Gompertz (OBXIIGo) distribution, which provides a more adaptable model for lifetime data by integrating the BXII family with the Gompertz distribution [36].

The work of Gad et al. [23] and other similar studies underscore a critical motivation in modern statistical modeling: the development of new distributions that can effectively capture complex data structures not adequately described by classical models. The literature reviewed herein indicates a clear trend toward creating higher-order BXII-based models that offer greater flexibility in their distributional shapes and hazard rate functions. These models consistently serve as a better alternative to their classical counterparts, a finding supported by statistical goodness-of-fit metrics.

One such area where the choice of distribution is paramount is in acceptance sampling, particularly within truncated life-testing experiments. In specialized manufacturing or medical diagnosis, using simple, stratified, or cluster sampling may be insufficient. Group acceptance sampling plans (GASPs) provide an efficient alternative by examining several units in batches, conserving resources. The efficacy of a GASP, however, hinges on the choice of the underlying product lifetime distribution. An ill-suited model can lead to misclassifying good or bad lots, resulting in either excessive consumer risk or uneconomical producer rejection. This has prompted researchers to propose GASPs under more generalized distributions to overcome the limitations of standard models. For example, GASPs have been developed for truncated life tests based on the inverse Rayleigh and log-logistic distributions by Aslam & Jun [37], the Marshall-Olkin extended Lomax distribution by Rao [38], and the inverse Weibull distribution with median lifetime as a quality measure by Singh & Tripathi [39]. Other recent studies have employed generalized distributions like the AG-transmuted exponential by Almarashi & Khan [40], transmuted exponential by Owoloko et al. [41], transmuted Rayleigh by Saha et al. [42], and exponential-logarithmic distributions by Ameeq et al. [43] in the GASP framework. This body of research demonstrates that superior modeling of the lifetime distribution directly leads to more efficient and risk-balanced sampling decisions. Most recently, Ekemezie et al. [44] presented the odd Perks-Lomax (OPL) distribution and used it to construct a GASP, showcasing how a novel, versatile distribution can enhance both statistical theory and practical application simultaneously.

Building on this demonstrated need for more intricate models, this research proposes and thoroughly investigates a novel statistical distribution: the Type-I Heavy-Tailed Burr Type XII (TIHTBXII) distribution. The primary motivation is to address the persistent challenges in adequately modeling real-world data with extreme skewness, heavy tails, and complex, non-monotonic hazard functions, for which existing generalized distributions still exhibit limitations. We not only meticulously derive the mathematical properties of the TIHTBXII distribution and examine various parameter estimation procedures but, more significantly, we demonstrate its practical utility by constructing a group acceptance sampling plan (GASP) specifically tailored to truncated life tests based on this new model. This new tool leverages the TIHTBXII distribution’s excellent capability to provide more accurate and reliable data, thereby enabling better decision-making in fields like quality control, reliability engineering, and lifetime data analysis, where current models may fail to capture the true underlying risk mechanisms.

The rest of this paper is organized as follows: In Section 2, we define the basic functions of the proposed TIHTBXII model, including its probability density function (PDF), cumulative distribution function (CDF), survival function, and hazard function. Section 3 is dedicated to the derivation of several key characteristics of the TIHTBXII model, such as its quantile function, moments, incomplete moments, moment-generating function, mean residual life function, entropy, extropy, and order statistics. In Section 4, we discuss the estimation of the model parameters using maximum likelihood, maximum product spacing, least squares, and weighted least squares methods. Section 5 presents a Monte Carlo Markov Chain (MCMC) simulation study conducted with four different parameter settings, with the bias and root mean square error (RMSE) also plotted to support the numerical results. In Section 6, we design a group acceptance sampling plan (GASP) to illustrate how the TIHTBXII model can optimize quality control decisions by minimizing the average sample number. Finally, Section 7 reports the applications of the proposed model to real datasets, including the duration of active repairs of airborne communication transceivers, groundwater contaminant measurements, and Dominica COVID-19 mortality rate data, while Section 8 provides the conclusion of this study.

2  Development of the Type-I Heavy-Tailed Burr XII Distribution

The Burr XII distribution has cumulative distribution function (CDF) and probability density function (PDF), respectively, given as

F(x)=1−1(1+xc)k;x>0,c,k>0,(1)

and

f(x)=ckxc−1(1+xc)k+1.(2)

Reference [45] proposed the Type I heavy-Tailed (TI-HT) family of distributions with CDF and PDF respectively given as

G(x;θ,ξ)=1−(1−F(x;ξ)1−(1−θ)F(x;ξ))θ;x∈ℜ,θ>0,(3)

and

g(x;θ,ξ)=θ2f(x;ξ){1−F(x;ξ)}θ−1{1−(1−θ)F(x;ξ)}θ+1,(4)

where F(x;ξ) is the baseline distribution which depends on ξ∈ℜ. Introducing Eqs. (1) and (2) into (3) and (4) produces a special distribution referred to as type-I heavy-tailed Burr XII (TIHTBXII) distribution with CDF and PDF respectively defined as

G(x;θ,c,k)=1−[11−θ{1−(1+xc)k}]θ;x>0,θ,c,k>0,(5)

and

g(x;θ,c,k)=θ2ckxc−1(1+xc)k−1[1−θ{1−(1+xc)k}]θ+1.(6)

The survival function is

s(x;θ,c,k)=[11−θ{1−(1+xc)k}]θ;x>0,θ,c,k>0.(7)

The hazard function is

h(x;θ,c,k)=θ2ckxc−1(1+xc)k−1[1−θ{1−(1+xc)k}]1;x>0,θ,c,k>0.(8)

The proposed TIHTBXII distribution aims to improve model flexibility, provide the best fit to real-world data, and accommodate heavy-tailed data in fields such as reliability engineering, medical and financial sciences, among others.

Fig. 1a and b represents the graph of the PDF and hazard function, respectively. These plots are at various combination of the parameter values. The plots reflect left-skewed behaviour, both low and high peaks, and bumped-shape. This demonstrate that the TIHTBXII can be utilized in modeling different lifetime events. The hazard rate plots depict bump-shape, bathtub, L-shape, J-shape and strictly non-decreasing shape. This again implies that TIHTBXII model can benefit from modeling different lifetime events.

Figure 1: Plots of (a) PDF, (b) Hazard function of the TIHTBXII distribution

3  Characterization of the TIHTBXII Model

In this section, we study some basic properties of the new distribution which include quantile function, moment and incomplete moment, moment generating function, order statistic, entropy and extropy, and the mean residual life function.

3.1 Quantile Function

The quantile function Q(u)=G−1(u), for 0<u<1, is derived as:

Q(u)={[1−(1−(1−u)−1θθ)]1k−1}1c.(9)

In inference, the advantage of Q(u) is in generating random samples and overall simulation studies. Further, various quantile measures can be obtained such as the mean, variance, skewness and kurtosis.

3.2 Moment

The r-th crude moment of a continuous random variable X which follows any distribution say TIHTBXII (θ,c,k) can be expressed as follows

μr′=E(Xr)=∫−∞∞xrg(x) dx=∫0∞θ2ckxr+c−1(1+xc)k−1[1−θ{1−(1+xc)k}]θ+1 dx=θ2ck∫0∞xr+c−1(1+xc)k−1[1−θ{1−(1+xc)k}]−(θ+1) dx.

Using series expansion

(1+xc)k−1=∑i=0∞(k−1i)xic;for|x|<1.

For θ∈R or C, a generalized binomial theorem provides that

[1−θ{1−(1+xc)k}]−(θ+1)=∑j=0∞(θ+jj)θj[1−(1+xc)k]j.

Also, for j∈R/N0 and |x|<1, standard binomial expansion provides that

[1−(1+xc)k]j=∑l=0∞(−1)l(jl)(1+xc)lk.

Similarly,

(1+xc)lk=∑m=0∞(lkm)xmc.

By back substitution,

μr′=θ2ck∑i=0∞∑j=0∞∑l=0∞∑m=0∞(−1)l(k−1i)(θ+jj)(jl)(lkm)θj∫0∞xr+c+ic+mc−1 dx.

Notice that

(jl)=j!l!(j−l)!=(j,l)l!;Recalle−x=∑l=0∞(−1)lxll!,

so that

μr′=θ2ck∑i=0∞∑j=0∞∑m=0∞(k−1i)(θ+jj)(lkm)θj(j,l)∫0∞∑l=0∞(−1)lxll!1xlxr+c+ic+mc−1 dx=θ2ck∑i=0∞∑j=0∞∑m=0∞(k−1i)(θ+jj)(lkm)θj(j,l)∫0∞xr+c+ic+mc−l−1e−x dx.

Therefore,

μr′=θ2ck∑i=0∞∑j=0∞∑m=0∞(k−1i)(θ+jj)(lkm)θj(j,l)Γ(r+c+ic+mc−l);r=1,2,⋯(10)

3.3 Incomplete Moment

For a continuous random variable X∼TIHTBXII(θ,c,k), the r-th incomplete moment about the origin is given by:

μr(x0)=E[Xr∣X≤x0]=∫0x0xrg(x;θ,c,k)dx.

Substitute into the incomplete moment:

μr(x0)=θ2ck∫0x0xr+c−1(1+xc)k−1[1−θ(1−(1+xc)k)]−(θ+1)dx.

Using the generalized binomial expansion:

(1+xc)k−1=∑i=0∞(k−1i)xic.

Expanding the second term yields

[1−θ(1−(1+xc)k)]−(θ+1)=∑j=0∞(θ+jj)θj[1−(1+xc)k]j.

and

[1−(1+xc)k]j=∑l=0∞(−1)l(jl)(1+xc)lk.

Also,

(1+xc)lk=∑m=0∞(lkm)xmc.

The integrand becomes a power series in x:

xr+c−1⋅∑i=0∞(k−1i)xic⋅∑j=0∞(θ+jj)θj∑l=0∞(−1)l(jl)∑m=0∞(lkm)xmc.

Now collect all powers of x, and integrate:

μr(x0)=θ2ck∑i,j,l,m(−1)l(k−1i)(θ+jj)θj(jl)(lkm)∫0x0xr+c−1+ic+mcdx.∫0x0xr+c−1+ic+mcdx=x0r+c+ic+mcr+c+ic+mc.

Hence,

μr(x0)=θ2ck∑i=0∞∑j=0∞∑l=0∞∑m=0∞(−1)l(k−1i)(θ+jj)θj(jl)(lkm)⋅x0r+c+ic+mcr+c+ic+mc.(11)

3.4 Moment Generating Function

Let t∈R+, the moment generating function (MGF) of any continuous random variable X given as MX(t) which assumes the TIHTBXII(θ,c,k) is defined as

MX(t)=E(etX)=∫−∞∞etxg(x) dx=∫0∞θ2ckxc−1(1+xc)k−1etx[1−θ{1−(1+xc)k}]θ+1 dx=θ2ck∫0∞xc−1(1+xc)k−1etx[1−θ{1−(1+xc)k}]−(θ+1) dx.

Using series expansion

(1+xc)k−1=∑i=0∞(k−1i)xic;for|x|<1.

For θ∈R or C, a generalized binomial theorem provides that

[1−θ{1−(1+xc)k}]−(θ+1)=∑j=0∞(θ+jj)θj[1−(1+xc)k]j.

Also, for j∈R/N0 and |x|<1, standard binomial expansion provides that

[1−(1+xc)k]j=∑l=0∞(−1)l(jl)(1+xc)lk.

Similarly,

(1+xc)lk=∑m=0∞(lkm)xmc;ex=∑n=0∞xnn!.

By back substitution,

MX(t)=θ2ck∑i=0∞∑j=0∞∑l=0∞∑m=0∞∑n=0∞(−1)ln!(k−1i)(θ+jj)(jl)(lkm)θj∫0∞xn+c+ic+mc−1 dx.

Notice that

(jl)=j!l!(j−l)!=(j,l)l!;Recalle−x=∑l=0∞(−1)lxll!,

so that

MX(t)=θ2ck∑i=0∞∑j=0∞∑m=0∞∑n=0∞(k−1i)(θ+jj)(lkm)θj(j,l)n!∫0∞∑l=0∞(−1)lxll!1xlxn+c+ic+mc−1 dx=θ2ck∑i=0∞∑j=0∞∑m=0∞∑n=0∞(k−1i)(θ+jj)(lkm)θj(j,l)n!∫0∞xn+c+ic+mc−l−1e−x dx.

Therefore,

MX(t)=θ2ck∑i=0∞∑j=0∞∑m=0∞∑n=0∞(k−1i)(θ+jj)(lkm)θj(j,l)Γ(n+c+ic+mc−l)n!.(12)

3.5 Mean Residual Life Function

The mean residual life (MRL) function of a non-negative continuous random variable X is defined as

m(t)=E[X−t∣X>t]=1G¯(t)∫t∞G¯(x)dx,t≥0,(13)

where G¯(t)=1−G(t) is the survival function. This function quantifies the expected remaining lifetime given survival up to time t. Suppose X∼ TIHTBXII(θ,c,k), the MRL is given as follows

m(t)=[1−θ+θ(1+tc)k]θ∫t∞[1−θ+θ(1+xc)k]−θ dx.

Recall that for j∈R/N0 and |x|<1, the standard binomial expansion provides that

(1+x)−n=∑r=0∞(−1)r(n+r−1r)xr,

hence

m(t)=[1−θ+θ(1+tc)k]θ∫t∞∑r=0∞(−1)r(θ+r−1r)θr[1−(1+xc)k]r dx.

Similarly, by expanding [1−(1+xc)k]r, we have

m(t)=[1−θ+θ(1+tc)k]θ∫t∞∑r=0∞∑s=0∞(−1)r+s(θ+r−1r)(rs)θr(1+xc)ks dx.

Also, by expanding (1+xc)ks, we realize

m(t)=[1−θ+θ(1+tc)k]θ∫t∞∑r=0∞∑s=0∞(−1)r+s(θ+r−1r)(rs)θr∑y=0∞(ksy)xyc dx=[1−θ+θ(1+tc)k]θ∑r=0∞∑s=0∞∑y=0∞(−1)r+s(θ+r−1r)(rs)(ksy)θr∫t∞xyc dx.

Recall

(rs)=r!s!(r−s)!=(r,s)s!;ande−x=∑s=0∞(−1)sxss!m(t)=[1−θ+θ(1+tc)k]θ∫t∞∑r=0∞∑y=0∞∑s=0∞(−1)sxss!1xs(−1)r(θ+r−1r)(ksy)(r,s)θrxyc dx=[1−θ+θ(1+tc)k]θ∑r=0∞∑y=0∞(−1)r(θ+r−1r)(ksy)(r,s)θr∫t∞xyc−se−x dx.

Therefore, the MRL function is

m(t)=[1−θ+θ(1+tc)k]θ∑r=0∞∑y=0∞(−1)r(θ+r−1r)(ksy)(r,s)θrΓ(yc−s+1,t).(14)

3.6 Entropy

Entropy is a fundamental concept in information theory and statistical mechanics, used to quantify the uncertainty or randomness associated with a probability distribution. For a given probability density function g(x), the Rényi entropy of order q>0, q≠1, is defined as:

Hq(x)=11−qlog⁡{∫−∞∞gq(x)dx}.

In this section, we derive the Rényi entropy of the proposed distribution. Due to the complexity of the density function, the integral involved is analytically intractable in closed form. However, by applying the generalized binomial expansion and power series techniques, we express the entropy in terms of infinite series and known special functions, particularly the Gamma function.

Hq(x)=11−qlog⁡{∫−∞∞gq(x) dx}=11−qlog⁡{(θ2ck)q∫0∞xq(c−1)(1+xc)q(k−1)[1−θ(1−(1+xc)k)]q(θ+1) dx}=11−qlog⁡{(θ2ck)q∫0∞xq(c−1)(1+xc)q(k−1)[1−θ(1−(1+xc)k)]−q(θ+1) dx}.

Using generalized binomial expansion when q(θ+1)∈R and |x|<1, it can be shown that

[1−θ(1−(1+xc)k)]−q(θ+1)=∑i=0∞∑j=0∞∑l=0∞(−1)j(q(θ+1)+i−1i)(ij)(ikl)θixlc.

Similarly,

(1+xc)q(k−1)=∑h=0∞(q(k−1)h)xhc.

Therefore,

Hq(x)=11−qlog⁡{(θ2ck)q∫0∞xq(c−1)∑h=0∞(q(k−1)h)xhc∑i=0∞∑j=0∞∑l=0∞(−1)j(q(θ+1)+i−1i)(ij)(ikl)θixlc}=11−qlog⁡{(θ2ck)q∑i=0∞∑j=0∞∑l=0∞∑h=0∞(−1)j(q(θ+1)+i−1i)(ij)(ikl)(q(k−1)h)θi∫0∞xqc−q+hc+lc dx}.

Recall that

(ij)=i!j!(i−j)!=P(i,j)j!;ande−x=∑j=0∞(−1)jxjj!,

Hq(x)=11−qlog⁡{(θ2ck)q∑i=0∞∑l=0∞∑h=0∞∑j=0∞(−1)j∫0∞xjj!1xj(i,j)θi(q(θ+1)+i−1i)(ikl)(q(k−1)h)xqc−q+hc+lc dx}=11−qlog⁡{(θ2ck)q∑i=0∞∑l=0∞∑h=0∞P(i,j)θi(q(θ+1)+i−1i)(ikl)(q(k−1)h)∫0∞xqc−q+hc+lc−je−x dx}=11−qlog⁡{(θ2ck)q∑i=0∞∑l=0∞∑h=0∞P(i,j)θi(q(θ+1)+i−1i)(ikl)(q(k−1)h)Γ(qc−q+hc+lc−j+1)}.

3.7 Extropy

This subsection defines extropy, a measure of order or information content within a system, serving as a dual to the more commonly known concept of entropy. We consider a specific mathematical formulation of extropy, denoted by the functional J(g). The initial representation of J(g) is an infinite series where each term encompasses an integral over a complex function g(x). To analytically evaluate this integral, we systematically apply the generalized binomial expansion to several nested terms, including [1−θ{1−(1+xc)k}]−n(θ+1), [1−(1+xc)k]i, and (1+xc)jk. This decomposition transforms the original integral into a convergent multifold summation. The final integration step is then performed, yielding an expression for J(g) in terms of the Gamma function, Γ(z), thus providing a complete analytical series expansion.

J(g)=∑n=2∞1n(n−2)∫0∞gn(x) dx=∑n=2∞1n(n−2)∫0∞(θ2ck)nxn(c−1)(1+xc)n(k−1)[1−θ{1−(1+xc)k}]n(θ+1) dx=∑n=2∞(θ2ck)nn(n−2)∫0∞xn(c−1)(1+xc)n(k−1)[1−θ{1−(1+xc)k}]−n(θ+1) dx.

Using generalized binomial expansion; for n(θ+1)∈R and |x|<1,

[1−θ{1−(1+xc)k}]−n(θ+1)=∑i=0∞(n(θ+1)+i−1i)θi[1−(1+xc)k]i.

Also,

[1−(1+xc)k]i=∑j=0∞(−1)j(ij)(1+xc)jk.

Similarly,

(1+xc)jk=∑h=0∞(jkh)xhc;and(1+xc)n(k−1)=∑l=0∞(n(k−1)l)xlc.

So that,

J(g)=∑n=2∞(θ2ck)2n(n−1)∑i=0∞∑j=0∞∑h=0∞∑l=0∞(−1)j(ij)(n(k−1)l)(jkh)(n(θ)+i−1i)θi∫0∞xnc−n+lc+hc dx.

Since,

(ij)=i!j!(i−j)!=P(i,j)j!;e−x=∑j=0∞(−1)jxjj!,

J(g)=∑n=2∞(θ2ck)nn(n−1)∑i=0∞∑l=0∞∑h=0∞(n(k−1)l)(jkh)(n(θ+1)+i−1i)∫0∞∑j=0∞(−1)jxjj!1xjP(i,j)θixnc−n+lc+hc dx=∑n=2∞(θ2ck)nn(n−1)∑i=0∞∑l=0∞∑h=0∞(n(k−1)l)(jkh)(n(θ+1)+i−1i)θiP(i,j)∫0∞xnc−n+lc+hc−je−x dx=∑n=2∞(θ2ck)nn(n−1)∑i=0∞∑l=0∞∑h=0∞(n(k−1)l)(jkh)(n(θ+1)+i−1i)θiP(i,j)Γ(nc−n+lc+hc−j+1).

3.8 Order Statistic

In the realm of probability theory and statistics, an order statistic is a fundamental concept representing the r-th smallest value among a set of n random variables. Given a set of n independent and identically distributed (i.i.d.) random variables, say X1,X2,…,Xn, arranged in ascending order as X(1)≤X(2)≤⋯≤X(n), X(r) is the r-th order statistic. The PDF of the r-th order statistic, fXr:n(x), is derived from the PDF of the individual random variables, g(x), and their CDF, G(x). For the proposed TIHTBXII model, it is given as

fXr:n(x)=n!(r−1)!(n−r)![G(x)]r−1[1−G(x)]n−rg(x)=n!(r−1)!(n−r)!{1−[1−θ(1−(1+xc)k)]−θ}r−1{1−θ(1−(1+xc)k)}−θ(n−r)×θ2ckxc−1(1+xc)k−1[1−θ{1−(1+xc)k}]θ+1=n!θ2ckxc−1(r−1)!(n−r)!∑i,j=0∞∑l,h=0∞∑m,p=0∞∑q,s=0∞∑t,w=0∞(−1)i+p+m+sθh+j+l(r−1i)(iθ+h−1h)(hp)(pkq)×(θ(n−r)+j−1j)(jm)(kmw)(θ+ll)(ls)(kst)xqc+wc+tc.

4  Estimation

We proceed with the estimation of the parameters of the TIHTBXII distribution using maximum likelihood, maximum product spacing, least squares and weighted least squares procedures.

4.1 Maximum Likelihood Estimation

Maximum Likelihood Estimation (MLE) involves finding the values of the parameters that maximize the likelihood of observing the given data. First, we formulate the likelihood function. Let X1,X2,…,Xn be a random sample of size n from the distribution with the given PDF g(x;θ,c,k). Assuming the observations are independent and identically distributed (i.i.d.), the likelihood function L(θ,c,k|x) is the product of the individual PDFs:

L(θ,c,k|x)=∏i=1ng(xi;θ,c,k)=∏i=1nθ2ckxic−1(1+xic)k−1[1−θ{1−(1+xic)k}]θ+1.

Formulate the log-likelihood function since it is almost always easier to work with the natural logarithm of the likelihood function, called the log-likelihood function, denoted as l(θ,c,k|x) or ln⁡L. This is because the logarithm is a monotonically increasing function, so maximizing ln⁡L is equivalent to maximizing L. Also, products become sums, which are much simpler to differentiate. So, the log-likelihood function is:

l(θ,c,k|x)=n(2ln⁡θ+ln⁡c+ln⁡k)+(c−1)∑i=1nln⁡xi+(k−1)∑i=1nln⁡(1+xic)−(θ+1)∑i=1nln⁡[1−θ{1−(1+xic)k}],