# Modeling Cancer Remission Time Data by Means of the Max Erlang Binomial Distribution

Authors: Bogdan Gheorghe Munteanu
Journal: Computational and Mathematical Methods in Medicine
 (2021)

---

## Abstract

In this paper, a statistical simulation algorithm for the power series distribution, called the Max Erlang Binomial distribution, is proposed, analyzed, and tested for bladder cancer remission time data. In order to present the simulation technique, the EM algorithm for statistical estimation aimed at estimating the model parameters is described.


---

## Body

## 1. Introduction

The introduction of this new (generalized) distribution addresses reliability problems when lifetime can be expressed as the maximum or minimum of a sequence of independent and identically distributed (iid) random variables, which represents the system components’ risk times. In recent years, some researchers have proposed a series of new distributions for the maximum and minimum of a sequence of iid random variables. For example, Adamidis and Loukas [1], Kus [2], Tahmasbi and Rezaei [3], Louzada et al. [4], and Cancho et al. [5] were interested in determining the maximum or minimum distribution when the components in a sequence of iid random variables are exponentially distributed, and a number of components are of a discrete type. Next, Flores et al. [6] treated the distribution of a vector’s maximum with components that are exponentially distributed in a random number of a power series distribution type. This type of distribution is called the complementary exponential power series (CEPS) distribution. Also, Morais and Barreto-Souza [7] considered analyzing the Weibull distribution class by means of the power series distribution class (WPS). Recently, Louzada et al. [8] have developed a mathematical model that unifies the procedure for obtaining a distribution of the maximum and minimum of a sequence of iid random variables of the absolutely continuous type in a random number N characterized by the generating function. But the problem of determining the general formula when the random variable N forms a part of a power series distributions remains unsolved.In this paper, the simulation algorithms for these family distributions are proposed. This study is intended as a completion of the research by Balkema and de Haan (1974), Bryson (1974), Ahsanullah (1991), Balakrishnan and Ahsanullah (1994), Childs and others (2001), Al Awadhi and Ghitany (2001, 2007), Zahrani and Harbi (2013), Al-Zahrani and Sagor (2014), Tahir and Cordeiro ([9], 2016), Hassan and Abd-Elfattah (2016), and Munteanu ([10], 2013). The above-mentioned algorithm was implemented by means of the Eclipse SDK programming environment.This work has the following structure: Section2 defines the mathematical properties of the Max Erlang Binomial power series distribution (the cumulative distribution function, the probability density function, the mean, and variance). The simulation techniques targeting the Max Erlang Binomial distribution are analyzed and formulated in Section 3, with results validation via the Pearson test. In Section 4, the simulation algorithm for the Max Erlang Binomial distribution parameters is proposed and tested using the method of the maximum likelihood estimation. Section 5 discusses an application of the proposed distribution using a real-life dataset. Lastly, in Section 6, some useful conclusions are drawn.


## 2. Development of the Mathematical Model

