Renewal and Renewal-Intensity Functions with Minimal Repair
The renewal and renewal-intensity functions with minimal repair are explored for the Normal, Gamma, Uniform, and Weibull underlying lifetime distributions. The Normal, Gamma, and Uniform renewal, and renewal-intensity functions are derived by the convolution method. Unlike these last three failure distributions, the Weibull except at shape does not have a closed-form function for the n-fold convolution. Since the Weibull is the most important failure distribution in reliability analyses, the approximate renewal and renewal-intensity functions of Weibull were obtained by the time-discretizing method using the Mean-Value Theorem for Integrals. A Matlab program outputs all reliability and renewal measures.
Renewal (RNL) functions give 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 . Renewal functions (RNFs) have particular importance in analysis of warranty [2–6]. They have wide variety of applications in decision making such as inventory theory , supply chain planning [8, 9], continuous sampling plans [10, 11], insurance application, and sequential analysis [2, 8, 12, 13].
Since RNFs play an important role in many applications, it is important to obtain them analytically, if possible. Based on analytical approach, the RNF is the inverse of its Laplace transform , where Laplace transforms will be defined later. Blischke and Murthy  state that “the advantage of analytical method is that one can carry out parametric studies of the RNF, that is, the behavior of as a function of the parameters of the distribution.” However, for most distribution functions, obtaining the RNF analytically is complicated and even impossible . Therefore, development of computational techniques and approximations for RNFs has attracted researchers .
One of the well-known approximations given by Täcklind  is , where and are the first and second raw moments, which is generally known as the asymptotic approximation and is also cited in numerous papers such as Smith . The asymptotic approximation has a closed-form expression and thus it is easy to apply to optimization problems that involve a RNL process . However since the asymptotic approximation is not accurate for small values of , Parsa and Jin  propose better approximation by keeping the positive features of asymptotic approximations such as simplicity, closed-form expression, and independence from higher moments of an underlying distribution. Jiang  proposes an approximation for the RNF with an increasing failure rate (IFR) which is also useful in areas such as optimization where a RNF needs to be evaluated.
There are series methods available in the literature to approximate RNFs such as Smith and Leadbetter  who developed a method to compute the RNF for the Weibull by using power series expansion of where is the shape parameter of the Weibull. On the other hand, instead of using power series expansion, Lomnicki  proposes another method by using the infinite series of appropriate Poissonian functions of . There are also many other approximations available such as Xie , Smeitink and Dekker , Baxter et al. , Garg and Kalagnanam , and From . There are usually three criteria: model simplicity, applicability, and accuracy to evaluate the value of RNF approximation . Increasing the complexity may lead to more accurate approximation but may make the process complicated and difficult to implement in practice .
Furthermore, in the literature, lower and upper bounds on RNFs have been discussed, which can be used to obtain upper and lower bounds on warranty costs such as Blischke and Murthy . Marshall  provides lower and upper linear bounds on the RNF of an ordinary RNL process. Ayhan et al.  provide tight lower and upper bounds for the RNF which are based on Riemann-Stieljes integration. There are also many other studies conducted about bounds on RNF such as Barlow , Leadbetter , Ozbaykal , Xie , Ran et al. , and Politis and Koutras .
Finally, simulation can be considered as an alternate approach to estimate the value of a RNF. Brown et al.  use the Monte Carlo simulation to estimate the RNF for a RNL process with known interarrival time distribution.
This paper proposes a convolution method to obtain the Gamma and Normal renewal and renewal-intensity functions; throughout, we are assuming that repair time is negligible (or minimal) relative to Time to Failure (TTF). As a result of this assumption, the replacement time of a failed component in a system is minimal. A further objective is to obtain the renewal and RNL-intensity functions of the uniform distribution by using to convolutions and applying the normal approximation for convolutions beyond 12. Because the Weibull distribution, except at shape , does not have a closed-form function for its -fold convolution, the last objective is to approximate its renewal and RNL-intensity functions by discretizing time using the Mean-Value Theorem for Integrals. We will also highlight the differences between the renewal-intensity function and the hazard function .
Suppose that component (or system) failures occur at times measured from zero and assuming that replacement (or restoration time) is negligible relative to operational time, then represents the operating time (measured from zero) until the failure. Because the probability density function (pdf) of (or ) may be different from the intervening times , , , we consider only the simpler case of pdf of time to first failure being identical to those of intervening times . That is, all ’s are assumed iid (independently and identically distributed). Therefore, sum of the times to the first failure from zero plus the intervening times of 2nd failure until the th failure. If , then the central limit theorem (CLT) states that the distribution of approaches normality with mean , where (the expected time between successive renewals) , and with variance of , denoted by , equal to . However, if the pdf of , , with being the parent variable, is highly skewed and/or is not sufficiently large, then the exact pdf of is given by the -fold convolution of with itself denoted by . Thus, by renewal we mean either the replacement of a failed component with a brand-new one, or the case when the failed component can be repaired in a negligibly short time to an almost brand-new condition. It has been proven by many authors both in stochastic processes and reliability engineering that the RNF for the duration is given bywhere the random variable represents the number (or count) of renewals that occur during the time interval and is the cumulative distribution function (cdf) of the -fold convolution of with itself; that is, , where is the -fold convolution density of . It is also very well known that the RNI (renewal intensity) of a distribution function by definition is given by .
Authors in stochastic processes refer to as the renewal density because describes the (unconditional) probability element of a renewal during the interval ; further explanation is forthcoming in Section 5.2, while a few authors in reliability engineering, such as Ebeling  and Leemis , refer to as intensity function. Because is never a pdf over the support set of for all failure distributions, throughout this paper we will refer to it as the renewal intensity function (RNIF).
The simplest and most common renewal process is the homogeneous Poisson process (HPP), where the intervening times are exponentially distributed at the constant interrenewal (or failure) rate . Because is a constant and intervening times are iid, a Poisson process is also referred to as a homogeneous renewal process. It is also very well known that for a HPP, the pdf of interarrivals ’s is given by , and that of the time to failure (or renewal), measured from zero, is given by (the Gamma density with shape and scale ). As a result, the use of (1a) for the interval leads to the RNF: , a fact that has been known for well more than a century. Further, the RNIF for a HPP is also a constant and is given by . Further, it is also widely known that for a HPP the .
3. The Renewal and RNL-Intensity Functions for a Normal Baseline Failure Distribution
Suppose that the time between failures, , , are NID , that is, normally and independently distributed with mean time between failures (MTBF) , and process variance . Then, statistical theory indicates that time to the failure (measured from zero) is the -fold convolution of with itself, that is, time-to-the--failure . Hence, in the case of minimal repair, from (1a), the RNF is given by where universally stands for the cumulative distribution function (cdf) of the standardized normal deviate , and gives the probability of at least renewals by time . Cui and Xie  gave the same exact expression for the Normal RNF as in the above equation (2), which they used as an approximation for the Weibull renewal with shape . It should be kept in mind that the normal failure law is approximately applicable in reliability analyses only if the coefficient of variation CVT ≤ 15% 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 CVT is not sufficiently small, then the truncated Normal can qualify as a failure distribution.  discussed the RNF for the truncated Normal. From a practical standpoint, the restriction on CVT can be relaxed to less than 20%.
As an example, suppose that a cutting tool’s TTF has the lifetime distribution ( operating hours, hours2) (extracted from Example 9.5 of [36, page 225]) with minimal repair (or replacement time), where . Then, (2) shows that the expected number of renewals (or replacements) during 42 hours of use is given by . Further, for a Laplace-Gaussian nonhomogeneous process, letting , the RNIF by direct differentiation of (2) is given by or where the symbol stands for the standard normal density . The value of RNIF at 42 hours for the baseline distribution, from (4), is (at hours) = 0.07883674 failures/hour. Note that the value of the hazard rate at 42 hours is given by failures/hour, where is the reliability function at time . Because the Normal failure probability (Pr) law always has an increasing failure rate (IFR), . Ross [39, pages 426-427] proves that (least upper bound) only if the hazard-function (HZF), , is a decreasing failure-rate (i.e., “the item is improving,” [37, page 165]), and Ross further proves when is a decreasing failure rate (DFR), then for all , and as a result , and we may add that equalities can occur only at .
4. The Renewal Function for a Gamma Baseline Failure Distribution
In order to obtain the RNF and RNIF of a Gamma nonhomogeneous Poisson process (NHPP), we first resort to Laplace transformation because it is the common procedure in stochastic processes. It has been proven (see [40, pages 277–280]) in theory of stochastic processes that the Laplace transforms (LTs) of and are, respectively, given byFurther, it is also widely known that the LT of the Gamma density is given by , , and the Gamma . However, these last closed-form expressions for LTs are valid only if the shape is an exact positive integer because the integration by parts does not terminate for noninteger values of . As a result, when the shape parameter is not a positive integer, there exists no closed-form inverse-Laplace transform for the Gamma ; 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 . In fact, for the specific Erlang density with shape and scale , it is very well known that . Upon inversion to the -space, we obtain the well-known . Hence, at , the RNIF is given by , which is quite different from the corresponding Gamma (at ) IFR for . We used Matlab’s ilaplace function as an aid in order to obtain the RNF and RNIF for the Erlang at shapes and 4 which are, respectively, given below The Gamma HZF at is given by , which exceeds the above at for all .
We now use the cdf , see the excellent text by Ebeling [36, page 226, Equation 9.10], directly in order to obtain the Gamma RNF from (1a): where = Matlab’s gammainc represents the incomplete-Gamma function at point and shape . In fact, gives the cdf of the standard Gamma density at point and shape .
4.1. The Gamma Renewal Intensity Function
Next, we directly differentiate (7) in order to obtain the Gamma RNIF. That is, However, the function under the summation in (8) is simply the Gamma density with shape and scale . Then, we used our Matlab program to obtain the value of (8) at years, at shape and scale years, which yielded renewals/year. The value of the HZF at 12 years is failures/year. In order to check the validity of , we resort to the limiting form . Because the expected TTF of the Gamma density is , then for, this example, years, which yields (12 years) .
Because the Gamma density is an IFR model if and only if the shape , then for and all . Only at , the Gamma baseline failure density reduces to the exponential with CFR, the only case for which . Note that since the well-known renewal-type equation for the RNIF is given by , this last equation clearly shows that ; further for certain yields , and hence for all baseline failure distributions. Moreover, if the minimum life , then for certain .
4.2. Examining Results of the Convolution Method
Jin and Gonigunta  propose the use of generalized exponential functions to approximate the underlying Weibull and Gamma distributions and solve RNFs using Laplace transforms. Table 1 shows a comparison of the RNF from (7) with their results. They refer to as the actual RNF [obtained from Xie’s  numerical integration] and is their approximated RNF. Table 1 shows that their becomes more accurate as compared to from (7) for larger values of .
Further, it is well known from statistical theory that the skewness of Gamma density is given by and its kurtosis is , both of these clearly showing that their limiting values, in terms of shape , is zero, which are those of the corresponding Laplace-Gaussian . We compared our results, using our Matlab program, for the Gamma at , , and which yielded (Matlab gives 15 decimal accuracy), while the corresponding Normal yielded expected renewals.
5. The Renewal and RNI Functions When the Underlying TTF Distribution Is Uniform
5.1. Renewal Function When TTF Is Uniform
Suppose that the TTF of a component or system (such as a network) is uniformly distributed over the real interval [, weeks]; then, , , , , and the cdf is , weeks. Further, succeeding failures have identical failure distributions as . From a practical standpoint, the common value of minimum life . Then, the fundamental renewal equation is given by [40, pages 277–280]. Since we are considering the simplercase of time to first failure distribution being identical to those of succeeding times to failure, then whereas represents the cdf of . However, Hildebrand  proves that , the integral inside brackets representing the convolution of with . Conversely we conclude that . Upon inversion of this last LT, we obtain the widely known . Hence, (9a) can also be represented as The renewal equation of the type (9b) has been given by many well-known authors such as Xie , Murthy et al. , Leemis , and other notables.
In order to obtain the RNF for the Uniform density, we substitute into (9a), for the specific Uniform baseline distribution, for which minimum life , in order to obtain ; letting in this last equation yields The above equation (10) shows that the RNIF is given by . Solving this last differential equation and applying the boundary condition , we obtain , where time must start at zero (or the last RNL); that is, this last expression is valid only for , . Note that Ross [39, Problem 3.7, page 154] gives the same identical only for the standard underlying failure density, whose base is equal to 1. (It is worth mentioning the fact that Ross [39, page 7] uses the common notation in stochastic processes of for the reliability function at lifetime ; the authors of this paper are indebted to Ross’s text that has provided so much help and guidance in writing this paper).
Next in order to obtain the RNF for the , we transform the origin from zero to minimum life = by letting in the RNF . This yields , and hence , , and is the Uniform-density base. The corresponding RNIF is given by , .
Because the Uniform renewal function is valid only for the interval , we will obtain the -fold convolution of the -density which in turn will enable us to obtain for by making use of (1a) that uses the infinite-sum of convolution cdfs, . As stated by Olds , the convolutions of Uniform density of equal bases, , have been known since Laplace. However there are other articles on this topic since 1950. The specific convolutions of the Uniform density with itself over the interval were obtained by Maghsoodloo and Hool  only for -fold, 3, 4, 5, 6, and 8-folds. There are other articles on the Uniform convolutions such as Shiu , Killmann and Collani  and Kang et al. . We used the procedure in Maghsoodloo and Hool  but re-developed each of the through convolutions of by a geometrical mathematical statistics method (the (1983) article is available from the authors, although the 8th convolution has some typos). Further, we programmed this last geometric method in Matlab in order to obtain the exact 9 through 12-fold convolutions of the with itself. For example, our 7-fold convolution-density of with itself is given below, where
5.2. The Renewal-Intensity Function Approximation When TTF Is Uniform
The Uniform RNIF is given by only for the interval . Bartholomew  describes as the (unconditional) probability element of a renewal during the interval , and in the case of negligible repair-time, also represents the instantaneous failure intensity function at . However, as described by nearly every author in reliability engineering, the HZF gives the conditional hazard-rate at time only amongst survivors of age ; that is, . The hazard function for the baseline distribution is given by , , , which is infinite at the end of life interval , as expected. Because the uniform HZF is an IFR, then, for the uniform density, it can be proven, using the infinite geometric series for and the Maclaurin series for ; that for all .
In order to compute the RNIF for , we used two different approximate procedures. First, by directly differentiating (1a) as follows: , where the exact , and uniform convolutions , for can be calculated by our Matlab program. For , the program uses the ordinate of approximation, where and . And the second method which uses the right hand and left hand derivatives will be explained in Section 6.3.
5.3. The Uniform Approximation Results
The method we propose uses the exact through 12 convolutions and then applies the Normal approximation for convolutions beyond 12. The question now arises how accurate is the normal approximation for ? We used our 12-fold convolution of the standard Uniform to determine the accuracy. Clearly, the partial sum , each and mutually independent, has a mean of 6 and variance , where , and maximum life . The summary in Table 2 shows the normal approximation to for intervals of length . Table 2 clearly shows that the worst relative error occurs at and that the normal approximation improves as moves toward the right tail. The accuracy is within 2 decimals up to one- and 3 decimals beyond . Therefore, we conclude that the normal approximation to each of Uniform , , due to the CLT, should not have a relative error at exceeding 0.002960.
It should also be noted that the normal approximation to the Uniform , must be very accurate from the standpoint of the first 4 moments. Because the skewness of -fold convolution of with itself is identically zero, which is identical to the Laplace-Gaussian , and hence a perfect match between the first 3 moments of with those of the corresponding Normal. It can be proven (the proof is available on request) that the skewness of the partial sum , ’s being iid like , is given by Further, the kurtosis of is given by where are the universal notation for the central moments. Equation (12b) clearly shows that the kurtosis of the uniform , for , respectively, is , with the amount −1.20 being the kurtosis of a underlying failure density. Thus, an is needed in order for the kurtosis of to be within 0.01 of Laplace-Gaussian . Fortunately, summary in Table 1 clearly indicates that the Normal approximation is superior at the tails, where kurtosis may play a more important role than the middle of density [48, page 86].
6. Approximating the Renewal Function for Unknown Convolutions
6.1. The Three-Parameter Weibull Renewal Function
Unlike the Gamma, Normal, and Uniform underlying failure densities, the Weibull baseline distribution (except when the shape parameter ) does not have a closed-form expression for its -fold cdf convolution , and hence (1a) cannot directly be used to obtain the renewal function for all . When minimum life = 0 and shape , the Weibull specifically is called the Rayleigh pdf; we do have a closed-form function for the Rayleigh but it cannot be inverted to yield a closed-form expression for its .
It must be highlighted that there have been several articles on approximating the RNF such as Deligonul , and also on numerical solutions of by Tortorella , From , and other notables. The online supplement by Tortorella provides extensive references on renewal theory and applications. Note that Murthy et al.  provide an extensive treatise on Weibull Models, referring to the Weibull with zero minimum life as the standard model. Murthy et al.  state, at the bottom of their page 37, the confusion and misconception that had resulted from the plethora of terminologies for intensity and hazard functions. We will use the time-discretizing method used by Elsayed , and others such as Xie , with the aid of Matlab to obtain another approximation for .
6.2. Discretizing Time in order to Approximate the Renewal Equation
Because and the underling failure distributions are herein specified, the first term on the right-hand side (RHS) of , , can be easily computed. However, the convolution integral on the RHS, , except for rare cases, cannot in general be computed and has to be approximated. The discretization method was first applied by Xie , where he called his procedure “THE RS-METHOD,” RS for Riemann-Stieltjes. However, Xie  used the renewal type (9b) in his RS-METHOD, while we are discretizing (9a).
The first step in the discretization is to divide the specified interval into equal-length subintervals, and only for the sake of illustration, we consider the interval and divide it into 10 subintervals (0, 0.50 week), (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 weeks. As a result, , where the index pertains to the subinterval (0, 0.50) and pertains to the last subinterval (4.5, 5). We now make use of the Mean-Value Theorem for Integrals, which states the following: if a function is continuous over the real closed interval , then for certain, there exists a real number such that , , being the ordinate of the integrand at . Because both and the density are continuous, then, for example, applying the above Mean-Value Theorem for Integrals to the 4th subinterval, there exists for certain a real number such that , . As a result, , where , , . Clearly the exact values of , cannot in general be determined, and because in this example , it follows that . As proposed by Elsayed , who used the end of each subinterval and a conditional probability approach to arrive at his equation (7.10), we will approximate this last function in the same manner by which results in The above equation (13) is similar to that of (7.10) of Elsayed , where his subintervals are of length . It should be highlighted that using the interval midpoints (instead of endpoints) inorder to attain more accuracy leads to more computational complexity. We first used the information at to calculate the last term of (13); further, at , (13) yields . However, represents the expected number of renewals during an interval of length . Assuming that is sufficiently small relative to such that is approximately Bernoulli, then . Hence, at , the value of the term before last in (13) reduces approximately to . At , the value of (13) is given by , where , where has been approximated. Continuing in this manner, we solved (13) backward-recursively in order to approximate . The smaller always leads to a better approximation of .
6.3. Renewal-Intensity Function Approximation for the Weibull Distribution
Next, after approximating the Weibull RNF, how do we use its approximate to obtain a fairly accurate value of Weibull RNIF ? Because , then, for sufficiently small the approximate , which uses the right-hand derivative, and , using the left-hand derivative. Because the RNF is not linear but strictly increasing, our Matlab program computes both the left- and right-hand expressions and approximates by averaging the two, where and are inputted by the user. It is recommended that the user inputs such that the probability of 2 or more failures during an interval of length is almost zero. Further, our program shows that if and only if the shape .
6.4. Time-Discretization Accuracy
The accuracy of above discretization method was checked in two different ways.
First, we verified that at , the Weibull mean, median, and mode become almost identical at which the Weibull skewness is and Weibull kurtosis is , with these last 2 standardized-moments being very close to those of Laplace-Gaussian of identically equal to zero. Our Matlab program at , , , and yielded the Weibull (with CPU-time = 240.4183 sec), while the corresponding Normal (with MTBF = 1797.84459964 and ) resulted in (CPU ). Secondly, we ran our program for the exponential (which is the Weibull with minimum-life , and ), at and for varying values of . Table 3 depicts the comparisons against the exact expected RNLs.
The above time-discretization-method, using the Mean-Value Theorem for Integrals, can similarly be applied to any lifetime distribution and should give fairly accurate results for sufficiently small subintervals. However, Table 3 clearly shows that even at , the computational time exceeds one minute and the corresponding error may not be acceptable; further, it is not a closed-form approximation.
This paper provided the exact RNF and RNIF for the Normal, Gamma, and Uniform underlying failure densities. We have devised a Matlab program that outputs nearly all the renewal and reliability measures of a 3-parameter Weibull, Normal, Gamma, and Uniform. We have quantified the differences between and for , except in the only case of CFR for which . Further, when is an IFR, then so that , , for the four baseline failure distributions that we have studied. However, this is not quite consistent with given by some authors in reliability engineering, such as Ebeling , for a NHPP. Further, some authors such as Leemis  use the same notation for the Weibull RNIF and also use for the Weibull HZF (see his Example 6.5, page 165).
Similar works for other underlying failure densities are in the immediate future. Work is also in progress that incorporates nonnegligible repair-time requiring some knowledge ofconvolution of TTF distribution with that of time-to-restore (TTR). Such a stochastic process is referred to as an alternating renewal process [51, page 350].
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
E. A. Elsayed, Reliability Engineering, John Wiley & Sons, Hoboken, NJ, USA, 2nd edition, 2012.
E. W. Frees, “Warranty analysis and renewal function estimation,” Naval Research Logistics Quarterly, vol. 33, pp. 361–372, 1986.View at: Google Scholar
W. Blischke and D. N. P. Murthy, Warranty Cost Analysis, CRC Press, New York, NY, USA, 1st edition, 1994.
S. Karlin and A. J. Fabens, “Generalized renewal functions and stationary inventory models,” Journal of Mathematical Analysis and Applications, vol. 5, no. 3, pp. 461–487, 1962.View at: Google Scholar
G. L. Yang, “Application of renewal theory to continous sampling sampling plans,” Technometrics, vol. 25, no. 1, pp. 59–67, 1983.View at: Google Scholar
G. L. Yang, “Application of renewal theory to continous sampling sampling plans,” Naval Research Logistics Qsuarterly, vol. 32, no. 1, pp. 45–51, 1983.View at: Google Scholar
S. G. From, “Some new approximations for the renewal function,” Communications in Statistics B, vol. 30, no. 1, pp. 113–128, 2001.View at: Google Scholar
K. Politis and M. V. Koutras, “Some new bounds for the renewal function,” Probability in the Engineering and Informational Sciences, vol. 20, no. 2, pp. 231–250, 2006.View at: Google Scholar
C. E. Ebeling, An Introduction To Reliability and Maintainability Engineering, Waveland Press, Long Grove, Ill, USA, 2nd edition, 2010.
L. M. Leemis, Reliability? Probabilistic Models and Statistical Methods, 2nd edition, 2009.
S. M. Ross, Stochastic Processes, John Wiley & Sons, New York, NY, USA, 2nd edition, 1996.
U. N. Bhat, Elements of applied stochastic processes, John Wiley & Sons, New York, NY, USA, 2nd edition, 1984.
F. Hildebrand, Advanced Calculus for Applications, Prentice Hall, New York, NY, USA, 1962.
D. N. P. Murthy, M. Xie, and R. Jiang, Weibull Models, John Wiley & Sons, Hoboken, NJ, USA, 2004.
E. G. Olds, “A note on the convolution of uniform distributions,” The Annals of Mathematical Statistics, vol. 23, no. 2, pp. 282–285, 1952, 10.2307/2236455.View at: Google Scholar
S. Maghsoodloo and J. H. Hool, “On fourth moment limits required to achieve. 01 accuracy of normal approximation for symetrical linear combinations,” Journal of Alabama Academy of Science, vol. 54, no. 1, pp. 1–13, 1983.View at: Google Scholar
F. Killmann and V. E. Collani, “A note on the convolution of the uniform and related distributions and their use in quality control,” Economic Quality Control, vol. 16, no. 1, pp. 17–41, 2001.View at: Google Scholar
J. S. Kang, S. L. Kim, Y. H. Kim, and Y. S. Jang, “Generalized convolution of uniform distributions,” Journal of Applied Mathematics & Informatics, vol. 28, no. 5-6, pp. 1573–1581, 2010.View at: Google Scholar
M. G. Kendall and A. Stuart, The Advanced Theory of Statistics, vol. 1, Charles Griffin and Company, London, UK, 1963.
D. R. Cox and H. D. Miller, The Theory of Stochastic Processes, John Wiley & Sons, New York, NY, USA, 1968.