Generating Renewal Functions of Uniform, Gamma, Normal and Weibull Distributions for Minimal and Non Negligible Repair by Using Convolutions and Approximation Methods by Dilcu Helvaci A dissertation submitted to the Graduate Faculty of Auburn University in partial fulfillment of the requirements for the Degree of Doctor of Philosophy Auburn, Alabama December 14, 2013 Keywords: Reliability, Renewal Function, Renewal Intensity, Convolutions Copyright 2013 by Dilcu Helvaci Approved by Saeed Maghsoodloo, Chair, Professor Emeritus of Industrial and System Engineering Alice Smith, Professor of Industrial and System Engineering Mark Carpenter, Professor of Statistics Fadel M. Megahed, Professor of Industrial and System Engineering ii Abstract This dissertation explores renewal functions for minimal repair and non-negligible repair for the most common reliability underlying distributions Weibull, gamma, normal, lognormal, logistic, loglogistic and the uniform. The normal, gamma and uniform renewal functions and the renewal intensities are obtained by the convolution method. In the uniform distribution case complexity becomes immense as the number of convolutions increases. Therefore, after obtaining twelve convolutions of the uniform distribution, we applied the normal approximation. The exact Weibull convolutions, except in the case of shape parameter , as far as we know are not attainable. Unlike the gamma and the normal underlying failure distributions, the Weibull base-line distribution does not have a closed-form expression for the n-fold convolution. Since the Weibull is the most important and common base-line distribution in reliability analyses and its renewal and intensity functions cannot be obtained analytically, we used the time-discretizing method. Most calculations have been done with the aid of MATLAB Programming Language. When MTTR (Mean Time to Repair) is not negligible and that TTR has a pdf denoted as r(t), the expected number of failures, expected number of cycles and the resulting availability were obtained by taking the Laplace transforms of renewal functions. Finally, the approximation method for obtaining the expected number of cycles, number of failures and availability using raw moments of failure and repair distribution is provided. Keywords: Reliability, Renewal Function, Renewal Intensity, Convolutions iii Acknowledgements The author is grateful to Dr. Saeed Maghsoodloo for his guidance, encouragement, meticulous attention to details, insight and agreement to oversee this work. The author owes this very dissertation effort and Ph.D. completion to Professor Emeritus Saeed Maghsoodloo?s generosity, motivation, patience and wisdom. Dr. Maghsoodloo has always been more than an advisor; he is a dear friend and the author is thankful to him for his efforts. The author must also recognize Dr. Alice Smith for her demanding, but continuous support for this effort. Dr. Alice Smith has contributed immensely in the academic growth of the author through the efforts, assignments, mentoring and opportunities that she provided to the author. The author also wants to thank Dr. Mark Carpenter for his agreement and cooperation on this research. In addition, the author must recognize Dr. Fadel Megahed for his positive support and guidance. This work is dedicated to Filiz Kaya (mother) and Dr. Jeremy L. Barnes (Husband). Without their support, this work would be more difficult. The author shares his deepest gratitude to all of these people. The author wants to send appreciation to the Turk Petrol Foundation and Professor Zikri Altun for supplying the author financial support that was critical for the completion of this degree. The author wants to add a debt of gratitude for all the friends and colleagues that encouraged the continuation of this research, despite all odds and setbacks. iv Table of Contents Acknowledgements ........................................................................................................................ iii List of Tables ............................................................................................................................... viii List of Figures ................................................................................................................................ ix 1 Introduction .................................................................................................................................. 1 1.1 History of Reliability ........................................................................................................ 1 1.2 Research Objectives ......................................................................................................... 2 1.3 Research Methods ............................................................................................................ 3 1.4 Dissertation Layout .......................................................................................................... 4 2 Literature Review......................................................................................................................... 6 2.1 Reliability ......................................................................................................................... 6 2.2 Availability and Maintainability ...................................................................................... 7 2.3 Reliability versus Quality ................................................................................................. 9 2.4 Reliability Applications.................................................................................................... 9 2.5 Renewal Processes ......................................................................................................... 10 2.6 Importance of Renewal Functions.................................................................................. 11 2.7 Previous Work on Renewal Functions ........................................................................... 12 2.8 Contribution of this work to the Literature .................................................................... 15 2.9 Literature Review Summary .......................................................................................... 17 3 Methodology .............................................................................................................................. 18 v 3.1 Reliability Measure ........................................................................................................ 18 3.1.1 The Reliability (or Survivor) Function R(t) ................................................................... 18 3.1.2 The Mean Time to Failure .............................................................................................. 20 3.1.3 The Hazard (or the Failure Rate) Function h(t).............................................................. 21 3.2 Renewal Intensity Function ............................................................................................ 23 3.3 Convolutions .................................................................................................................. 24 3.4 The Weibull Distribution ............................................................................................... 25 3.5 The Normal Distribution ................................................................................................ 29 3.6 The Gamma Distribution ................................................................................................ 31 3.7 The Uniform Distribution............................................................................................... 34 3.8 Laplace Transforms ........................................................................................................ 36 3.9 Chapter Summary ........................................................................................................... 36 4 Renewal Processes with Minimal Repair .................................................................................. 38 4.1 The Renewal Function M(t) during [0, t] with Minimal Repair..................................... 40 4.2 The Renewal Functions for Known Convolutions ......................................................... 42 4.2.1 The Renewal Function for a Normal Distribution ......................................................... 42 4.2.2 The Renewal Function for a Gamma Baseline Distribution .......................................... 44 4.3 The Renewal Function when the Underlying TTF Distribution is Uniform .................. 49 4.4 Approximating the Renewal Function with Unknown Convolutions ............................ 55 4.5 Discretizing Time in Order to Approximate the Renewal Equation .............................. 56 4.6 Chapter Summary ........................................................................................................... 62 5 Expected Number of Renewals for Non-Negligible Repair ...................................................... 63 5.1 Expected Number of Failures and Cycles ...................................................................... 63 vi 5.2 Availability ..................................................................................................................... 66 5.3 Markov Analysis ............................................................................................................ 69 5.4 Renewal and Availability Functions when TTF is Gamma and TTR is Exponential .... 71 6 The Approximate Expected Number of Renewals for Non-Negligible Repair ......................... 75 6.1 Weibull TBF and Uniform TTR .................................................................................... 76 6.2 Uniform TBF and Weibull TTR..................................................................................... 78 6.3 Gamma TBF and Uniform TTR ..................................................................................... 79 6.4 Uniform TTF and Gamma TTR ..................................................................................... 80 6.5 Intractable Convolutions of f(t) with r(t) ....................................................................... 80 6.6 Approximation Results ................................................................................................... 88 7 MATLAB Program .................................................................................................................... 91 7.1 Minimal Repair .............................................................................................................. 91 7.1.1 Gamma Distribution ....................................................................................................... 92 7.1.2 Normal Distribution ....................................................................................................... 93 7.1.3 Weibull Distribution ....................................................................................................... 93 7.1.4 Uniform Distribution ...................................................................................................... 95 7.2 Non-Minimal Repair ...................................................................................................... 95 7.2.1 The Exact Non-Minimal Repair ..................................................................................... 96 7.2.2 Approximate Non-Minimal Repair ................................................................................ 96 8 Conclusions and Proposed Future Research .............................................................................. 99 8.1 Summary and Conclusions ............................................................................................. 99 8.2 Future Work ................................................................................................................. 100 Appendix ..................................................................................................................................... 108 vii Appendix A: Uniform Convolution Equations, where c = b-a > 0....................................... 108 Appendix B: 3rd and 4th Moment of Sums of iid rvs ............................................................. 116 Appendix C: Moments of the Most Common Base-line Distributions in Reliability........... 118 viii List of Tables Table 1: Gaps and Contributions .................................................................................................. 16 Table 2: Reliability Measures of Most Common Base-line Distributions .................................... 37 Table 3: Normal Approximation to the 12-fold convolution of U(0,1) ........................................ 54 Table 4: Time Discretizing Approximation Method Results ........................................................ 60 Table 5: Moment Based Approximation Results for M1(t) ........................................................... 88 Table 6: Moment Based Approximation Results for M2(t) ........................................................... 89 Table 7: Moment Based Approximation Results for A(t) ............................................................. 89 Table 8: Parameters and Density Functions of Most Common Baseline Distributions in Reliability .................................................................................................................................... 118 ix List of Figures Figure 1: Weibull Graph Based on Different Beta Values ........................................................... 27 Figure 2: The Graph of Weibull Reliability Function ................................................................. 29 Figure 3: Normal Distribution Graph When R(t0.05) = 0.95 .......................................................... 30 Figure 4: Uniform Distribution Graph .......................................................................................... 35 Figure 5: The Intervening Times of Two Successive Renewals ................................................... 38 Figure 6: Relative Error versus t? ............................................................................................... 61 Figure 7: On & off system ............................................................................................................ 69 Figure 8: Relative Error versus Time Graph for M1(t), M2(t) and A(t) ......................................... 90 Figure 9: MATLAB Based Minimal Repair GUI ......................................................................... 92 Figure 10: MATLAB Based Non Minimal Repair GUI ............................................................... 98 x List of Acronyms TTF = T Time to Failure MTTF Mean Time to Failure MTTR Mean Time to Repair or to Restore MTBF Mean Time Between Failures CDF Cumulative Distribution Function FR Failure Rate pdf Probability Density Function pmf Probability Mass Function rv Random Variable HZF Hazard Function CHF Cumulative Hazard Function Pr Probability RHS Right Hand Side CPU Central Processing Unit IFR Increasing Failure Rate xi CFR Constant Failure Rate DFR Decreasing Failure Rate MO Modal Point NID Normally Independently Distributed RF Renewal Function RNIF Renewal Intensity Function xii Nomenclature t A specified value in the range space of T ()Ht Cumulative hazard function ()ET Expected Value of TTF ()ft Probability density function (pdf) ()TQt Unreliability function at time t ()Ft Cumulative distribution function ()ht Instantaneous hazard function ()Rt Reliability function ()sNt Number of survivors at time t ()fNt Number of failed items by time t ? Population standard deviation ? Population mean 2? Population variance 0 t?? Minimum or guaranteed life ? The Exponential failure rate ? Weibull shape parameter (or slope) xiii ? Weibull characteristic life ()t? Renewal intensity function at time t ()Mt Renewal function within the interval [0, t] t? Time Increment ()At Availability Function iTTF Time to the i th failure ( , )Uab Uniform distributions over the real interval [a, b] ()n? Gamma function at n > 0 ? The Gamma distribution shape parameter ? The Gamma distribution scale parameter ()()nft n-fold convolution at time t L Laplace transformation symbol 1 CHAPTER 1 1 Introduction This chapter introduces the history of reliability, outlines the objectives that the research intends to achieve, and then gives the research methods adopted. The introduction provides the framework for the research that follows. The chapter concludes by describing the layout of the dissertation. 1.1 History of Reliability ?During the expansion after World War I, the aircraft industry was the first to use reliability concepts. Initially everything was qualitative. As the number of aircraft grew during the 1930?s, reliability was slowly being quantified as function of mean failure rate and average number of failures of an airship or airplane? [1]. Before World War II, reliability studies were mostly intuitive, qualitative and subjective [2]. The early development of mathematical reliability models began in Germany during World War II, where ?a group led by Wernher Von Braun was developing the V-1 missile? [3]. During the 1950?s, nuclear industry started to develop and to use reliability concepts in nuclear power plants and control systems [1]. Further, this last decade witnessed the initial stages of the use of component reliability in terms of ?failure rate, life expectancy, design adequacy and success prediction? [3]. 2 Since the turn of the 20th century, much of the reliability research results started to transfer to industry and academia. In the past eight decades universities have been teaching reliability theory and applications. There has also been steady growth on reliability publications. Furthermore, in order to compete in today?s global economy, manufacturing and other industries should consider reliability as a primary concern [4]. 1.2 Research Objectives The research objectives are a combination of mathematical methods along with the some approximations. This work focuses on the renewal function (RF), whose definition will follow, and has mainly four objectives. The RF is simply the mathematical expectation of number of renewals in a stochastic process. The first objective explores the RF,Mt() , for minimal repair by using convolutions of gamma, normal and uniform distributions. The closed-form expressions for the RFs are provided in the case of normal, gamma (which includes the exponential) and the uniform distributions. Both renewal functions and the renewal intensities )(t=? dM t d t( )/ ( ) are provided for these distributions. Further, the exact uniform convolutions through the eighth have been known since 1983[5]. The second objective is discretizing time in order to approximate the fundamental renewal equation. Unlike the gamma and normal underlying failure distributions, the Weibull base-line distribution (except when the shape parameter 1?? ) does not have a closed form expression for the n-fold convolution nft()() . Therefore, we cannot obtain the Weibull renewal 3 and intensity functions by using the convolution method. Since the Weibull distribution is the most important of all underlying distributions in reliability analyses, this dissertation uses an approximation method by discretizing time in order to estimate the RF, which can be applied to any baseline distribution [see also E. A. Elsayed, (pp. 428-432)]. The third objective is obtaining renewal and availability functions for some time to failure and time to repair distributions by using Laplace transforms when repair time is not negligible. The fourth objective is obtaining an approximation for expected number of failures, number of cycles and availability when repair is not negligible for some common reliability distributions by using the first four raw moments of failure and repair distributions. 1.3 Research Methods The dissertation presented so far generates the renewal functions for some common distributions under the case of minimal repair and non-minimal repair. It develops and enhances an overarching process based on multiple mathematical and statistical theories that will generate renewal and intensity functions, and will also provide approximation methods. For the exponential distribution obtaining a closed-form expression for the RF was documented for well over 100 years ago. Unfortunately, for other distributions such as Weibull and uniform this is very difficult. We used the convolution method for the uniform distribution in order to obtain its renewal and intensity functions. For the uniform distribution, as the number of convolutions increases the problem becomes more complex both geometrically and 4 mathematically. Therefore, after convoluting twelve uniforms we applied the normal approximation. Unlike the gamma and the normal underlying failure distributions, the Weibull base-line distribution (except when the shape parameter 1?? ) does not have a closed-form expression for the n-fold cumulative convolution nFt() () , and hence nn1M t E N t F t ? ??? ? ()( ) [ ( )] ( ) cannot be used to obtain the value of the RFMt() . Since the Weibull distribution is the most important baseline distribution in reliability analyses, we used the time-discretizing approximation method [6] described in Chapter 4. Further, most calculations are done with the aid of MATLAB Programming Language. We also used moment based method and Laplace transforms method for non-negligible repair. 1.4 Dissertation Layout The dissertation is divided into eight chapters including this first chapter entitled ?Introduction?. The layout and organization of remaining chapters are as follows. Chapter 2 presents the literature review on reliability. This includes the concept of reliability in mathematical details. It also explains the terms availability and maintainability and underlines the differences of these two concepts. Moreover, it identifies the differences between quality and reliability. Chapter 2 continues on renewal stochastic processes and it discusses reliability applications. Finally, it explains the importance of renewal functions, contributions to the literature and the general previous works that have been done. 5 In Chapter 3, the reliability methodology is presented. The principles, mathematical development, and structure are identified and introduced. Reliability measures like, MTTF (Mean Time to Failure), hazard function are explained both mathematically and conceptually. It explains the mathematical concepts of convolution method. Chapter 3 continues by some distribution functions such as uniform, gamma, normal and Weibull that was studied in this work. Chapter 3 concludes with a brief summary table. Chapter 4 describes the renewal processes for minimal repair in detail. It gives the renewal and intensity functions for the normal and gamma distributions. It also provides convolutions of the uniform through order n = 12 for the case of minimal repair. By using these convolutions it calculates the renewal and intensity functions for the interval [0, t] when the underlying distribution is uniform. Chapter 4 continues by explaining the time-discretizing method in order to approximate the fundamental renewal equation, which can be applied to any base-line distribution. It concludes with approximating the Weibull distribution?s RF. In Chapters 5 and 6, unlike the previous chapters, we assume that MTTR (Mean Time to Repair) is not negligible and that TTR has a pdf denoted as r(t). Chapter 5 gives the expected number of failures, expected number of cycles and availability by taking the Laplace transforms of renewal functions, and Chapter 6 gives the approximate number of renewals for most common distributions by using raw moments of failure and repair distributions. Chapter 7 introduces the Matlab program, describes the inputs and outputs of the program for minimal and non-minimal repair cases. Finally, Chapter 8 gives the conclusion and possible future work. 6 CHAPTER 2 2 Literature Review This chapter introduces the literature review on reliability. It defines reliability, availability and maintainability. It also discusses differences between quality and reliability. Furthermore, it briefly discusses the applications of reliability which abound. Then, it continues with importance of renewal functions and previous work on renewal functions. Finally, it summarizes the contribution of this study. 2.1 Reliability The concept of reliability is not new. Both manufacturers and customers have long been concerned with the reliability of products they produce and use [7]. Basically the general perception about reliability is functioning without any problem. Stating something is reliable implies it can be depended on to work satisfactorily. However, the real definition of reliability involves quantifying measures. ?The reliability of an item is the probability that it will adequately perform its specific purpose for a given period of time under specified stressed conditions? [8]. Another definition is as follows: ?Reliability is defined to be the probability that a component or system will perform a required function for a given period of time when used under stated operating conditions? [4]. As it is seen from the two above definitions, reliability is a probability. Therefore, this implies reliability can never be negative or greater than one. Since reliability is a probability, probability 7 axiom results are used in reliability [8]. Reliability studies deal with different complexity levels of units or systems [9]. Unreliability induces heavy losses to organizations. It also causes them to lose reputation. Failure of products before reaching their warranty period is costly. Therefore, today we have reached a point where reliability is considered as a major performance measure [10]. In summary, reliability is the survival probability of an item or system beyond a certain point in time. Since the cost of unreliability is high, reliability is considered as a major performance parameter. Unlike most classical statistical parameters, reliability is always a function of time (t). Although in cases of very short mission time, reliability can be considered as static merely as an approximation. 2.2 Availability and Maintainability Availability includes both failure and repair rates of a system, and therefore, it is considered to be one of the most important reliability performance measures [6]. ?Availability is the probability that a system or component is performing its required function at a given point in time or over a stated period of time when operated and maintained in a prescribed manner? [4]. There are generally four types of availability measures [8]. These are: 1) Point (or instantaneous) availability ()At : Instantaneous availability at time t is the probability that a repairable unit is functioning reliably at time t. Therefore, if there is no repair, the availability, ()At , is equal to the reliability function Rt() . 8 2) Limiting availability: As time increases, instantaneous availability will approach a constant value, once an item has stationary times to failure Ti, stationary times to repair TTRi, and this value is known as the limiting availability [8]. 3) Average availability on (0, t]: This can be explained as the ?expected fraction of time? that a system is operational during (0, t], [8]. 4) Limiting average availability on (0, t]: This is also known as ?steady-state availability?. It is the availability of a component or system when the time interval is very large [6]. Maintainability is the probability that a failed item can be repaired or restored to become operationally effective within a specified period of time when repair action is performed in accordance with prescribed procedures [11]. Maintainability is important for eliminating defects, correcting defects and their causes, meet new requirements and adapt to a changing world. A design process should start with defining system maintainability objectives. Both maintainability and availability are very important measures in reliability theory and have a wide range of applicability. However, in the context of this work the brief amount of information given above on maintainability and availability should perhaps be satisfactory. Maintainability directly affects availability because time for repairs or preventive maintenance can change a system from available to unavailable state. Therefore, there is a close relationship between reliability and maintainability. Reliability affects maintainability, maintainability affects reliability and they both impact availability [12]. Further, like reliability, availability and maintainability are also probabilities. Therefore, probability theory rules can be applied to both availability and maintainability. 9 2.3 Reliability versus Quality Quality may be defined in many different ways. It can be conceptual, perceptual or conditional and it may be interpreted differently. Basically quality means ?fitness to use? and it has several dimensions such as ?performance?, ?reliability?, ?durability?, ?serviceability?, ?conformance to standards?, etc. [13], [14]. Therefore, reliability which is closely associated with product quality is one of several quality dimensions. Quality can be defined qualitatively and can be achieved through a satisfactory quality assurance program [4]. ?Quality assurance is the set of activities that ensures quality levels of products and services properly maintained and that supplier and customer quality issues are properly resolved? [14] . On the other hand, reliability is largely concerned with how long a product can continue to operate under specified conditions once it becomes functional. ?Reliability may be viewed as the quality of a product?s operational performance over time, and as such it extends quality into the time domain? [4]. Improving reliability is an important part of improving product quality [15]. 2.4 Reliability Applications Examples of high-reliability systems abound worldwide, such as aircraft systems([16], [17] etc.), electric power generating stations ([18], [19] etc.), chemical plants ([20] etc.), power systems to telephone and communication systems ([21], [22], [23], [24], etc.), computer systems and networks ([25] etc.) [1]. Reliability studies are conducted at either the component or system level. Generally, reliability calculations are easier at the component than system level. If it is the system level, 10 then there are two different classifications. A system can be static or dynamic. Further, in both cases it can be serial, parallel or mixed. Explanations of each follow. In static reliability component or subsystem reliabilities are considered to be approximately constant for a specified duration of time [26]. In this case, the mission time is sufficiently short so that the assumption of constant reliability is almost tenable. Whereas, in dynamic models there will be continual reliability degradation of subcomponents with respect to time. In series system-models all subsystems must operate reliably in order for the system to function properly. As soon as one subsystem fails, the system fails [27]. There is also another model type called ?chain model? or ?weakest link model? that is described in the literature under the title of series systems [26]. Based on this model, a system will fail as soon as the weakest component (or link) fails. Therefore, system reliability is equal to reliability of the weakest among all components or subsystems [26]. However, a pure parallel system is composed of subsystems or components that the success of any one of which results in system success [28]. Therefore, such a system is reliable if at least one component is reliable. There are many studies in the literature that deal with system reliabilities such as [29], [30], [31], [32], etc. 2.5 Renewal Processes Renewal processes are stochastic events such that their nth stage value is the sum of n independent random variables of common distributions with nonnegative ranges [33]. As an example, consider a machine component that is replaced as soon as it fails with a new one. Let 11 Nt() be the number of replacements during the interval [0, t] of length t. Then,Nt() is called a renewal counting process. The study of renewal processes focus on the following topics: (1) The pmf (Pr mass function) of Nt() , (2) The expected number of renewals during [0, t] or [t0, t0+t], ENt[ ()] , denoted by ()M t E N t ? [ ( )], M for mean, is called the RF. Henceforth, the symbol E will represent the Expected-Value operator. (Note that this case also includes the negligible repair-time.), (3) The occurrence Pr mass or density function of a renewal at specific epochs of time, and (4) The time needed for the occurrence of n events (such as failures that are followed by a replacements) to occur [34]. [For more details see U. N. Bhat (1984), Elements of Applied Stochastic Processes, 2nd Ed., Chapter 8.]. 2.6 Importance of Renewal Functions Renewal functions, gives the expected number of failures of a system or a component during a time interval and this is used to determine the optimal preventive maintenance schedule of a system [35]. Renewal functions are quantities that have particular importance in analysis of warranty ([36], [37], [38], [39] etc.) [40]. Expected cost of warranty estimation is closely related to the RF estimation [36]. For example, consider the case that cell phone manufacturer has a two year free replacement warranty which means that if a cell phone fails manufacturer agrees to replace it with a brand new one without any charge. Then , suppose M(2) is the expected number of failures(replacements) during two years warranty period and C(2) expected warranty cost is C(2)= c*M(2), where c is assumed as the fixed cost per replacement [36]. It is 12 obvious that the cost of warranty is greatly affected by the number of replacements [35]. In today`s competitive environment product`s warranty policy is important to attract customer. Offering longer warranty terms usually attract more customers but it means more cost [41]. So, warranty has two important roles protection and promotion [42]. Therefore, it is very important to determine an optimal warranty time and this means obtaining M(t) with greater accuracy is very essential especially if manufacturer produces large number of units or very expensive items. Renewal functions play an important role and have wide variety of applications in decision making such as inventory theory ([33]), supply chain planning ([43], [44]), continuous sampling plans ([45],[46]), insurance application and sequential analysis ([47][48]) [36],[43]. 2.7 Previous Work on Renewal Functions As we have seen in the previous section renewal functions play an important role in many applications. Therefore it is important to obtain renewal functions analytically. Based on analytic method, M(t) is the inverse Laplace transform of ??Mswhere ? ? f(s)M s = s[1 f(s)]? ([35]), where Laplace transforms will be defined later. ?The advantage of analytical method is one can carry out parametric studies of the RF, i.e., the behavior of M(t) as a function of the parameters of the distribution? [40]. However, for most distribution functions obtaining the RF analytically is complicated and even impossible [43]. Therefore development of computational techniques and approximations for renewal functions has attracted researchers [49]. 13 One of the well-known approximations is 21 2 1( ) / / 2 1ttM ? ? ?? ? ??? which is generally known as asymptotic approximation and was generated by Tacklind,S (1945)[50] and also cited in numerous papers such as [51]. The asymptotic expression has a closed-form expression thus it is easy to apply optimization problems that involve renewal process [43]. However since asymptotic expansion is not accurate for small values of t, Parsa&Jin (2013) [43] propose better approximation by keeping the positive features of asymptotic approximations such as simplicity, closed-form expression, and independence from the distribution. Jiang (2010) [52], proposes an approximation for the RF with an IFR which is also useful in areas such as optimization where renewal function needs to be evaluated. There are series methods available in the literature to approximate renewal functions such as, Smith and Leadbetter (1963) [53] who developed a method to compute the RF for Weibull by using power series expansion of t? where ? is the shape parameter of the Weibull. On the other hand instead of using power series expansion, Lomnicki (1966) [54] proposes another method by using the infinite series of appropriate Poissonian functions of t?. There are also many other approximations methods available such as Xie (1989) [55], Smeiticnk & Dekker (1990) [56], Baxter et al (1982) [57], Gang&Kalagnaman (1998) [58], From ( 2001) [59]etc. For example Xie (1989) [55] proposed RS-method for solving renewal-type integral equations based on direct numerical Riemann-Stieltjes integration. There are usually three criteria: model simplicity, applicability and approximation accuracy to evaluate the value of the analytical RF approximation [52]. Increasing the complexity may lead more accurate approximation but may make the process complicated and difficult to implement in practice [60]. 14 Studies on renewal functions for some particular underlying renewal functions such as Weibull have been done. Jiang (2009) [61], proposes an approximation for the RF of Weibull distribution with an IFR which is accurate for time t up to a certain value of larger than the characteristic life. On the other hand, Jin & Gonigunta (2010) [62], proposes an approximation method for Weibull RF with DFR. Sinha (1985) [63] obtains Bayes estimation of the survivor function of the s-normal distribution. Papdopoulos &Tsokos (1975) [64] obtain confidence bounds for the Weibull failure model. Many others are also available like [65], [66], [67], etc. Furthermore, in the literature bounds on renewal functions have been discussed. Since they provide upper and lower bounds on warranty costs bounds on M(t) are very useful for many warranty models [40]. Ross (1996) [68], shows that if a distribution has DFR then the RF is bounded as ? ? 2 211 1 12?ttMt?? ??? ? ? ??? ? where 1? and 2?? are the first and second raw moments. Marshall(1973) [69] provides lower and upper linear bounds on the RF of an ordinary renewal process. Ayhan et.al. (1999) [70] provide tight lower and upper bounds for the RF which are based on Riemann-Stieljes integration. There are also many other studies available about bounds on RF such as [71], [72], [73], [74], [75], [76] and etc. Finally, simulation can be considered as an alternate approach to estimate the value of renewal function. Brown et al (1981) [77] use the Monte Carlo simulation to estimate the RF for a renewal process with known interarrival time distribution. Papadopoulos & Tsokos (1975) [64], perform Monte Carlo simulation to obtain 90% and 95% Bayes confidence bounds for the random scale parameter and reliability function to illustrate their results. ?The simulation 15 approach offers an alternate method for obtaining the solution to the problem without the need of solving complicated mathematical formulations? [40]. 2.8 Contribution of this work to the Literature The renewal and intensity functions with minimal repair for the most common lifetime underlying distributions normal, gamma, uniform and Weibull are explored. The exact normal, gamma, and uniform renewal and intensity functions are derived by the convolution method. Unlike these last three failure distributions, the Weibull distribution, except at shape ? = 1, does not have a closed-form function for the n-fold convolution. Since the Weibull is the most important failure distribution in reliability analyses, its approximate renewal and intensity functions were obtained by the time-discretizing method. And also for non-minimal repair moment based approximation method was generated. Table 1 summarizes the contributions of this dissertation. 16 Table 1: Gaps and Contributions Gaps in the Literature Contributions Obtaining the renewal functions for the most baseline distributions is not possible. Obtaining renewal functions and renewal intensities for normal and gamma distribution by using convolution method. Approximating Weibull renewal function and renewal intensity by using time discretizing method. ?The renewal function of uniform distribution using the convolution method is not available. We used geometrical mathematical statistics method to obtain uniform convolutions from n= 2 through n = 12 and then applied the normal approximation for convolutions beyond 12 to obtain renewal function of uniform distribution. For the case of non-negligible repair only the closed-form renewal functions exist in the case of exponential TTF and TTR. We obtained the closed form renewal function for gamma TTF and exponential TTR when ? (shape parameter of gamma) is an exact positive integer from 2 to 7. We obtained approximations for expected number of cycles ? When TTF is Weibull and TTR is uniform ? When TTR is Weibull and TTF is Weibull ? When TTF is gamma and TTR is uniform ? When TTR is gamma and TTF is uniform by first obtaining the convolution density functions and then using the time discretizing method. We obtained approximations for two parameter exponential, three parameter Weibull, gamma, normal, lognormal, logistic and loglogistic distributions by using the first four moments of failure and repair distributions. 17 2.9 Literature Review Summary The literature describes a multitude of research in the area of reliability. Reliability has a wide variety of applications from aircraft to power systems. The literature supports increasing applications of reliability in today?s competitive global economy. In this chapter, we explained some important concepts such as availability and maintainability and described the differences and association between quality and reliability. This chapter also explains the importance of renewal functions and the general previous work that has been done about them. 18 CHAPTER 3 3 Methodology This chapter describes the mathematical concepts and relationships that are needed in reliability and renewal theory presented throughout this dissertation. 3.1 Reliability Measure There are three measures of reliability: (1) The reliability functionRt() , (2) The mean time to failure (MTTF), and (3) The hazard (rate) function ht(). If either function (1) or (3) is known, then the other 2 measures can be uniquely determined, but the knowledge of (2), i.e., ? = MTTF, is not sufficient to obtain unique functions forRt() andht() . In fact, for the same MTTF of an underlying distribution, there are uncountably infinite other base-line distributions that have identically the same MTTF. 3.1.1 The Reliability (or Survivor) Function R(t) The reliability of a component is the probability (Pr) that the device will perform without failure during the mission time t, under specified stress conditions. For example, R(of a new passenger tire for t = 500 interstate miles) ?100% = 1. However, the reliability of the same passenger tire under racing conditions at Indianapolis 500 would be almost zero. Note that the terminology survivor function, ()St is also used for non- 19 repairable items, such as light bulbs, transistors, or rocket-motor of an unmanned spacecraft. Let T = the random variable lifetime, or time to failure (TTF), with Pr density function . Then, the reliability function at time t, or the survival Pr for a mission of length t, is given by (the Pr that a component lifetime exceeds time t) ( ) ( ) d 1 1 1 TT t R t T t f x x P r T t F tt QP r ( ) ( ) ( ) ( ) ? ? ? ? ? ? ? ? ? ? ??, (3.1) where ()Ft = the cdf of T at t, and ( ) ( )?TTQ t F t represents the unreliability function at time t, or the cumulative failure Pr by time t. The pdf (Pr density function), ft() , is also referred to as the failure (or mortality) density function. Some authors use the notation ()tS for the reliability or survivor function at time t to imply survival probability beyond t; however, the notation Rt() is a bit more prevalent in engineering applications. We now obtain the relationship between Rt() and by differentiating equation (3.1) with respect to t, recalling that ??TTf t d F t d t d Q t d t( ) ( ) / ( ) / because ? ?() ? ???t F t f T d TT . Due to the fact that time cannot be negative, and thus the lower limit in this last integral must be zero instead of ? ?, i.e. 0 ( ) ( )? ?t F t f T dTT = Failure Probability by time t. The developments below show the relationship between Rt() and ()f t 20 1??? ? ? ? ? ? ? ? ???d R t d d F tF t f t d R t d F td t d t d t( ) ( )( ) ( ) ( ) ( ) and ? ? ? ? ?f t d R t d t f t d t d R t( ) ( ) / ( ) ( ) . 3.1.2 The Mean Time to Failure We are now in a position to obtain the relationship betweenRt() and the mean time to failure (MTTF=?), which is defined as the mathematical expectation of T = TTF. Henceforth, the symbol V will represent the variance operator, and the reader must be cognizant of the fact that anytime the operator E or V is applied to any rv (random variable), the end-result will always be a population parameter. The Mean Time to Failure (MTTF) is given by; 00 M T T F ?? ? ? ? ???E T tf t d t t d R( ) ( ) ( ) M TT F 0( ) ( ) 0 ( ) ( )0 0 0 ? ? ? ?? ? ? ? ? ?? ? ?tR t R t d t R t d t R t d t] . (3.2) The above equation clearly shows that the unconditional mean-life, ( ) MTTF?E T , starting at age zero, of any device or system is given by the integral of its reliability function evaluated always from zero to infinity so that the MTTF is the total area under Rt() and the abscissa t (i.e., from 0 to ? even if the minimum life is greater than zero). 21 Note that if a component or system is repairable (or renewable, i.e., failed units are almost immediately replaced), then ET()is called the mean time between failures (MTBF). The MTBF of any system also can be found by integrating its reliability function from zero to infinity [78] (a simple technique of obtaining MTBF for complex systems). Again, the lower limit of the integral must always be zero even if the minimum life t0 = ? > 0. 3.1.3 The Hazard (or the Failure Rate) Function h(t) By definition the failure rate (FR) of any device is defined as the failure rate during ?t t ?t(, ) given that the device age is t, i.e., ? ? ? ? ? ?PrP t T t ? t T t P r t T t ?t R t R t ?tF R t ? t T t ? t R t ? t? ? ? ? ? ? ? ??? ? ?? ( ) ( )( ) ( ) (3.3) This implies that the failure rate of a component at time t is the probability that it will fail in the interval ?t t ?t(, ) given that its life has exceeded time t (i.e., given that the age of the component is t). Put differently, the failure rate of a population of identical items at time t is the proportion of the units failing per unit of time in the interval ?t t ?t(, ) amongst all survivors at time t. The hazard function (HZF), (),ht is simply the instantaneous FR, i.e., ? ? 0 1 ][ ? ??????? ? ? ? ? ???t Rt ?t R t d R t d t f tht Rt ?t R t R t () ( ) / ( )( ) l im ( ) ( ) ( ) (3.4) In other words, if we have a system with N0 identical items on test at time 0, ??sNt survivors at time t and ??fNt failed items by time t, then by the above definition ht() is the rate 22 of failure, fdN t dt( )/ (this derivative is not quite appropriate because ??fNt is discrete), only amongst the survivors ??sNtbeyond time t, i.e., 0 0 0 () () () () ()( ) ( ) ( ) s f s s ss s s Nt ddN d N N t Nd N t d R t d t f td t d t d t d tht NtN t N t N t R t R t N ?? ????? ??? ? ? ??? ?? ? ? ? ? ?( ) / ( ) ( ) ( ) as before. The quantity ? ? ? ?( t ) ( t ) ( t ) ( t ) |h d t d R R d F R P t T t d t R P t T t d t T t/ / / ,? ? ? ? ? ? ? ? ? ? ? ? gives the proportion of items that will fail within (t, t dt)? amongst those that are still functioning at time t. From the cumulative of hazard function,CHF H t? () , we can obtain the relationship amongst the reliability measures ()ft , Rt() , and ht()as shown below: t 000( ) ( ) ( ) ?? ? ? ?? ? ? ? ? ? ? ? ?? ? ? ??? ()/ l n ( ) l n ( )tt HtH t h x d x d R R R x R t R t e t 0 x dx ( ) e ?? ? h R t () ; ?? ? ? ? Htfth t f t h t R t h t e Rt ()()( ) ( ) ( ) ( ) ( ) () (3.5) Properties of the Hazard function (a) () 0?ht for all t. 23 (b) 00 0 0 0fh f hR? ? ?()( ) ( ) ( )() must be finite at t = 0 unless f(0) = ? (c) () ?? ??t htlim iff h(t) is an IFR, simply implying that any man-made system must have a finite life (or a finite TTF), i.e., no man-made system can last forever! (d) () 00 ( ) ( ) 1 ?? ? ???? Htf t d t h t e d t. Note that the above relationships imply that the assumption of almost constant reliability during a short interval (0, t) leads to an almost zero HZF value. 3.2 Renewal Intensity Function The renewal intensity, ()t? , gives the instantaneous renewal rate at time t, i.e., ? ? 0 () ?t Mt ? t M t d M tt ? t d t? ()( ) lim , ? ???? (3.6) so that ()?t ?t? gives the unconditional probability element of a renewal during the interval t t ?t?(, ) , and in the case of negligible repair time, ()t? also represents the instantaneous failure intensity function; hence, 0 ( ) d()??t xM t x? . Note that nearly authors in Stochastic Processes and some in reliability literature refer to ()t? as the renewal density because ()?t ?t? gives the probability element of a renewal during the interval t t ?t?(, ) ; however ()t? is never a pdf. Hence, we have chosen to refer to ()t? as the renewal intensity in lieu of renewal density. The renewal intensity, ()t? , should not be confused with the hazard rate function h(t) [79], because ()?ht?t is the conditional probability of a failure during time interval (t, t+?t ) given that the unit?s age is t, whereas ()t ?t? ? is the unconditional 24 probability of a failure during ?t . The hazard rate function is a relative rate pertaining only to the first failure, whereas the intensity function is an absolute rate of failure for also repairable systems [4], including minimal repair. Only in the case of exponential base-line distribution (CFR) both ht() and ()t? are identically equal to the constant failure rate ?. 3.3 Convolutions The n-fold convolution of a statistical distribution arises in a wide variety of applications of probabilistic models such as reliability theory, renewal theory, inventory theory, queuing theory, continuous sampling plans, insurance risk analysis, and sequential analyses [80]. Mathematically, a convolution of two density functions ??1ft and ??2ft denoted 12 f*f , gives the density of sum of two variates 12TT? . It can be proven that the convolution of with is given by [81], 1 2 2 1 00 12 tt f t * f t f t u f u d u = f t u f u d u? ? ???( ) ( ) ( ) ( ) ( ) ( ) (3.7) Note that in general the lower limit would be ?? but in reliability theory, the lower limit is always zero. In the literature review there are some studies about uniform convolutions [82], [83], [81], but we used the geometrical approach to obtain the precise uniform convolutions through order 12 that are presented in Chapter 4. Maghsoodloo and Hool [5] obtained the uniform convolutions for orders 2, 3, 4, 5, 6 and 8. 25 3.4 The Weibull Distribution It is well known that the underlying distribution of almost any manufactured dimension by man can be approximately modeled by a normal (or Laplace-Gaussian) pdf. The Weibull pdf plays the exact same important role for the underlying distribution of TTF (or lifetime) of most mechanical and electrical components or systems. The key events in the derivation of what is now known as the Weibull distribution took place between the years 1922 and 1943. There were three groups of scientists working independently for different aims. Waloddi Weibull (1887- 1979) was one of the three working on this distribution. The reason that the distribution bears his name is the fact that he propagated it internationally and interdisciplinary [84]. To arrive at a Weibull pdf, consider an exponential pdf at the constant failure rate ? = 1 failure per unit of time. Note that the symbol ? is used throughout this dissertation if and only if ht() is a constant failure rate (CFR). Clearly, 0 1xe dx . ? ? ?? (3.8) Then, for convenience letting t0 = ?, we make the transformation 0?t ?x ? ????????????, . As a result, 1 1?d x t ?? dt ? ? ? ? ?? ? ? ??? ? ? ? ???? ? ? ?, and substitution into (3.8) yields 11 1 ??t ? t ??? ? ? ? ? ?? t ? d t ? t ?e ? e d t ? ? ? ? ? ? ? ? ? ? ? ???????? ? ? ? ? ? ? ?? ? ? ?????? ? ? ?? ? ? ? ? ? ? ??? . (3.9) 26 Since the value of the integral in Eq. (3.9) is equal to 1 (or 100%), the integrand 1 ?t ?? ??? t ? e? ? ? ? ???? ??? ???? ?? ??? ?? ??? ? ?? ? must be a probability density function (pdf) over the range [? ? . The pdf, ? ? 1 0 ?t ?? ??? t ?f t e t ?t? ? ? ? ???? ? ?? ???? ? ? ??? ? ?? ? ? ??? ? ? ? ? , (3.10) and () 0ft? for 00 t ?t? ? ? , is called the Weibull model, denoted ? ? ? , with minimum (or guaranteed) life ? (the location parameter), the characteristic life ?, and slope (or the shape parameter) ?; ? ? ? is called the scaling parameter. Different authors tend to use different symbols for the three parameters of a Weibull pdf, but ? is the most common symbol for the slope. Figure 1 shows the Weibull density based on different shape (? ) values. The reliability function for Weibull distribution is, ? ? 1 () ? ? x ? ? ?? u t t ? ?? ? x ?R t T t e d x e d u . ? ? ? ?Pr ??? ???? ?? ? ??? ????? ??? ???? ? ? ? ???? ???? ? ? 1 , 0 , t t Rt et ?? ?? ? ? ???? ?? ??? ???? ? ?? ? ? ?? (3.11) 27 Figure 1: Weibull Graph Based on Different Beta Values The MTTF of Weibull distribution, 00 ( ) ( ) ????? ? ?? ???? ? ?? ? ? ?t ? ? ?? ? E T R t d t d t e d t; letting ?t ?x ????????? ?? in the second integral results in, ? ? ? ?1 111 00 () ??? ????? ? ? ? ???? ?x ?x? ? ? ?E T ? e x d x ? x e d x?? //. Since by definition, 1() 0 ? ??? ? nx ? n x e dx, and ( ) ( 1)?n? n ? n+ we obtain ? ? ? ?1M TTF ( ) 1 1??? ? ? ? ? ? ? ???? ?ET ? ? ? ? ? ? ? ? ?/( ) ( ) / (3.12) The above MTTF ?? iff 1?? (i.e., only for the CFR = constant failure rate, and IFR = 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0.0 5 0.3 0.5 5 0.8 1.0 5 1.3 1.5 5 1.8 2.0 5 2.3 2.5 5 2.8 3.0 5 3.3 3.5 5 3.8 4.0 5 4.3 4.5 5 4.8 Beta=1 Beta=2 Beta=3 28 increasing failure rate cases). The modal point of Weibull is given by 11 ?MO ? ? ? ? ???? ? ? ??? /( ) ( ) / (3.13) The variance of Weibull is given by ? ? 222 2111??? ? ? ?? ? ? ? ? ???? ? ? ? ? ? ? ???T? V T ? ? ????() (3.14) Special cases of the Weibull model (i) When ? , the Weibull reduces to the exponential pdf with minimum life t0 = ? and mean-life (or characteristic life tc) equal to ?. The Weibull pdf with slope ? can be used to model the TTF of a component during its useful life (constant failure rate = CFR). (ii) When ? and ? , the Weibull becomes the Rayleigh density function ? ? 2 2?22 ?? ???( ) ( ) ( ) t ? tft ? t ? e te/ ///, where 2 ???c?t / . (iii) When 0 < ? < 1, it can be shown that the hazard function of the Weibull is a decreasing function of time (DFR = decreasing failure rate) so that the Weibull may be used to model the TTF during the burn-in (or debugging, or infant-mortality) period of a component (see Ebeling pp. 31-32). (iv) When ? > 1, the hazard function is increasing (i.e., increasing failure rate = IFR) and the Weibull density can be used to model the TTF during the wear-out period of a component. 29 (v) When 1 < ? < 2, h(t) is an IFR concave function, and for ? > 2, h(t) is an IFR convex function because 2 2d ()htdt > 0 when ? > 2. Figure 2: The Graph of Weibull Reliability Function 3.5 The Normal Distribution The normal distribution which is also called Laplace-Gaussian is the most commonly used distribution in the field of statistics [85]. The reasons are related to its mathematical properties, central limit theorem, and various experimental responses often have distributions that are approximately normal [86]. The normal density was first discovered by Abraham de Moivre (a French mathematician) in 1738 as the limiting distribution of the binomial pmf and T h e W e i b u l l R E f u n c t i o n f o r d e l t a = 3 0 0 0 0 , t h e t a = 6 0 0 0 0 , a n d s l o p e = 4 0 0 . 2 0 . 4 0 . 6 0 . 8 1 1 . 2 0 20000 40000 60000 80000 100000 120000 Ti m e i n K m R ( t ) 30 was named the ?normal? by Karl Pearson in 1894 [87]. It can be proven that the points of inflection occur at ? ? ?. Figure 3: Normal Distribution Graph When R(t0.05) = 0.95 A continues random variable T is said to have a normal distribution with parameters ? and ? (or ? and ?2), where M T T F??? ? ? ? ? and , iff the pdf of T is [88] ? ? ? ? ? ?2 221 2 t? ?f t ? ? e ????? t ?? ??? ?? ? ? ?/;, (3.15) The reliability function is derived from ? ? 2 2 11 1 ? 22t x? t?R t e x p d x ??? ? ? ??? ????? ? ? ? ? ???? ?? ?() , (3.16) 31 where ? represents the cdf of the N(0,1). There is no closed-form solution to the above integral, and it must be evaluated numerically. The hazard function cannot be written in closed form either and is provided below. ? ? ? ?2 221 2 1 ? t? ?e ft ??ht Rt t? ? ?? ?? ?? ?? ???? / ()() () (3.17) It must be noted that the normal failure law is applicable only if its MTTF = ? is more than 10 standard deviations to the right of the origin (or zero) because ?(?10) < 7.62?10?24. 3.6 The Gamma Distribution By definition 1x 0 ( ) x xn? n e d? ??? ? . Dividing both sides of this last definition we obtain: 1x 0 11 x x()? ??? ? n ed?n , and per force the integrand ? ? 11 ??? ntf t n t e?n(); must be a density function called the standard gamma pdf. When n is not an integer, the common notation for the shape parameter n is ?, ?, or ?. A single integration by parts of ()?n will show that ( ) ( 1 ) ( 1 )? ? ?? n n ? n. After 1?n() integration by parts, we obtain ? ?? ? ? ?( ) 1 2 1? ? ?? n n n ...?. Inserting 1n= into the definition of ()?n yields 11 00 ( 1 ) 1? ? ?? ? ???tt? ? t e d t e d t ??; therefore, 32 ? ? ? ? ? ? ? ?( ) 1 2 1 1 1 0 1? ? ? ? ? ? ? ?? n n n n ?!!.... Further, this last result also implies that ?? n+ 1 n ? n( ) ( ) ( ); it can be proven that (1/ 2)? ?? and hence (3 / 2 ) (1 / 2 ) (1 / 2 )??? , 5 3 3 3 12 2 2 2??? ? ? ? ? ?? ? ?? ? ? ? ? ?? ? ? ? ? ?? , etc. For example, ? ? 9 36288010 ??? ! . To obtain the gamma pdf, we make the transformation x t?? in the definition of gamma function: 11 ? 1 ? 0 0 0 1() ? ? 1 ? ? () n x n t n t? n x e d x ? t e d t ? t e d t ?n ? ? ? ? ? ?? ? ? ?? ? ?? ? ?( ) ( ), or 1 ? 0 1 ? ? 1() ntt e d t?n ?? ??? () . (3.18) The above equation clearly shows that the integrand must be a pdf over the range [0, ?) because its integral over [0, ?) yields 100%. The function under the integral is called the gamma pdf, as shown below, in statistical literature with rate ? (or ? = 1/? the scale parameter) and the shape parameter n. 1 () ntf ?nt t e? ?????( ) ( ) . (3.19) When n is a positive integer the above density is called Erlang and has extensive applications in Queuing Processes. When n is not a positive integer, the most common notion for shape is ?; thus, a better representation for the general gamma pdf is 1 () tf t?te? ????? ??( ) ( ) . 33 The meaning of the parameters n and ? will be made clear in the application example provided below, where 1/? is also called the scale parameter. The major application of the gamma pdf occurs from the fact that the sum of n independent and identical exponential rvs, TSystem = T1 + T2 +...+ Tn, has a gamma time to failure density function with parameters n and ?, where each Ti is distributed according to e??? t . Therefore, an n-unit standby system with quiescent failure rates of almost zero for standby units and perfect switching, has a lifetime that has the gamma density with parameters n and ?. Further, the gamma density has applications in maintenance scheduling where the amount of deterioration during an interval [t1, t2] has a gamma pdf with scale parameter 1/? and shape 21n ? t -t? (), where ?>0 is a constant of proportionality. The first four moments of the gamma pdf are given by ? ? ? ?12 ?? ? ? ? ? ?S y s t e m nE T E T T T n / ? ? ? ? 212 ?? ? ? ? ? ?S y s t e m nV T V T T T n / the standardized third moment 3 2?n? / , and 4 36?n??( / ); hence the kurtosis is equal to 4 6 ? ? n? / , and it can be verified that 1 ?MO n??( )/ . Note that the values of 3 2?n? / and 4 36?n??( / ), clearly show that the limiting distribution (i.e., as n ? ?) of the gamma density is the Gaussian N(n/?, n/?2). Unfortunately, n must exceed 200 (because the exponential is highly skewed) before the normal approximation to gamma becomes fairly adequate. 34 To obtain the reliability function for the gamma pdf, we make use of its relationship with the Poisson pmf as described below, where X(t) describes the number of failures during a mission time of length t and Tn represents time to the nth failure. 1 0 ( ) = 1) ? ?? ? ??? ? ? ? ?? ?? ? kn t k R t T t P r X t n fa il u r ts ekenP r ( ) ( ) () !( . So far, we know the expressions for the gamma pdf, the reliability function, and the MTTF ?? /n . The hazard function can be obtained from f t R t( )/ ( ) . 3.7 The Uniform Distribution If a < b, the random variable T is said to have a continuous uniform probability density on the interval [a, b] if and only if the density function of T is given by [89] ? ? 1 / ( ) , 0, b t be ls e wa h e ra eft ? ? ??? ?? (3.20) where a is the minimum-life and b is the maximum-life. The graph of the above pdf is shown atop the next page. Note that all values of t from a to b are equally likely in the sense that the probability that t lies in an interval of width entirely contained in the interval from a to b is equal to a ? regardless of the exact location of the interval [90]. Its reliability measures are listed below. The MTTF is: ? ?12? b a??. (3.21) Median of the uniform distribution is: ? ? 05 12t b a??. . (3.22) 35 Figure 4: Uniform Distribution Graph Variance of the distribution is: ? ?22 112? b a?? (3.23) Because 11? ??? ? ? ??? b tt btd x d xb a b a b a, the reliability function is given by: ???Rt 1 , 0 , 0 a bt a b t tb , b < t a ??? ?? ??? ?? ? ?? (3.24) The hazard function is: ? ?? ? ? ?1() b b t b aftht Rt a?? ? ? ??/ /() () 0 , 01 , a a t b bt t ???? ? ?? ??? , which is an IFR. Hence the CHF is given by H(t) = 0 n 0 l , , a ba a b tb b t t ,t ??? ????? ??? ?? ??? ? ???? (3.25) ? b 1/(b- ) 36 3.8 Laplace Transforms If ()ft is a known function of for values of , its Laplace transform ? is defined by the equation, ? ? -s t 0 ( ) (t) t ? ? ?f t e f d ,L (3.26) and abbreviated as, ? ? ? ?f s ( )ft?L (3.27) Thus, the above operation transforms the function ()ft of the real variable t into a new function ? of the subsidiary variable s, which may be real or complex [91]. Often, it is easier to work in the s-space than the t-space to obtain solutions to mathematical problems. 3.9 Chapter Summary This chapter introduced the mathematical concept of the dissertation. The table below shows the summary of reliability measures of most important base-line distributions. Only the Logistic and Loglogistic density reliability measures, both of which have also reliability applications, are not yet provided. 37 Table 2: Reliability Measures of Most Common Base-line Distributions Lifetime Distribution Failure Density f(t) Survival Function R(t) Hazard Function MTTF Exponential te??? ??te ? 1/? Weibull ()1() tt e ??? ????? ? ? ? ??? ???? ()te ?????? ? 1()t ???? ? ? ? ???? ()[(1/ ]1) x ? ? ?????? Gamma 1()() ntten ?? ?? ?? 1 0 ()!kn t k t ek ??? ? ?? , when n is a pos. integer ? ? ? ?f t tR/ n/? Lognormal 2ln()1 2 te t ?? ?? ?? ln( )[]?? t? ? ? ? ? ?f t tR/ 2/2e??? Beta, t= ? + (U ? ?) x , 0 ? x ?1 11() (1 )( ) ( ) abab xx??? ??? ? Standard Beta pdf 1 , ,tb e ta d is t a bU ?????? ????? ? ? ? ?f t tR/ ? ?()Uaab?? ??? Normal ? ? ? ?2 2/21 2 x?e? ? ??? 1 ????? ????t?? ? ? ? ?f t tR/ ? Uniform 1 / ( ), 0, b t belsewheraae? ? ???? 1 , 0 , 0 a bt a b t tb , b < t a ???? ? ??? ? ?? ?? 0 , 0 1 , a a t b bt t ???? ? ?? ??? ? ?12b?? 38 CHAPTER 4 4 Renewal Processes with Minimal Repair Suppose that failures occur at times Tn (n = 1, 2, 3, 4, ?) measured from zero and assuming that replacement (or restoration time) is negligible relative to operational time, then Tn represents the operating time (measured from zero) until the nth failure, where T0=0 . Because the pdf of T1 may be different from the intervening times X2 = (T2 ? T1), X3 = (T3 ? T2), X4 = (T4 ? T3), ..., we consider only the simpler case of probability density function (pdf ) of time to first failure f1(t) being identical to those of intervening times X2, X3, X4, ? as depicted below. Figure 5: The Intervening Times of Two Successive Renewals Note that X1, X2, X3 ? represent intervening times between failures, while Ti represents time to the ith renewal measured from zero. Further, all Xi?s are assumed iid (independently and identically distributed). The above figure clearly shows that Tn (Time to the nth renewal) = ii1X?? n = sum of the times to the 1st failure from zero plus the intervening times of 2nd failure until the nth failure. If n > 60, then the Central Limit Theorem (CLT) states that the distribution of Tn 0 X 1 T 1 T 2 X 2 T 2 T 3 X 3 T ime 39 approaches normality with mean n?, where ? = E(Xi) = the mean time between successive renewals, i = 1, 2, 3, 4, ? and with variance 2n? , where ?2 = V(Xi). However, if the pdf of Xi, f(x), X being the parent variable, is highly skewed and/or n is not sufficiently large, then the exact pdf of Tn = ii1X?? n is given by the n-fold convolution of f(x) with itself denoted as ? ? ? ? ? ? ? ?1 1*nT nf t f t f t?? , where ? ???1nft? is the pdf of the sum X2 +X3 +?+X n, or the (n ? 1) convolution of f(t) with itself. Note that Figure 5 is also approximately valid for the case of Minimal-Repair, i.e., the case when MTTR ? 0, or MTTR = Mean Time to Repair is negligible relative to MTTF. Therefore, in this section by renewal we mean either the replacement of a failed component with a brand-new one, or the case when the failed component can almost immediately be repaired and consequently be put back on-line. The simplest and most common renewal process is the homogeneous Poisson process (HPP), where the intervening times are exponentially distributed at the constant inter-renewal (or failure) rate ?. Because ? is a constant and intervening times are iid, a Poisson process is also referred to as a homogeneous renewal process. Throughout this dissertation, we will establish that only in the case of exponential failure Pr (Probability) law with CFR (Constant Failure- Rate), the RNIF is identical to the constant instantaneous hazard function h(t) = ?. Further, it is also well-known that for a HPP the ? ?V tNt ???? ?? , and hence the coefficient variation of ()Nt is given by () 1/Nt tCV ??. 40 4.1 The Renewal Function M(t) during [0, t] with Minimal Repair Because the two events ? ? ( ) N t n? and ? ? nTt? are equivalent, it follows that ? ? ? ? ? ? ( )n nP N t n P T t F t? ? ? ????? , where ? ? ? ? ? ? ? ? ? ?11( ) *nnF t F t F t? is the n-fold convolution representing the cdf of ii1X n nT ??? . Thus, ? ? ? ? ? ? ? ? ? ? ? ? ? ?1 ( ) 1 nnP N t n P N t n P N t n F t F t?? ? ? ? ? ? ? ?? ? ? ?? ? ? ? It has been proven by many authors both in Stochastic Processes and Reliability Engineering that the RF for the duration [0, t] is given by n1 ()( ) [ ( ) ] ( ) ? ??? ? nM t E N t F t , (4.1a) and ? ? ? ? ()1 2( 2 1 ) ( )nn nV FN t M tt ? ?? ? ?? ? ? ?? ? ? ?? (4.1b) where the random variable N(t) represents the number of renewals that occur during the time interval [0, t], and ()()nFt is the cdf (cumulative distribution function) of the n-fold convolution of f(t) with itself, i.e., ( ) ( )0( ) ( ) t nnF t f d? ? ?? , where ()()nft is the n-fold convolution density of f(t). It is also widely known that the RNI (Renewal Intensity) of Ft() by definition is given by )( ( ) /t=? dM t dt. Authors in Stochastic Processes refer to )(t? as the renewal density, while some authors in Reliability Engineering refer to )(t? as the RNIF (Renewal Intensity Function); 41 because it is also well known that )(t? is never a pdf, throughout this dissertation we refer to it as the RNIF. The reader should distinguish between the RNIF )(t? and the generally constant and unit-less traffic intensity parameter in Queuing Theory defined as ? = (average arrival rate of customers)/ (average service rate) = expected fraction of time a single server is busy. Further, an approximate expression for the third raw moment of N(t) is given by Kambo et al (2012) [92]. Their expression for E[N(t)3] can be used to approximate the skewness of N(t). It is also well-known that for a homogeneous renewal process the pdf of interarrivals Xi?s is given by () ???? tf t e , and that of the time to nth failure (or arrival, or renewal) is given by ? ? ? ? ? ? ? ? 1 ?() ?n n tTnf t f t t e?n ? ???? (the gamma density with shape n and scale ? = 1/ ?). As a result, the use of Eq. (4.1a) for the interval [0, t] leads to the RF ? ?()M t E N t?????? ()=1 () ?? nn Ft = 1 1 0 ()() t nx n x x e d xn ? ? ?? ? ? ? ??? ? = 1x n1x0 ( x )e d x( 1) !???? ?? ?? ???t nn= x 0x0 ( x )e d x!??? ?? ?? ??t n n n = x0 dx ? ??t =?t , a fact that has been known for more than 100 years. Further, the RNIF for a HPP is a constant and is given by )(t=? ( ) / /??d M t d t d ( t) d t= ?. It has also been proven in the theory of stochastic processes [see D. R. Cox and H. D. Miller (1968), pp. 340-347][93] that the ? ?lim = / ??t M tt ? , where ? = MTTF. In the case of CFR (Constant Failure-Rate), because the mean of the exponential base-line distribution is ? = 1/?, then it follows that this last limiting result is an exact identity only for the exponential density ??? te , 0 ? t < ?. Further, when MTTR (Mean 42 Time-to-Restore) is almost zero, then the renewal intensity is practically the same as the failure intensity. In the case of exponential base-line distribution, with minimal repair, the point availability function is simply equal to its reliability function ? ? ??? tR t e . 4.2 The Renewal Functions for Known Convolutions Unfortunately, obtaining a closed-form expression for the RF for all distributions is not as simple as the case of exponential interarrival times. This is due to the fact that the n-fold convolution of most baseline distributions used in reliability analyses is not either known or attainable. 4.2.1 The Renewal Function for a Normal Distribution For the sake of illustration, suppose that the time between failures, Xi , i = 1,2, 3, 4, ?, are NID(? = MTBF, ?2), i.e., normally & independently distributed with MTBF (Mean Time Between Failures), and process variance ?2. Then, Statistical Theory dictates that time to the nth failure (measured from zero) is the n-fold convolution of N(?, ?2) with itself, i.e., Time-to-the- nth-Failure Tn = TTFn ~ N(n?, n?2). Hence, in the case of minimal-repair, from Eq. (4.1a) the RF is given by ()11( ) ( ) () ?? ?? ??? ???n nn tnt F t nM ?? (4.2) where ? universally stands for the cdf (cumulative distribution function) of the standardized normal deviate N(0, 1), and ()()nFt gives the Pr of at least n renewals by time t. Xie et al (2003) 43 [94] gives the same exact expression for the normal RF as in Eq. (4.2), which they used as an approximation for the Weibull renewal with shape ? ? 3. It should be born in mind that the normal failure law is approximately applicable in reliability analyses only if the coefficient variation of T, denoted CVT, ? 0.15625 = 15.625% because the support for the normal density is (? ?, ?), while TTF can never be negative (this assures that the size of left-tail below zero is less than 1?10?10). If the CV is not sufficiently small, then the truncated-normal can qualify as a failure distribution; From (2001) [59] discusses the RF for the truncated-normal. As an example, suppose a cutting tool?s TTF has the lifetime distribution N(? = MTBF = 15 operating hours, ?2 = 2.25) with minimal-repair (or replacement-time), where CVT = 0.10. Then, Eq. (4.2) shows that the expected number of renewals (or replacements) during 42 hours of use is given by 1 4 2 1 542 15 2 1 2 4 1 0 7( ) ). .(? ? ????? nM nn while M(62 hours) = 3.747561 expected renewals. Using the limiting result ()Mt ? /t ? , we obtain M(42 hours) ? 42/15 = 2.80 (% relative-error = 31.82, while M(62) ? 62/15 = 4.133333 with % relative-error = 10.29). There is a more accurate approximation for ()Mt , through the second moment, given by 2 2 2( ) / + 2?? ? ? ?M tt ( ) / ( )= 2/ 0 .5 0(C V 1)Tt ? ?? [51]. For the above normal nonhomogeneous Process, M(62 hours) = 4.1333333 + (2.25?225)/450 = 3.638333 (a much closer approximation). We attempted to obtain a more accurate approximation through the third moment for RF M(t), but the corresponding approximate Laplace transform had complex roots, 44 consistent with other findings. Further, for a Laplace-Gaussian process, the renewal intensity by direct differentiation of (4.2) is given by 2222 11 1 2 ?? ?? ?? ? ? ??? ( ) / /[]()( ) ( , ) tn ? n ? nn d M t? ? n ? n ? e dt ?t n ? (4.3) where the symbol ? stands for the standard normal density. The value of renewal intensity at 42 hours for the N(15, 2.25) baseline distribution, from Eq. (4.3), is ?(at 42 hours) = 0.07883674 failures/hour. Note that the value of the hazard-rate at 42 hours is given by h(42) = f(42)/R(42) = 10.358138 failures/hour, where R(t) is the reliability function at time t. Because the normal failure Pr law always has an IFR (Increasing Failure-Rate) h(t), then h(t) > ?(t). (Sheldon M. Ross, 1996, pp.426-427) [68] proves that / ( )??t M t? / ?t ? 20.50(CV 1)T ? if h(t) is a DFR (Decreasing Failure-Rate), and he further proves when h(t) is a DFR, then h(t) ? ?(t) for all t ? 0, and as a result R(t) ? ()eMt? , and we add that equalities can occur only at t = 0. 4.2.2 The Renewal Function for a Gamma Baseline Distribution Suppose that the TTF of a hot-water heater has a gamma failure density with shape parameter ? = 1.5 and scale ? = 1/? = 3.5 years [this is quite similar to the Example 9.6 on p. 226 of Ebeling (2010) [4]]. When the heater fails, it is replaced with a new one with the same identical shape and scale (i.e., minimal replacement-time). Our objective is to obtain the RF M(t) for t years of operation. In order to obtain the RF and RNIF of a gamma NHPP, we first resort to Laplace-transforms (LTs). 45 It has been proven in theory of stochastic processes that the LTs of ()t? and M(t) are, respectively, given by [for a proof see (Bhat,1984) [34], pp. 277-280)] L{?(t)} = s 0 f ( s )() 1 f ( s )? (s) te t d t? ??? ?? ? , s > 0, (4.4a) and M()s = L{M(t)} = s 0 f ( s ) se ( ) d s [1 f ( s ) ] s? ( )t M t t? ? ???? (4.4b) Further, it is also widely known that the LT of the gamma density ? ? 1( ) e () tt tf ???? ?? ??? is given by f(s) = ? (?+s) ? ? , ?(s) =L{?(t)} = ? (?+s) ?? ? ?? , and the gamma M(s) = L{M(t)}= ?s[(? +s) ? ]???? . We used Matlab?s ilaplace at ? = 1/3 and ? =1.5, but Matlab(R2012a,64bit) could not invert ?(s) to the t-space at ? = 1.5. It seems that when the shape parameter ? is not a positive integer, there exists no closed-form inverse-Laplace transform for the gamma density; however, when the underlying failure distribution is Erlang (i.e., gamma with positive integer shape), then there exists a closed-form inverse Laplace-transform for ? = 2, 3, and 4; for positive integers beyond 4 there does exist complicated closed-form expressions. In fact, for the specific Erlang density with shape ? ? 2 and scale ? = 1/?, it is well known that M(s) = st 0 e (t)dtM? ?? = 46 2 22? 1 14 s 4 ( s 2 )s ( s 2 ? ) 2 s??? ? ? ??? = 2114s 4 (s 2 )2s???? ?? . Upon inversion to the t-space, we obtain the well-known M(t) = E[N(t)] =L? 1{M(s) } = L? 1 2114 s 4 (s 2 )2s???????????? = 21 ?1e4 2 4 tt+ ???? . Hence, at ? ? 2, the RNIF is given by 2?22 t?? ? ? ?() ?() d M t?edtt , which is quite different from the corresponding gamma (at ? ? 2) IFR HZF ? ? 10 ()( ) e ! k tkt ttekht ?? ? ??? ?? ? ? ??= ( ) / (1 )tt? ? ?? > ()?t for t > 0. We used Matlab?s ilaplace function as an aid in order to obtain the RF and RNIF for the Erlang at shapes ? ? 3 and 4 which are, respectively, given below. ?()Mt L?1 333?s[(? + s) ? ]??????? ?? = ( 3 / 2 )e c o s ( 3 / 4 ) s i n ( 3 / 4 ) / 3 / 3 1 / 3[]tt t t? ?? ? ? ? ? ? ? ?, ??t() L?1 333?(? + ) ?s????????? = ( 3 / 2 )1 e [ c o s ( 3 / 4 ) 3 s i n ( 3 / 4 ]3 {}? ?? ? ? ? ?t tt, ?()Mt 2 /8e / 8 e c o s ( ) s i n( ) / 4 3[]ttt t t????? ? ? ? ? ? ? ? ?, ??t() L?1 444?(?+s) ???????? ?? = 21 e 2 e s in ( )4 []????? ? ? ?tt t The gamma HZF at ? ? 3 is given by h(t) = 2 2 2( ) / ( 3 ) ] / (1 / 2 )t t t??? ? ? ? ? ?, which is not the same function as the first ()?t above. At ? = 5, 6, 7 ? Matlab(R2012a) provides an expression 47 for ()?t only in terms of roots of a polynomial of at least order 5. The user has to find the roots in order to obtain ()?t . Referring back to our example where ? = 1/3 and ? =1.5, we use the cdf F(n)(t) = Pr(Tn ? t) in order to obtain the RF directly from Eq. (4.1a): 1 11 0 ?() ?( ) ( ) ( ?) ( ?? ?? ?? ?? ??? ? n ?? t x nnnF t x eM nt dx? 1 110 ( ? , )( ? ?????? ?? ?????? ?t n u nn u e du ? t n ?n , (4.5) where ( , )?tn?? = Matlab?s gammainc 1 0 1( , ) () nutt n u e d un ?? ?? ????? ? represents the incomplete-gamma function at point ?t and shape n?. In fact, ( , )?tn?? gives the cdf of the standard gamma density at point ?t and shape n?. Thus, for the Water-heater example ? ?12M t years? = 1 (12 / 3.5, 1.5 )n n? ? ?? = 2.11934672 expected failures. Using the 2nd-order approximation we obtain ? ? 2 2 2 2/ ( / ) ( / / ) / [ 2 ( / ) ]? ? ? ? ? ?? ? ? ? ?tM t = / 0 .5 0 (1 ) /t ??? ? ? ? = 2.1190476 (which yields a ?0.0141% relative error). We next directly differentiate Eq. (4.5) in order to obtain the gamma RNIF. That is, ??()() dM t? dtt 11 0 ()( ? ? ? ? ? ?? ???? ?t nx n x e d xtn ??? = 1 1 0 ()( ? ? ? ? ? ?? ???? ?t nx n x e d xtn ??? = 1 1 ()( ? ? ?? ? ? ??? nt n ten ??? (4.6) 48 However, the function under the summation in Eq. (4.6) is simply the gamma density with shape n? and scale ? = 1/?. Then, we used our Matlab program, to obtain the value of Eq. (4.6) at t = 12 years, at shape ? = 1.5 and scale ? = 3.5 years, which yielded ?(t =12) = 0.190348 renewals/ year. The value of the hazard-function at 12 years is h(12) = f(12)/R(12) = 0.0193613/0.0765932 = 0.252781 failures/year. Because the gamma density is an IFR model iff (if and only if) the shape ? > 1, then ?(t) < h (t) for t > 0. Only at ? = 1, the gamma baseline failure distribution reduces to the exponential with CFR, the only case for which ?(t) ? h (t) = ?. In order to check the validity of ?(t = 12) = 0.190348, we resort to the limiting form ? ? 1 1 M T B Ftlim / /t ???? ??. Because the expected TTF of the gamma density is ? =???, then for the Water-heater example ? = MTBF = 1.5?3.5 = 5.25, which yields ?(12 years) ? 1/5.25 = 0.190476/year. Note that since the renewal-type equation for the RNIF ?(t) is given by ? ? ? ? 0 ( x ) ( x ) d x ,t tft f t ?? ?? ?? this last equation clearly shows that ?(0) = f(0); further h(t) = f(t)/R(t) for certain yields h(0) = f(0), and hence ?(0) = f(0) = h(0) for all baseline failure distributions. Moreover, if the minimum-life ? > 0, then ?(?) = f(?) = h(?). Some authors in Reliability Engineering, such as Ebeling (2010) [4], use the expression 0 ()(0, )ee t x dx Mt ??? ?? to represent the reliability function R(t) for a NHPP. Clearly, for the case of gamma baseline distribution, which represents a NHPP, the above Ebeling?s expression is an approximation because the exact unconditional reliability is always given by 49 0 ( x ) d x ()( ) e e t h HtRt ? ????, H(t) being the CHF (Cumulative Hazard Function). Leemis (2009) [8] defines the RF as the cumulative intensity function using his notation ? 0 ( ) ( )dtt? ? ? ? ?? ?, which is identical to the RF M(t), where ?(t) is his notation for the RNIF. It seems that he is also using, atop p. 146, the notation ?(t) as the hazard function for the Weibull. In section 4.3 it will also be established that the HZF h(t) and the RNIF ?(t) of the Weibull are not the same, except at t = 0. It is well known from statistical theory that the skewness of gamma density is given by ?3 = 2/ ? and its kurtosis is ?4 = 6/?, both of these clearly showing that their limiting values, in terms of shape ?, is zero, which are those of Laplace-Gaussian N(?/?, ?/?2). We compared our gamma program at ? = 70, ? = 15, and t = 5000 which yielded M(5000) ? 4.186155, while the corresponding normal program yielded M(5000) ? 4.185793 expected renewals. 4.3 The Renewal Function when the Underlying TTF Distribution is Uniform ?Uniform distribution is used to model the time of occurrence of events that are equally likely to occur at any time during an interval? [95]. Kececioglu (2002)[95] states that ?the most frequently used distributions in Reliability Engineering are exponential, Weibull, normal, lognormal, extreme value, Rayleigh (the Rayleigh being a special case of the Weibull with minimum-life ? = 0 and shape ? = 2), and uniform?. Electrical bulbs, stress of mechanical component which has lower and upper psi limits, network systems are some examples in which uniform distribution is used in reliability models. Zhao and Duan [96] propose a reliability 50 estimation model of IC`s interconnect based on uniform distribution of defects on a chip. And also ?the uniform distribution is used in Bayesian estimation as a prior reliability distribution? ([97], [95]) as an example please see [64]. Accordingly, suppose the TTF of a component or system (such as a network) is uniformly distributed over the real interval [a, b weeks]; then ? ? 1/f t c? , a ? 0, b > a, c = b?a > 0, and the cdf is F(t) = (t?a )/c, a ? t < b weeks. Further, succeeding failures have identical failure distributions as U(a, b). From a practical standpoint, the common value of minimum-life a = 0. Then, the fundamental renewal equation is given by ? ? ? ? 1 0 ( ) ( ) d t MtM t F t f??? ? ? ? ? [[34], pp.277-280]. Since we are considering the simpler case of time to first failure distribution being identical to those of succeeding times to failure, then ? ? ? ? 0 ( x ) (x )d x ,? ?? ?tM t F t M t f (4.7a) whereas before F(t) represents the cdf of T = TTF. However, Hildebrand (1962) [98] proves that f(s)M(s)? = L 0 ( x ) (x ) xt M t f d?????????? = L 0 ( x) (x)t M t dF?????????? , the integral inside the first brackets representing the convolution of M(t) with f(t). Conversely we can conclude that f(s)?M(s) = 51 L 0 ( x) (x)t F t dM?????????? . Upon inversion of this last LT we obtain L ?1{ f(s) M(s)? } = 0 ( x) (x)t F t dM?? = 0 ( x) (x ) (x )t F t d??? . Hence, Eq. (4.7a) can also be represented as ? ? ? ? ? ? 00 ( x ) ( x ) ( x ) ( x ) xttF t d M F tt t F t dMF ?? ? ????? ? (4.7b) The renewal-equation of the type (4.7b) has been given by many authors such as [55], [99], [8], and other notables. In order to obtain the RF for the uniform density, we substitute into Eq. (4.7a), for the specific uniform U(0, b) baseline distribution, for which a = 0, in order to obtain ? ? ? ? 0 t ( x ) d x/???????? ?tM t E N t M t bb; letting t? x= ? in this last equation yields M(t) = 01 ( )( d ) t tb Mb??? ??= 0 1 ( )d .ttb Mb? ? ?? (4.8) The above Eq. (4.8) shows that the RNIF is given by ( ) 1 ( )() ? ? ?d M t M tt d t b b? ? d ( ) ( )dtM t M tb? = 1b ? ttd ( ) ( )eedt bbM t M tb//??? = t1e b/? ? /d [ ( ) ]dt tbM t e? = t1e ba /? ? /t() bMt e? = t/1 e d t Cbb ? ?? = t/eCb???, where C is the constant of integration. Applying the boundary condition M(t = 0) = 0, we obtain ? ? ? ? = e 1bM t E N t ????? ? t/, where time must start at zero, i.e., this last expression is valid only for 0 ? t ? b, b > 0. Note that Ross (1996) [68] gives, 52 without proof, the same identical ??Mt only for the standard U(0, 1) underlying failure density. For example, if the time to down-state of a network is U(0, 8 weeks), then the expected number of down-states during the interval (0, 8 weeks) is given by ()M 8weeks = 8/8e ? 1 = 1.718282. In order to calculate the RF for the same uniform distribution during the interval (2, 4) we may use the above result: 4 2 44 / 8 / 8 22 1( ) ] 0 . 3 6 4 6 9 68() ttt d t e d t eM 2 , 4 w e e k s = ? ? ? ??? Bartholomew (1963) [100] describes ()t? ??t as the (unconditional) Pr element of a renewal during the interval (t, t+?t), and in the case of negligible repair-time,?t() also represents the instantaneous failure intensity function. However, as described by nearly every author in Reliability Engineering, the HZF h(t) gives the instantaneous conditional hazard-rate at time t only amongst survivors of age t, i.e., ()ht ??t = Pr(t ? T ? t +?t)/R(t). The hazard function for the U(0, b) baseline distribution is given by h(t) = 1bt? , 0 ? t ? b, b > 0, which is infinite at the end of life-interval b, as expected. Because the uniform HZF is an IFR, then for the uniform density it can be proven, using the infinite series for 1/1/btb? and the Maclaurin series for ebt/ , that h(t) >?t()for all 0 < t ? b. Next in order to obtain the RF for the U(a, b), we transform the origin from zero to minimum-life = a > 0 by letting ? = a + (b?a)t/b = a + ct/b in the RF ? ? ? ? =M t E N t????? 53 e1b?t/ . This yields, ()M? = ? ?e 1ac?? ? , and hence ? ?( ) = [ ( ) ] e 1t a cM t E N t ???, 0 ? a ? t < b, and c = b ? a > 0. The corresponding RNIF is given by ? ?( ) e t a ctc? ???? ??, 0 ? a ? t < b. Because the uniform renewal function is valid only for the interval [a, b], we will obtain the n-fold convolution of the U(a, b)-distribution which in turn will enable us to obtain M(t) for t > b by making use of Eq. (4.1a) that uses the infinite-sum of convolution cdfs, F(n)(t). As stated by Olds (1952) [101], the convolutions of uniform density of equal bases, c, have been known since Laplace. The specific convolutions of the uniform density with itself over the interval [?1/2, 1/2] were obtained in [5] only for n = 2-fold, 3, 4, 5, 6, and 8-folds. There are other articles on the uniform convolutions such as [102], [103], and [81]. We used the procedure in Maghsoodloo & Hool (1983) [5] but re-developed each of the n = 2 through n = 8 convolutions of U(a, b) by a geometrical mathematical statistics method. Further, we programmed this last geometric method in Matlab in order to obtain the exact 9 through 12-fold convolutions of the U(a, b) with itself. Convolution-densities of U(a, b) are given at the appendix. The Matlab program, uses the exact n = 2 through 12 convolutions F(n)(t) and then applies the normal approximation for convolutions beyond 12. The question now arises how accurate is the normal approximation to F(n)(t), n = 13, 14, 15, ?? We used our 12 -fold convolution of the standard uniform U(0, 1) to determine the accuracy. Clearly, the partial sum T12 = 12 ii1X?? , each Xi ~ U(0, 1) and mutually independent, has a mean of 6 and variance 12(b?a)2/12 = 1, where a = 0, and maximum-life b = 1. Table 3 shows the normal approximation to F(12)(t) for intervals of 54 0.50-Stdev. The table clearly shows that the worst relative-error occurs at ? StDev, and that the normal approximation improves as Z moves toward the right-tail. The accuracy is within 2 decimals up to one StDev and 3 decimals beyond 1.49 StDevs. Therefore, we conclude that the normal approximation to each of F(n)(t), n = 13, 14, 15, ? , due to the CLT , should not have a relative error at Z = 0.50 exceeding 0.002960. Table 3: Normal Approximation to the 12-fold convolution of U(0,1) It should also be noted that the normal approximation to the uniform F(n)(t), n = 13, 14, 15, ? must be very accurate from the standpoint of the first 4 moments. Because the skewness of n-fold convolution of U(a, b) with itself is identically zero, which is identical to the Laplace- Gaussian N(n(a+b)/2, n(b?a)2/12), and hence a perfect match between the first 3 moments of ? ? ? ?? ?nT nf t f t? with those of normal. It can be proven (the proof is at the appendix) that the skewness of the partial-sum ii1 Xn nT ??? , Xi?s being iid like X, is given by 3 / 2 2 3 / 2 33 3 2 3 3 3[ X ) ( X (( ) ( ) ( ) ] ( ( ))( X))? ? ???nn nT ? ? n ? n ? ? ? nTT n = (4.9a) Further, the kurtosis of Tn is given by ? ? ? ? ? ?4 4 444 33 X /? ( X ) / + 3 ( 1 ) ( X )/ [ ] /? ? ??? ? ? ? ? ?nn=T T nn n n n, (4.9b) Z 0.5 1 1.5 2 2.5 3 3.5 4 F(12)(z) 0.689422 0.839273 0.932553 0.977724 0.994421 0.998993 0.999879 0.999991 Normal Approx. 0.691462 0.841345 0.933193 0.97725 0.99379 0.99865 0.999767 0.999968 Rel-Error 0.002960 0.002469 0.000686 -0.00049 -0.00063 -0.00034 -0.00011 -2.3E-05 55 where ?i ( i = 2, 3, 4) are the universal notation for the ith central moments. Eq. (4.9b) clearly shows that the kurtosis of the uniform ?? )( nft , for n = 13, 14, and 15, respectively, are 4? = ?1.20/13 = ?0.09231, ?1.20/14 = ?0.08571, and ?1.20/15 = ?0.08000, the amount ?1.20 being the kurtosis of a U(a, b) underlying failure density. Thus, an n = 120 is needed in order for the kurtosis of ? ??? nft to be within 0.01 of Laplace-Gaussian N(n(a+b)/2, n(b?a)2/12). Fortunately, the previous summary table clearly indicates that the normal approximation is superior at the tails, where kurtosis plays a more important role, than the middle of ? ??? nft density. In order to compute the RNIF ?t() for t > b, we used two different approximate procedures, one of which will be detailed in section 4.5. Eq. (4.1a) clearly shows that nnn 1 n 1dt F t f tdt ?? ??????? ( ) ( )( ) ( ) ( ) ; the Matlab program knows the exact 1()()ft= f(t) and uniform convolutions nft()() , for n = 2 , 3, ?, 12. For n > 12, it uses the ordinate of normal density N(n?, n?2) approximation, where ? = (a+b)/2 and ?2 = (b?a)2/12. 4.4 Approximating the Renewal Function with Unknown Convolutions Unlike the gamma, normal and uniform underlying failure distributions, the Weibull base-line distribution (except when the shape parameter ? = 1) does not have a closed-form expression for the n-fold cdf convolution F(n)(t), and hence Eq. (4.1a) cannot directly be used to calculate the renewal function M(t) for all ? > 0. When minimum-life = 0 and shape ? = 2, the Weibull specifically is called the Rayleigh pdf; we do have a closed-form function for the Rayleigh ?(s) but it cannot be inverted to yield a closed-form expression for its ?(t). Because the 56 Weibull is the most important underlying mortality density in reliability analyses, we will use the discretizing method to approximate the renewal function for three parameter Weibull distribution. The parameter ? is the minimum-life or threshold (a location parameter), ? is called the characteristic-life (? = ? ?? being the scale parameter), and ? is the slope (or shape) parameter. Further, it is well-known that the Weibull?s HZF 1() () ? 0 , t < ? ht ? t ? ,t ?? ?? ? ?? ?? ??? is an IFR (with CV < 100%) iff the slope ? > 1, h(t) = ? = 1/? is a CFR iff ? ? 1 (with CV = 100%), and it is a DFR iff 0 < ? < 1 (with CV > 100%). It must be highlighted that there have been many articles on approximating the Weibull RF such as Jiang (2007)[66], From (2001) [59], and other notables. Note that Murthy et al (2004) [99] provide an extensive treatise on Weibull Models, referring to the Weibull with zero minimum-life as the standard model. These last three authors also highlight the confusion and misconception resulting from the terminologies of intensity and hazard function for the Weibull. Jin & Gonigunta (2008) [60] first approximated the cdf of the 2-parameter Weibull (i.e., threshold ? = 0) by an optimum generalized exponential function; then they obtained the LT of the corresponding generalized exponential, which could be inverted to yield their actual Weibull RF. 4.5 Discretizing Time in Order to Approximate the Renewal Equation Because M (t) = F(t) + 0 ( ) ( )dt M t f? ? ??? and the underling distributions are herein specified, the first term on the RHS (Right-Hand Side) of M (t), F(t), can be easily computed. 57 However, the convolution integral on the RHS, 0 ( ) ( )d t M t f? ? ??? , except for rare cases, cannot in general be computed and has to be approximated. The discretization method was first applied by Xie [55], where he called his procedure ?THE RS-METHOD?, RS for Riemann-Stieltjes. However, Xie [55] used renewal-type Eq. (7b) in his RS-METHOD. The first step in the discretization is to divide the specified interval (0, t) into equal-length subintervals, and only for the sake of illustration we consider the interval (0, t = 5 weeks) and divide it into 10 subintervals (0, 0.50), (0.50, 1), ? , (4.5, 5). Note that Xie?s method does not require equal-length subintervals. Thus, the length of each subinterval in this example is 0.50?t? weeks. As a result, 5 0 (5 ) ( )dMf? ? ??? ? /210 i = 1 ( 1 ) / 2 (5 ) ( )di i Mf? ? ? ? ?? ? , where the index i = 1 pertains to the subinterval (0, 0.50), and i = 10 pertains to the last subinterval (4.5, 5). We now make use of the Mean-value Theorem for Integrals, which states: if a function f(x) is continuous over the real closed interval [a, b], then for certain there exists a real number x0 such that b a f(x)dx? ? f(x0)?(b?a), a ? x0 ? b, f(x0) being the ordinate of the integrand at x0. Because both M(t??) and the density f(t) are continuous, applying the above Mean-value Theorem for Integrals to the 4th subinterval, there exists for certain a real number ?4 such that 2 3 / 2 (5 ) ( )dMf? ? ??? ? 44( 5 ) ( ) ( 2 3 / 2 )Mf????, 3/2 ? ?4 ? 2. As a result, 5 0 (5 ) ( )dMf? ? ??? ? 58 10 =1 (5 ) ( )(1 / 2 )iii Mf???? , where 0 ? ?1 ? 0.50, 0.5 ? ?2 ? 1, ?, 4.5 ? ?10 ? 5. Clearly the exact values of ? ?5 iM ?? , i = 1, ?, 10 cannot in general be determined, and because in this example ( )(1/2)if ? = /2 ( 1)/2 ( )di i f ?? ?? = Pr[(i?1)/2 ? TTF ? i/2], it follows that 5 0 (5 ) ( )dMf? ? ??? ? /210 =1 ( 1 ) / 2 (5 ) ( )dii i i Mf? ? ? ? ?? ?. As proposed by Elsayed [35] who used the end of each subinterval, we will approximate this function in the same manner by (5 0.50 )?Miwhich results in M (5) ? F (5) + 5 0 (5 ) ( )d?? Mf? ? ? ? F (5) + /210 i = 1 ( 1 ) / 2 (5 0 .5 0 ) ( )di i M i f ?? ? ?? ?= /210 i=1 ( 1 /2 ( )di i f ?? ?? ? + /210 i = 1 ( 1 ) / 2 (5 0 .5 0 ) ( )di i M i f ?? ? ?? ?= /210 i = 1 ( 1 ) / 2 1 (5 0 . 5 0 ) ( )d[]{} ? ? ?? ?i i + M i f ?? (4.10) The above Eq. (10) is similar to that of (7.10) of Elsayed [35], where his subintervals are of length 1?t? . We first used the information M (0) ? 0 at i = 10 to calculate the last term of Eq. (10); further, at i = 9, Eq. (10) yields ( 8 1 ) / 2 8 / 2 1 (5 4 .5 ) ( )d[] ?? ?+ M f ??= 4 .5 41 (0 .5 0 ) ( )d[]?+ M f ?? . However, M (0.50) represents the expected number of renewals during an interval of length t?? 0.50. Assuming that t? is sufficiently small relative to t such that N(t) is approximately Bernoulli, then M( t?? 0.50) ? 1 (0.50)+?F 0 (0.50)R? . Hence, at i = 9, the value of the term 59 before last in Eq. (10) reduces approximately to 1 (0.50)[]+F ?4.5 4 ( )df ??? . At i = 8, the value of Eq. (10) is given by 4 3.51 (1 ) ( )d[]?+ M f ?? , where M (1) = 2 i = 1 1 (1 0 .5 0 )[]{ ?? + M i /2 ( 1)/2 ( )d } ? ? ?i i f ??, where M(0.50) has been approximated. Continuing in this manner, we backward recursively solved Eq. (4.10) to approximate M (t). The smaller ? always leads to a better approximation of M (t). In order to check the accuracy of this approximation method, we first used it to approximate M(t) at ? =1 (which is the exponential failure law with M(t) = ?t), t = 10000, minimum-life ? = 0, ? = ? = 1000 = 1/?, and t? = 50 = 0.005t, approximation yielded M(10000) ? 9.754115099857199 compared to the exact ?t = 0.001?10000 = 10, a percent relative error of ?2.459 with cpu-time = 65.112829 seconds. While, at the same exact parameters, our Matlab function at t? = 40 = 0.004t, M(10000) ? 9.802640211919197 ( a % relative error of -1.97360) with cpu-time = 567.432046. We ran the same program with same parameters by just changing t? . We observed that for smaller values of t? we really approach the exact value. The table below depicts the results. 60 Table 4: Time Discretizing Approximation Method Results t? M(t) App_M(t) Relative Error Elapsed Time (seconds) 10 10.0000000000 9.9501662508319 -0.49834% 40238.61714 25 10.0000000000 9.8760351886669 -1.23965% 2018.707614 40 10.0000000000 9.8026402119192 -1.97360% 567.432046 50 10.0000000000 9.7541150998572 -2.45885% 82.738732 100 10.0000000000 9.5162581964040 -4.83742% 65.112829 200 10.0000000000 9.0634623461009 -9.36538% 32.320588 250 10.0000000000 8.8479686771438 -11.52031% 21.066944 500 10.0000000000 7.8693868057473 -21.30613% 16.510933 1000 10.0000000000 6.3212055882856 -36.78794% 16.347242 Figure 6 shows relationship between t? and Relative Error. As t? increases Relative Error increases. This implies that smaller t? leads to more accurate result. 61 Figure 6: Relative Error versus t? Next, after approximating the Weibull RF, how do we use its M(t) to obtain a fairly accurate value of Weibull RNIF ?(t)? Because 0 ( ) ( )()( ) l i mt M t t M ttd M tt dt ? ???? ? ??? , then for sufficiently small t? > 0 the approximate ( ) ( )() M t t M t tt ??? ??? , which uses the right-hand derivative, and ( ) ( )() M t M t t tt ??? ??? , using the left-hand derivative. Because the RF is not linear but strictly increasing, our Matlab program computes both the left-and right-hand expressions and approximates ?t() by averaging the two, where t and t? are inputted by the user. It is recommended that the user inputs 0 < t? ? 0.10t. -80.00000% -70.00000% -60.00000% -50.00000% -40.00000% -30.00000% -20.00000% -10.00000% 0.00000% 10.00000% 20.00000% 0 200 400 600 800 1000 1200 ?t ve rsu s R el_ Err ?t Relative Error Linear (Relative Error) 62 As a result, all available approaches have merits and demerits in terms of complexity, computational time and accuracy. The method here is easy to implement basically for any lifetime distribution and it gives fairly accurate results for smaller subintervals. However, the downside is for smaller subintervals the computational time increases. And also it is not a closed- form approximation. Based on given parameters the program calculates both renewal and reliability measures. Therefore, this approach cannot be used in some maintenance optimization models if is desired to have closed-form expression. 4.6 Chapter Summary This chapter provided the RF and RNIF for the gamma and uniform underlying failure densities. We also devised Matlab programs that output all the renewal and reliability measures of a 3-parameter Weibull, normal, gamma, and uniform. We have highlighted that the RNIF ?t()is different from the HZF h(t) for t > 0, except in the case of CFR. 63 CHAPTER 5 5 Expected Number of Renewals for Non-Negligible Repair Unlike the previous chapters, we now assume that MTTR (Mean Time to Repair) is not negligible and that TTR has a pdf denoted as r(t). This chapter gives expected number of failures, number of cycles and availability by taking the Laplace transforms of renewal functions. 5.1 Expected Number of Failures and Cycles Let the variates X1, X2, X3,? represent TTF i be iid with the underlying failure density f(x) having means MTBF = ?x and variance 2x? ; further, let Y1, Y2, Y3 , ? represent the ith Time-to-Repair (TTRi), i = 1, 2, 3, 4,? with the same pdf r(y) having means MTTR = ?y and variance 2y? . Then, Ti = Xi + Yi represents the time between cycles (TBCs) which are also iid whose density is given by the convolution ? ? ( ) * ( )t tg f r t? , and whose Laplace transform (LT) is given by g(s) f (s) r (s)??. Clearly the mean and variance of the cycle-times Ti?s are ?x + ?y and 22xy???. As described by U. N. Bhat (1984) [34] there will be two types of renewals: (1) A transition from a Y-state (i.e., when system is under repair) to an X-state (at which the system is operating reliably), (2) A transition from an X-state (or operating-reliably-state) to a Y-state (where system will go under repair). 64 Let M1(t) represent the expected number of cycles (or number of renewals of type 1), and M2(t) represent the expected number of failures (or renewals of type 2). Then, as stated by Bhat (1984) [34] and E. A. Elsayed (2012), the LTs (Laplace-Transforms) of the two renewal functions, respectively, are given by 1 g (s ) f (s ) r (s )M (s ) s [1 g (s ) ] s [1 f (s ) r (s ) ]???? ? ? (5.1a) 2 f(s) f (s )M (s ) s [ 1 g (s ) ] s [ 1 f (s ) r (s ) ]??? ? ? (5.1b) The corresponding LTs of RNIFs (Renewal-Intensity Functions) are given by 1 f (s ) r (s )(s ) 1 f (s ) r (s )? ?? ?? , and 2 f (s )(s ) 1 f (s ) r (s )? ? ?? (5.2) As an example, suppose TTFi ~ Exp(?) and TTRi ~ Exp(r); then as has been documented by numerous other authors, s 0 f ( s ) e (/e) stt dt? ? ? ? ??? ???? and ? ?rs 0 r ( s ) r e e r / r stt dt? ???? ?? . On substituting these last 2 LTs into Eq. (5.1a), we obtain 1 rM ( s ) s [ ( s ) ( r s ) r ]?? ? ? ? ? ? = 2 2 2r r rs s (s )?? ? ???? ? ? ? ? , where r???? M1(t) = ? ?1 1M (s)?L = 1 2 2 2r r rs s ( s )? ???? ? ???????? ? ? ? ???L = 22r r e tt ???? ? ?????? , which gives the expected number of transitions from a repair-state to an operational-state (or expected number of cycles). Similarly, 65 2 f (s )M (s ) s[1 f (s ) r (s ) ]? ?? = 22 2 2 2rs s (s )? ? ???? ? ? ? ? , which upon inversion yields ? ? 22 222r e tt tM ??? ? ???? ?? ? , representing the expected number of failures during an interval of length t. Note that the limit of both renewal functions M1(t) and M2(t) as r ? ? (i.e., MTTR ? 0) is exactly equal to ?t, as expected. Further, a comparison of M2(t) with M1(t) reveals that M2(t) > M1(t) for all t > 0, which is intuitively meaningful because the expected number of failures must exceed the expected number of cycles for all t > 0. As an example, if ? = 0.0005/hour and repair-rate = 0.05, then M1(t =1000 hours) = 0.485246544456426, while M2(1000) = 0.495147534555436, and hence the availability will be shown below that A(at t = 1000 hours) = 1+0.485246544456426?0.495147534555436 = 0.990099009900990. We now obtain a general expression for the RF s M1 and M2 by inverting equations (5.1). Eq. (5.1a) shows that 1 g (s)M (s)[ 1 g (s) ] s?? 11g (s)M (s) M (s) g (s)s?? ? ? ? ? 11 t 0 M ( x ) g ( x ) dG x,M ttt ??? ?where G(t) is the cdf of f(x)*r(x) = g(x). Eq. (5.1b) now shows that 2 f(s)M (s )[ 1 g (s ) ] s?? 22f(s)M (s ) M (s ) g (s )s?? ? ? ? ? 02 t 2M ( x )F g ( x ) d xtM t t ??? ? ; thus, in general the expected number of cycles is given by 66 ? ? ? ? 01 t 1M ( x )G g ( x ) d xtM t t ??? ? (5.3a) While the expected number of failures ? ? ? ? 02 t 2M ( x )F g ( x ) d xtM t t ??? ? (5.3b) For example, suppose TBFs ~ N(?x = MTBF, 2x? ) and TTR is also N(?y, 2y? ); then TBCs ~ N(?x+?y, 22xy???). Then, M1(t) = 1 ()n tnn? ? ??? ?? , where ? =?x+?y, 22xy? ??? ? , and M1(t) gives the expected number of renewals of the first type, i.e., the expected number of cycles. However, because the system is under repair a fraction of the times, then ? ? 12 ()xnxtn nMt ?? ?? ?? ??. In order to obtain a good approximation for M2(t) and the resulting A(t), we may argue that the expected duration of time the system is under repair is given by M1(t) ?MTTR; letting t2 = t ? M1(t) ?MTTR, then Eq. (5.3b) shows that the expected number of failures is approximately given by ? ? 2 2 2( ) ( ) ? ? ??? ? ? ????? ?xx x xn tn nt tM . Clearly, M2(t) ? M1(t) for all t ? 0. 5.2 Availability Because we are assuming that a system can be either in an operational-state, or under repair, then the reliability function must be replaced by the instantaneous (or point) availability 67 function at time t, denoted A(t), which represents the Pr that a repairable unit or system is functioning reliably at time t. Thus, if there is no repair, the availability function is simply A(t) = R(t), the reliability function . However, if the component (or system) is repairable, then there are two mutually exclusive possibilities: (1) The system is reliable at t, in which case A1(t) = R(t), (2) The system fails at time x, 0 < x < t, gets renewed (or restored to almost as-good-as- new) in the interval (x, x+?x) with Pr element ?(x) dx, and then is reliable from time x to time t Trivedi (1982) [104]. This second Pr is given by ? ? 02 ( x ) d x ( x)? ??t RA tt ? . Because the above two cases are mutually exclusive, then ? ? ? ? ? ? ? ? 1 02 ( x ) ( x ) d x ? ? ? ? ?? t RA t A t A t tRt ? (5.4). Taking Laplace transform of the above Eq. (5.4) [and observing that the integral is the convolution of R(t) with ?(t)] yields A(s) = R(s) + R(s) (s)? = R(s) [1+ (s)? ]= R(s) [1 + f (s) r (s) 1 f (s) r (s)??? ] = R(s) 1 f (s) r (s)?? , (5.5) where r(t) is the density of repair-time. For the case when the TTF (of a component or a system) has a constant failure rate ? and time to repair is also exponential at the rate r, 68 s 0 R ( s ) e e ( )/stt dt? ? ? ? ?? ??? ? ; hence, the Laplace-transform of availability from Eq. (5.5) is given byA(s) = 1 / ( s ) 1 [ / ( s )][ r / ( r s )]??? ? ? ? ? = rs s[s ( r)]?? ?? = r/ ss???? ?? ? ? ? ? ?1trA ( s ) eAt ? ? ??? ? ???L where ? = ? + r, which is provided by many other authors in Reliability Engineering such as C. E. Ebeling (2010) [4], E. A. Elsayed (2012) [6], etc. For example, given that ? = ?failure-rate = 0.0005 and r = ?repair-rate = 0.05 per hour, then ? = ? + r = 0.0505 and the Pr that the unit is available (i.e., not under repair) at t = 1000 hours is given by ? ? 0 . 0 5 0 5 ( 1 0 0 0 )0 . 0 5 0 . 0 0 0 5 e10 0 . 0 5 0 5 0 .0 0 0 . 9 9 0 005 990099005 1A ????, while R(1000 hours W/O Repair) = 0.5e? = 0.60653066 < A(1000) = 0.9901. Note that in the exponential case, we can also obtain the availability function A(t) directly from Eq. (5.4) as follows: ? ? ? ? ( x )11 00 ( x ) ( x ) d x e e ( x ) d x? ? ? ? ??? ? ? ???tt ttRAR t ?tt ?, where xx11 22d r r r r r( x ) d ( x ) / d x x e edx?M ? ? ? ????? ? ? ? ?? ? ? ? ? ???? ? ?????. Upon substitution of this RNIF into the expression for A(t), we obtain ? ? ( x ) x 0 rre e ( 1 e ) d x e? ? ? ? ? ? ? ? ???? ? ? ? ?? ? ??tt t tAt , as before. As pointed out by E. A. Elsayed (2012, pp. 466-467) [6], we also observe that 69 s 0 R (s) ( )te R t dt? ?? ? = s 0 e [(1 ( )]t F t dt? ? ?? = s 0 1 e ( )s t F t dt? ?? ? =1 F(s)s? . Hildebrand (1962) [98] proves that F(s) f(s) / s? so that 1 f(s)R(s)= s? ; on substitution into Eq. (5.5) we obtain 1 f (s )A (s ) s [1 f (s ) r (s )]?? ?? = 1 f ( s )s [1 f ( s ) r ( s ) ] s [1 f ( s ) r ( s ) ]?? ? ? ? = 1 f ( s ) r ( s ) f ( s )s s [1 f ( s ) r ( s ) ] s [1 f ( s ) r ( s ) ]???? ? ? ?. Inverting these 3 LTs, we obtain ? ? ? ? ? ?12 1 A t M t M t? ? ? (5.6) for all underlying failure distribution f(t) and TTR-distribution r(t). Eq. (5.6) is the same as that of E.A. Elsayed (2012) [6] on his page 467. 5.3 Markov Analysis Note that we can also use Markov analysis, as has been done by many authors in stochastic processes, in the case of constant failure and repair rates to obtain the availability of a simple on & off- system as depicted in the following Figure: 0 1 ? r Figure 7: On & off system 70 where state ?0? represents a system in the reliable-state and ?1? represents the same system under repair. The above figure clearly shows that 0 0 1d P ( ) / P ( ) r P ( ) ,t d t t t? ?? ? where P0(t) represents the unconditional Pr of finding the systems in the operational state ?0? at time t, and similarly for P1(t). Because P1(t) = 1?P0(t) for all t, we obtain 0 0 0d P ( ) / P ( ) r (1 P )t d t t? ?? ? ? and hence 00d P ( ) / ( r ) P ( ) rt d t t? ? ? ?, or 00d P ( ) / P ( ) rt d t t? ? ?. This is a simple differential equation with the integrating factoret? . Solving and applying the boundary condition 0P ( 0) 1t ?? results in 0 rP ( ) e tt ???????, which is identical to the A(t) obtained above in 2 other different methods. It should be noted that, although the solution 0 rP ( ) e tt ??????? is valid exactly iff both failure and repair rates are constants, it can be used to obtain a rough approximate solution for the A(t) when only the MTBF and MTTR are available. For example, suppose that a system?s TTF ~ U(0, 2000 hours), while TTR is also U(10, 30 hours); then, the MFR (Mean Failure Rate) ? 1/2000 = 0.0005, and the MRR = 1/20 = 0.05. Thus, a rough approximate solution for A(t) in this case of uniform TTF and TTR is also given by 0 0 . 9 9 0 0 9rP ( ) e 901tt ???? ? ??? , while this last availability value is exact only for the exponential cases of TTF and TTR. 71 5.4 Renewal and Availability Functions when TTF is Gamma and TTR is Exponential It is well known that the LT of an underlying gamma failure density with shape ? and scale ? = 1/? is given by f (s) ? / (? + s)? ??; note that only when ? is a positive integer this last closed-form is valid. When ? is not an exact positive integer, there is no closed-form solution for the LT of a gamma density because the integration-by-parts never terminates. Thus, in the case of shape being an exact positive integer, i.e., Erlang underlying failure-density, we have: A(s) = 1 f (s)s[1 f (s) r (s)]???= ?1 (? + s) ?rs[ 1 ] rs(? + s) ? ??? ? ? ? ? = (? + s) (s r) ? (s r) s [ (? + s) (r s) ? r]? ? ??? ?? ?? . At ? = 2 this last LT reduces to A(s) = 2 22s ( 2 r ) s 2 ?rs [ s ( 2 r ) s 2 ? r]? ? ? ?? ? ? ? ? ? = 12 2 3 1r cccs s s r????, where 1r & 2r are the roots of the polynomial 22s ( 2 r ) s 2 ? r 0=? ? ? ? ? ?. Thus, 21 ( ) ( r / 2 ) rr r / 2? ? ? ? ???, 22r r / 2 ( ) ( r / 2 ) r?? ? ? ? ?? , 1 2r /c (2r? ??, 1 2 2( 2 r )c ( 2 r ) r 4 rr?? ? ? ?? ? ? ? ? , and 23 2( 2 r )c ( 2 r ) r 4 rr? ? ? ?? ? ? ? ? . Inverting back to the t-space we obtain ? ? 12rr232 r / 2 r c ec( ) e ttAt ?? ? ? ?. This last availability function clearly shows that as t ? ?, ??At ? 2r/ 2( r )?? = /2(rr )/?? = MTBF/(MTBF+MTTR), and further, A(0) ? 1, as expected. For example, if a system has an underlying gamma failure distribution with shape ? = 2, scale 72 ? = 1/? = 1000 hours and TTR has a constant repair-rate r = 0.05, then the availability at 500 hours is given by A(500) = 0.993855509027565; while the same system with minimal repair has an ? ? ? ? ( ) e5 0 0 5 0 0 0 . 9 0 9 7 9 5 9 8 9 5 6 8 9 5 0( 1 + )xt t x e dA t R xt? ? ? ? ?? ? ? ? ?? ?? ?. That is, repair will improve availability by 9.24%. The same system has an A(1000 hours) = 0.991466622031406, and R(1000 hours, no repair) = 0.735758882342885; now repair will improve availability by 34.75%. Thus, the steady-state (or long-term) availability of such a system as discussed by many other authors is A = 0.05/(0.05+0.0005) = 0.99009901. At ? = 2 the LT of expected number of cycles reduces to 2 1 2 2 2r ?M ( s ) s [ s ( 2 r ) s 2 ? r ]? ? ? ? ? ? ? = 5 6 74 2 12c c ccs s s r s r? ? ??? , where 1r & 2r are the same roots, 24 r(2 r)2)c (r? ????? , 5 (c r / 2r)?? ?? , 4 6 522c c rc r 4 r?? ?? , and 5 7 412c r cc r 4 r?? ?? . Upon inversion, we obtain ? ? 12rr671 4 5 cc e c ec ttM t t? ???. For the same parameters as above, we obtain M1(t = 10,000 hours) = 4.695618077086354 expected cycles. Similarly, it can be shown that the LT of the expected number of failures is given by 2 2 2 2 2? ( r s )M ( s ) s [ s ( 2 r ) s 2 ? r ]+? ? ? ? ? ? ? = 8 9 1 0 112 12c c c cs s r s rs? ? ???, where 228 2( r ))c ( 2r????? , 9 (c r / 2r)?? ?? , 810 192( 2 r r )rc cc4r??? ????, and 281 921 ( 2 r r ) cc cr 4 r? ? ?? ?? ?? . Upon inversion to the t-space we obtain 73 ? ? 12rr1028 119 c e c c c ettM t t ?? ? ?. The value of expected number of failures during a mission of length 10,000 hours is M2(t =10000) = 4.705519067168098, which exceeds M1(10000) = 4.695618077086354, as expected. Further, M1(10000) ?M2(10000) + 1 = 0.990099009918256, which is identical to the value availability function obtained from ???At 12r t r t( ) c 22 r / ee2r c3???? at t = 10000. Unfortunately, when TTF is Erlang at ? = 3, 4, 5 & 6 and a specified constant repair rate r, the corresponding denominators ? ? s [1 f (s)Ds r (s) ]??? has at least 2 complex roots, which are generally complex conjugate pairs. Yet, after partial-fractioning, the LT?s can be inverted to yield real-valued M1(t) and M2(t), as demonstrated below. At ? = 3, 1 f (s ) r (s )M (s ) s[1 f (s ) r (s ) ]?? ?? = 3 3 3 3 ?r rs(? + s) ?rs[ 1 ] rs(? + s) ? ? ??? = 3 33?rs [(? + s) r s) r? ]( ?? = 3 2 3 2 2 3 2?rs [ + ( 3 ? + r ) s ( 3 ? r 3 ) s ? 3 r ]s+ ? ? ? ? ? = 351 2 4 2 1 2 3ccc c cs s s r s r s r? ? ? ?? ? ? , where the root 1r will be real, while 2r and 3r will be complex conjugates, i.e., both 23rr? and 23rr? will be real numbers. In order to maintain equality in the above PFRAC (Partial Fraction), it can be shown that 3 2 1 2 3rc rr r??? , 1 2 1 3 2 2 31 2 31c r r r r + r rrc rr? ?? ; further, letting the constants ? ?1 2 1 2 3 1 1 2 1 3 2 3 c ( r r r r + r r )r r r ca ? ? ? ? ?, 21 1 2 32 ( r r +c )c ra ?? ?, then 3c , 4c , and 5c 74 are the unique solution given by C = ? ?34 15c c c Ab?? ??, where C is the 3?1 solution vector, b is a 3?1 vector b= 1 2 1 a a c ???? ????? ?? and the 3?3 matrix A = 2 3 1 3 1 2 2 3 1 3 1 2 r r r r r r r r r r r r 1 1 1 ???? ? ? ? ?? . A Matlab program was devised to obtain the expected number of cycles M1(t) as outlined above. The program also uses similar procedure as above to compute M2(t) and the resulting A(t). The Matlab program has the capability to compute the 3 renewal measures M1(t), M2(t), and A(t) for ? = 2, 3, 4, 5, 6 and 7. 75 CHAPTER 6 6 The Approximate Expected Number of Renewals for Non-Negligible Repair As in chapter 5, we assume that MTTR (Mean Time to Repair) is not negligible and that TTR (Time to Restore, or repair) has a pdf denoted as r(t) but this chapter gives the approximate number of cycles, number of failures and the resulting availability for particular distributions. Availability was explained in the previous chapter. The inverse Laplace transform of Equation 5.5 results in the point availability A(t). If the underlying distributions are not exponential, problems arise in inverting the Laplace transform [105]. Therefore numerical solutions and approximations become the only alternatives for obtaining A(t) [35]. There are numerous approximation techniques in the literature such as Sarkar & Chaudhuri (1999) [105] uses Fourier transform technique to determine the availability of a maintained system under continuous monitoring and with perfect repair policy. They also obtain closed-form expressions when the system has gamma life distribution and exponential repair time. Ananda and Gamage (2004) [106] consider statistical inference for the steady state availability of a system when repair distribution is two-parameter lognormal and failure distributions are Weibull, gamma and lognormal. There are also other papers in the literature that work on confidence limits for steady state availability of a system like [107], [108] etc. In this chapter in order to approximate availability and renewal functions two different approximation techniques are discussed. First for some cases like Weibull TTF and uniform TTR we managed to obtain the convolution of failure density f(t) and repair density r(t). Then 76 we used these convolution densities to approximate M1(t), M2(t) and A(t) by using time discretizing approximation method that was discussed in Chapter 4. However, obtaining the convolutions of f(t) with r(t) for the general classes of failure and repair distributions is not always tractable, such is the case of both TTF and TTR being Weibull. In these cases we used moment based approximation which only requires knowing the first four row moments of failure and repair distributions. ?There are a number of cases where the moments of a distribution easily obtained, but theoretical distributions are not available in closed form? [109]. And also, efficient estimators for the various moments of the underlying distribution could be calculated from the observed sample data [92]. Kambo et. al. (2012) [92], uses first three moments of failure distribution in order to approximate the renewal function for negligible repair and they conclude that the method produces exact results of the renewal function for certain important distributions like mixture of two exponential and Coxian-2. In this chapter, we propose an approximation for the evaluation of expected number of cycles, number of failures and availability based on first four row moments of failure and repair distributions where convolution of f(t) and r(t) is intractable. We conclude that the method produces very accurate results for especially large values of time t. 6.1 Weibull TBF and Uniform TTR Let the variates X1, X2, X3,? represent TTF i be iid with the underlying failure density f(x) having means MTBF = ?x and variance 2x? ; further, let Y1, Y2, Y3 , ? repres ent the ith Time-to-Repair (TTRi), i = 1, 2, 3, 4,? with the same pdf r(y) having means MTTR = ?y and 77 variance 2y? . Then, Ti = Xi + Yi represents the time between cycles (TBCs) which are also iid whose density is given by the convolution ? ? ( ) * ( )t tg f r t? , and whose Laplace transform (LT) is given by g(s) f (s) r (s)??. Clearly the mean and variance of the cycle-times Ti?s are ?x + ?y and 22xy???. Suppose the TBFs of a component (or a system) has the Weibull distribution with minimum life ? = t0 ? 0, characteristic life ? > ? = t0, and shape (or slope) ? > 0, i.e., TTF ~ W(?, ?, ?). Letting ? = 1/(? ? t0), the density of X = TBFs (Time Between Failures) is given by 000[ ( x t ) ]1( x ) [ ( x t ) ] e , t x 0. We are considering only the simpler case of the TTFF (Time to first Failure) and TTFi, i = 2, 3, 4, ? having the same Weibull distributions, and also succeeding repairs have the same identical U(a, b) distributions. Then, the Time-Between-Cycles is given by TBCs = TBF + TTR; we used a geometrical mathematical statistics method to obtain the exact convolution of f(x) with r(y), denoted g(t). The corresponding pdf of TBCs, g(t), is given below. ? ? 0 00 [ ( t ) ] 00 [ ( t ) ] [ ( t ) ] 0 1 e c t t( ) * ( ) e e c , t < { } / { } / ta t b t a , a +gt t b +f t r t b + t ? ? ? ? ? ? ? ? ? ? ? ? ? ??? ??? ?? ? ? ? ? ? ?? (6.1a) 78 The above density has no closed-form (or explicit) antiderivative, except when ? ? 1, but Matlab can integrate g(t) within any desired limits (t1, t2), a+t0 ? t1 < t2 < ?. 6.2 Uniform TBF and Weibull TTR Conversely, suppose that the TBFs of a component or system has the U(a, b) density function and its TTR has the W(?, ?, ?) density. Thus, the repair-rate function is given by ?r[r(t??]?, where ? represents the minimum repair-time. Only when ? =1, the repair-rate is constant and is denoted by r, and at ? =1 the TTR has the exponential distribution. Because most of TTR distributions in Reliability Engineering are positively-skewed, it is recommended that the value of shape ? not to exceed 3. Then, the failure density is f(x) =1/c, a ? x < b, c = b?a, and repair density is given by ? ? 1 [ r ( y ) ]y r [ r ( y ) ] e , y 0 and scale ? =1/?. The density of X = TBF (or uptime) is given by 1 ?( ()?() ? ) e 0xf x = x , x 0. We are considering only the simpler case of the TTFF (Time to first Failure) and TTFi, i = 2, 3, 4, ? having the same gamma distributions, and also succeeding repairs have the same identical U(a, b) distributions. Then, the Time-Between-Cycles is given by TBCs = TBF + TTR, and it can be proven that TBCs has the following density, which is the convolution of f(t) with r(t), and is denoted by g(t). ? ? ( ) * ( ) , <[ ( ) , ] /[ ( ) , ] [ ( ) , ] /{} c , a tgt bf t r t c b ttat a t b ???? ? ???? ??? ? ? ? ???? ? ? ? (6.2a) where () 1x 0 1 x e d x()[ ( ), ] ? ????? ?tata ? ??? ? ? is the cdf of the standard gamma density, at ?(t??). The above density has no closed-form antiderivative but Matlab can integrate g(t) for any interval within a ? t < ?. 80 6.4 Uniform TTF and Gamma TTR Conversely, suppose that the TTF of a component or system has the U(a, b) density function and its TTR has the gamma density. Then, the failure density is f(x) =1/c, a ? x < b, c = b?a, and repair density is given by ? ? 1rrx ( r () ) e 0xr = x , x 0 2 ( 2 ) 2 2() 2( 2 ) /( 2 ) / ata t cb abft a b ttc b ? ? ? ???? ? ?? ????? ? ? ? ? ? ?? ? 23 2 23 23 ( 3 ) ( 3 ) / 2 3 3 3 3 / 2 / ( 3 ) / 32 ( ) 2 2 2( 2 ) 3 a t a b f t a b t a t c a t c a t c ab ab c t tbbc ? ? ? ? ? ? ?? ? ? ? ? ?? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? 3 4 32 2 3 4 ( 4 ) 32 2 3 4 34 4 / 6 ( 3 4 1 2 4 1 2 4 4 ) / 6 ( 3 4 2 4 4 6 0 4 4 4 ) / 6 (4 43 3 2 2 () 2 ) / ( 6 ) 23 34 a t c a t c a t c a a t a b a b t a b ft a b t a b t c c a t c a t c a t c bb c bt tc a ?? ??? ? ? ? ? ? ? ? ? ? ? ? ? ? ?? ? ? ? ? ? ? ? ? ?? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? 4 5 4 3 22 3 4 5 4 3 22 3 4 5 4 3 22 3 4 ( 5 5) 5 / 2 4 ( 4 5 2 0 5 3 0 5 2 0 5 5 ) 54 4 3 2/ 2 4 ( 6 5 6 0 5 2 1 0 5 3 0 0 5 1 5 5 ) / 2 4 ( 4 5 6 0 5 3 3 0 5 7 8 0 5 () 3 2 2 3 236 5 5 ) / 2 4 a t c a t c a t c a t c a t c c a t c a t c a t a b a b t a b ft a b a ba t c a t c c a t c a t c a t c a c c bt a ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?? ? ? ? ? ?? ? ? ? ? ? 45( 5 ) / ( 2 4 ) 4 45 t a b a b tb t c b ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?? ? 109 ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? 5 6 5 4 3 22 3 4 5 6 5 4 3 223 ( 6 ) 4 5 6 5 4 32 6 1 2 0 [ 5 6 3 0 6 6 0 6 6 0 6 3 0 6 6 1 2 0 [ 5 6 6 0 6 2 7 0 6 5 65 ] 5 4 2 ]47 0 6 5 8 5 6 2 3 7 2 3 360 [ 5 6 9 0 6 6 3 0 6( 2) / / / a t c a t c a t c a t c a t c a t c c a t a b a b t a b a b ta t c a t c a t c a t c a t c c a t c a ab ft t c a t ? ? ? ? ? ? ? ? ? ? ?? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? 23 4 5 6 5 4 3 22 3 4 56 56 1 3 0 6 3 4 6 5 6 2 1 9 3 6 0 [ 5 6 1 2 0 6 1 1 4 0 6 5 3 4 0 6 ] 3 3 1 2 2 7 0 6 1 0 9 7 4 1 2 0 24 24 ( 6 ) ( 1 2 0 ) 5 ] 56 / / / c a t c a t c c a t c a t a b t a b a b t a b ab c a t c a t t c a t cc b t c b ? ? ? ? ? ? ? ? ? ? ?? ? ? ? ? ? ? ? ? ? ? ? ?? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? 6 7 6 5 4 323 24 5 6 6 5 4 323 24 7 ( 7 ) 5 76 76 6 5 2 ) 7 72 0 6 7 42 7 10 5 7 14 0 7 10 5 7 42 7 7 15 7 21 0 7 11 55 7 32 20 7 49 35 7 39 90 7 13 3 () 7 52 , ] ( 72 0 , ] / ( 72 / 0, [ ) [ / ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?? ? ? ? ? ? ? ? ? a t c a t c a t c a t c a t c a t c a t c a t c a t c a t a b a b t a b a t c a t c a t c a t a t c c c f ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? 7 6 5 4 323 24 5 6 6 5 4 323 2 56 74 10 7 21 0 7 17 85 7 78 40 7 18 79 5 7 23 52 0 7 12 08 9 15 7 42 0 7 48 30 7 29 12 0 7 96 81 0 7 16 80 00 7 11 91 8 43 4 3 3 4 34 2 ] / ( 36 0 ) , ] / ( 7 [ 20 ) , [ ? ? ? ?? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?? ? ? ?? ? a t c a t c a t c a t c a t c a t c a t c a t c a t c a t c a t b c a t c t a b a b t a b a c c ? ? ? ? ? ? ? ? ? ? ? ? 6 5 4 32 7 3 24 5 6 67 6 7 21 0 7 30 45 7 23 38 0 7 10 00 65 7 22 57 50 7 20 89 43 ( 7 ) 25 ( 72 0 ) 2 5 6 67 ] / ( 72 0 ) [ / , , ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? a t c a t c a t c a t c a t c a t b t a b c b a b t a b btt bc a c 110 ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? 7 8 7 6 5 423 324 5 6 7 8 7 6 5 423 3245 (8 6 ) 8 / 504 0 [ 7 8 56 8 168 8 280 8 280 8 168 8 56 8 8 ] 504 0 [ 21 8 336 8 218 4 8 756 0 8 1540 87 7 6 2 ( 0 8 186 48 8 1 ) 2488 / a t c a t c a t c a t c a t c a t c a t c a t c c a t c a t c a t c a t c a t c a t a t a b a b t a b ft c ?? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?? ? ? ? ? ?? ? ? ? ? ?? ? ? ? ? ? ? ? ? ? ? ? ? ?? ? ? ? ? ?? ? 78 7 6 5 423 324 5 6 7 8 7 6 5 423 3245 8 357 6 504 0 [ 8 / 144 8 / 6 5 8 / 3 9 8 256 8 / 9 53 8 488 8 / 9 247 7 / 105 6 2 5 3 ] 53 ] [ 8 / 144 2 8 / 9 3 8 199 8 / 9 96 8 73 44 7 8 / 3 / / a b t a b ab a t c c a t c a t c a t c a t c a t c a t c a t a b t c c a t c a t c a t c a t c a t c a t ? ? ? ? ? ? ? ? ?? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?? ? ? ?? ? ? ?? ? ? ?? ? ? ?? ? ? ? ? ? ? ? ? ? ? ? ? ? 6 7 8 7 6 5 423 324 5 6 7 8 7 6 5 423 344 8 642 49 / 315 ] [ 8 / 240 8 / 6 17 8 / 6 53 8 / 2 264 7 8 / 18 967 4 4 3 5 35 8 / 2 156 83 8 / 18 139 459 / 210 ] [ 7 8 336 8 6 2 8 6 88 8 781 20 8 5 / / c a t c c a t c a t c a t c a t c a t c a t c a t c c a t c a t c a b t a b a a t c a t b t a b ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?? ? ? ? ? ? ? ? ? ? ? ? ? ? ? 324 5 6 7 8 78 289 20 8 213 544 8 8 475 333 6 8 449 119 2 ] 504 0 ( 8 ) / ( 5 2 040 ) 67 78 / a b t a b ab c a t c a t c a t c c b t c tb ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?? ? ? ? ? ? ? ? ? ? ? ?? 111 ? ? ? ? ? ? ? ? ? ? ? ?? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? 8 9 4 3 2 2 3 4 4 3 2 2 3 4 9 8 7 6 5 23 4 4 ( 5 9) 9 / 403 20 [ 2 9 12 9 18 9 12 9 3 * ( 4 9 12 9 18 9 12 9 3 ) ] 403 20 [ 28 9 504 9 780 9 315 624 9 396 90 9 640 08 98 8 7 2 () 9 / a t c a t c a t c a t c a t c a t c a t c a t c a t c c a t c a t c a t c a t ca a t a b a b t a b f t t t c a ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? 32 6 7 8 9 8 7 6 2 5 4 3 3 4 5 2 6 7 8 9 8 642 60 9 367 92 9 920 7 ] 403 20 [ 56 9 151 2 9 173 88 9 111 384 9 43 7 2 6 3 6 3 5659 0 9 107 906 4 9 165 034 8 9 143 287 2 9 541 917 ] 403 20 [ 70 4 9 / / c a t c a t c c a t c a t c a t c a t c a t c a t c a t c a a b t a b a b t a t c c a b t ? ? ? ?? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?? ? ?? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? 76 2 5 4 3 3 4 5 2 6 7 8 9 8 7 6 2 54 34 252 0 9 390 60 9 340 200 9 182 133 0 9 614 628 0 9 12800340 9 15082200 9 7715619 ] 40320 [ 56 9 252 0 9 491 40 9 541 800 9 369 117 0 9 159 037 5 2 5 4 4 / c a t c a t c a t c a t c a t c a t c a t c c a t c a t c a t c a t c a a b t a b t ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? 3 5 2 6 7 8 9 8 7 6 2 5 4 3 3 4 5 2 6 7 8 09 42324660 9 63667800 9 41503131 ] 40320 [ 28 9 151 2 9 355 32 9 474 264 9 392 931 0 9 206 745 84 9 67410252 9 124449192 9 99584613 ] 40320 4 5 3 6 / / c a t c a t c a t c c a t c a t c a t c a t c a t c a t c a t c a t a cc b t a b?? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? 9 4 3 2 2 3 4 4 3 2 2 3 4 9 8 9 [ 2 9 60 9 666 9 322 8 9 572 7 * ( 4 9 132 9 163 8 9 906 0 9 188 49 ) ] 403 20 9 40 3 6 2 7 27 32 90 8 8 / / a t c a t c a t c a t c a t c a t c a t c a t a b t a b a b t a b a b t b cc b t c ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?? ?? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? 112 ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? 9 10 9 8 7 6 23 5 4 3 4 5 6 2 7 8 9 10 7 ) 9 (1 2 0 8 10 / 362 880 [ 9 10 90 10 360 10 840 10 126 0 10 126 0 10 10 9 9840 10 360 10 90 10 10 ] 362 880 [ 18 10 360 10 306 0 0 28 () 1 / a t c a t c a t c a t c a t c a t c a a t a b abt c a t c a t c a t c c a t c a t c a t t a b ft ?? ? ? ? ? ? ? ? ? ? ? ? ? ?? ? ? ? ? ? ? ? ? ?? ? ?? ? ? ?? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? 6 3 5 4 3 4 5 6 2 7 8 9 10 9 8 7 6 23 54 45 147 00 10 447 30 10 900 90 10 120 540 10 103 500 10 517 95 10 115 15 ] 181 440 [ 42 10 126 0 10 163 80 10 121 380 10 567 630 10 174 699 0 28 3 37 10 / c a t c a t c a t c a t c a t c a t c c a t c a t c a t c a t ca ab t c a t t a b ?? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?? ?? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? 3 6 2 7 8 9 10 9 8 7 6 23 5 4 3 4 5 6 2 7 553 620 10 462 042 0 10 349 114 5 10 116 946 5 ] 181 440 [ 63 10 252 0 10 441 00 10 443 100 10 281 925 0 10 118 005 30 10 325 731 00 10 57311100 10 5 3 7 4 8 6 / c a t c a t c a t c c a t c a t c a t c a t c a t c a t c a t a b b at a c t? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?? ? ? ?? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? 8 9 10 9 8 7 6 23 5 4 3 4 5 6 2 7 8 9 440 375 10 263 556 55 ] 181 440 [ 63 10 315 0 10 693 00 10 879 900 10 710 325 0 10 378 119 70 10 132 801 900 10 297063900 10 384528375 10 219738095 ] 1 4 8 6 5 5 / / c a t c c a t c a t c a t c a t c a t c a t c a t c a t c a a b t a t b c ?? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? 10 9 8 7 6 23 5 4 3 4 5 6 2 7 8 9 10 1440 [ 42 10 252 0 10 667 80 10 102 522 0 10 10042830 10 65064510 10 278704020 10 761094180 10 1202708745 10 838419985 ] 181440 5 5 4 6 46 81 37 [1 / ab c a t c a t c a t c a t c a t c a t a b t c a t c a t c a b t tc a c b a ? ? ? ? ? ? ? ?? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? 9 8 7 6 23 5 4 3 4 5 6 2 7 8 9 10 98 2 0 126 0 10 390 60 10 703 500 10 810 873 0 10 619 964 10 10 314 246 940 10 1017758700 10 1910283795 10 1582796435 ] 18144 0 [ 9 10 72 37 0 10 255 60 1 28 / a t c a t c a t c a t c a t c a t c a t c a t c a t c c at a b t a ac b ct ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? 76 3 5 4 3 4 5 6 2 7 8 9 10 9 10 0 528 360 10 700 686 0 10 618 017 40 10 362 410 440 10 1361868840 10 2974204890 10 2874204890 ] 36288 0 10 / 362 88 2 8 9 9 100 / a t c a t c a t c a t c a t c a t c a t c c bt a b t a b a b t bc ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? 113 ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? 10 11 9 10 8 92 7 6 5 3 4 5 4 3 2 6 7 8 10 1 (1 1 1 1) 11 3628800 [ 110 11 110 11 10 11 495 11 1320 11 2310 11 11 2772 11 2310 11 1320 11 495 11 11 ] 3628800 [ 45 11 10 10 9 2 () / / a t c c a t c a t a t c a t c a t c a t c a t c a t c a t c a a t a b a b t a b ft t c c at ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? 0 9 8 2 7 6 5 4 3 4 5 6 32 7 8 9 10 11 10 9 8 2 990 11 9405 11 51480 11 182490 11 440748 11 736890 11 843480 11 633105 11 281490 11 56309 ] 3628800 [ 30 11 990 11 14355 11 2 9 2 8 3 / c a t c a t c a t c a t c a t c a t c a t c a t c a t c c a a b t a t c a t c a b t ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? 76 34 5 4 3 5 6 7 2 8 9 10 11 10 9 8 2 7 3 0780 11 1656040 11 2415798 11 6130740 11 10614780 11 12020580 11 8048865 11 2421694 ] 907200 [ 105 11 4620 11 90090 11 1025640 11 7 8 3 4 5 7 / c a t c a t c a t c a t c a t c a t c a t c c at a b t a c a t c a t c b at ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? 65 45 43 67 2 8 9 10 11 10 9 8 2 7 3 58320 11 37746324 11 129664920 11 303173640 11 462563640 11 416439870 11 168171652 ] 1814400 [ 126 11 6930 11 169785 11 2439360 11 227 7 4 6 5 60 / c a t c a t c a t c a t c a t c a t c c a t c a t c a a b t b t ca a t ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? 65 45 43 67 2 8 9 10 11 10 9 8 2 7 3 430 11 144166176 11 628303830 11 1862451360 11 3597983235 11 4095278880 11 2087687723 ] 18144 00 [ 105 11 6930 11 204435 11 3548160 1 4 5 0 6 11 65 / c a t c a t c a t c a t c a t c a t c c a t c a t c a t a b t a b c a t ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? 65 45 4 3 2 6 7 8 9 10 11 10 9 8 2 3 08530 11 308490336 11 1634978730 11 5897374560 11 13861625085 11 19184198880 5 6 4 11 11879998933 ] 1814400 [ 30 11 2310 11 79695 11 1 7 621620 11 / c a t c a t c a t c a t c a t c a t c c a t c a t c a t ab ct ta a b? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? 7 6 5 45 4 3 2 6 7 8 10 11 10 9 8 2 7 34 21543060 11 195172362 11 1220779560 11 5204388420 11 14471011170 11 173642 47 08576 ] 907200 [ 45 11 3960 11 156420 11 3651120 11 5575416 38 0 11 / c a t c a t c a t c a t c a t cc a a b t t c a t c a t c a t c a ab ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? 65 5 4 3 2 6 7 8 9 10 11 10 9 8 2 7 3 581803992 11 4200171360 11 20706055920 11 66686784120 11 126660745860 11 107710566656 ] 3628800 [ 10 11 990 11 44055 11 1160280 11 20 3 0 8 2 9 253 / t c a t c a t c a t c a t c ab a t c c a t c a t c a a t at ct b ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? 65 45 4 3 2 6 7 8 9 10 11 10 11 90 11 236615148 11 1937972190 11 10861539480 11 39853850355 11 86420523090 11 84062575399 ] 3628800 11 11 36 2 9 10 12 10 18 08 / / c a t c a t c a t c a t c a t c a t a b t a b ab cc a c t tbc ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? 114 ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? 11 12 11 10 9 823 7 6 5 44 5 6 7 3289 ( 12 11 ) 0 1 12 12 399 168 00 [ 11 12 132 12 660 12 198 0 12 3960 1 12 554 4 12 554 4 12 396 0 12 198 0 12 660 12 132 12 12 ] 399 168 00 2 11 11 () / / a t c a t c a t c a t c a t c a t c a t c a t c a t ca a t a b ab ft t c a t c a t c c ?? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? 11 10 92 8 7 63 4 5 5 4 36 7 8 29 10 11 12 [ 55 12 132 0 12 138 60 12 851 40 12 344 520 12 970 200 12 194 594 4 12 278 388 0 12 278 586 0 12 185 790 0 12 743 292 12 135 156 ] 399 168 0 10 2 2 0 10 / t a t c a t c a t c a t c a t c a t c a t c a t c a t c a t c a c ab tc ab ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? 11 10 92 8 7 63 4 5 5 4 36 7 8 29 10 11 1 [ 55 12 198 0 12 316 80 12 298 320 12 184 536 0 12 790 944 0 12 24049872 12 51997440 12 78459480 12 78768800 12 47385096 12 12945728 ] 13305600 93 / a t c a t c a t c a t c a t c a t c a t c a t c a t c a t c a t ta cc b ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? 2 11 10 92 8 7 63 4 5 5 4 36 7 8 29 10 [ 55 12 264 0 12 567 60 12 722 040 12 604 692 0 12 350 750 40 12 144094104 12 420055680 12 852879060 12 1150094000 12 927890172 12 339 93 1 8 5 72 4 5 a t c a t c a t c a t c a t c a t c a t c a t c a t c a t c a b t a b at ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? 11 12 11 10 92 8 7 63 4 5 5 4 36 7 8 29 6 ] 665 280 0 [ 77 12 462 0 12 124 740 12 200 046 0 12 211 780 80 12 155 499 960 12 808780896 12 2983069320 12 7654933440 12 13029593500 12 1325 8 4 7 5 1 /cc a t c a t c a t c a t c a t c a t c a t c a t c a t c a t a b t a b ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? 10 11 12 7 5 6 6 797 328 12 610 575 528 4 ] 665 280 0/ca a c a tc b t b ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?? ? ? ? ? ? ? ? ? 115 ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? 11 10 9 2 8 7 6 3 4 5 5 ( 12 43 6 7 8 10 ) 2 9 [ 77 12 554 4 12 180 180 12 348 810 0 12 446 846 40 12 397 746 888 12 2510700192 12 11243278200 12 35024109120 12 72328491620 12 8917 () 790 481 6 12 4 a t c a t c a t c a t c a t c a t c a t c a t c a t c a t f c a t t ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? 11 12 11 10 9 8 23 7 6 5 4 5 6 43 78 9764991340 ] 6652800 [ 55 12 462 0 12 175 560 12 398 244 0 12 59902920 12 627211200 12 4664006424 12 24630254880 12 90533256660 12 220638 6 695 7 200 65 /cc a t c a t c a t c a t c a t c a t c a t c a t c a b c a t a b t ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? 2 9 10 11 12 11 10 9 2 8 7 6 3 4 5 54 67 12 320976156732 12 211242138736 ] 6652800 [ 55 12 528 0 12 229 680 12 597 432 0 12 103 221 360 12 124 348 224 0 12 10655224272 12 64929416640 12 5 7 4 8 / at c a t c c a t c a t c a t c a t c a t c a t c a t c a ab t t a b ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?? ? ? ? ? ?? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? 32 89 10 11 12 11 10 9 2 8 7 6 3 4 5 275693192280 12 7767484491200 12 1306889097096 12 994854930208 ] 13305600 [ 55 12 594 0 12 291 060 12 853 974 0 12 166 664 520 12 227 129 4 8 3 9 364 0 12 220 / c a t c a t c a t c c a t c a t c a t c a t c a t c a tb t a b a? ? ? ? ? ? ? ? ?? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? 54 67 32 89 10 11 12 11 10 9 8 23 49990424 12 152455299480 12 735516395460 12 2357542443300 12 4517350959132 12 3919268323356 ] 3 39916800 [ 11 12 132 0 12 719 40 12 235 026 8 2 10 0 12 / c a t c a t c a t c a t c a t c c a t c a t c a a b t a t c a t b? ? ? ? ? ? ? ?? ? ? ? ? ? ? ?? ? ? ?? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? 76 45 54 67 32 89 10 11 12 11 51135480 12 777906360 12 8442009576 12 65344700520 12 353483604540 12 1272457556700 12 2742649040868 12 2680731676644 ] 39916800 2 10 1 12 12 399 168 0 1 / / c a t c a t c a t c a t c a t c a t c a t c c a b t a c a t b ?? ? ? ? ? ?? ? ? ? ? ? ? ? ?? ?? ? ? 12 0 11 12a b t bc ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?? ? 116 Appendix B: 3rd and 4th Moment of Sums of iid rvs The Skewness of the Sum of n independent and Identically Distributed (iid) Variates Suppose X1, X2, ?, Xn are independently and identically distributed random variables each with means ? and variances ? ? 22 2Xi V X ? ?? ? ? ?, where each Xi is identically distributed like X. Let the nth partial sum n ii1nS X??? ; then, clearly ? ?=nE S n? and ? ? n2S 2== nV S n?? . It is well known that the ? ? ? ? ? ?33 3 3 / () /,nia S X x n a X n? ? ?? where ? ? 33 3 ,aX ?? ? and 33 , [( ) ]EX? ? ? ? the 3th central moment of X. The proof follows. 3 3 3 113 ( ) [ ( ) ] [ ( ) ] [ ]? () nnn n i i iiS E S n E X n E X? ? ???? ? ? ? ? ??? ? ? ? ? ? ?321 1 16n n ni i i ji i jEX ? X ? X ?? ? ?????? ? ? ?????? ? ?= ? ? ? ? ? ? ? ? 3 3 3 3 11[ ] [ ] [ nn i i i i i i iiiEX ? E X ? n E X ? n ? X??? ? ? ? ? ? ??? The corresponding skewness of Sn is given by; ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?3 / 23 / 2 233 3 2 3 3 3() nn i i? S n ? X ? n ? X n ? ? X ? n ? X nS ?? The Kurtosis of the Sum of n iid Variants Suppose X1, X2, ?, Xn are independently and identically distributed random variables with means ? and variances ? ? 22 2Xi V X ? ?? ? ? ?, where each Xi is identically distributed like X. It is well known that the kurtosis of each xi is equal to ? ?44 3,X???? where ? ? 444, /X? ? ? ? and 44 , [( ) ]EX? ? ? ? the 4th central moment of X. 117 Now consider the partial sum ii1S n n X??? ; our objective is to compute the 4th central moment of Sn from the known central moment of each identical Xi, i =1, 2,?, n. Clearly, the mean of Sn is given by ? ? nE S n??, the variance is given by ? ? ? ? 2 niV S n V X n? ? ?, and thus ? ? ? ? 44 ii4 ii 1 1E X E ( )S nn n nX???? ??? ? ? ???? ? ? ? ? ? ? ???? ? ? ?? ? ? ? ? ? ??? 14 2 2 i i ji 1 i 1 j42 1E ( ) ( ) ( C ) n n nX X X? ? ?? ? ? ? ??? ? ? ? ? ??? ??? ? ? ? ?4i ijn2 6 C V( ) V( )XnX X?? ??? Note that in the binomial expansion of 4n ii1(X )????????? ? , the expectation of odd products such as 312[( )( ) ]E X X???? vanish due to mutual independence of Xi and Xj for all i ? j. Hence, ? ? ? ? 444 3 ()1nS n X n n? ? ? ? ? ? Thus, ? ? 4 4 44 22( S ) ( X ) 3 ( 1 )S V ( S ) () ( X ) / + 3 ( 1 ) /nn n 4n n nn ? n n n??? ? ? ???? ? ? The corresponding kurtosis of Sn is given by, ? ? ? ? ? ?4 4 4 / .( X ) / + 3 ( 1 ) / 3 ( X ) 3 /[]nn44 Sn=SXn? n n ? n? ? ??? ? ? ? ? 118 Appendix C: Moments of the Most Common Base-line Distributions in Reliability Table 8: Parameters and Density Functions of Most Common Baseline Distributions in Reliability Lifetime Distribution Failure Density f(t) Threshold Or: Minimum-life Location Or: Shape Scale Exponential ()e t ??? ?? ? 1/? Three Parameter Weibull ()1( ) ( ) ttf t e ??? ????? ? ? ? ??? ??? ?? ? ? ??? Gamma ? ? 1 ( )[ ( ) ]() tteft ? ? ?? ???? ? ? ??? ? ? ? =1/? Lognormal 2l n ( )[]1() (): , , 2 tt tf e ????? ??? ? ???? ? ? ? ? Normal ? ? ? ? ? ?2 2/21 2 t?f t e? ? ???? ? ? Logistic ( ) / ( ) / 2e[1 e ]g; 1),( t tt ?? ??? ?? ?? ????? ? ? Loglogistic [ ln ( t ) ] / [ ln ( t ) ) ] / 2e[1 e ]1() ( t )ft ? ? ? ? ? ? ? ? ? ? ? ?????? ? ? ? ? 119 The Normal N(?, ?2); 1?? = ?; 2?? = ?2 + ?2; 3?? = 3?2?+ ?3; 4?? = 3?4 + 6?2?2 + ?4 ; 5?? = 15?4?+ 5?2?3 + ?5 . Two-Parameter Exponential pdf: 1?? = MTTF = ? + 1/?; 2?? = 2/?2 + 2?/? + ?2 The skewness ?3 = 2, while the kurtosis ?4 = 6. 3?? = 6/?3 + 6?/?2 + 3?2/? +?3 4?? = 9/? 4 + 4 31???? ?6 221???? + 3 41?? . Three-Parameter Weibull pdf: Shape =?; 1?? = MTTF = ? + (? ? ?)??[(1/?) + 1]; Note that when ? = 1, then (? ? ?) = 1/?; the 2nd raw moment is given by 2?? = (? ? ?) 2??[(2/?) + 1] + 2?(? ? ?)??[(1/?) + 1] + ?2 ; 2? = V= (? ? ?) 2?[?[(2/?) + 1] ? ?2[(1/?) + 1)] = ?2 3 3 2 3/ 2 3 2 1 1? ( 1 + ) - 3 ? ( 1 + ) ? ( 1 + ) + 2 ? ( 1 + ) 21[ ? ( 1 + ) - ? ( 1 + ) ] ? ? ? ? ?? ?? , and the kurtosis is 4 4 3 1 2 1 124( 1 ) 4 ( 1 ) ( 1 ) 6 ( 1 ) ( 1 ) 3 ( 1 ) 321 22 [ ( 1 ) ( 1 ) ] ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?? ? ? ? ? ? ? ? ? ? ? ? ? ?? = ?4-3 The 3rd raw moment is 3?? = (? ? ?) 3??[(3/?) + 1] + 3?(? ? ?)2??[(2/?) + 1] +3?2(? ? ?)??[(1/?) + 1] + ?3 4?? = ? 4?4 +4 31???? ?6 221???? + 3 41?? 120 Lognormal pdf: 2l n ( )[]1() (): , , 2 tt tf e ????? ??? ? ???? ? , where ? =Threshold, or Min-Life, ? is a location parameter, and ? is scale. The MTTF = 1?? = ? +exp(?+?2/2) and ? ? 222 2 2eT eV ? ? ? ?????= 22ex ( ) ( )p 2 e x p 2? ? ?? ? ? ? = 222e (e 1)? ? ?? ? , and 222 2 / 2 22 e 2 e? ? ? ??????? ? ? ? , the skewness is ?3 = ?3 /?3 = 22 2 3 1.5 e 3e 2 (e 1) ?? ? ?? ? > 0. ?4 = ?4 /?4 = 2 2 2 2 63 2 e 4 e 6 e 3 ( e 1) ? ? ? ? ? ? ? ? , and the kurtosis is 2 2 2 2 2 6 3 24 4 24 e 4 e 3 e/ 3 1 2 e 6( e 1 ) ? ? ? ? ?? ? ? ? ??? ?? ? ? , and the 3rd origin-moment is given by 3?? = 2 2 23 4 . 5 2 2 2 / 2 3e 3 e 3 e? ? ? ? ? ?? ? ?? ? ?? ? ? 4?? = 243 1 2 1 14 4 4 6 3? ? ? ? ?? ? ? ? ? ?? ??a Gamma pdf: ? ? 1 ( )[ ( ) ] () tteft ????? ? ? ? ??? ? ? , ? = Threshold, ? = Shape, and ? =1/? = scale. 1?? = ? +?/?; V(T) = ?/? 2 =?2 22 2 / ( / )? ? ??? ? ? ? ? ? ; ?3 =2/? , ?4 =3+6/? and the kurtosis is ?4 = 6/?, 3 3 2 2 3 3 2 2 3 23 2 / 3 / / 3 / 3 / 3 /? ? ? ? ? ? ? ? ? ??? ? ? ? ? ? ? ? ? ? ? ? ? ?, and 4?? = 243 1 2 1 14 4 4 6 3? ? ? ? ?? ? ? ? ? ?? ??a . 121 Logistic pdf: ( ) / ( ) / 2e[1 e ]g; 1),( t tt ?? ??? ?? ?? ????? ; ? ? 1E T ???? ? and ? ? 2 2( ) / 3VT ?? ??? ? ? 1ET? ?? = location = ?, ? = Scale; 2221( ) / 3 ( )???? ? ? ? ?, ? ? < t < ?; therefore, for RE- analyses ? must exceed 10(??/ 3 ); the cdf is G(t; ?, ?) = ( ) / ( ) / 11[1 e ] [1 e ]t t?? ???? ? ? ?? ?? ? ? < t < ?. The skewness is zero due to symmetry about ? and the kurtosis is 4 224 7 / 1 5( / 3 3 1 .2 0 0 0 0)? ? ? ??? , i.e., the Logistic distribution has thicker tails than the corresponding normal with mean ? and variance 2 2 2 /3?? ?? , i.e., the N(?, 2 2 2 /3???? ). For example, if the location of the mean is ? = 200 hours and scale ? = 4.5 hours, then the cdf G(t = 185; 200, 4.5) = 0.034445195666211, while ?[(185-200)/ 8.162097139054] = Pr(Z? ?1.837762984739307) = 0.033048669103432; so, it seems the left-tail of the Logistic is a bit heavier than that of the corresponding normal. At t = 180, the respective cdfs are 0.011607316445305 (Logistic), and 0.007135857827108 for the corresponding normal. So, as we move further on the tail, the Logistic seems to become even heavier than the normal. The 3rd raw moment is given by: 233 ()? ? ??? ? ? ?, and 4 2 2 44 7 ( ) / 1 5 2 ( )? ? ? ??? ? ? ? ? ?. Loglogistic pdf: A rv T has a Loglogistic pdf with threshold ? iff X = ln(T??) has logistic pdf . T = eX + ?, where X is logistic. The cdf of Loglogistic is obtained as follows: FT(t) = Pr(T ? t) = Pr(eX + ? ? t) = Pr(eX ? t??) = Pr[X ? ln(t ? ?)] = FX[ln(t ? ?)]; thus, 122 f(t) = dFT(t)/dt = dFX[ln(t ? ?)]/dt = dFX[ln(t ? ?)]/dx]?(dx/dt) = fX[at ln(t ? ?)]?(t??)?1 = = -[ ln ( - ) - ]/ -[ ln ( - ) - )] / 2 -1e[ 1 + e ]1 ? ? ( )- t ? ? ? t ? ? ? t ?? = [ ln ( ) ] / [ ln ( ) ) ] / 2e[1 e ]1() t tt ? ? ? ? ? ??? ? ? ? ? ? ???? , t > ?, ? > 0. Thus, cdf is F(t) = [ln( ) )]/1[1 e ]t ? ? ?? ? ?? = [ln ( ) )]/ 1[1 e ]t ? ? ?? ? ? ?? , ? ? t < ?, where ? is the Location-parameter and ? is the scale of the Loglogistic pdf. It seems that the 1st two moments do not generally exist from Johnson & et al. on p. 152, where the rth raw moment is given by E(Tr) = r/e ( r / ) c o s e c ( r / )?? ??? ??, where on a comparison of their Eq. (23. 89) atop their p. 152 with the above cdf of Loglogistic we must have 1/? = ?, and ?? = ?/?; then csc(r / )?? = 1/sin(r / )?? = 1/sin(r )? ; for the 1st raw moment r = 1, and 1/sin( )? = 1/0 , which does not exist. Similarly none of the raw moments for r = 2, 3, 4, 5, ? , n exist. Note that Johnson, Kotz & Balakrishnan define on their p. 151 the log-logistic variate, X, as Y = ? + ?ln(X), where Y is the standard logistic, and X is Log-logistic according to their definition. The pdf of their standard logistic Y is given by: ? ? y y2Y e(1 ef )y ? ??? , ? ? y y2Y e(1 ef )y ?? ? ? < y < ? (Eq. 23.8 J. K.,& B [110]) In our above notation, we have X = ? +? ln(T), where now X is logistic and T is Loglogistic. Then, T = e(X??)/?; this last expression clearly shows that the 3 authors (J. K. and B) [110] are using ? as location and ? as scale, and no minimum-life. Proceeding as before, we have: 123 FT(t) = Pr(T ? t) = Pr[e(X??)/? ? t ] = Pr[(X??)/? ? ln(t) = Pr[X ? ?+ ?ln(t)] = FX[?ln(t) +? ]; thus, f(t) = dFT(t)/dt = dFX[?ln(t)+ ? ]/dt = dFX[?ln(t) + ? ]/dx]?(dx/dt) = fX[at ?ln(t) + ?]?( ?/ T) = [ ln(t) ] [ ln(t) ] 2et{1 e } ? ? ?? ? ? ???? ; However, [ ln(t) ]e?? ?? = ln(t)ee?? ??? = ln(t )ee?? ??? = et?? ??? ; Thus, f(t) = 2e[1 e ]ttt ? ? ? ?? ? ? ? ?? = 2 2 2(e[1 e) ]ttttt ? ? ?? ??? ? ? ? ??? = 2e[e]ttt ? ? ? ?? ? ??? = 2e[e 1]ttt ? ? ? ???? f(t) = 2 1e (e 1) t t ? ? ? ?? ? ? t ? 0 and scale ? > 0. (Eq. 23.88, p.151 of J. K. &B) [110] The cdf of T is given by F(t) = 1 1et ?? ??? = 1(1 e )t ??? ? ?? . The raw moments of T are given atop p. 152 of [110] as follows: E(Tr) = r/ rre csc( )?? ??? ???? , where rcsc( )?? = rcosec( )?? = r1/sin( )?? . (Eq. 23.90 p.152 of [110]) From the above Eq. of [110] we obtain ? ? /e / s n )T i(E ?? ??? ???? ; this last clearly shows that the mean exists iff sin( ) 0?? ? . Further, as ? ? ? , sin()?? ? 0 and then E(T) does not exist. Further, for some 0 < ? < 1, sin()?? < 0 leading to a negative MTTF, which is not admissible. 124 We now compare the cdf of Loglogistic from Minitab, F(t) = [ln( ) )]/1[1 e ]t ? ? ?? ? ?? [ ln ( ) ) ] / 1[1 e ]t ? ? ?? ? ? ?? , ? ? t < ?, with the form provided by [110] as F(t) = 1 1et ?? ??? = ln( ) 1 1 e et ? ?? ?? = ln( ) 1 1et ? ?? ?? = ln( )11e t??? ? ?? = [ ln( ) ]11e t??? ? ?? . Because, [110] do not have minimum life, then comparing ln( )t???? against [ln(t?0) ??]/? shows that ? =1/?, and ?/? = ??. Hence, the raw moments of Minitab?s Loglogistic at zero-min-life are E(Tr) = r/e ( r / ) c sc ( r / )?? ??? ??= re (r ) cosec(r )? ????= re (r ) / sin(r )? ????; this again shows that when scale ? =1, 2, 3, 4, 5, ?. E(T) and E(T 2) do not exist. When ? = 0.50, E(T) exists but E(T2) does not. We now derive the first four moments of Loglogistic as follows. E(T) = [ l n ( ) ] / [ l n ( ) ) ] / 2 e [1 e ]() t t t d tt ? ? ? ? ? ?? ?? ? ? ? ? ? ? ? ?? ?? ; letting z= [ln( ) )] /t ? ? ??? results in dz = 1[( ) / ]t dt???? E(T) = z z2 e [1 e ] dzt ? ? ? ?? ?? ; however, ?z + ? = ln( )t ?? zet ??? ??? ; thus E(T) = zz z2 ( e e [1 e ] ) d z??? ?? ? ? ?? ? ?? = z (1 ) z2 e [1 e ] dze ??? ?? ? ? ?? ? ? ? . We now let u = ze? in the above last integral; then, du = ze dz?? , and du ( ze? ) = dz 1zue?? and 1dz u du??? 125 E(T) = (1 ) 2 0 zu [1 u ] d u ( e )e ??? ? ?? ? ?? ? = (1 ) 20 u [1 u ] d u (1 / u )e ??? ??? ?? ? = 20 u [1 u ] due ??? ??? ?? ? . Hildebrand (1962, p. 91) [98] proves that c 20 x (1 x ) dx (1 c) (1 c) ,?? ? ? ? ? ? ? ?? for all c within the open interval (?1, +1), i.e., c must lie within ?1 < c < +1, or else the integral does not exist. Therefore, from Hildebrand?s formula, we obtain 1?? = E(T) = e (1 ) (1 )?? ? ?? ? ? ? ? ? = e B e ta [ (1 ) , (1 ) ]?? ? ? ? ? ? = e [(1 ), (1 )]B?? ? ?? ? ? , where 0 < ? < 1, where B represents the Beta-function. The above Eq. clearly shows that the 1?? = E(T) of a Loglogistic exists iff 0 < ? < 1. The Beta-function is defined as Beta( , )ab = ( , )Bab = 1 11 0 x (1 x ) dxab???? = ( ) ( ) / ( )a b a b? ?? ? ?, a & b > 0. Next we derive the 2nd raw moment of the Loglogistic density. E(T2) = [ l n ( ) ] / [ l n ( ) ) ] / 2 2 e [1 e ]() t t t d tt ? ? ? ? ? ?? ? ? ? ? ? ? ? ? ? ?? ?? ; The same z-transformation yields 2?? = E(T 2) = z z22 e [1 e ] dzt ? ? ? ?? ? ?? = z z2z2 e[1 e ]dz( e )??? ? ?? ? ?? ? ??? = zz z 2 z 22 z 2 z 2 ee [ 1 e ] [ 1 e ] d z d z2 e e? ? ? ??? ?? ?? ???? ? ? ? ??? ? ? ? ??? 126 2?? = z ( 1 ) z ( 1 2 ) z 2 z 222 ee [ 1 e ] [ 1 e ] d z d z2 e e?????? ? ? ? ? ?? ?? ? ? ? ??? ???? Again, letting u = ze? results in u?1 = ez , ze dzdu ??? , or z1d z e d u u d u?? ? ? ?, we obtain z(1 ) z2 e [1 e ] dz??? ? ? ?? ?? = (1 ) 2 0 1 (1 u ) u ( u du )?? ? ? ? ?? = 20 (1 u) u du??? ?? = (1 ) (1 ) (1 ,1 )B? ? ? ?? ? ? ? ? ? ? ?. Similarly, z(1 2 ) z2 e [1 e ] dz??? ? ? ?? ?? = (1 2 ,1 2 )B ????, 0 < ? < 0.50. Hence, 2?? = E(T2) = 222 e ( 1 , 1 ) e ( 1 2 , 1 2 )BB??? ? ? ? ? ?? ? ? ? ? ?, 0 < ? < 0.50 The variance of the Loglogistic is given by 2 2 2V ( T) 2 e ( 1 , 1 ) e ( 1 2 , 1 2 ) { e [ ( 1 ) , ( 1 ) ] }B B B? ? ?? ? ? ? ? ? ? ?? ? ? ? ? ? ? ? ? ? ? = 2 2 2e ( 1 2 , 1 2 ) e [ ( 1 ) , ( 1 ) ]BB??? ? ? ?? ? ? ? ? ?2= 22e [ (1 2 , 1 2 ) (1 , 1 ) ]BB? ? ? ? ?? ? ? ? ?, 0 ? ? < 0.50; thus the variance does not exist outside the range 0 ? ? < 0.50; it does not exist at ? = 0.50, and is identically equal to zero at ? = 0, which is not permissible. On a comparison with J. K. B. Eq.(23.90, p. 152), we may deduce that V(T) = 2 2 2 2e [ 2 c s c ( 2 ) c s c ( ) ]? ? ? ? ?? ? ? ? ? ?, where csc( ) 1 / sin( )??? ? ?. To determine the skewness, we proceed as follows: 3?? = E(T 3) = z z23 e [1 e ] dzt ? ? ? ?? ? ?? ; because we know that minimum life does not impact the variance, then for simplicity we obtain the 3rd raw moment for the case of ? = 0; thus zet ???? , 127 and hence E(T3) = z z23 z 3 e [1 e ] dze ?? ? ? ? ? ?? ? ?? = z (1 3 )z23 e (1 e dze )?? ???? ?? ?? = 3e (1 3 ,1 3 )B? ????, 0 < ? < 1/3. The 3rd central moment is given by ? ?23 33e ( 1 3 , 1 3 ) 3 e ( 1 2 , 1 2 ) e ( 1 , 1 ) 2 e ( 1 , 1 )B B B B? ? ? ?? ? ? ? ? ? ? ?? ? ? ? ? ? ? ? ? ? ??? = 3 3 3 3e ( 1 3 , 1 3 ) e ( 1 2 , 1 2 ) ( 1 , 1 ) 2 e ( 1 , 1 )3B B B B? ? ?? ? ? ? ? ? ? ?? ? ? ? ? ? ? ? ? ? ?Hence, the skewness is 3 2 1 . 5 03 ( 1 3 , 1 3 ) ( 1 2 , 1 2 ) ( 1 , 1 ) 2 ( 1 , 1 ) ,[ ( 1 2 ,3 011 2 ) ( 1 31 , ) ] /B B B BBB? ? ? ? ? ? ? ??? ???? ? ? ? ? ? ? ? ? ? ?? ? ? ?? ? ? ?? The 3rd raw moment is given by 3?? = E(T 2) = 3 2 2 33 e ( 1 , 1 ) 3 e ( 1 2 , 1 2 ) e ( 1 3 , 1 3 )B B B? ? ?? ? ? ? ? ? ? ? ?? ? ? ? ? ? ? ? ? Similarly, the 4th standardized moment is ?4 = 2 24 2( 1 4 , 1 4 ) 4 ( 1 3 , 1 3 ) ( 1 , 1 ) 6 ( 1 2 , 1 2 ) ( 1 , 1 ) 3 ( 1 , 1 )( 1 2 , 1 2 ) ( 1 , 1 )[]B B B B B BBB? ? ? ? ? ? ? ? ? ? ? ?? ? ? ?? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?? ? ? ? ? , Letting (1 ,1 ) B e t a (1 , )C 1B ? ? ? ?? ? ? ? ? ?, the above reduces to ?4 = 2 24 2( 1 4 , 1 4 ) 4 C ( 1 3 , 1 3 ) 6 C ( 1 2 , 1 2 ) 3 C( 1 2 , 1 2 ) C[]B B BB? ? ? ? ? ???? ? ? ? ? ? ? ? ?? ? ? ; thus, the kurtosis is ?4 = ?4?3, only for 0 < ? < 1/4 = 0.25.