The Development of Parameter Estimation on Hazard Rate of Trivariate Weibull Distribution

In this study, the interrelation concepts of trivariate distribution function, trivariate survival function, trivariate probability density function and trivariate hazard rate function of trivariate Weibull distribution are presented. The goal of this contribution is to estimate the trivariate Weibull hazard rate parameters. To reach this goal, we will use an analitical approach in estimating called the Maximum Likelihood Estimation (MLE) method. Using numerical iterative procedure the scale parameters, the shape parameters and the power parameter estimators on trivariate hazard rate of trivariate Weibull distribution must be obtained. The MLE technique estimates accurately the trivariate Weibull hazard rate parameters.


INTRODUCTION
The Weibull distribution is characterized in the class of absolutely continuous distributions by two parameters, one is the shape parameter, say γ and the other is the scale parameter, say η. The Weibull distribution plays a central role in the analysis of survival or failure time data. For survival data, the proportional hazard model is the most popular model. Moreover, the Weibull distribution is good to be used in parameterized the proportional hazard model because of its flexibility, allowing increasing, decreasing or constant hazard rate function.
The literature on parameter estimation for the Weibull distribution is vast. So, we will select only some key papers. The first study, Tuerlinckx (2004), he derived the Laplace tranform of the positive stable density. His proposition had been applied as a useful vehicle to derive a multivariate survival function of Weibull distribution.
The other papers based on analytical methods of parameter estimation are Hanagal (2005), he propose a new bivariate Weibull regression model on censored samples with common covariates and obtaining the maximum likelihood estimators for parameters, Hanagal (2006) developed a maximum likelihood estimation procedure for a bivariate frailty regression model. Lee and Wen (2009) propose a multivariate Weibull model and derived the explicit form of PDF, CDF and general moment.
In this study, we propose a new parameter estimation on trivariate Weibull distribution related to the consideration of hazard rates. The plan of this study is organized as follows. We will introduce the Weibull distribution and its related functions, next we derive estimators for the parameters of the proposed model and the lastly, we present the results, discussion and deduce a few conclusion. Unfortunately, at this stage of the investigations there are no real data for iterating our proposed parameter estimation.