In [11], the properties of a new power distribution type series, called the Max Erlang Binomial (MaxErlB), are introduced and researched. As a mathematical model, this distribution describes the probabilistic behavior of lifetimes, widely used in researching the reliability of systems. In [11], this distribution is presented as the distribution of the maximum value in a random volume sample Z from a statistical population, Erlang distributed, where Z is a binomially distributed, zero-truncated random variable. Formally, things are presented as follows.Let us consider the random variableZ such that ℙZ∈1,2,⋯=1.Definition 1 ([12]).
We say that the random variableZ has a power series distribution if
(1)ℙZ=z=azΘzAΘ,z=1,2,⋯;Θ∈0,τ,where a1,a2,⋯ are nonnegative real numbers, τ is a positive number bounded by the convergence radius of power series (series function) AΘ=∑z≥1azΘz,∀Θ∈0,τ, and Θ is the power parameter of the distribution (Table 1).Table 1
The representative elements of the PSD families for various truncated distributions.
DistributionazΘAΘτBinom∗n,pnzp1−p1+Θn−1∞Poisson∗α1z!αeΘ−1∞Logp1zp−ln1−Θ1Geom∗p11−pΘ1−Θ1Pascalk,pz−1k−11−pΘ1−Θk1Bineg∗k,pz+k−1zp1−Θ−k−11PSD denotes the power series distribution function families. If the random variableZ has the distribution in Equation (1), then we write Z∈PSD.We consider thatXi~Erlangk,λ,k∈ℕ,k≥1, λ>0, where Xii≥1 are iid random variables with the distribution function FXix≡FErlx=1−∑i=0k−1λxi/i!e−λx,x>0, and the probability density function fXix≡fErlx=λkxk−1e−λx/k−1!,x>0.We note thatUErl=maxX1,X2,⋯,XZ.The results in this section are obtained using the general framework in [13], for which reason some proofs are not presented.Proposition 1 (see [11]).
If the random variableUErl=maxX1,X2,⋯,XZ, where Xii≥1 are nonnegative iid random variables, Xi~Erlangk,λ,k∈ℕ,k≥1,λ>0 and Z~Binom∗n,p, Z∈PSD, n∈1,2⋯, with AΘ=Θ+1n−1, Θ=p/1−p, p∈0,1, Θ∈0,τ, τ>0, the random variables Xii≥1 and Z being independent; then, the cumulative distribution functions and the probability density function of the random variable UErlB are the following:
(2)UErlBx=1−pe−λx∑i=0k−1λxi/i!n−1−pn1−1−pn,x>0,(3)uErlBx=npλkxk−1e−λx1−pe−λx∑i=0k−1λxi/i!n−11−1−pn,x>0.Definition 2 (see [11]).
We say that the random variableUErlB has a Max Erlang Binomial power series distribution with parameters k,λ, n, and p (UErlB~MaxErlBk,λ,n,p), if it has the cumulative distribution function (cdf) defined by Equation (2) and probability density function (pdf) defined by Equation (3).The numerical characteristics (mean, variance) of a random variable with a MaxErlB distribution, in a particular case (k=2), are presented in the following result:Proposition 2.
The mean and variance of the random variableUErlB~MaxErlB2,λ,n,p, λ>0, n∈1,2⋯, p∈0,1, are characterized by the following relations:
(4)EUErlB=nλ1−1−pn∑z=1n−1z−1nz−1z+1!pzzz+2,(5)VarUErlB=nλ21−1−pn∑z=1n−1z−1nz−1z+2!pzzz+3−n1−1−pn∑z=1n−1z−1nz−1z+1!pzzz+22.Proof.
After Equation (3) and the definition of the mean, we obtain
(6)EUErlB=npλ21−1−pn∫0∞x2e−λx1−pe−λx1+λxn−1dx,where 1−pe−λx1+λxn−1=∑z=1n−1z−1nz−1pz−11+λxz−1e−z−1λx, as developed by Newton’s binomial. A sum of n-integrals then can be solved with elementary methods (method of integration by parts), which leads to Equation (4).
Similarly, evaluating the second-order moment(7)EUErlB2=npλ21−1−pn∫0∞x3e−λx1−pe−λx1+λxn−1dx,together with the definition of variance, leads us to Equation (5).Remark 1.
We notice that fork=1, we obtain the complementary exponential distribution introduced by Flores et al. [6].


## 3. Statistical Simulation for the MaxErlB Distribution

Taking advantage of the fact that the random variableUErlB~MaxErlBk,λ,n,p,λ>0, k,n∈1,2⋯, p∈0,1, has the same distribution as the random variable max1≤i≤ZXi, where Xii≥1 are iid random variables, Xi~Erlangk,λ, k∈ℕ,k≥1, λ>0, and the value of the random variable Z~Binom∗n,p, p∈0,1, n∈1,2,⋯, coincide with the value of the random variable zero-truncated binomial distributed with the same parameters, but provided this is a nonzero value, we can briefly describe the following algorithm.


### 3.1. Statistical Simulation Algorithm for the MaxErlB Distribution

