Parameter estimation of the Weibull probability distribution

ATGEMATICS AND COMPUTERS N SIMULATION Mathematics and Computers in Simulation 37 (1994) 47-55

Parameter

estimation of the Weibull probability distribution Hongzhu

Qiao ‘, Chris P. Tsokos

*

Department of Mathematics, University of South Florida 4202 East Fowler Avenue, PHY 114, Tampa FL 33620-5700 USA

Abstract Newton-Raphson’s method plays a fundamental role in the maximum likelihood estimation of the two parameters of the Weibull probability distribution. It is well known that the method depends on the initial point of the iterative process and the iteration does not always converge. In the present article, we introduce an effective iterative procedure for estimating the subject parameter. Our proposed method always converges and does not depend on the initial point of the iteration. It will also be shown by numerical examples that our method converges faster than the popular Newton-Raphson method. Keywords: Two parameter Simple Iterative Procedure

Weibull distribution; Three parameter (SIP); Maximum likelihood estimation

Weibull (MLE)

distribution;

Newton-Raphson

method;

1. Introduction

The Weibull probability density function is an important probability distribution in characterizing the probabilistic behavior of a large number of real world phenomena. It is especially useful as a failure model in analyzing the reliability of different types of systems. The complete Weibull probability distribution for a random variable X contains three parameters, given by