MATERIALS AND METHODS
Let T be a continuous nonnegative random variable representing the failure time of an individual from a homogeneous population with probability density function f T (t) and cummulative distribution functon F T (t) = P(T ≤ t). The complementary cumulative distribution AJBS Science Publications function S T (t) = 1 -F T (t) = P(T > t) is called survival function. S T (t) is the probability of surviving an age of t.
The probability density fuction f T (t) of random variable T is then equal to dF T (t)/dt or equivalenly-dS T (t)/dt.
The univariate hazard rate h T (t) of failure time T is defined as the probability of failure during a very small time interval, assuming that the individual has survived to the beginning of the interval: for all t such that S T (t)>0. Another expression of h T (t) as defined above is h T (t) = f T (t)/S T (t) (Navarro, 2007). This hazard rate plays an important role in survival data analysis. Correspondingly, the cumulative hazard rate function of failure time T is defined as: hence, the link between the cumulative hazard rate H T (t) and the survival function S T (t) is done by the following relationship: this equation has been noted by (Singpurwalla, 2006). Unlike the univariate setup, there are various definitions of multivariate hazard (failure) rate functions. One can define the multivariate hazard rate in more than one way. In the case of bivariate hazard rate function (for overviews, see Kundu and Gupta, 2010a).
Let (T 1 ,T 2 ) T be a continuous nonnegative random vector with the joint survival function S T1 , T2 (t 1 ,t 2 ) = P(T 1 >t 1 ,T 2 >t 2 ), the joint distribution has the form F T1 , T2 (t 1 ,t 2 ) = P(T 1 ≤t 1 ,T 2 ≤t 2 ). This distribution can be defined in terms of survival functions as: If we view the density as the unconditional failure rate, we can define a conditional failure rate as being the same quantity after having accounted for the fact that the individual has already survived until the time point t 1 and t 2 . Then, the bivariate hazard rate function of the random variable T 1 ,T 2 can be written as a scalar quantity, given by: The other approach in defining bivariate hazard rate is based on a concept of vector-valued bivariate hazard rate as was stated in Kundu and Gupta (2010b.), that is: AJBS Similarly, in the case of trivariate, the cumulative distribution function of a continuous nonnegative random vector (T 1 ,T 2 ,T 3 ) T is defined by Eq. 1: is the joint survival function of the two failure time random variables T i ,T j with 1≤i<j≤3 and S Ti (t i )=P(T i >t i ) is the marginal survival function of T i .
The relationship between the joint survival function S T1,T2,T3 (t 1 ,t 2 ,t 3 ) and the joint cumulative distribution function F T1,T2,T3 (t 1 ,t 2 ,t 3 ) can also be written as: Which implies that if the cumulative distribution function (1) is absolutely continuous, then the joint probability density function f T1,T2,T3 (t 1 ,t 2 ,t 3 ) has the following expression Eq. 2:
The joint survival function for the trivariate Weibull distribution of random variable T 1 , T 2 and T 3 proposed by Lee and Wen (2009) for 0 < α ≤ 1; 0 ≤ t 1 , t 2 , t 3 < ∞ is given by Eq. 3: where T 1 , T 2 ,T 3 are three failur times with univariate survival functions S Ti (t i )=P(T i >t i ). Each univariate survival function has a shape parameter γ i and a scale parameter η i . The parameter α represents the degree of dependence in the association of T 1 , T 2 and T 3 . The case α = 0 and α = 1 correspond to maximal positive dependence and independence, respectively. It follows directly from (2) and (3) we have the joint density function of trivariate Weibull distribution and it takes the form Eq. 4: where,C(m,k,α) is the generalized factorial coefficient, which for all positif integer m, k with k ≤ m and for real number α is defined by: Furthermore, the summation is extended over all partitions of m into k parts, that is over all nonnegative integer solutions of the equations k 1 +2k 2 +...+mk m = m, k 1 + k 2 +…+ k m = k. This implies that expanding C(3,k,α) causes the expression (4) may be written, equivalently, as Eq. 5: By introducing the relationship of the trivariate Weibull survival function and the trivariate Weibull density function above we have the trivariate expression of the Weibull hazard rate Eq. 6: for all (t 1 ,t 2 ,t 3 ) such that S T (t)>0, or in the concept of vector-valued trivariate hazard rate, it is defined as Eq. 7: where for ℓ = 1,2,3 Suppose that there are n individuals under investigation and for each individual there are L failure types. Let T jℓ denote the exact failure times for the ℓ-th failure type of the j-th individual (j = 1 ,…, n; ℓ = 1 ,…, L). If their failure times follow the trivariate Weibull distribution and we assume that the observations are n independent and identically distributed (i.i.d), the likelihood function is Eq. 8: where, θ = (γ 1 , γ 2 , γ 3 , η 1 , η 2 , η 3 , α) T denotes a vector of unknown parameters of trivariate Weibull distribution,  (8). Taking natural logarithm on both sides of (8) and denote it by L (θ), we have Eq. 9: Note that the equation (9) is known as the log likelihood function of θ.
Maximization of (8) is rarely done by some procedure of direct optimization, but usually by some methods based on derivatives of (8). The process of forming these derivatives is made easier by departing from the log likelihood function (9) which is a sum instead of a product, as the logarithmic transformation is isotonic, the extremal points î θ of L (θ) and L(θ) will be the same.
It is clear that if vector θ is a solution of the following simultaneous equations, which are obtained by taking the partial derivative of (9) with respect to each elements of vector θ and setting the result to zero Eq. 10: Then θ is a maximizer of (9). Note that the equation (10) above is known as the log likelihood equation and the 7×1 vector ∂L(θ)/∂θ in the left hand side of (10) is called the score vector and denoted by Let L(θ) be a scalar-valued function of a vector θ as in Eq. 9. The second derivative of L(θ) with respect to vector θ is a matrix of the partial derivatives of L(θ) with respect to the elements of the vector θ. This matrix is called the Hessian and is denoted by H(θ): where the elements of H(θ) in (11) are ∂ 2 L(θ)/∂θ u ∂θ v , u,v = 1,...,7. The Hessian matrix (11) plays a crucial role in solving estimating Eq. 10 with a Newton Raphson algorithm. The MLEs of the trivariate Weibull parameters cannot be given in closed form since Eq. 10 is a system of interdependent non linear equation. Hence, we have to apply iterative methods, carried through to convergence or terminated after reaching some given stopping criterion, to calculate or approximate the MLE resulting in a so called iterated MLE. The MLEθ is the root to the estimating (10), which can be solved by newton raphson algoritm. AJBS

Science Publications
Maximization of (9) directly is consuming time. For simplicity in computation equation (10), without loss of generality, we have to split (9) into four parts: The results of analytical calculation for the four functions L i (θ) ( i = 1,…, 4) above can be described as follows.