Step 1.
We generate a valuez⋆ of the random variable Z⋆~Binomn,p, p∈0,1, n∈1,2,⋯Step 2.
Ifz⋆=0, then GO TO Step 1; otherwise, z=z∗Step 3.
For the valuez of the random variable Z (generated in Steps 1 and 2), simulate the values xi,i=1,2,⋯ as a values of z-iid random variables with distribution Erlangk,λ, k∈ℕ,k≥1, λ>0Step 4.
It is considereduErlB=max1≤i≤zxi, STOP.Following the simulation, we can apply the Chi-square test of concordance. Based on a test, based on the resultsuErlB1,uErlB2,⋯,uErlBm, the Chi-square criterion (Pearson’s criterion) is applied, and the basic and alternative hypotheses are verified, respectively:H0: sample values uErlB1,uErlB2,⋯,uErlBm are values of the random variable distributed MaxErlB2,λ,n,pH1: sample values uErlB1,uErlB2,⋯,uErlBm are not the values of the random variable distributed MaxErlB2,λ,n,p.The test is considered valid if the empirical value ofχc2 is less than the upper critical value of the Chi-square with r−1−L=12−1−4=7 freedom degrees (χ0.05;72=14.067). The statistics of Pearson’s test is calculated using the following relation:
(8)χc2=∑j=1rnj−n0pj2n0pj,where nj,j=1,r¯ represents the number of observed values in the interval tj−1,tj, n0=∑j=1rnj.The probabilitiespj that the random variable UErlB takes the values in the interval tj−1,tj are calculated using the following relation:
(9)pj=UErlBtj−UErlBtj−1=211−1−pn1−pe−λtj∑i=0k−1λtjii!n−1−pe−λtj−1∑i=0k−1λtj−1ii!n,where tj,j=0,r−1¯ represent the ends of each interval after they have been merged.Based on the algorithm presented above, we can notice (see Table2) that the mean and the empirical variance of the simulation results are well approximated by the mean and the theoretical variance of the random variable UErlB~MaxErlB2,10,3,0.2 (Proposition 2), and the Chi-square criterion validates each time the basic hypothesis according to which the simulated values are indeed governed by this distribution.Table 2
The validation of the MaxErlB simulation results with the application of the Chi-square test.
Sample sizeMeanVarianceChi-squareTheoreticalEmpiricalTheoreticalEmpirical1000.21970.20800.02020.01822.06110000.21410.02286.952100000.21510.02175.3171000000.21700.02164.83310000000.21680.02157.636100000000.21670.02142.901Moreover, the validation is confirmed for samples valuesm∈100,1000,10000,100000,1000000,10000000.The histogram of the simulated data and the plot of the probability density function of the simulated distribution (Figure1) also confirm the validity of the basic hypothesis, but visually.Figure 1
Histograms of relative frequencies of samples sizem=1000,10000,100000,1000000 and probability density function of the simulated values that are governed by the MaxErlB2,10,3,0.2 distribution.


## 3.1. Statistical Simulation Algorithm for the MaxErlB Distribution

