A Bounded Derivative Method for the Maximum Likelihood Estimation on Weibull Parameters

JOURNAL OF L A T E X CL ASS FILES, VOL. 6, NO. 1, J ANUAR Y 2007 1 A Bounded Deri v ation Method for the Maximum Like lihood Estimation on the P arame ters of W eib ull Distrib ution DeT ao Mao, and W enyuan Li, F ell ow , IEEE Abstract —Fo r the basic maximum likelih ood estimating func- tion of th e two parameters W eibull distributi on, a simpl e proof on i ts global monotonicity is giv en to en sure the existence and uniquen ess of its soluti on. The boundary of the function’s first- order deriva tive is defined based on its scale-free pr operty . With a bounded derivativ e, th e possible range of the root of this function can be d etermined. A nov el root-finding algorithm employing th ese established results is proposed accordingly , its con ver gence is pro ved a nalytically as well. Compared with other typical algorithms for this problem, the effici ency of th e proposed algorithm i s also d emonstrated b y n umerical experiments. Index T erms —T wo parameter W eibull Distribution; Maximum Likelihood Estimation; Global Monotonicity; Scale-free Pr operty . I . I N T R O D U C T I O N The W eibull distribution [1] is an impo rtant distribution in reliability and m aintainability analysis?. The estimatio n o n its parameters has been widely discussed, and there are several methodo logical categories for this pa rameter estimation issue [2] [3] [4]. For example , th e graph ic m ethods [5], tran scen- dental equation -solving method applying bifurc ation algorithm [4], maximum likelihood estimation (ML E) meth od [6] [7] [8]. The graph ic meth ods, su ch as W eib ull p robability plotting (WPP) [9], ar e straigh t for ward, but can not give a precise estimation. The transcend ental equatio n-solving method has a closed form, it ca n avoid the com puting expense in the iterativ e compu tation o f the r aw data, but in o rder to solve th e transcenden tal equatio n, advanced mathem atical techniq ues are required . In the MLE-b ased methods, as the b asic estimating equation is not in closed for m, it can b e solved only numer ically . There a re several typic al MLE-based metho ds for solving th is equation, suc h as th e secant m ethod, the b isection meth od and the Newton-Raphson method. Howev er , in both the secant method and the b isection meth od, th e convergence rates are very low; in the Newton-Raphson method , it has to compute both the b asic estimating fu nction an d its d eriv ati ve [10] at each iterative step. In some cases, the Newton-Raphson method cannot en sure convergence [1 1]. Furth ermore, these above MLE-b ased m ethods require either initial values or trial computatio n of the estima ted parameters. W enyuan Li is with Grid Operations, BC Hydro, Burnaby , BC V3N 4X8, Canada . e-mail: weny uan.li@bct c.com. DeT ao Mao is with ECE Department, Univ ersity of British Columbia, V anc ouver , B.C., Canada V6T 1Z4 e-mail: detao m@ece.ubc.ca. In [1 2], the auth or claims that th e existence and u niqueness of the solutio n of th e basic estimating eq uation und er MLE method cannot be assured. Howe ver , in [13] [14], based on Cauchy-Sch warz inequality , the author s h av e given similar proof s of the existence and uniq ueness on th e solution of the MLE-based estimator . In th is paper, w e present a straight forward proof by mathem atics in duction, which is much simpler than both the pr oofs g iv en in [13] [14]. For the basic estima ting f unction, we ha ve p roved th e scale free pro perty of its first or der derivati ve, the b oundar y of this deriv ative can thus be defined. Mo reover , with a bou nded deriv ativ e, a t each iterative step, the p ossible range of th e root of th e basic e stimating functio n can be determined . W e thus propo se a n ovel MLE -based root-find ing algo rithm based on these proper ties, its co mputationa l efficiency and a dvantages are well demon strated b y numerical experiments. The rem ainder o f this paper is organized as follows. I n Section II, a pro of on the mo notonicity of th e b asic estimating function is given. In Section III, the scale-fr ee pr operty of the first or der der i vati ve of this function is p roved, the bo undaries of the function itself a nd its fir st order deriv ativ e are defined . The feasible ran ge of its root th us can be determin ed.In Section IV, by employing these p roved results, a n ovel root- finding algorith m is hence design ed. T he conv ergence of the propo sed algor ithm is p roved. In Sectio n I V, the per forman ce of this p roposed root-fin ding algorithm are dem onstrated by numerical experiments. Conclusions are giv en in VI. I I . T H E G L O BA L M O N OT O N I C I T Y O F T H E B A S I C E S T I M A T I N G F U N C T I O N A. The B asic Estimatin g Function for the T wo P arameter W eibull Distrib ution The density function of the two par ameter W eibull dis- tribution is: f ( x ) = ( k x )( x λ ) k e − ( x λ ) k ( x ≥ 0 , k > 0 , λ > 0) (1) For a sampled data of n observations with the above Eqn .1 as the app licable d ensity function, its likelihood function is L ( x 1 , · · · , x n ; k , λ ) = n Y i =1 ( k x i )( x i λ ) k e − ( x i λ ) k (2) JOURNAL OF L A T E X CL ASS FILES, VOL. 6, NO. 1, J ANUAR Y 2007 2 By MLE method , the following eq uations can be obtained : ∂ ln L ∂ k = n k + n X i =1 ln x i − 1 λ k n X i =1 x k i ln x i = 0 , (3) ∂ ln L ∂ λ = k λ ( − n + 1 λ k n X i =1 x k i ) = 0 . (4) by eliminating λ , we get P n i =1 x k i ln x i P n i =1 x k i − 1 k = 1 n n X i =1 ln x i (5) by th e above Eqn.5, we can get the value of k by related numerical alg orithms. With k determined, λ can be calculated by Eqn.4 as ˆ λ = ( P n i =1 x ˆ k i n ) 1 ˆ k = ˆ k v u u t 1 n n X i =1 x ˆ k i (6) we can therefo re calculate bo th k and λ . Here ˆ k ( ˆ λ ) r efers to maximum likelihood estimators for param eter k ( λ ). By Eqn.5, we can define the basic e stimating function F ( k ) as F ( k ) = P n i =1 x k i ln x i P n i =1 x k i − 1 n n X i =1 ln x i − 1 k (7) Thus F ( k ) = 0 here is defin ed as the basic estimating equation. By mathematics induction , with n ≥ 1 and k > 0 , a simple proo f on the glob al mo notonicity o f F ( k ) , i.e., ∂ F ( k ) ∂ k > 0 can be given in the f ollowing section. B. Pr oof on the glo bal mo notonicity of F ( k ) Since ∂ F ( k ) ∂ k = 1 k 2 +  n X i =1 x k i ln 2 x i n X i =1 x k i − ( n X i =1 x k i ln x i ) 2  ( n X i =1 x k i ) − 2 (8) to p rove ∂ F ( k ) ∂ k > 0 , we n eed on ly to p rove that fo r k > 0 a nd any x i ∈ R + , P 1 ( x i , n, k ) = n X i =1 x k i ln 2 x i n X i =1 x k i − ( n X i =1 x k i ln x i ) 2 ≥ 0 for n = 1 P 1 ( x i , 1 , k ) = x k 1 ln 2 x 1 · x k 1 − x 2 k 1 · ln 2 x 1 = 0 (9) for n = 2 P 1 ( x i , 2 , k ) = x k 1 x k 2 (ln x 2 − ln x 1 ) 2 ≥ 0 (10) for n = m − 1 ( m ≥ 3 , m ∈ N + ), supp ose P 1 ( x i , m − 1 , k ) ≥ 0 (11) then for n = m P 1 ( x i , m, k ) = P 1 ( x i , m − 1 , k ) + x k m m − 1 X i =1 x k i (ln x m − ln x i ) 2 ≥ 0 (12) Therefo r ∂ F ( k ) ∂ k > 0 , and f unction F ( k ) is global monoton ic. W ith the global monoto nicity of F ( k ) , the existence and uniquen ess of the r oot of F ( k ) = 0 can be assured. C 1 -C 1 F(k) F U (k)=C 1 -1/k x F L (k)=-1/k k C 1 Fig. 1. The boundarie s of F ( k ) : − 1 k < F ( k ) < − 1 k + C 1 . Here C 1 = 1 n P n − 1 i =1 ln x max x i , k > 0 . I I I . B O U N D A R I E S O F T H E B A S I C E S T I M AT I N G F U N C T I O N A N D I T S F I R S T - O R D E R D E R I V A T I V E A. Bou ndaries of th e Basic Estimating Func tion F ( k ) Let G ( k ) = P n i =1 x k i ln x i P n i =1 x k i − 1 n n X i =1 ln x i (13) as    G (0) = 0 G (+ ∞ ) = 1 n P n − 1 i =1 ln x max x i > 0 also because ∂ G ( k ) ∂ k = P 1 ( x i , n, k )( n X i =1 x k i ) − 2 ≥ 0 , (14) thus ∀ k > 0 − 1 k ≤ F ( k ) ≤ C 1 − 1 k (15) here C 1 = G (+ ∞ ) = 1 n n − 1 X i =1 ln x max x i (16) we define the lower bo undary of F ( k ) is curve: F 2 ( k ) = − 1 k , th e up per boun dary of F ( k ) is cu rve: F 1 ( k ) = − 1 k + C 1 . Then the bou ndary c urves of F ( k ) can be seen in Fig. 1. B. Bou ndaries of th e F irs t Order derivative: ∂ F ( k ) ∂ k Let H ( n, x i , k ) =  n X i =1 x k i ln 2 x i n X i =1 x k i − n X i =1 x k i ln x i n X i =1 x k i ln x i  ( n X i =1 x k i ) − 2 (17) Theorem A : ∀ n ∈ N + , ∀ x i > 0 , and ∀ k > 0 , ∂ F ( k ) ∂ k ∈ [ k − 2 , k − 2 + ln 2 ( x max x min )] . JOURNAL OF L A T E X CL ASS FILES, VOL. 6, NO. 1, J ANUAR Y 2007 3 PR O OF: For a certain λ j > 0 satisfying ln( λ j · x i ) ≥ 0 , H ( n, x i , k ) = H ( n, λ j · x i , k ) ≤ max { H ( n, λ j · x i , k ) } = ln 2 ( λ j · x max ) − ln 2 ( λ j · x min ) = ln x max x min ln( λ 2 j · x max · x min ) (18) Since to satisfy    ln( λ j · x i ) ≥ 0 ln( λ 2 j · x max · x min ) ≥ 0 the minimu m value of λ j is x − 1 min , i.e., λ j ∈ [ x − 1 min , + ∞ ] . let λ = γ · x − 1 min ( γ ≥ 1 ) , max n ∈ N + ,x i > 0 , k> 0 ,γ ≥ 1 { H ( n, λ · x i , k ) } = ln( x max x min ) ln( γ 2 x max x min ) (19) therefor e the super mum value of H ( n, x i , k ) is sup { H ( n, x i , k ) } = min { max x i > 0 ,n ∈ N + , k> 0 ,γ ≥ 1 { H ( n, λ · x i , k ) }} = min γ ≥ 1  ln( x max x min )[ln γ 2 + ln( x max x min )]  = ln 2 ( x max x min ) Hence H ( n, x i , k ) ≤ ln 2 ( x max x min ) (20) It is proved in Section II-B, the infimum value of H ( n, λx i , k ) is: inf { H ( n, x i , k ) } = 0 (21) with ∂ F ( k ) ∂ k = H ( n, x i , k ) + 1 k 2 , therefo re 1 k 2 ≤ ∂ F ( k ) ∂ k ≤ C 2 + 1 k 2 (22) here C 2 = ln 2 ( x max x min ) .  W ith the boun dary of ∂ F ( k ) ∂ k defined, the feasible solu tion range of F ( k ) = 0 can b e determined. C. P ossible Ran ge of the F in al Solutio n of Equatio n F ( k ) = 0 The basic estimating fu nction F ( k ) c an be written as F ( k ) = Z k k 0 F ′ ( τ ) dτ According to Theor em A,      R k k 0 τ − 2 dτ < R k k 0 F ′ ( τ ) dτ < R k k 0 [ C 2 + τ − 2 ] dτ ( k > k 0 ) R k k 0 τ − 2 dτ > R k k 0 F ′ ( τ ) dτ > R k k 0 [ C 2 + τ − 2 ] dτ ( k < k 0 ) let denote      F L ( k , k 0 ) = R k k 0 τ − 2 dτ F U ( k , k 0 ) = R k k 0 [ C 2 + τ − 2 ] dτ (0,0) (k 0 ,F(k 0 )) k F(k) F L (k,k 0 ) F U (k,k 0 ) k L 0 k U 0 Fig. 2. The two boundary curves of F ( k ) : F L ( k, k 0 ) and F U ( k, k 0 ) , which can dete rmine the feasible range of the solutio n of F ( k ) = 0 . thus    F L ( k , k 0 ) < F ( k ) < F U ( k , k 0 ) ( k > k 0 ) F L ( k , k 0 ) > F ( k ) > F U ( k , k 0 ) ( k < k 0 ) Here F U ( k , k 0 ) and F L ( k , k 0 ) can be seen the bo undary curves of F ( k ) at poin t ( k 0 , F ( k 0 )) ( see Fig.2). W ith knowing F ( k 0 ) , th e two bound ary curves of F ( k ) can be d etermined as,    F L ( k , k 0 ) = k − 1 0 + F ( k 0 ) − 1 k F U ( k , k 0 ) = C 2 k − 1 k + ( F ( k 0 ) − C 2 k 0 + 1 k 0 ) W e can see that both F U ( k , k 0 ) and F L ( k , k 0 ) are alg ebraic function s with simp le m athematical forms. Suppose F U ( k , k 0 ) cuts k -axe at po int k 0 U , F L ( k , k 0 ) cuts k -axe at point k 0 L , in Fig .2, it is obvious that the root o f equation F ( k , k 0 ) = 0 must exist in [ k 0 L , k 0 U ] (if k 0 L < k 0 U ) o r in [ k 0 U , k 0 L ] (if k 0 L > k 0 U ). A deta iled algorith m employing th e prope rties in Sectio n III-A, III-B a nd III-C will be proposed in the following section , and a p roof on its con vergence will be demo nstrated as well. I V . R O OT - FI N D I N G A L G O R I T H M D E S I G N T o nu merically get the solution of F ( k ) = 0 , an in tuitiv e idea is to c alculate F ( k 1 ) at k 1 = 1 2 ( k 0 L + k 0 U ) , to f urther find out a n arrower interval [ k 1 U , k 1 L ] or [ k 1 L , k 1 U ] . By the same proced ure, at k i +1 = 1 2 ( k i L + k i U )( i → + ∞ ) , the fin al solution of F ( k ) = 0 can be iteratively app roximated . From the m athematical form of the basic estimating function F ( k ) , we can see that for large n, k and x i , the c omputation al complexity for F ( k ) is very high, thu s the calc ulation time of F ( k ) is a n imp ortant index in measurin g the efficiency of a r oot-findin g algorithm for F ( k ) = 0 . A novel root-fin ding algorithm h as been d esigned based on these prope rties proved in Section II , II I and III-C. A. Conver gence of the Boun ded deriva tive Algorithm As stated befo re, a straigh t forward m ethod employing th is idea is to calculate the final solution iteratively . The feasibility of this method can be e nsured by the following theorems. JOURNAL OF L A T E X CL ASS FILES, VOL. 6, NO. 1, J ANUAR Y 2007 4 k i F( K i ) i+1 ( k i ,F( k i )) F U ( K i+1 ) k U i Case I: F( k i+1 )> 0 Case II: F( k i+1 )< ! k L i+1 k U i+1 k L i+1 k U i+1 F L ( K i+1 ) k L i F L ( K, K i ) F U ( K, K i ) Fig. 3. Estimation on the con vergenc e rate of the proposed algorithm with kno wing F ( k i ) and F ( k i +1 ) , here F L ( k, k i ) and F U ( k, k i ) refer to the boundary functi ons of F ( k ) at point ( k i , F ( k i )) , k i +1 = 1 2 ( k i L + k i U ) . Lemma . For any point ( k i , F ( k i )) of func tion F ( k ) , the existence an d u niqueness of the r oots of F U ( k , k i ) = 0 and F L ( k , k i ) = 0 can be assured , th erefore the interval [ k i L , k i U ] ( i ≥ 1) th at covers the ro ot of F ( k ) = 0 alw ays exists. PR O OF: (1 ) By Inequ ality 15, we know that ∀ k > 0 , F ( k ) > − 1 k . Th e roo t of equation F L ( k , k 0 ) = k − 1 0 + F ( k 0 ) − 1 k = 0 is k 0 L = 1 F ( k 0 ) − ( − k − 1 0 ) > 0 (23) which is reasonab le. (2) Since F U ( k , k 0 ) = 0 can be rewritten as C 2 k 2 + ( F ( k 0 ) − C 2 k 0 + 1 k 0 ) k − 1 = 0 (24) Since here  ∆ = B 2 − 4 AC = ( F ( k 0 ) − C 2 k 0 + 1 k 0 ) 2 + 4 C 2 > 0 C = − 1 < 0 it is obvio us th at F U ( k , k 0 ) = 0 always has two roo ts with different signatu res. T o satisfy k > 0 , on ly the positive ro ot should be preser ved. Therefo re, the existence o f the inter val [ k i L , k i U ] ( i ≥ 1 ) can be assured.  Theorem B . Let [ k i L , k i U ] d enote th e feasible interval of the solution of F ( k ) = 0 at point ( k i , F ( k i )) , [ k i +1 L , k i +1 U ] denote the feasible interval at point ( k i +1 , F ( k i +1 )) , with k i +1 = 1 2 ( k i L + k i U ) , then the con vergence rate γ = || k i +1 L − k i +1 U || || k i L − k i U || < 1 2 . Here k i +1 L = max( k i L , k i +1 L ) , k i +1 U = min ( k i U , k i +1 U ) , i ∈ N + . PR O OF: W ith kn owing p oint ( k i , F ( k i )) ( i ≥ 1 ) , we can deduce the bou ndary cu rves: F L ( k , k i ) and F U ( k , k i ) , an d the related fe asible interval [ k i +1 L , k i +1 U ] can b e determin ed accordin gly (as shown in Fig.3). Since if F ( k i +1 ) = 0 , then K is the roo t we are loo king for . Thus here there are only two situations f or F ( k i +1 ) : F ( k i +1 ) > 0 and F ( k i +1 ) < 0 . Case I : F ( k i +1 ) > 0 Fig. 4. The flow chart for the combined algorithm. Here δ 1 is a scale threshold , δ 2 is the approximat ing precision or the computin g halt crite rion. W ith kn owing F ( k i +1 ) ( i ≥ 1) , the bou ndary curves F L ( k , k i +1 ) and F U ( k , k i +1 ) at point ( k i +1 , F ( k i +1 )) ca n be determined , and so is the related f easible interval [ k i +1 L , k i +1 U ] . As in this case, by Fig. 3, it is obvious that k i +1 U < K i +1 , then the conv ergence rate γ = k i +1 U − k i +1 L k i U − k i L ≤ k i +1 U − k i L k i U − k i L < k i +1 − k i L k i U − k i L = 1 2 (25) Case II : F ( k i +1 ) < 0 Similarly , as shown in Fig.3, in th is ca se, it is obvious that k i +1 L > K i +1 , then the convergence rate γ = k i +1 U − k i +1 L k i U − k i L ≤ k i U − k i +1 L k i U − k i L < k i U − k i +1 k i U − k i L = 1 2 (26) Thus the prop osed algo rithm is co n vergence, a nd its con- vergence rate γ < 1 2 .  B. Impr oved Algo rithm Combining the Seca nt Method an d the Bound ed derivative Method It is obviou s that at large scale the bo unded d eriv ati ve algorithm can co n verge rap idly , at least at th e same r ate as the bisectio n method. While in a small scale, e specially in the linearizab le n eighbor hood aroun d the fin al solution poin t of the basic e stimating equa tion, the secant m ethod method has a better conver gence r ate. Th erefore it is in tuiti ve to combine both the two me thods tog ether . The flow chart for the combined algor ithm ca n be seen in Fig.4. As shown in the flow chart in Fig.4, in the comb ined algorithm , at a large scale, at ea ch itera ti ve step, w e calculate F ( k i ) o ne tim e, then with o nly b asic algeb raic calcu lation, we find the possible range of the solution, and then nar row down it to a smaller rang e at next step. When the p ossible rang e of the roo t comes to a small scale, where th e linearization part of the basic estimating fu nction d ominates, th e co mbined algorithm will switch to the secant method, since it has a better perfor mance in this case. JOURNAL OF L A T E X CL ASS FILES, VOL. 6, NO. 1, J ANUAR Y 2007 5 Moreover , c ompared with the secan t metho d, the bisection method and th e Newton-Raphson method , which r equires initial values of F ( k ) , no initial value is r equired in the combined method. V . N U M E R I C A L E X A M P L E S A. General Cases W ith MA TLAB, 1000 times of numer ical experiments have been taken to study the p erform ance of this p roposed algo- rithm. In these simulations, the shape parameter of W eibull distribution k ∼ U (0 , 40] , the scale p arameter λ ∼ U (0 , 40] , the nu mber of data N ∼ U [2 , 1000] , here U refers to Un iform Distribution. Compared with the secan t me thod, the bisection method as well as the bi-secan t metho d ( which has combine d the secant and the bisection meth od), the compu ting complexity (here refers to the calculation time of F ( k ) ) of each meth od under various appro ximating prec ision ǫ can be seen in the following T able I. T able I: A verag ed c alculation times of F ( k ) Appro- Secant Bisection Bi-secant Pro posed precision ǫ method method[ 6 ] Method method 10 − 1 9.42 5.21 3.12 1.6 8 10 − 2 30.22 6.33 4.29 2.5 8 10 − 3 88.5 11.22 9.39 4.2 8 10 − 4 130.6 14.19 11.51 5.05 B. A Con structed Case The advantage of the prop osed algor ithm can also be verified b y a c oncrete case, which is gener ated from a W eibull Distribution. The samp led data is in T ab le II, the comp arison of the performan ce of different method s can b e seen in T able III. T able II :Sampled Data from a W eibull Distribution 2.614 4 4 .1834 4.325 8 4 .3496 4.374 0 4.40 06 3.207 3 4 .2573 4.327 3 4 .3544 4.382 8 4.40 51 3.980 0 4 .2884 4.333 4 4 .3646 4.387 3 4.41 23 4.176 7 4 .3150 4.340 3 4 .3698 4.395 9 4.41 94 4.431 7 4 .4919 4.444 8 4 .5082 4.462 3 4.54 39 4.475 6 4 .5715 T able II I: Calculation times of F ( k ) o f various method s Appro- Secant Bisection Bi-Secant Propo sed precision ǫ method method Method method 10 − 1 4 4 3 1 10 − 2 79 6 3 1 10 − 3 178 10 3 2 10 − 4 271 13 21 3 10 − 6 459 20 22 4 10 − 10 833 33 23 5 10 − 14 1206 46 24 6 V I . C O N C L U S I O N In this pa per , to assure th e un iqueness and existence of the solution of the basic estimating fu nction, we hav e proved its global monoto nicity in a very simp le way . W ith provin g the scale-free prop erty of this function ’ s fir st orde r de riv ati ve, the possible range of its solution at ea ch iterativ e step can b e determined . Based o n these pr operties, a novel roo t-finding algorithm is pr oposed in this pape r . Its efficiency h as b een demonstra ted by num erical exper iments. For fu ture work, this method can b e extended to type I an d type II data cen soring cases. R E F E R E N C E S [1] W . W eib ull, “ A statistic al distrib ution function of wide appli cabilit y , ” J ournal of Mecha nics , vol. 18, pp. 293–297., 1951. [2] A. Cohen, T runcated and Censor ed Samples: Theory and Applicat ions . Dekke r , N e w Y ork, 1991. [3] L. Meeker , W .Q. Escobar , Statist ical Methods for Reliabil ity Data . W ile y , New Y ork, 1998. [4] W . Li, Risk Assessment Of P ower Systems: Models, Methods, and Applicat ions . W ile y-IEEE Press, April 2005. [5] N. R. Mann, R. E . Schafer , and N. D. Singpurwalla , Methods for statisti cal analysi s of reliabil ity and life data . John W ile y and Sons, Ne w Y ork, 1974. [6] A. C. Cohen, “Maximum like lihood estimation in the weibull distrib u- tion based on complete and on censored samples, ” T echnomet rics , vol. 7(4), pp. 579–588., 1965. [7] H. L . Harter and A. H. Moore, “Point and interv al estimators, based on m order statisti cs, for the scale paramete r of a weibull populati on with kno wn shape parameter , ” T ec hnometrics , vol. 7(3), pp. 405–422, Aug, 1965. [8] ——, “Maximum-likel ihood estimation of the parameters of gamma and weib ull populations from complete and from censored data, ” T echn omet- rics , vol. 7(4), pp. 639–643, 1965. [9] R. Ross, “Graphica l m ethods for plotting and e va luating weibull dis- trib uted data, ” IEEE T rans Dielec tr Electr Insul , vol. 1(2), pp. 247–253, 1994. [10] H. Qiao and C. P . Tsokos, “Parameter estimat ion of the weibu ll proba- bility distributio n, ” Mathematics and Compute rs in Simulat ion , vol. 37, pp. 47–55, 1994. [11] L. P . Gupta, R. C. Gupta, and S. J. Lvin, “Numerical methods for the maximum likeli hood estimati on of weibull parameters, ” Journal of Statist ical Computati on and Simulation , vol. 62(1), pp. 1–7, 1998. [12] B. Dobson, The W eibull Analysis Handbook, 2nd ed. ASQ Quality Press, Mil wauk ee, 2006. [13] N. Balakrishna n and M. Kateri, “On the maximum likel ihood estimati on of parameters of weibull distrib ution based on complete and censore d data, ” Statistic s and Probab ility Letters , vol. 78, pp. 2971–2975, 2008. [14] N. Farnu m and P . Booth, “Uniquene ss of maximum lik elihood estimators of the 2-parameter weibull distributio n, ” IEEE T ransactions on Relia- bility , vol. 46(4), pp. 523–525, 1997.