CC BY
24
LIFETIME ANALYSIS OF INCANDESCENT LAMPS: THE MENON-AGRAWAL MODEL REVISITED
Mark Bebbington, Chin-Diew Lai
Institute of Information Sciences and Technology, Massey University, Private Bag 11222, Palmerston North, New Zealand.
e-mail: M.Bebbington@massey.ac.nz, C.Lai@massey.ac.nz
Ricardas Zitikis
Department of Statistical and Actuarial Sciences, University of Western Ontario, London, Ontario, Canada, N6A 5B7.
e-mail: zitikis@stats.uwo.ca
Abstract
The use of the Weibull distribution to model lifetimes of incandescent lamps was originally suggested by Leff (1990). Following this suggestion, Agrawal and Menon have offered and investigated, in a series of papers, an improved model constructed from physical considerations and laws of mathematical statistics. In the present paper we offer supplementary thoughts concerning the Agrawal-Menon model and its several modifications. In addition, we discuss the use of Pinelis's l'Hospital-type calculus rules in the analysis of ageing properties of lifetime distributions.
Keywords: Survival function, hazard rate function, mean residual life function, Weibull distribution, normal distribution, truncated normal distribution, lognormal distribution.
1. Introduction
The laws of physics are commonly taught using incandescent lamps (see, e.g., Evans, 1978; Leff, 1990; Maclsaac et al., 1999; Menon and Agrawal, 2003). Interestingly, statistical analysis of the lifetime of incandescent lamps does not appear to be an old science despite the fact that lamps have been around for more than two centuries: H. Davy created the first incandescent lamp in 1802, and T. Edison created the first practical incandescent lamp in 1879 (see, e.g., Wikipedia, 2007). Recently, Agrawal and Menon (1998), and Menon and Agrawal (2003, 2006, 2007, 2008) have analyzed their reliability characteristics based on theoretical models and experimental data.
Leff (1990) argues that since the hazard rate (HR) function h (t) = -S '(t)/ S (t) of the exponential survival function,
S^ (t) = S^ (t\P) = e-/P,
is constant (i.e., h(t) = 1/P), the lifetimes of incandescent lamps cannot follow the exponential law, unlike radioactive decay. To include the necessary dependence on history and thus improve upon the model's fit to experimental data, Leff (1990) therefore suggested using the Weibull survival function
Sw(t) = Sw(t \ a,P) = e-(/P)a,
where a, J3> 0 are unknown parameters. The Weibull HR function h(t) = (a/P)(t/P)a 1 is increasing for every a > 1. Leff (1990) notes that a = 5 has given a good fit to his data. It is interesting to note that when a = 5 the Weibull survival function is close to the normal survival function (see, e.g., Johnson et al, 1994, p. 632), which hints that the latter may be the basis for an alternative hazard formulation.
Among other things, Leff (1990) also notes that the 'average life' indicated on bulb's package is actually the 'median life'. From the mathematical point of view, the mean and median lifetimes are different, respectively:
f S(t)dt and t^ = F"1(^> :
av j0
where F l(u) is the inverse of the cumulative distribution functionF(t) = 1 - S(t). Leff (1990) observes that
despite being mathematically different, tv and tmed are nearly equal in practice, thus hinting at the symmetric
nature of the lifetime distributions of incandescent lamps. Menon and Agrawal (2003) corroborate this observation.
In a series of papers, Menon and Agrawal (2006, 2007, 2008) suggest and investigate an improved model for the survival function based on laws of physics and the normal approximation to the binomial distribution. Specifically, Menon and Agrawal (2007) argue that on the unit-less scale of an argument t (see below) the survival function is
(1) S (t) = 1 + erf(f-T)) with erf(t) = A.J0V dy.
where y is a parameter associated with variability (see below). We note that the survival function S(t) can be written as O 0 1(y(1 - t)>/2) , and we thus have the equation
(2) S(t) = 1 -O ^(r), where / = 1 and ct2 = 1/(2Y),
where O ^ denotes the normal distribution function. Thus, the unit-less t scale has been chosen in such a way that the mean lifetime / is equal to 1, and thus we have the equation
t
(3) t =
t
av
Hence, on the t-scale we have the following representation for the Agrawal-Menon survival function: (4) S(t) = 1 -O^t), where / = tm and CT = t2av/(2r2).
We shall find it convenient to use the notation
SN (t) = SN (t | UCT) = [-T^je-(^^dy
V2nr
for the normal survival function 1 -O 2 (t).
U,CT v '
Since (4) is the normal survival function, it is strictly smaller than 1 for all t e even though we expect
the survival function to be exactly 1 for all t < 0. When the mean lifetime is notably larger than the variance, the survival function is close to 1 for all t < 0. In practical terms this justifies the use of the normal distribution for modeling the lifetime of lamps under the aforementioned caveat. Nevertheless, from the rigorous point of view we expect lamp lifetimes to follow distributions whose survival functions are exactly 1 at t = 0. Menon and Agrawal (2006) scaled the distribution to have unit mass on the positive half-line, which is the truncated normal survival function
where the normalizing constant is a = 001(j/a). Note also that the constant a is practically equal to 1 when the mean j is larger than, say, 3a (see Table 1 below) and so we have thatSTN(t) « SN(t). The latter observation and
our numerical findings in Table 1 below do indeed justify the use of the normal distribution in the current context, as is done by Menon and Agrawal (2007, 2008).
Given the practical performance of the Menon and Agrawal (2006, 2007) models, we expect that any candidate survival function should be close to the normal survival function. For this reason we suggest considering the lognormal survival function
This is defined for t > 0 and is equal to 1 at t = 0. Limpert et al. (2001) provide a discussion on which distribution - normal or lognormal - should be preferred in various situations, accompanied with numerous illustrative examples. The fundamental difference between the two is that, while both are based on a variety of forces acting independently, in the former the effects are additive, while in the lognormal case they are multiplicative. Lamps can fail from a variety of causes, any one of which is sufficient, even though the other factors may not impede the lamps functionality. Thus lamps can be modeled as a series system, where the reliability function is a product of individual factors. Hence the lognormal with its multiplicative interpretation is an attractive alternative. When the coefficient of variation is small, it is difficult to distinguish the lognormal distribution from the normal distribution (Limpert et al, 2001). The major observable difference between them is that the lognormal is non-symmetric, i.e., the median and mean may differ.
Interestingly, in a survey of published data sets, Limpert et al. (2001) found that the only ones that were not fitted satisfactorily by the lognormal consisted of differences, sums, means or other functions of original measurements. However, for many data sets where a lognormal distribution was acceptable, a normal distribution was statistically rejected. We note also that Xie and Pecht (2003) selected a lognormal distribution to model the reliability of semiconductor light emitting devices.
2. Analysis of the data set of Menon and Agrawal (2008)
Menon and Agrawal (2008) provide the data of the failure times of 50 new Phillips (India) lamps, which we will use to examine the fit of the following four survival functions: SW (t), SN (t), STN (t) and SIN (t). The lamps were
monitored at regular time intervals of twelve hours to count the fused lamps. The instants when at least one fused lamp was found were recorded and there were thirty-two such instances. The minimal recorded value was 840 and the maximal one was 2568. Naturally, several fused lamps were found at some instances.
Hence, we have 'grouped data' with each failure time that has occurred during a twelve-hour period (ti-1, tt ]
recorded as ti . To simplify the estimation procedure, we follow the obvious course, and instead of randomly
'dispersing' the observations throughout the corresponding time periods (ti-1, tt ] we replace them by the mid-values
ti-1 + (t - t-1)/2. Hence, the fifty failure times have been reassigned one of the values 6 +12k hours, for
k = 0,1,2,... Denoting these fifty 'observations' by t*,..., t5*0, we fit the survival functions using the maximum
likelihood method. The numerical results are presented in Table 1 with the corresponding survival functions shown in Figure 1.
Table 1. Fitted distributions for the data set from Menon and Agrawal (2008)
Distributions Parameters LL t av ^ med
Weibull a = = 4.2556 ß = 1541.4 -364.25 1402.1 1414.3
Normal = 1407.8 G = 343.10 -362.84 1407.8 1407.8
Truncated normal = 1407.8 G = 343.10 -362.84 1407.8 1407.8
Lognormal = 7.2198 G= 0.24968 -362.06 1409.5 1366.3
Note that the Weibull distribution has estimated shape parameter a =4.2556, less than the value a = 5 suggested by Leff (1990), but still making the Weibull distribution fairly close to the normal (see, e.g., Johnson et al., 1994, p. 632). The log-likelihood values in Table 1 show that the Weibull distribution has poorest fit of the four, the normal and truncated normal are tied for second place, and the lognormal is slightly superior. The difference between the mean tav and median tmed is largest for the lognormal distribution. Not surprisingly, the mean and median lifetimes of the normal and truncated normal are same.
0.9
\
CD
0.2
0.7
0.6
0.5
0.4
0.1
0
0
500
1000
1500
2500
3000
Figure 1. The empirical (stepwise) and four fitted survival functions: Weibull (dotted), normal (solid), truncated normal (dot-dashed, coinciding with solid), and lognormal (dashed).
Here we see that the reason the lognormal fits better is that it models the skewness, in particular the right-hand tail. May et al. (2000) note that for samples with n > 30, data fit well by the normal has significantly smaller skewness than that not well fit by the normal, and suggest the use of the Shapiro-Wilk (or Ryan-Joiner) test for normality to identify the better of the normal and lognormal distributions. The lamp failure data has skewness of 0.51, compared to -0.21 after taking logs. The Ryan-Joiner test for normality provides P-values of 0.075 and " > 0.1" for the raw and logged data, respectively. Both of these results are evidence in favor of the lognormal over the normal distribution. We note that it is common when data is skewed to reject some observations as outliers, leaving a symmetric distribution. The question, which can only be resolved by the collection of more data, is whether this tail is a real phenomenon.
We will now look at some other characteristics of the four candidate distributions. In the two panels of Figure 2 we show graphs of the estimated hazard rate (HR) functions (top panel) and the mean residual life (MRL) functions
ju(t) = J S(x)dx / S(t)
(bottom panel).
Figure 2. The fitted HR (top panel) and MRL (bottom panel) functions: Weibull (dotted), normal (solid), truncated normal (dot-dashed, coinciding with solid), and lognormal (dashed). The saw-like function in the bottom panel is the empirical MRL function.
Note that ju(0) is the mean tav whose numerical values are recorded in Table 1. The set of Menon and Agrawal
(2008) data shows that the last fused lamp was found at 2568 hours. This explains our choice of the range [0, 3000] in the plots of Figures 1 and 2. A visual assessment of Figure 2 suggests that all the fitted HR functions are increasing and the fitted MRL functions are decreasing. Theoretical results (see Section 3 below) confirm these observations for the Weibull, normal, and truncated normal distributions. In the case of the lognormal distribution, however, the HR function is upside-down bathtub-shaped and the MRL function is bathtub-shaped. (We refer to Section 3 for the definition of UBT and BT shapes.) Thus we have plotted the lognormal HR and MRL functions on the interval [0, 10000] in Figure 3 using the same parameters as those in the bottom line of Table 1. Naturally, as in classical regression analysis, fitting distributions well beyond the data range - the interval [840, 2568] for the Menon and Agrawal (2008) data set - can be, and indeed frequently is, parlous. Hence, although they may be real, rather than artificial, we should not be too disturbed by the decreasing nature of the HR function and the increasing nature of the MRL function outside the interval [0, 3000]. In fact, even the monotonically increasing and decreasing, respectively, HR and MRL functions beyond the interval [0, 3000] may not be natural even for the Weibull, normal, and truncated normal distributions. Indeed, the discussion in Agrawal and Menon (1998) indicates that more realistic monotonicity patterns of the HR and MRL functions should likely be more pronounced on the right-hand tails, due to physical properties such as nearly instantaneous fusing of light bulbs at the end of their lifetimes.
x 10"5
Figure 3. The HR (solid) and MRL (dashed) functions corresponding to the lognormal distribution with the parameters as in Table 1
3. Shapes of HR and MRL functions and Pinelis's calculus rules
We have described above the shapes of the HR and MRL functions of the Weibull, normal, truncated normal, and lognormal distributions. These shapes have been known for a long time (see, e.g., Lai and Xie, 2006, and references therein). The shapes of HR functions are usually identified using the result of Glaser (1980) concerning the similarity of the shapes of the ratios f (t)/ S (t) and f'(t)/ S'(t). In the case of MRL function, the results of Gupta and Akman (1995) have frequently been employed. In this section we will describe how to achieve the same goals using general results, reported by I. Pinelis in a series of papers, concerning the similarity of the shapes of generic ratio-functions u(t)/ v(t) and u' (t)/ v'(t). We believe that this is the first instance of Pinelis's calculus rules being utilized in the context of reliability engineering and, specifically, for analyzing the ageing properties of lifetime distributions.
Determining the shape of the Weibull HR function does not actually require sophisticated tools of analysis, since the function is easy to calculate and possesses the simple form:
hW (t) = — ta-1.
W pa
Hence, we have that
(i) the HR function hW (t) is decreasing when a < 1, constant when a = 1 and increasing when a > 1.
The corresponding MRL function, on the other hand, is difficult to derive explicitly and is therefore challenging to analyze. We shall therefore employ one of Pinelis's results:
Theorem 1. (Pinelis, 2001; see Proposition 1.1). Let u(t) and v(t) be differentiable functions on the interval (a, b) where -ro < a < b <ro. Assume that v(t) and its derivative v'(t) are non-zero on the interval (a, b) and v'(t) does not change its sign on (a, b). Furthermore, assume that u(b-) = 0 = v(b-). The following two statements hold:
(1) If u'(t) / v' (t) is increasing on (a, b), then u(t) / v(t) is also increasing on (a, b).
(2) If u'(t)/ v'(t) is decreasing on (a, b), then u(t)/ v(t) is also decreasing on (a, b).
We can now finish our discussion of the Weibull distribution by determining the shape of its MRL function. Note that the condition u(b-) = 0 = v(b-) is satisfied for the functions £ S(x)dx and S(t) when t t b = ro . Let us write the equation
u (t) pro
Mw (t) = -^r, where uw (t) = I Sw (x)dx and vw (t) = Sw (t).
vw(t) h
The derivative v'W (t) = -fW (t) is always negative on (0, ro) and thus neither vanishes nor changes its sign on(0, ro). Hence, according to Theorem 1, whatever monotonicity we have for the ratio u'W (t)/ v'W (t), the same monotonicity holds for the ratio uW (t)/ vW (t) as well. Writing the former ratio as follows:
(5)
uW (t) = -Sw (t) = _±_ vW (t) SW (t) V (tr
we see that monotonicity of u'W (t)/ v'W (t) is just the opposite to that of the HR function h(t) which has already been determined in (i) above. We thus have that
(ii) the MRL function uW (t) is increasing when a< 1 and decreasing when a > 1. When a = 1, then the function uW (t) is constant; easy to check.
We shall also use statement (1) of Theorem 1 to investigate the HR functions of the normal and truncated normal distributions. Note that the condition u(b-) = 0 = v(b-) is satisfied for the functions f (t)and S(t) when
1T b = ro . Furthermore, we shall use statement (2) of Theorem 1 to investigate the MRL functions of the normal and truncated normal distributions. Note that the condition u(b-) = 0 = v(b-) is satisfied for the functions
| S(x)dx and S(t) when t T b = ro .
Consider now the normal distribution. We start with the equation (6) hN (t) = ^^, where uw (t) = fN (t) and vn (t) = Sn (t).
vN (t)
The derivative v'N(t) = - fN(t) is negative on (-ro, ro) and thus neither vanishes nor changes its sign on (-ro, ro). Hence, according to Theorem 1, whatever monotonicity we have for the ratio u'N (t)/ v'N (t), the same monotonicity holds for the ratio uN (t)/ vN (t) as well. Note that
u'n (t) = t-¡u
vN (t) '
which is an increasing function of t . Hence, we have that (iii) the HR function hN (t) is increasing.
As to the shape of the normal MRL function, we proceed analogously to the Weibull case (see equation (5) in particular). Since we have already established that hN (t) is increasing, we therefore have that
(iv) the MRL function ¡uN (t) is decreasing.
From the mathematical point of view, the only difference between the normal and truncated normal survival functions on [0, ro) is the constant a, which does not influence the shapes of the HR and MRL functions. Hence, from (iii) and (iv), we conclude that
(v) the HR function hTN (t) is increasing, and
(vi) the MRL function /uTN (t) is decreasing.
Finally, we will check the shapes of the lognormal HR and MRL functions. With uLN (t) = fLN (t) and vLN (t) = SLN (t), we have that
(7) uN(t) = logt+&2 - u
) v'N (t) ta2 •
The derivative of ratio (7) is (1 + /u-a2 - log t)/(ta2) and so the ratio must be increasing on (0, c) and decreasing on (c,ro), where c = exp{1 + /- a2}. Such functions are called upside-down bathtub (UBT) shaped. Similarly, if a function is first decreasing and then increasing, it is called bathtub (BT) shaped. (We understand 'increasing' and 'decreasing' in the strict sense throughout the current paper.) The following result of Pinelis tells us what to expect from the ratio u (t)/ v (t) when u'(t)/ v'(t) is either UBT or BT.
Theorem 2. (Pinelis, 2006; third and fourth lines of Table 4.1). Let u(t) and v(t) satisfy same assumptions as in Theorem 1, including u(b-) = 0 = v(b-). Then the following two statements hold:
(1) If u'(t)/ v' (t) is UBT, then u(t)/ v(t) is either decreasing or UBT.
(2) u' (t)/ v'(t) is BT, then u(t)/ v(t) is either increasing or BT.
Since ratio (7) is UBT, part (1) of Theorem 2 implies that the lognormal HR function must be either decreasing or UBT. In turn, the latter statement together with part (2) of Theorem 2 (see also equations (5) above) imply that the lognormal MRL function must be either increasing or BT. To find out which of these alternatives are true, we appeal to Proposition 4.4 in Pinelis (2006), which requires us to calculate the ratio
A = limv2(t)(u/v) '(t) .
'¿0 v | v'(t) |
Theorem 3. (Pinelis, 2006; lower half of Table 4.2). Let u(t) and v(t) satisfy same assumptions as in Theorem 1, including u(b-) = 0 = v(b-). Then the following two statements hold:
(1) If A < 0, then the ratio u(t)/ v(t) in part (1) of Theorem 2 is decreasing, and if A > 0, then the ratio is UBT.
(2) If A > 0, then the ratio u(t) / v(t) in part (2) of Theorem 2 is increasing, and if A < 0, then the ratio is BT. Consider the lognormal HR function hN (t). The corresponding ratio A is:
h'N(t) fLN(t)_ logt + a2-u
2
A = lim v2N (t) = lim^^ = - lim
40 Wn (t)l 40 fLN (t) 40 ta2
Hence, according to part (1) of Theorem 3, we have that
= ro
(vii) the HR function hLN (t) is UBT.
Consider next the lognormal MRL function ¡iLN(t). The corresponding ratio A is:
A = lim vLN (t) ¡ML = lim = - lim -Sln (t) = -ro .
WLN (t)l ^l^ (t)| ^ fLN (t)
Hence, according to part (2) of Theorem 3, we have that
(viii) the MRL function ¡uLN (t) is BT.
4. Concluding remarks
We have shown that the (truncated) normal and lognormal distributions satisfactorily describe the lifetimes of incandescent lamps. The latter is slightly superior statistically and easier to interpret in terms of failure being due to any one of a number of possible causes. The main difference resulting from the use of the lognormal is that the mean and median lifetimes diverge, but the available data set appears to be too small to distinguish this or, equivalently, the tail behavior of the distribution.
The shapes of the HR and MRL functions corresponding to the four distributions employed above to analyze the lifetime of incandescent lamps were analyzed using general calculus results for the ratios of functions and derivatives. Reflecting on these derivations in the context of the existing literature (see, e.g., Lai and Xie, 2006, and references therein), we note that except in some trivial cases such as the Weibull HR function, monotonicity of the HR and MRL functions has been investigated mainly using Glaser's (1980) n-function n(t) = - f '(t)/ f (t). (It is of course the l'Hospital-type ratio f'(t)/ S' (t) corresponding to the HR function f (t)/ S (t).) In this sense, therefore, Glaser (1980) is a precursor to Pinelis's research of the past decade, and so are Gupta and Akman (1995) in the context of determining MRL shapes in terms of the corresponding HR shapes. We again refer to Lai and Xie (2006) for various uses of, and further theoretical developments related to, Glaser's n -function and Gupta and Akman's (1995) results; see also Bebbington et al. (2008) for a recent application of n(t) to the Birnbaum-Saunders distribution. We conclude the discussion with a passage from Pinelis (2004) which refers to the shapes of the generic ratio-functions u(t)/ v(t) and u'(t)/ v'(t):
"...In contrast, the argument just presented is straightforward and rather mechanical. This is exactly the point that we wish to make in this paper. Now a wide class of inequalities become almost trivial in that ad hoc creativity is no longer needed for many such problems. But then is there any excitement left? Yes, what is exciting now is to have such general rules for monotonicity!" (Pinelis, 2004, p. 908)
We concur.
Acknowledgements
We sincerely thank Professor D.C. Agrawal for sharing with us the findings of Menon and Agrawal (2008), and the data set presented therein, prior to its publication. Research of the third author (R.Z.) has been supported by the Natural Sciences and Engineering Research Council (NSERC) of Canada._
References
[1] D.C. Agrawal and V.J. Menon, 1998. Lifetime and temperature of incandescent lamps. Physics Education, 33, 55-58.
[2] M. Bebbington, C.D. Lai and R. Zitikis, 2008. A proof of the shape of the Birnbaum-Saunders hazard rate function. The Mathematical Scientist (to appear).
[3] J. Evans, 1978. Teaching electricity with batteries and bulbs. The Physics Teacher, 16, 15-22.
[4] R.E. Glaser, 1980. Bathtub and related failure rate characterizations. Journal of the American Statistical Association, 75, 667-672.
[5] R.C. Gupta and H.O. Akman, 1995. Mean residual life function for certain types of non-monotonic ageing. Communications in Statistics - Stochastic Models, 11, 219-225.
[6] N.L. Johnson, S. Kotz and N. Balakrishnan, 1994. Continuous Univariate Distributions, Vol 1. (second edition), Wiley, New York.
[7] C.D. Lai and M. Xie, 2006. Stochastic Ageing and Dependence for Reliability. Springer, New York.
[8] H.S. Leff, 1990. Illuminating physics with light bulbs. The Physics Teacher, 28, 30-35.
[9] E. Limpert, W.A. Stahel and M. Abbt, 2001. Log-normal distributions across the sciences: keys and clues. Bioscience, 51, 341-352.
[10] D. MacIsaac, G. Kanner and G. Anderson, 1999. Basic physics of the incandescent lamp (light bulb). The Physics Teacher, 37, 520-525.
[11] J.H. May, D.P Strum and L.G. Vargas, 2000. Fitting the lognormal distribution to surgical procedure times. Decision Sciences, 31, 129-148.
[12] V.J. Menon and D.C. Agrawal, 2003. Lifetimes of incandescent bulbs. The Physics Teacher, 41, 100-101.
[13] V.J. Menon and D.C. Agrawal, 2006. A theory of filament lamp's failure statistics. The European Physical Journal: Applied Physics, 34, 117-121.
[14] V.J. Menon and D.C. Agrawal, 2007. A theory for the mortality curve of filament lamps. Journal of Materials Engineering and Performance, 16, 1-6.
[15] V.J. Menon and D.C. Agrawal, 2008. Renewal rate of filament lamps: theory and experiment. Journal of Failure Analysis and Prevention (to appear).
[16] I. Pinelis, 2001. L'Hospital type rules for oscillation, with applications. Journal of Inequalities in Pure and Applied Mathematics, Volume 2, Issue 3, Article 33.
[20] I. Pinelis, 2004. L'Hospital rules for monotonicity and the Wilker-Anglesio inequality. American Mathematical Monthly, 111, 905-909.
[22] I. Pinelis, 2006. On L'Hospital-type rules for monotonicity. Journal of Inequalities in Pure and Applied Mathematics, Volume 7, Issue 2, Article 40.
[24] Wikipedia, 2007. Incandescent Light Bulb . Accessed on 8 October 2007 at the address http://en.wikipedia.org/wiki/Incandescent light bulb
[25] J. Xie and M. Pecht, 2003. Reliability prediction modeling of semiconductor light emitting device. IEEE Transactions on Device and Materials Reliability, 3, 218-222.