Step 1.
We generate a valuez⋆ of the random variable Z⋆~Binomn,p, p∈0,1, n∈1,2,⋯Step 2.
Ifz⋆=0, then GO TO Step 1; otherwise, z=z∗Step 3.
For the valuez of the random variable Z (generated in Steps 1 and 2), simulate the values xi,i=1,2,⋯ as a values of z-iid random variables with distribution Erlangk,λ, k∈ℕ,k≥1, λ>0Step 4.
It is considereduErlB=max1≤i≤zxi, STOP.Following the simulation, we can apply the Chi-square test of concordance. Based on a test, based on the resultsuErlB1,uErlB2,⋯,uErlBm, the Chi-square criterion (Pearson’s criterion) is applied, and the basic and alternative hypotheses are verified, respectively:H0: sample values uErlB1,uErlB2,⋯,uErlBm are values of the random variable distributed MaxErlB2,λ,n,pH1: sample values uErlB1,uErlB2,⋯,uErlBm are not the values of the random variable distributed MaxErlB2,λ,n,p.The test is considered valid if the empirical value ofχc2 is less than the upper critical value of the Chi-square with r−1−L=12−1−4=7 freedom degrees (χ0.05;72=14.067). The statistics of Pearson’s test is calculated using the following relation:
(8)χc2=∑j=1rnj−n0pj2n0pj,where nj,j=1,r¯ represents the number of observed values in the interval tj−1,tj, n0=∑j=1rnj.The probabilitiespj that the random variable UErlB takes the values in the interval tj−1,tj are calculated using the following relation:
(9)pj=UErlBtj−UErlBtj−1=211−1−pn1−pe−λtj∑i=0k−1λtjii!n−1−pe−λtj−1∑i=0k−1λtj−1ii!n,where tj,j=0,r−1¯ represent the ends of each interval after they have been merged.Based on the algorithm presented above, we can notice (see Table2) that the mean and the empirical variance of the simulation results are well approximated by the mean and the theoretical variance of the random variable UErlB~MaxErlB2,10,3,0.2 (Proposition 2), and the Chi-square criterion validates each time the basic hypothesis according to which the simulated values are indeed governed by this distribution.Table 2
The validation of the MaxErlB simulation results with the application of the Chi-square test.
Sample sizeMeanVarianceChi-squareTheoreticalEmpiricalTheoreticalEmpirical1000.21970.20800.02020.01822.06110000.21410.02286.952100000.21510.02175.3171000000.21700.02164.83310000000.21680.02157.636100000000.21670.02142.901Moreover, the validation is confirmed for samples valuesm∈100,1000,10000,100000,1000000,10000000.The histogram of the simulated data and the plot of the probability density function of the simulated distribution (Figure1) also confirm the validity of the basic hypothesis, but visually.Figure 1
Histograms of relative frequencies of samples sizem=1000,10000,100000,1000000 and probability density function of the simulated values that are governed by the MaxErlB2,10,3,0.2 distribution.


## 4. EM Algorithm for the MaxErlB Distribution