f(x; a, b,

c>=

c(x -

a)C-l

bc

exp{ -( 7)‘)

(1.1)

where x > a, b > 0, c > 0. The three parameters a, b and c are the location, scale, and shape parameters, respectively. We shall denote the three parameter Weibull distribution by W(a, b, c>. 1 Present address: Dept. of Mathematics * Corresponding author.

and Physics, Fort Valley State College.

0378-4754/94/$07.00 0 1994 Elsevier Science B.V. All rights reserved 0378-4754(94)00055-O

SSDI

48

H. Qiao, C.P. Tsokos /Mathematics

and Computers in Simulation 37 (1994) 47-55

Successful application of the Weibull distribution, such as in reliability analysis, hazard rate, etc., depends primarily on having good estimates of its parameters. There are several studies that have investigated the analytical and numerical aspects of estimating the parameters, see Rockette [3], McCool [2], Bain and Engelhardt [l], among others. In the present study, this part of our research, we will focus our attention on the two parameter Weibull model. That is, a will be assumed to be known. In this case, in order to find the maximum likelihood estimation of the two parameters, we have to solve the nonlinear likelihood equations. Most studies simply use the popular Newton-Raphson method as the numerical procedure to estimate the parameters, ignoring the fact that the Newton-Raphson procedure is unstable and may diverge when one chooses inappropriate initial points. In view of this difficulty, we have developed a new method to solve the two parameter Weibull problem. It takes less computation time for than the fast speed Newton-Raphson method. In addition our proposed procedure converges for any choice of the initial points. A numerical example is performed to illustrate the effectiveness of our proposed iterative procedure.

2. Two parameter Weibull model In (l.l), if we let a = 0 (consider random variable T =X - a in case a # 0. similar results follow), the three parameter Weibull distribution will be reduced to C-l

f(x;

0, b, C) =

G exp( -( 5)‘)

(24

where x 2 0, b > 0 and c > 0. Let (xJ~=~ > 0 be the independent identically distributed samples from W(0, b, c). To avoid triviality, from now on we will assume that {Xi, i = 1,. . . , n} are not all the same, that is, at least two of them are distinct. The log likelihood function of W(0, b, c) can be written as log L(0, b,

c) = n log c - ~tc log b + (c - 1) i

log xi - 2

i=l

(2)’

i=l

where xi > 0, b > 0, c > 0. The maximum likelihood estimates of the scale and shape parameters solving the following system of equations:

a log L(0,

b, c)

ab a

log L(0, b, c)

ac

=-

-?lC

b

+ 2 cx;b-c-l

(2.2) are obtained

by

(2.3)

= 0,

i=l

C

=;-nlogb+k

log; i=l

= 0.

P-4)

H. Qiao, C.P. Tsokos /Mathematics

and Computers in Simulation 37 (1994) 47-55

49

Or equivalently,

(2.5)

(2.6)

and t xi” log xi i=l

1 _-=-

1” y1 .C log

c

kx:

(2.7)

‘i.

l-l

i=l

The above Eq. (2.7) is a nonlinear equation with unknown c. The existence of solution of Eq. (2.7) has been proved by McCool [2] in the following theorem. Theorem 1 (McCool). Likelihood equation (2.7) has a unique positive solution. From McCool’s theorem, we know the positive solution of (2.7) exists and is unique. We will use c* to denote the unique positive solution of Eq. (2.7). To our knowledge, there are no methods to solve Eq. (2.7) analytically. Thus, we have to resort to numerical methods to obtain an approximate estimate to the unique solution. The famous Newton-Raphson method has been quoted extensively in literature for numerically solving Eq. (2.7). Using Newton-Raphson method, we can write Eq. (2.7) as

” 2 log

Xi 2 Xf - c i xf log Xi + 2 Xi”= 0.

n j=i

(2.8)

i=l

i=l

i=l

Assume g(c)

= ;

,$log Xi & 1-l

i=l

-c 2X:_ log Xi + &. i=l

(2.9)

i=l

It follows that g’(c)

= ;

,$log I-1

Xi 2 Xi”+ c t X’ log Xi -c k X’ log2 xj. 1 i i=l i=l i=l

Once an initial estimate for c, say, c0 is given, we can set up the Newton-Raphson [4] as follows. Let ‘kfl

= ‘k

-~ d4 g’@k)

for k = 0, 1, 2,. . . . ’

(2.10) iteration

(2.11)

50

H. Qiao, C.P. Tsokos /Mathematics

and Computers in Simulation 37 (1994) 47-55

To simplify our notations, let q(c) =

k log Xi,

s*(c) = 2 Xf,

i=l

i=l n

n

s3( c) = c xi” log xi

and

sq( c) = c xi” log2 xi.

(2.12)

i=l

i=l

We can rewrite g and g’ as g(c) = ;s,s,

(2.13)

- cs3 + s2

and 1 g’(c) = -sl( s* + CSJ - csq.

(2.14)

n

Then the Newton-Raphson

iteration procedure

can be rewritten as

(2.15)

The Newton-Raphson method is a popular method to use because it converges fast, at a quadratic rate, if we can select a good acceptable starting point c,, ([4]). However, it is also well known that Newton-Raphson method which relies on a starting point does not always converge. That is, if applying Newton-Raphson method one would experience the difficulty in locating an initial point for the Newton-Raphson method to assure convergence. This instability of the Newton-Raphson method is part of our motivation in pursuing the development of a more effective and stable procedure. Thus, we proceed as follows. We can rewrite Eq. (2.7) using our notations of Eq. (2.121, 1

$3 --s2

=-

c

Sl n

Solving for c, we obtain c=

ns2

.

Pm3- SlS2

If we let ns2 4(c)ns3 -

s1s2

,

then the Eq. (2.7) is equivalent to c=

q(c).

(2.16)

H. Qiao, C.P. Tsokos / Mathematics and Computers in Simulation 37 (1994) 47-55

51

Thus, the problem of solving Eq. (2.7) has been reduced to finding the fixed point of q(c). By McCool’s theorem, we know c” = q(c*) and c* is the unique positive fixed point of q(c). We shall now proceed to develop a new iterative procedure, which we shall call it “Simple Iterative Procedure”, or SIP. The SIP begins with ‘k ‘k+l

+ &k)

=

2

(2.17)

’

starting at an initial point c = cO for k = 0, 1, 2,. . . . In developing our iteration procedure, the key questions that we need to address are: (i) Does the proposed SIP always converge? (ii) How fast does the proposed SIP converge? (iii) Does the convergence of the SIP depend on what condition, whatever? Consistent with the aim of the present study we introduce the following theorem answers the above questions.

that

Theorem 2. The Simple Iteration Procedure always converges for any positive starting point. That is, ck generated by SIP (2.17) always converges to the unique faedpoint c” of q(c) as k + ~0, for any starting point c0 > 0. And the convergence is at least at a geometric rate of 4.

To prove this theorem, we need the following lemmas. Lemma

1. sz - s2s4 < 0.

Proof. We start with

s*3 =

( 2 xf

xi)* =( i

log

(x;‘*)(x;‘*

log 2

([$l(xF’2y]‘[ if(‘:/’log xi)2II t

s

li))j

i=l

i=l

= t

i=l

xi” f: Xf log2 xi = s*s4. i=l

where the equality holds only for the case that all xi, i = 1,. . . , n, are the same. In our case, because of the assumption we made in the beginning of this section, the inequality is sharp. [7 Lemma 2. ns3 < s1s2. Proof. We begin by verifying ns3 - s1s2 =

C 14i
(Xf-Xi”)(lOg

Xi-log

Xj)20

H. Qiao, C.P. Tsokos /Mathematics

52

and Computers in Simulation 37 (1994) 47-55

where the equality holds only when all X, i = 1,. . . , n are identical. For the same reason given in Lemma 1, thus the inequality is sharp. q Lemma 3. q(c) > 0 and q’(c) < 0 for any c > 0. That is, q(c) is a decreasing function for-positive C.

Proof. From Eq. (2.131, we have n

n 4(c)

xx; i=l

=

nkxf

=

log xi-

klog

XitX;

i=l

i=l

ns2

ns3 - s1s2

> 0.

i=l

Simplifying the expression and applying Lemma 1 and Lemma 2, we have q,(c) = n2(si (ns3

-s2s4) -

q

co.

v2)2

Now we are ready to prove our main Theorem 2. Proof of Theorem 2. From Lemma 3, we know if c > 0, then q(c) > 0, therefore, (c + q(c)/2

> 0.

Thus, from Eq. (2.171, if the starting point c,, > 0, then all ck > 0. The starting point c,, may have three options: (1) cg = c*, then q(c,) = c* and cr = c2 = * * * = c*, obviously ck converges to c*. (2) If c,, > c”, then from Lemma 3, q(c,) < q(c*) = c*. Namely, c* E [q(c,), cJ. (31 If cg < c*, then from Lemma 3, q(c,) > q(c*) = c*. Namely, c* E [c,, q(c,)]. Let us use Int Cc, dl to denote the closed interval with two end points c and d. Thus, Int Cc, d) implies [c, d] if c < d or [d, c] if c > d. In all three options, we can conclude that c* E Int(c,, q(c,)). Similarly, we can conclude that c* E I&Cc,, q(c,)). for k = 1, 2.. . . Note that cr generated from Eq. (2.17) is in the middle of ca and q(c,). In view of Lemma 3, q(c,) must be in the middle of cg and q(c,), too. It follows clearly that I Cl - 4(c,)

I 5 + I CrJ- 4(c,)

I

and Int(c,,

44)

c Int(c,,

4(c,)).

Similarly, we can show that

Ic/c+1 -4(c,+,)I

I3lc~-q(ck)l

and Int(++r,

&+J>

c In+,,

By induction we conclude that

Ic/y&)l

+c+IJ(co)l.

44).

H. Qiao, C.P. Tsokos /Mathematics

and Computers in Simulation 37 (1994) 47-55

53

Also since we have proved that c* E Int(c,,

q(4),

we have, Ic~--c*l”~lc0-4(c,)l. Thus, it is clear that ck converges to c* for k + 03. C* -

a

Ck+l

Ck

t

ck+2

&k+l)

‘dck)

dck+d

Fig. 1 illustrates the behavior of the proposed Simple Iteration Procedure. Now we proceed to discuss the rate of convergence of the proposed SIP, that is, how fast does the proposed SIP converge? It is clear that we can write

1ck

- CO 15 ;

1CO - q(c,,) I*

which implies that the proposed SIP will converge at least at a geometric rate of 4.

q

Furthermore, as it will be seen from our numerical results, SIP converges much faster than the geometric rate of $. In fact, it is even faster than the famous Newton-Raphson method when, of course, this method converges.

3. Numerical examples We have simulated many sets of data {xi}FZl from the Weibull probability distribution W(0, b, c), for different parameter configurations of b, c and different sample sizes 12. We proceed with the popular Newton-Raphson method and the proposed SIP to estimate the two parameters of the Weibull model by solving Eq. (2.7). For Newton-Raphson Procedure, we use Eq. (2.15); for the method SIP, we use Eq. (2.17). The simulation was conducted for many different starting points co for each procedure and we used the same stopping criterion for both methods, that is, we stop if I+-c&i1


where we chose tolerance numerical results.

to1 = 0.001 for our numerical illustrations.

(34 Table 1 summarizes our

Analysis. (i) In about one third of the cases that we simulated, the Newton-Raphson Procedure failed to converge whereas all the cases using the proposed Simple Iterative Procedure converges. In fact, by our Theorem 2, the Simple Iterative Procedure always converges for any starting point co.

Results of the Numerical

54

H. Qiao, C.P. Tsokos /Mathematics

Table 1 Comparison

of Newton-Raphson

b

procedure

and the proposed

n

C

and Computers

in Simulation

SIP

CO

Convergence Newton’s

3 3 3 3 5 5 5 5 100 100 7 7

2 2 2 2 .5 .5 .5 .5 8 8 12 12

25 25 50 50 25 25 50 50 75 75 100 100

37 (1994) 47-55

10 5 10 5 10 5 10 5 10 5 10 5

Steps

Method

SIP 8 7 8 7 10 8 9 8 5 6 6 7

16 9 14 8 failed 31 failed 29 5 failed 4 failed

(ii) In the cases where the Newton-Raphson Procedure converged, the proposed Simple Iterative Procedure converged faster in all but one case, in which the initial point was quite close to the true value of c. (iii) In all cases the computational time for the Simple Iterative Procedure was much less than the Newton-Raphson Procedure. In fact, by inspecting (2.15) and (2.171, we can clearly see that in each step of the iteration, the Newton-Raphson Procedure needs more computation than the proposed Simple Iterative Procedure. Considering that the Newton-Raphson Procedure took more steps to converge, it took even more time for Newton-Raphson Procedure to obtain approximate estimates of the parameters within the specified tolerance.

4. Conclusion In the present study, we have developed an effective Simple Iterative Procedure (SIP), to obtain numerical estimates of the shape and scale parameters of the Weibull probability distribution. The method has the following features. It always Converges. It does not depend on the initial starting point of the iteration. It converges at least at a geometric rate of +. Thus these findings correct the shortcomings of the famous Newton-Raphson method in obtaining numerical estimates of the parameters of the Weibull probability distribution. In our numerical simulations, one third of the time the Newton-Raphson method fails to converge.

References 111 L.J. Bain and M. Engelhardt, York, 1991 2nd Ed.).

Statistical

Analysis

of Reliability

and Life Testing

Models

(Marcel

Dekker,

New

H. Qiao, C.P. Tsokos /Mathematics [2]

and Computers in Simulation 37 (1994) 47-55

55

J.T. McCool, Inference on Weibull percentiles and shape parameter from maximum likelihood estimates, IEEE Trans. Reliability 19 (1970) 2-9. [3] H. Rockette, C.E. Antle and L.A. Klimko, Maximum likelihood estimation with the Weibull model, J. Amer. Stat. Assoc. 69 (1974) 246-249. [4] K.E. Atkinson, An introduction to Numerical Analysis (Wiley, New York, 1978).