For the First Part
First of all, by recalling that θ is a sevendimensional vector of unknown parameters consisting of three shape parameters, three scale parameters and one power parameter, we see that each of Eq.12a-12d is a function of the seven unknown parameters. Moreover, notice that throughout this study, the subscript k or ℓ indicates the k or ℓ -th failure type and therefore the range of k and ℓ are always {1, 2, 3}. Now, let g 1 (θ) be defined as in the left hand side of (10). The first partial derivative of log likelihood L 1 (θ) in (12a) with respect to θ, termed score vector g 1 (θ), is: where, g 1 (θ) in (13) is 7×1 score vector for θ, which elements of the score vector may be obtained as ∂L 1 (θ)/∂θ = (∂L 1 (θ)/∂γ T , ∂L 1 (θ)/∂η T , ∂L 1 (θ)/∂α) T . By applying the chain rule in deriving (12a) with respect to θ, we have the following explicit expression: Let H 1 (θ) be the 7×7 Hessian matrix as given in (11), the second partial derivative of the log likelihood function L 1 (θ) in (12a) with respect to θ is a matrix of the form: where the elements of H 1 (θ) in (14) are ∂ 2 L 1 (θ) /∂θ u ∂θ v , u,v = 1,...,7. The k ℓ th entries of Hessian matrix (14) , for k, ℓ = 1, 2, 3 with k ≠ ℓ , can be derived easily as: Furthermore, for k = ℓ , the diagonal entries of Hessian matrix (14) can be displayed as: In addition, some other off-diagonal entries of Hessian matrix (14) can also be obtained, i.e., for k, ℓ = 1, 2 and 3:
By combining the derivative quotient rule and repeated use of the chain rule on (12b) with respect to θ, we have all elements of the mixed second order partial derivatives symmetric matrix (∂ 2 Q j /∂θ∂θ T ) or (∂Q j /∂θ T ) in (18) Where:

Science Publications
The diagonal elements of a square matrix (∂Q j /∂θ T ) in (18) which have been formulated more explicitly as in Eq. 20 consist of ∂ 2 Q j /∂α 2 itself and the diagonal elements of the following matrices: (∂ 2 Q j /∂γ∂γ T ) and (∂ 2 Q j /∂η∂η T ).

For the Third Part
Considering Eq. 12c that the first order partial derivatives of L 3 (θ) with respect to θ can be expressed as: Note that g 3 (θ) in (23) is the 7×1 score vector for θ, which elements of the score vector may be obtained as ∂L 3 (θ)/∂θ = (∂L 3 (θ)/∂γ T , ∂L 3 (θ)/∂η T , ∂L 3 (θ)/∂α) T . Therefore by applying the chain rule in deriving (12c) with respect to θ, the elements of the score vector g 3 (θ) above can be done easily: where, Vj * = Uj −K j ln (K j ). Using similar arguments to (14) and (17), symmetric Hesian matrix H 3 (θ) of second order partial derivatives of L 3 (θ) can be represented as: The diagonal entries of Hessian matrix (24) are found by deriving the right-hand side of Eq. 24 and take the form: Similarly, the k ℓ th entries of Hessian matrix (24), for k, ℓ = 1, 2, 3 with k ≠ ℓ , can be derived to be: After some straightforward calculation, the other entries of Hessian matrix (24) are accomplished, i.e., for ℓ = 1, 2, 3 satisfy:

For the Fourth Part
Let g 4 (θ) be defined in accordance with (13), (15) and (23). Then, for the L 4 (θ) as given in (12d), we get the score vector, which is written as: By employing the right-hand side of (25), the elements of score vector g 4 (θ) above can be easily calculated:

Science Publications
As before, we let H 4 (θ) be defined similar to (14), (17) and (24). Then, for the L 4 (θ) as already mentioned in (12d), we have: By differentiating L 4 (θ) in Eq. 12d two times with respect to the same elements of vector θ, the diagonal entries of Hessian matrix (26) can be found immediately and we have the following results: In a similar manner to those above, we can treat the off-diagonal entries of Hessian matrix (26), for k, ℓ = 1, 2, 3 with k ≠ ℓ , to that off-diagonal entries and write them in the form of: Analogously, for k = ℓ = 1, 2, 3, it follows that with a simple algebraic manipulation yields the rest entries of the Hessian matrix (26) (13), (15), (23) and (25) we get the score vector of log likelihood function (9) as below Eq. 27: Since the score vector and the Hessian matrx of the log likelihood are now known, the Newton Raphson or other gradient-based methods can be applied to calculate θ that optimizes L(θ) in Eq. 9.

RESULTS
Let θ be the maximizer of (9). When we insert the MLE θ into the survival function  (3), then we will get the MLE of S T (t j ) and an estimator of K j of the j-th individual, because the MLEs are functional invariant: i.e., whenθ is the MLE of θ, h(θ )will also be the MLE of h(θ) if h(θ) is a finite function and need not be one to one. When we substitute θ to the Hessian matrix (11) we will get the estimated covariance matrix of the MLE θ , that is given by Putting the estimatorsθ and j ⌢ K into (6), we obtain the scalar quantity trivariate Weibull hazard rate estimator, that is: This hazard rate estimator is useful for characterizing trivariate exponential distributions as the special cases of trivariate Weibull distributions.