The EM algorithm introduced in 1977 in the paper [14] comes to perfect the maximum likelihood method which, in the case of processing incomplete statistical data, becomes practically unusable. Next, the algorithm is implemented for the MaxErlB2,λ,3,p,λ>0,p∈0,1 distribution.We consider the values of a samplex1,x2,⋯,xm of size m a statistical population govorned by a MaxErlB distribution with the probability density function uErlBx,Ψ, x>0, which depends on the parameter vector Ψ=λ,p, given that the parameter n of the zero-truncated binomial distribution and the parameter k of the Erlang distribution are given. According to the definition of the maximum likelihood function and Equation (3), we have
(10)Lx1,x2,⋯,xm;Ψ=∏j=1m3pλ2xje−λxj1−pe−λxj−pλxje−λxj21−1−p3=3pλ2me−λ∑j=1mxj1−1−p3m∏j=1mxj1−pe−λxj−pλxje−λxj2.To obtain the maximum likelihood equations for the MaxErlB distribution regarding the estimationλ^,p^ for the parameters λ,p, we consider
(11)lnLx1,x2,⋯,xm;λ,p=mln3+lnp+2lnλ−λ∑i=0k−1xj−mln1−1−p3+∑j=1mlnxj+21−pe−λxj−pλxje−λxj.The parametersn and k being considered known, then the equations of the method for the maximum likelihood estimation function (MLE) are characterized by the nonlinear system SΨ=0, where SΨ=∂lnL/∂λ,∂lnL/∂p. Developing the system of equations SΨ=0, we notice that it becomes difficult to solve in relation to the unknowns λ and p. We are thus in the situation in which the application of the EM algorithm explained and analyzed by Dempster et al. [14], then expanded by McLachlan and Krishnan [15] is required. In this algorithm, the random variable Z is considered a random variable latency, that is, the random variable which cannot be observed directly.For this, we consider, formally, the following sample:(12)x1,z1,x2,z2,⋯,xm,zm,by m observations of the random variable UErlB,Z.This shows thatx1,z1,x2,z2,⋯,xm,zm can be interpreted as a complete set of statistics, being, in this case, a sample of incomplete data. The description of the EM algorithm supposes a known conditional mean EZUErlB;Ψ, where Ψ=λ,p.The probability density functionuErlBx,Ψ, x>0, of the random variable UErlB wich corresponds to a complete set of data, is defined by the following relation according to the definition of probability density in the case of the maximum (see [13], Consequence 2.2):
(13)uErlBx=Θλ2xe−λxA′Θ1−e−λx−λxe−λxAΘ,x>0,In these conditions, the probability density functionuErlBx,z of the random variable UErlB,Z which corresponds to a complete set of data is given by
(14)uErlBx,z;Ψ=zfErlxFErlxz−1ℙZ=z=zazΘzλ2xe−λx1−e−λx−λxe−λxz−1AΘ,whereAΘ=1+Θn−1, Θ=p/1−p, p∈0,1, az=nz, z≤n, fErlx, and FErlx,x>0 are the probability density function, respectively, the cumulative distribution function which has the Erlangk,λ,k∈ℕ,k≥1,λ>0 distribution.Then, the probability density function of the random variableZ conditioned by the random variable UErlB has the following expression:
(15)uErlBzx=uErlBx,zuErlBx=zazΘz−11−e−λx−λxe−λxz−1A′Θ1−e−λx−λxe−λx.Therefore, considering the obvious relation∑z≥1z2azΘz−2=A′′Θ+1/Θ·A′Θ, the conditional mean becomes
(16)EZUErlB;Ψ=∑z=1nz·uErlBzx;Ψ=∑z=1nz2azΘz−11−e−λx−λxe−λxz−1A′Θ1−e−λx−λxe−λx=Θ1−e−λx−λxe−λxA′Θ1−e−λx−λxe−λx∑z=1nz2azΘz−21−e−λx−λxe−λxz−2=Θ1−e−λx−λxe−λx·A′′Θ1−e−λx−λxe−λxA′Θ1−e−λx−λxe−λx+1.SinceZ~Binom⋆n,p∈PSD, k,n∈1,2,⋯ with AΘ=1+Θn−1, Θ∈0,+∞, Θ=p/1−p, p∈0,1, we have
(17)EZUErlB;Ψ=2p1−e−λx−λxe−λx1−pe−λx−pλxe−λx+1.We describe the EM algorithm for theMaxErlB2,λ,3,p distribution as an iterative process of estimating the unknown parameter Ψ=λ,p through Ψh=λh,ph calculated for a few steps h≥1 such that the following condition is satisfied:
(18)maxλh−λh−1,ph−ph−1≤ε,or h=K be accomplished when ε>0 and K represents the number of preset iterations.The steps of the EM algorithm for MaxErlB distribution are the following:Step 1.
We takeλ=λ0, p=p0, λ0>0, p0∈0,1Step 2. (Expectation).
To iterateh, h≥1, we calculate the mean value of zjh−1, j=1,m¯ according to Equation (17) for k=2:
(19)zjh−1=2ph−11−e−λh−1xj1+λh−1xj1−ph−1e−λh−1xj1+λh−1xj+1Step 3. (Maximization).
Through the maximum likelihood estimation (MLE) method, we take into consideration the following sample:(20)x1,z1h−1,x2,z2h−1,⋯,xm,zmh−1,with the maximum likelihood function:
(21)Lx1,⋯,xm,z1h−1,⋯,zmh−1;Ψh−1=∏j=1muErlBxj,zjh−1;Ψh−1=∏j=1mzjh−1xj3zjh−1ph−1zjh−11−ph−13−zjh−11··λh−12e−λh−1xj1−e−λh−1xj1+λh−1xjzjh−1−11−1−ph−13.=λh−121−1−ph−13m∏j=1mzjh−1xj3zjh−1ph−1zjh−11−ph−13−zjh−1··e−λh−1xj1−e−λh−1xj1+λh−1xjzjh−1−1.
Thus, we can find iterationΨh=λh,ph which estimates the parameters Ψ=λ,pStep 4.
We examine Equation (18). If NOT, then GO TO Step 2; otherwise, Ψ≔Ψh, STOP.Given the function(22)lnLx1,x2,⋯,xm,z1h−1,z2h−1,⋯,zmh−1;Ψh−1==2mlnλh−1−mln1−1−ph−13+∑j=1mln3zjh−1+lnzjh−1+zjh−1lnph−1+3−zjh−1ln1−ph−1++lnxj−λh−1xj+zjh−1−1ln1−e−λh−1xj1+λh−1xj,the maximum likelihood equations are characterized by the nonlinear system SΨh−1=∂lnL/∂λh−1,∂lnL/∂ph−1, namely
(23)SΨh−1:2mλh−1+∑j=1mxj2λh−1zjh−1−1e−λh−1xj1−e−λh−1xj1+λh−1xj−xj=0−3m1−ph−121−1−ph−13−3m1−ph−1+∑j=1mzjh−1ph−11−ph−1=0.Table3 shows the results obtained from the implementation of the EM algorithm (described above), in the Octave 1.5.4 GUI programming environment. We must also emphasize that for different sample sizes (m∈100,1000,10000,100000,1000000), we obtain very good approximations of the parameters λ and p that characterize the MaxErlB distribution, when the parameters k and n are known.Table 3
The estimate of the parameter vectorΨ=λ,p of MaxErlB2,λ,3,p distribution by Ψ^=λ^,p^.
Sample size(λ,p)λ^p^h100(1, 0.5)1.0170.49413410001.0300.518144100000.9990.4961521000000.9990.50315010000000.9980.497152


## 5. Application

We will now consider a dataset which represents the remission times (in months) of a random sample of 128 bladder cancer patients. The dataset itself has previously been used in [16–18]. It is summarized as follows: 0.08, 2.09, 3.48, 4.87, 6.94, 8.66, 13.11, 23.63, 0.20, 2.23, 3.52, 4.98, 6.97, 9.02, 13.29, 0.40, 2.26, 3.57, 5.06, 7.09, 9.22, 13.80, 25.74, 0.50, 2.46, 3.64, 5.09, 7.26, 9.47, 14.24, 25.82, 0.51, 2.54, 3.70, 5.17, 7.28, 9.74, 14.76, 26.31, 0.81, 2.62, 3.82, 5.32, 7.32, 10.06, 14.77, 32.15, 2.64, 3.88, 5.32, 7.39, 10.34, 14.83, 34.26, 0.90, 2.69, 4.18, 5.34, 7.59, 10.66, 15.96, 36.66, 1.05, 2.69, 4.23, 5.41, 7.62, 10.75, 16.62, 43.01, 1.19, 2.75, 4.26, 5.41, 7.63, 17.12, 46.12, 1.26, 2.83, 4.33, 5.49, 7.66, 11.25, 17.14, 79.05, 1.35, 2.87, 5.62, 7.87, 11.64, 17.36, 1.40, 3.02, 4.34, 5.71, 7.93, 11.79, 18.10, 1.46, 4.40, 5.85, 8.26, 11.98, 19.13, 1.76, 3.25, 4.50, 6.25, 8.37, 12.02, 2.02, 3.31, 4.51, 6.54, 8.53, 12.03, 20.28, 2.02, 3.36, 6.76, 12.07, 21.73, 2.07, 3.36, 6.93, 8.65, 12.63, and 22.69.Figure2 provides the histogram of relative frequencies of a sample size which characterizes the remission times of bladder cancer, where the curve represents the pdf of the random variable UErlB~MaxErlB2,10,3,0.2 distribution defined by Equation (3).Figure 2
A histogram and probability density function plot for remission times of bladder cancer.


## 6. Conclusion

The conclusions revealed by the present research are related to the study of power series distributions type of a maximum of a sequence of iid random variables which are found in a random number.Also, the distribution of a maximum number of iid random variables through the PSD family, characterized by the number of the random variable in the sequence, was presented in a compact, coherent approach.For this purpose, programs for the statistical simulation of the MaxErlB power series distributions type were developed. The validity of the maximum distributions was performed using Pearson’s test of consistency and is reflected in Table2. Describing the EM algorithm implemented in the GUI Octave 1.5.4 programming environment to estimate the parameters of the MaxErlB distribution is presented in Table 3.A real data sequence on bladder cancer remission times was used to illustrate and compare the histogram of the relative frequencies of remission times and the probability density function plot of the remission time values that are governed by the MaxErlB distribution (Figure2).



---

*Source: 9932729-2021-09-25.xml*