# Renewal and Renewal-Intensity Functions with Minimal Repair

Authors: Saeed Maghsoodloo; Dilcu Helvaci
Journal: Journal of Quality and Reliability Engineering
 (2014)

---

## Abstract

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β
=
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, 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.


---

## Body

## 1. Introduction

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 [1]. 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 [7], 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 RNFM
(
t
) is the inverse of its Laplace transform M
-
(
s
), where Laplace transforms will be defined later. Blischke and Murthy [6] state that “the advantage of analytical method is that one can carry out parametric studies of the RNF, that is, the behavior of M
(
t
) as a function of the parameters of the distribution.” However, for most distribution functions, obtaining the RNF analytically is complicated and even impossible [8]. Therefore, development of computational techniques and approximations for RNFs has attracted researchers [14].One of the well-known approximations given by Täcklind [15] is M
(
t
)
≃
t
/
μ
1
′
+
μ
2
′
/
2
μ
1
′
2
-
1, where μ
1
′ and μ
1
′ 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 [16]. The asymptotic approximation has a closed-form expression and thus it is easy to apply to optimization problems that involve a RNL process [8]. However since the asymptotic approximation is not accurate for small values of t, Parsa and Jin [8] 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 [17] 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 [18] who developed a method to compute the RNF for the 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 [19] proposes another method by using the infinite series of appropriate Poissonian functions of t
β. There are also many other approximations available such as Xie [20], Smeitink and Dekker [21], Baxter et al. [22], Garg and Kalagnanam [23], and From [24]. There are usually three criteria: model simplicity, applicability, and accuracy to evaluate the value of RNF approximation [17]. Increasing the complexity may lead to more accurate approximation but may make the process complicated and difficult to implement in practice [25].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 [6]. Marshall [26] provides lower and upper linear bounds on the RNF of an ordinary RNL process. Ayhan et al. [27] 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 [28], Leadbetter [29], Ozbaykal [30], Xie [31], Ran et al. [32], and Politis and Koutras [33].Finally, simulation can be considered as an alternate approach to estimate the value of a RNF. Brown et al. [34] 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 usingn
=
2 to n
=
12 convolutions and applying the normal approximation for convolutions beyond 12. Because the Weibull distribution, except at shape β
=
1, does not have a closed-form function for its n-fold convolution, the last objective is to approximate its renewal and RNL-intensity functions by discretizing time using theMean-Value Theorem for Integrals. We will also highlight the differences between the renewal-intensity function ρ
(
t
) and the hazard function h
(
t
).


## 2. Preliminaries

Suppose that component (or system) failures occur at timesT
n
(
n
=
1,2
,
3,4
,
…
) measured from zero and assuming that replacement (or restoration time) is negligible relative to operational time, then T
n represents the operating time (measured from zero) until the n
th failure. Because the probability density function (pdf) of T
1 (or X
1) may be different from the intervening times X
2
=
(
T
2
-
T
1
), X
3
=
(
T
3
-
T
2
), X
4
=
(
T
4
-
T
3
)
,
…, we consider only the simpler case of pdf of time to first failure f
1
(
t
) being identical to those of intervening times X
2
,
X
3
,
X
4
,
…. That is, all X
i’s are assumed iid (independently and identically distributed). Therefore, T
n
=
∑
i
=
1
n
X
i
= sum of the times to the first 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 T
n approaches normality with mean n
μ, where μ (the expected time between successive renewals) =
E
(
X
i
), i
=
1,2
,
3,4
,
… and with variance of T
n, denoted by V
(
T
n
), equal to n
σ
2. However, if the pdf of X
i, f
(
x
), with X being the parent variable, is highly skewed and/or n is not sufficiently large, then the exact pdf of T
n
=
∑
i
=
1
n
X
i is given by the n-fold convolution of f
(
x
) with itself denoted by f
T
n
(
t
)
=
f
(
n
)
(
t
). 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 [
0
,
t
] is given by(1a)
M
(
t
)
=
E
[
N
(
t
)
]
=
∑
n
=
1
∞
F
(
n
)
(
t
)
,
(1b)
V
[
N
(
t
)
]
=
∑
n
=
1
∞
(
2
n
-
1
)
F
(
n
)
(
t
)
-
[
M
(
t
)
]
2
,where the random variable N
(
t
) represents the number (or count) of renewals that occur during the time interval [
0
,
t
] and F
(
n
)
(
t
) is the cumulative distribution function (cdf) of the n-fold convolution of f
(
t
) with itself; that is, F
(
n
)
(
t
)
=
∫
0
t
f
(
n
)
(
τ
)
d
τ, where f
(
n
)
(
t
) is the n-fold convolution density of f
(
t
). It is also very well known that the RNI (renewal intensity) of a distribution function F
(
t
) by definition is given by ρ
(
t
)
=
d
M
(
t
)
/
d
t.Authors in stochastic processes refer toρ
(
t
) as the renewal density because ρ
(
t
)
×
Δ
t describes the (unconditional) probability element of a renewal during the interval (
t
,
t
+
Δ
t
) [35]; further explanation is forthcoming in Section 5.2, while a few authors in reliability engineering, such as Ebeling [36] and Leemis [37], refer to ρ
(
t
) as intensity function. Because ρ
(
t
) is never a pdf over the support set of f
(
t
) 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 X
i’s is given by f
(
t
)
=
λ
e
-
λ
t, and that of the time to n
th failure (or renewal), measured from zero, is given by f
T
n
(
t
)
=
f
(
n
)
(
t
)
=
(
λ
/
Γ
(
n
)
)
(
λ
t
)
n
-
1
e
-
λ
t (the Gamma density with shape n and scale β
=
1
/
λ). As a result, the use of (1a) for the interval [
0
,
t
] leads to the RNF: M
(
t
)
=
E
[
N
(
t
)
]
=
∑
n
=
1
∞
F
(
n
)
(
t
)
=
∑
n
=
1
∞
∫
x
=
0
t
(
λ
/
Γ
(
n
)
)
(
λ
x
)
n
-
1
e
-
λ
x
d
x
=
∫
x
=
0
t
λ
e
-
λ
x
∑
n
=
1
∞
(
(
λ
x
)
n
-
1
/
(
n
-
1
)
!
)
d
x
=
∫
x
=
0
t
λ
e
-
λ
x
∑
n
=
0
∞
(
(
λ
x
)
n
/
n
!
)
d
x
=
∫
x
=
0
t
λ
d
x
=
λ
t, 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 ρ
(
t
)
=
d
M
(
t
)
/
d
t
=
d
(
λ
t
)
/
d
t
=
λ. Further, it is also widely known that for a HPP the V
[
N
(
t
)
]
=
λ
t.


## 3. The Renewal and RNL-Intensity Functions for a Normal Baseline Failure Distribution

Suppose that the time between failures,X
i, i
=
1,2
,
3,4
,
…, are NID (
μ
,
σ
2
), that is, normally and independently distributed with mean time between failures (MTBF) =
μ, and process variance σ
2. Then, statistical theory indicates that time to the n
th failure (measured from zero) is the n-fold convolution of N
(
μ
,
σ
2
) with itself, that is, time-to-the-n
th-failure T
n
=
TT
F
n
~
N
(
n
μ
,
n
σ
2
). Hence, in the case of minimal repair, from (1a), the RNF is given by
(2)
M
(
t
)
=
∑
n
=
1
∞
F
(
n
)
(
t
)
=
∑
n
=
1
∞
Φ
(
t
-
n
μ
σ
n
)
,
where Φ universally stands for the cumulative distribution function (cdf) of the standardized normal deviate N
(
0,1
), and F
(
n
)
(
t
) gives the probability of at least n renewals by time t. Cui and Xie [38] 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 β
≥
3. 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. [24] 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 distributionN (μ
=
MTBF
=
15 operating hours, σ
2
=
2.25 hours2) (extracted from Example 9.5 of [36, page 225]) with minimal repair (or replacement time), where C
V
T
=
0.10. Then, (2) shows that the expected number of renewals (or replacements) during 42 hours of use is given by M
(
42
)
=
∑
n
=
1
∞
Φ
(
(
42
-
15
n
)
/
1.5
n
)
=
2.124107. Further, for a Laplace-Gaussian nonhomogeneous process, letting z
n
=
(
t
-
n
μ
)
/
σ
n, the RNIF by direct differentiation of (2) is given by
(3)
ρ
(
t
)
=
d
d
t
∑
n
=
1
∞
Φ
(
z
n
)
=
∑
n
=
1
∞
d
d
z
n
Φ
(
z
n
)
d
z
n
d
t
=
∑
n
=
1
∞
φ
(
z
n
)
1
σ
n
=
∑
n
=
1
∞
1
2
π
e
-
z
n
2
/
2
1
σ
n
;
or
(4)
ρ
(
t
)
=
∑
n
=
1
∞
1
2
π
n
σ
2
e
-
[
(
t
-
n
μ
)
/
n
σ
2
]
2
/
2
,
where the symbol φ stands for the standard normal density N
(
0,1
). The value of RNIF at 42 hours for the N
(
15,2.25
) baseline distribution, from (4), is ρ (at t
=
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 probability (Pr) law always has an increasing failure rate (IFR), ρ
(
t
)
<
h
(
t
). Ross [39, pages 426-427] proves that t
/
μ
=
the
glb
(
greatest
lower
bound
)
≤
M
(
t
)
≤
t
/
μ
+
0.50
(
C
V
T
2
-
1
)
=
the
lub (least upper bound) only if the hazard-function (HZF), h
(
t
), is a decreasing failure-rate (i.e., “the item is improving,” [37, page 165]), and Ross further proves when h
(
t
) is a decreasing failure rate (DFR), then ρ
(
t
)
≥
h
(
t
) for all t
≥
0, and as a result R
(
t
)
≥
e
-
M
(
t
), and we may add that equalities can occur only at t
=
0.


## 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 ρ
(
t
) and M
(
t
) are, respectively, given by(5a)
ℒ
{
ρ
(
t
)
}
=
ρ
-
(
s
)
=
∫
0
∞
e
-
s
t
ρ
(
t
)
d
t
=
f
-
(
s
)
1
-
f
-
(
s
)
,
s
>
0
,
(5b)
M
-
(
s
)

=

ℒ
{
M
(
t
)
}

=

∫
0
∞
e
-
s
t
M
(
t
)
d
t

=

f
-
(
s
)
s
[
1
-
f
-
(
s
)
]

=

ρ
-
(
s
)
s
.Further, it is also widely known that the LT of the Gamma density f
(
t
)
=
(
λ
/
Γ
(
α
)
)
(
λ
t
)
α
-
1
e
-
λ
t is given by f
-
(
s
)
=
λ
α
/
(
λ
+
s
)
α, ρ
-
(
s
)
=
ℒ
{
ρ
(
t
)
}
=
λ
α
/
(
(
λ
+
s
)
α
-
λ
α
), and the Gamma M
-
(
s
)
=
ℒ
{
M
(
t
)
}
=
λ
α
/
s
[
(
λ
+
s
)
α
-
λ
α
]. 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 M
-
(
s
); 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
,
4
,
…. In fact, for the specific Erlang density with shape α
≡
2 and scale β
=
1
/
λ, it is very well known that M
-
(
s
)
=
∫
0
∞
e
-
s
t
M
(
t
)
d
t
=
λ
2
/
s
2
(
s
+
2
λ
)
=
-
1
/
4
s
+
λ
/
2
s
2
+
1
/
4
(
s
+
2
λ
). Upon inversion to the t-space, we obtain the well-known M
(
t
)
=
E
[
N
(
t
)
]
=
ℒ
-
1
{
M
-
(
s
)
}
=
ℒ
-
1
{
-
1
/
4
s
+
λ
/
2
s
2
+
1
/
4
(
s
+
2
λ
)
}
=
-
1
/
4
+
(
λ
/
2
)
t
+
(
1
/
4
)
e
-
2
λ
t. Hence, at α
≡
2, the RNIF is given by ρ
(
t
)
=
d
M
(
t
)
/
d
t
=
λ
/
2
-
(
λ
/
2
)
e
-
2
λ
t, which is quite different from the corresponding Gamma (at α
≡
2) IFR h
(
t
)
=
λ
(
λ
t
)
e
-
λ
t
/
∑
k
=
0
α
-
1
(
(
λ
t
)
k
/
k
!
)
e
-
λ
t
=
λ
(
λ
t
)
/
(
1
+
λ
t
)
>
ρ
(
t
) for t
>
0. We used Matlab’silaplace function as an aid in order to obtain the RNF and RNIF for the Erlang at shapes α
≡
3 and 4 which are, respectively, given below
(6)
M
(
t
;
α
=
3
)
=
ℒ
-
1
{
λ
3
s
[
(
λ
+
s
)
3
-
λ
3
]
}
=
λ
t
3
+
e
-
(
3
λ
/
2
)
t
[
cos
⁡
(
λ
t
3
4
)
n
n
n
n
n
n
n
n
n
n
n
n
+
sin
(
λ
t
3
/
4
)
3
]
×
3
-
1
-
1
3
ρ
(
t
;
α
=
3
)
=
ℒ
-
1
{
λ
3
[
(
λ
+
s
)
3
-
λ
3
]
}
=
λ
3
{
3
4
1
-
e
-
(
3
λ
/
2
)
t
[
cos
⁡
(
λ
t
3
4
)
n
n
n
n
n
n
n
n
n
n
n
n
n
l
l
+
3
sin
⁡
(
λ
t
3
4
)
]
}
,
M
(
t
;
α
=
4
)
=
λ
t
4
+
e
-
2
λ
t
8
+
e
-
λ
t
[
cos
⁡
(
λ
t
)
+
sin
(
λ
t
)
]
4
-
3
8
,
ρ
(
t
;
α
=
4
)
=
ℒ
-
1
{
λ
4
(
λ
+
s
)
4
-
λ
4
}
=
λ
4
[
1
-
e
-
2
λ
t
-
2
e
-
λ
t
sin
(
λ
t
)
]
.
The Gamma HZF at α
≡
3 is given by h
(
t
)
=
[
λ
(
λ
t
)
2
/
Γ
(
3
)
]
/
(
1
+
λ
t
+
λ
2
t
2
/
2
)
=
λ
(
λ
t
)
2
/
(
2
+
2
λ
t
+
λ
2
t
2
), which exceeds the above ρ
(
t
) at α
≡
3 for all t
>
0.We now use the cdfF
(
n
)
(
t
)
=
Pr
(
T
n
≤
t
)
=
Pr
[
N
(
t
)
≥
n
)
], see the excellent text by Ebeling [36, page 226, Equation 9.10], directly in order to obtain the Gamma RNF from (1a):
(7)
M
(
t
)
=
∑
n
=
1
∞
F
(
n
)
(
t
)
=
∑
n
=
1
∞
∫
0
t
λ
Γ
(
n
α
)
(
λ
x
)
n
α
-
1
e
-
λ
x
d
x
=
∑
n
=
1
∞
∫
0
λ
t
u
n
α
-
1
Γ
(
n
α
)
e
-
u
d
u
=
∑
n
=
1
∞
I
(
λ
t
,
n
α
)
,
where I
(
λ
t
,
n
α
) = Matlab’s gammainc (
λ
t
,
n
α
)
=
(
1
/
Γ
(
n
α
)
)
∫
0
λ
t
u
n
α
-
1
e
-
u
d
u represents the incomplete-Gamma function at point λ
t and shape n
α
. In fact, I
(
λ
t
,
n
α
) gives the cdf of the standard Gamma density at point λ
t and shape n
α.


### 4.1. The Gamma Renewal Intensity Function

Next, we directly differentiate (7) in order to obtain the Gamma RNIF. That is,
(8)
ρ
(
t
)
=
d
M
(
t
)
d
t
=
∑
n
=
1
∞
f
(
n
)
(
t
)
=
∂
∂
t
∑
n
=
1
∞
∫
0
t
λ
Γ
(
n
α
)
(
λ
x
)
n
α
-
1
e
-
λ
x
d
x
=
∑
n
=
1
∞
∂
∂
t
∫
0
t
λ
Γ
(
n
α
)
(
λ
x
)
n
α
-
1
e
-
λ
x
d
x
=
∑
n
=
1
∞
λ
Γ
(
n
α
)
(
λ
t
)
n
α
-
1
e
-
λ
t
.
However, the function under the summation in (8) is simply the Gamma density with shape n
α and scale β
=
1
/
λ. Then, we used our Matlab program to obtain the value of (8) 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 HZF at 12 years is h
(
12
)
=
f
(
12
)
/
R
(
12
)
=
0.0193613
/
0.0765932
=
0.252781 failures/year. In order to check the validity of ρ
(
t
=
12
)
=
0.190348, we resort to the limiting form lim
⁡
t
→
∞
ρ
(
t
)
=
1
/
μ
=
1
/
MTBF. Because the expected TTF of the Gamma density is μ
=
α
×
β, then for, this example, μ
=
MTBF
=
1.5
×
3.5
=
5.25 years, which yields ρ (12 years) ≅
1
/
5.25
=
0.190476
/
year.Because the Gamma density is an IFR model if and only if the shapeα
>
1, then ρ
(
t
)
<
h
(
t
) for t
>
0 and all α
>
1. Only at α
=
1, the Gamma baseline failure density reduces to the exponential with CFR, the only case for which ρ
(
t
)
≡
h
(
t
)
≡
λ. Note that since the well-known renewal-type equation for the RNIF ρ
(
t
) is given by ρ
(
t
)
=
f
(
t
)
+
∫
0
t
ρ
(
t
-
x
)
f
(
x
)
d
x, 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 for certain ρ
(
δ
)
=
f
(
δ
)
=
h
(
δ
).


### 4.2. Examining Results of the Convolution Method

Jin and Gonigunta [25] 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 H
(
t
) as the actual RNF [obtained from Xie’s [20] numerical integration] and H
a
(
t
) is their approximated RNF. Table 1 shows that their H
(
t
) becomes more accurate as compared to M
(
t
) from (7) for larger values of t.Table 1
Comparison ofH
a
(
t
) with M
(
t
) for the Gamma renewal function.
t
β
=
1.5
β
=
3
β
=
5
H
(
t
)
H
a
(
t
)
M
(
t
)
H
(
t
)
H
a
(
t
)
M
(
t
)
H
(
t
)
H
a
(
t
)
M
(
t
)
0
0.000
0.000
0.000
0.000
0.000
0.000
0.000
0.000
0.000
1
0.517
0.508
0.517
0.081
0.065
0.081
0.004
0.001
0.004
2
1.170
1.167
1.170
0.340
0.331
0.340
0.053
0.042
0.053
3
1.834
1.833
1.834
0.665
0.664
0.665
0.186
0.179
0.186
4
2.500
2.500
2.500
0.999
1.000
0.999
0.379
0.383
0.379
5
3.167
3.167
3.167
1.333
1.333
1.333
0.592
0.597
0.592Further, 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 the corresponding Laplace-Gaussian N
(
α
/
λ
,
α
/
λ
2
). We compared our results, using our Matlab program, for the Gamma at α
=
70, β
=
15, and t
=
5000 which yielded M
(
5000
)
=
4.186155 (Matlab gives 15 decimal accuracy), while the corresponding Normal yielded M
(
5000
)
≅
4.185793 expected renewals.


## 4.1. The Gamma Renewal Intensity Function

Next, we directly differentiate (7) in order to obtain the Gamma RNIF. That is,
(8)
ρ
(
t
)
=
d
M
(
t
)
d
t
=
∑
n
=
1
∞
f
(
n
)
(
t
)
=
∂
∂
t
∑
n
=
1
∞
∫
0
t
λ
Γ
(
n
α
)
(
λ
x
)
n
α
-
1
e
-
λ
x
d
x
=
∑
n
=
1
∞
∂
∂
t
∫
0
t
λ
Γ
(
n
α
)
(
λ
x
)
n
α
-
1
e
-
λ
x
d
x
=
∑
n
=
1
∞
λ
Γ
(
n
α
)
(
λ
t
)
n
α
-
1
e
-
λ
t
.
However, the function under the summation in (8) is simply the Gamma density with shape n
α and scale β
=
1
/
λ. Then, we used our Matlab program to obtain the value of (8) 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 HZF at 12 years is h
(
12
)
=
f
(
12
)
/
R
(
12
)
=
0.0193613
/
0.0765932
=
0.252781 failures/year. In order to check the validity of ρ
(
t
=
12
)
=
0.190348, we resort to the limiting form lim
⁡
t
→
∞
ρ
(
t
)
=
1
/
μ
=
1
/
MTBF. Because the expected TTF of the Gamma density is μ
=
α
×
β, then for, this example, μ
=
MTBF
=
1.5
×
3.5
=
5.25 years, which yields ρ (12 years) ≅
1
/
5.25
=
0.190476
/
year.Because the Gamma density is an IFR model if and only if the shapeα
>
1, then ρ
(
t
)
<
h
(
t
) for t
>
0 and all α
>
1. Only at α
=
1, the Gamma baseline failure density reduces to the exponential with CFR, the only case for which ρ
(
t
)
≡
h
(
t
)
≡
λ. Note that since the well-known renewal-type equation for the RNIF ρ
(
t
) is given by ρ
(
t
)
=
f
(
t
)
+
∫
0
t
ρ
(
t
-
x
)
f
(
x
)
d
x, 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 for certain ρ
(
δ
)
=
f
(
δ
)
=
h
(
δ
).


## 4.2. Examining Results of the Convolution Method

Jin and Gonigunta [25] 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 H
(
t
) as the actual RNF [obtained from Xie’s [20] numerical integration] and H
a
(
t
) is their approximated RNF. Table 1 shows that their H
(
t
) becomes more accurate as compared to M
(
t
) from (7) for larger values of t.Table 1
Comparison ofH
a
(
t
) with M
(
t
) for the Gamma renewal function.
t
β
=
1.5
β
=
3
β
=
5
H
(
t
)
H
a
(
t
)
M
(
t
)
H
(
t
)
H
a
(
t
)
M
(
t
)
H
(
t
)
H
a
(
t
)
M
(
t
)
0
0.000
0.000
0.000
0.000
0.000
0.000
0.000
0.000
0.000
1
0.517
0.508
0.517
0.081
0.065
0.081
0.004
0.001
0.004
2
1.170
1.167
1.170
0.340
0.331
0.340
0.053
0.042
0.053
3
1.834
1.833
1.834
0.665
0.664
0.665
0.186
0.179
0.186
4
2.500
2.500
2.500
0.999
1.000
0.999
0.379
0.383
0.379
5
3.167
3.167
3.167
1.333
1.333
1.333
0.592
0.597
0.592Further, 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 the corresponding Laplace-Gaussian N
(
α
/
λ
,
α
/
λ
2
). We compared our results, using our Matlab program, for the Gamma at α
=
70, β
=
15, and t
=
5000 which yielded M
(
5000
)
=
4.186155 (Matlab gives 15 decimal accuracy), while the corresponding Normal yielded M
(
5000
)
≅
4.185793 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 [a, b weeks]; then, f
(
t
)
=
1
/
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 M
(
t
)
=
F
1
(
t
)
+
∫
0
t
M
(
t
-
τ
)
f
(
τ
)
d
τ [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(9a)
M
(
t
)
=
F
(
t
)
+
∫
0
t
M
(
t
-
x
)
f
(
x
)
d
x
,
whereas F
(
t
) represents the cdf of T
=
TTF. However, Hildebrand [41] proves that M
-
(
s
)
×
f
-
(
s
)

=

ℒ
{
∫
0
t
M
(
t
-
x
)
f
(
x
)
d
x
}
=
ℒ
{
∫
0
t
M
(
t
-
x
)
d
F
(
x
)
}, the integral inside brackets representing the convolution of M
(
t
) with f
(
t
). Conversely we conclude that f
-
(
s
)
×
M
-
(
s
)
=
ℒ
{
∫
0
t
F
(
t
-
x
)
d
M
(
x
)
}. Upon inversion of this last LT, we obtain the widely known ℒ
-
1
{
f
-
(
s
)
×
M
-
(
s
)
}
=
∫
0
t
F
(
t
-
x
)
d
M
(
x
)
=
∫
0
t
F
(
t
-
x
)
ρ
(
x
)
d
(
x
). Hence, (9a) can also be represented as
(9b)
M
(
t
)
=
F
(
t
)
+
∫
0
t
F
(
t
-
x
)
d
M
(
x
)
=
F
(
t
)
+
∫
0
t
F
(
t
-
x
)
ρ
(
x
)
d
x
.The renewal equation of the type (9b) has been given by many well-known authors such as Xie [20], Murthy et al. [42], Leemis [37], and other notables.In order to obtain the RNF for the Uniform density, we substitute into (9a), for the specific Uniform U
(
0
,
b
) baseline distribution, for which minimum life a
=
0, in order to obtain M
(
t
)
=
t
/
b
+
∫
0
t
M
(
t
-
x
)
d
x
/
b; letting t
-
x
=
τ in this last equation yields
(10)
M
(
t
)
=
t
b
+
1
b
∫
t
0
M
(
τ
)
(
-
d
τ
)
=
t
b
+
1
b
∫
0
t
M
(
τ
)
d
τ
.
The above equation (10) shows that the RNIF is given by ρ
(
t
)
=
d
M
(
t
)
/
d
t
=
1
/
b
+
M
(
t
)
/
b. Solving this last differential equation and applying the boundary condition M
(
t
=
0
)
=
0, we obtain M
(
t
)
=
e
t
/
b
-
1, where time must start at zero (or the last RNL); that is, this last expression is valid only for 0
≤
t
≤
b, b
>
0. Note that Ross [39, Problem 3.7, page 154] gives the same identical M
(
t
) only for the standard U
(
0,1
) 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 F
-
(
x
)
=
1
-
F
(
x
) for the reliability function at lifetime x; 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 theU
(
a
,
b
), we transform the origin from zero to minimum life = a
>
0 by letting τ
=
a
+
(
b
-
a
)
t
/
b
=
a
+
c
t
/
b in the RNF M
(
t
)
=
e
t
/
b
-
1. This yields M
(
τ
)
=
e
(
τ
-
a
)
/
c
-
1, and hence M
(
t
)
=
e
(
t
-
a
)
/
c
-
1, 0
≤
a
≤
t
<
b, and c
=
b
-
a
>
0 is the Uniform-density base. The corresponding RNIF is given by ρ
(
t
)
=
[
e
(
t
-
a
)
/
c
]
/
c, 0
≤
a
≤
t
<
b.Because the Uniform renewal functionM
(
t
)
=
e
(
t
-
a
)
/
c
-
1 is valid only for the interval [
a
,
b
], we will obtain the n-fold convolution of the U
(
a
,
b
)-density which in turn will enable us to obtain M
(
t
) for t
>
b by making use of (1a) that uses the infinite-sum of convolution cdfs, F
(
n
)
(
t
). As stated by Olds [43], the convolutions of Uniform density of equal bases, c, 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 [
-
1
/
2,1
/
2
] were obtained by Maghsoodloo and Hool [44] only for n
=
2-fold, 3, 4, 5, 6, and 8-folds. There are other articles on the Uniform convolutions such as Shiu [45], Killmann and Collani [46] and Kang et al. [47]. We used the procedure in Maghsoodloo and Hool [44] but re-developed each of the n
=
2 through n
=
8 convolutions of U
(
a
,
b
) 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 U
(
a
,
b
) with itself. For example, our 7-fold convolution-density of U
(
a
,
b
) with itself is given below, where c
=
b
-
a
>
0
(11)
f
(
7
)
(
t
)
=
{
(
7
a
-
t
)
6
(
720
c
7
)
,
7
a
≤
t
≤
6
a
+
b
-
[
6
(
7
a
-
t
)
6
+
42
c
(
7
a
-
t
)
5
+
105
c
2
(
7
a
-
t
)
4
+
140
c
3
(
7
a
-
t
)
3
6
a
+
b
≤
t
≤
5
a
+
2
b
+
105
c
4
(
7
a
-
t
)
2
+
42
c
5
(
7
a
-
t
)
+
7
c
6
]
×
(
720
c
7
)
-
1
,
[
15
(
7
a
-
t
)
6
+
210
c
(
7
a
-
t
)
5
+
1155
c
2
(
7
a
-
t
)
4
+
3220
c
3
(
7
a
-
t
)
3
5
a
+
2
b
≤
t
≤
4
a
+
3
b
+
4935
c
4
(
7
a
-
t
)
2
+
3990
c
5
(
7
a
-
t
)
+
1337
c
6
]
×
(
720
c
7
)
-
1
,
-
[
10
(
7
a
-
t
)
6
+
210
c
(
7
a
-
t
)
5
+
1785
c
2
(
7
a
-
t
)
4
+
7840
c
3
(
7
a
-
t
)
3
4
a
+
3
b
≤
t
≤
3
a
+
4
b
+
18795
c
4
(
7
a
-
t
)
2
+
23520
c
5
(
7
a
-
t
)
+
12089
c
6
]
×
(
360
c
7
)
-
1
,
[
15
(
7
a
-
t
)
6
+
420
c
(
7
a
-
t
)
5
+
4830
c
2
(
7
a
-
t
)
4
+
29120
c
3
(
7
a
-
t
)
3
3
a
+
4
b
≤
t
≤
2
a
+
5
b
+
96810
c
4
(
7
a
-
t
)
2
+
168000
c
5
(
7
a
-
t
)
+
119182
c
6
]
×
(
720
c
7
)
-
1
,
-
[
6
(
7
a
-
t
)
6
+
210
c
(
7
a
-
t
)
5
+
3045
c
2
(
7
a
-
t
)
4
+
23380
c
3
(
7
a
-
t
)
3
2
a
+
5
b
≤
t
≤
a
+
6
b
+
100065
c
4
(
7
a
-
t
)
2
+
225750
c
5
(
7
a
-
t
)
+
208943
c
6
]
×
(
720
c
7
)
-
1
,
(
7
b
-
t
)
6
(
720
c
7
)
,
a
+
6
b
≤
t
≤
7
b
.


### 5.2. The Renewal-Intensity Function Approximation When TTF Is Uniform

The Uniform RNIF is given byρ
(
t
)
=
[
e
(
t
-
a
)
/
c
]
/
c only for the interval 0
≤
a
≤
t
<
b. Bartholomew [35] describes ρ
(
t
)
×
Δ
t as 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 at t. However, as described by nearly every author in reliability engineering, the HZF h
(
t
) gives the conditional hazard-rate at time t only amongst survivors of age t; that is, h
(
t
)
×
Δ
t
=
Pr
(
t
≤
T
≤
t
+
Δ
t
)
/
R
(
t
). The hazard function for the U
(
a
,
b
) baseline distribution is given by h
(
t
)
=
1
/
(
b
-
t
), a
≤
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 geometric series for h
(
t
)
=
(
1
/
b
)
/
(
1
-
t
/
b
) and the Maclaurin series for e
t
/
b; that ρ
(
t
)
<
h
(
t
) for all a
<
t
≤
b.In order to compute the RNIFρ
(
t
) for t
>
b, we used two different approximate procedures. First, by directly differentiating (1a) as follows: ρ
(
t
)
=
(
d
/
d
t
)
∑
n
=
1
∞
F
(
n
)
(
t
)
=
∑
n
=
1
∞
f
(
n
)
(
t
), where the exact f
(
1
)
(
t
)
=
f
(
t
), and uniform convolutions f
(
n
)
(
t
), for n
=
2,3
,
…
,
12 can be calculated by our Matlab program. For n
>
12, the program uses the ordinate of N
(
n
μ
,
n
σ
2
) approximation, where μ
=
(
a
+
b
)
/
2 and σ
2
=
(
b
-
a
)
2
/
12. 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 exactn
=
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 for 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 T
12
=
∑
i
=
1
12
X
i, each X
i
~
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. The summary in Table 2 shows the normal approximation to F
(
12
)
(
t
) for intervals of length 0.50
-
σ. Table 2 clearly shows that the worst relative error occurs at σ
/
2 and that the normal approximation improves as Z moves toward the right tail. The accuracy is within 2 decimals up to one-σ and 3 decimals beyond 1.49
σ. Therefore, we conclude that the normal approximation to each of Uniform 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 2
Normal approximation to the 12-fold convolution ofU
(
0,1
).
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.3
E
-
05It should also be noted that the normal approximation to the UniformF
(
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 f
T
n
(
t
)
=
f
(
n
)
(
t
) with those of the corresponding Normal. It can be proven (the proof is available on request) that the skewness of the partial sum T
n
=
∑
i
=
1
n
X
i, X
i’s being iid like X, is given by(12a)
β
3
(
T
n
)
=
μ
3
(
T
n
)
[
μ
2
(
T
n
)
]
3
/
2
=
n
μ
3
(
X
)
(
n
σ
2
)
3
/
2
=
μ
3
(
X
)
(
σ
3
n
)
=
β
3
(
X
)
n
.
Further, the kurtosis of T
n is given by
(12b)
β
4
(
T
n
)
=
α
4
(
X
)
n
+
3
(
n
-
1
)
n
-
3
⟶
β
4
(
T
n
)
=
[
α
4
(
X
)
-
3
]
n
=
β
4
(
X
)
n
,where μ
i
(
i
=
2,3
,
4
) are the universal notation for the i
th central moments. Equation (12b) clearly shows that the kurtosis of the uniform f
(
n
)
(
t
), for n
=
13,14
,
and
15, respectively, is β
4
=
-
1.20
/
13
=
-
0.09231
,
-
1.20
/
14
=
-
0.08571
,
and
-
1.20
/
15
=
-
0.08000, with 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 f
(
n
)
(
t
) to be within 0.01 of Laplace-Gaussian N
(
n
(
a
+
b
)
/
2
,
n
(
b
-
a
)
2
/
12
). 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 f
(
n
)
(
t
) density [48, page 86].


## 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 [a, b weeks]; then, f
(
t
)
=
1
/
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 M
(
t
)
=
F
1
(
t
)
+
∫
0
t
M
(
t
-
τ
)
f
(
τ
)
d
τ [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(9a)
M
(
t
)
=
F
(
t
)
+
∫
0
t
M
(
t
-
x
)
f
(
x
)
d
x
,
whereas F
(
t
) represents the cdf of T
=
TTF. However, Hildebrand [41] proves that M
-
(
s
)
×
f
-
(
s
)

=

ℒ
{
∫
0
t
M
(
t
-
x
)
f
(
x
)
d
x
}
=
ℒ
{
∫
0
t
M
(
t
-
x
)
d
F
(
x
)
}, the integral inside brackets representing the convolution of M
(
t
) with f
(
t
). Conversely we conclude that f
-
(
s
)
×
M
-
(
s
)
=
ℒ
{
∫
0
t
F
(
t
-
x
)
d
M
(
x
)
}. Upon inversion of this last LT, we obtain the widely known ℒ
-
1
{
f
-
(
s
)
×
M
-
(
s
)
}
=
∫
0
t
F
(
t
-
x
)
d
M
(
x
)
=
∫
0
t
F
(
t
-
x
)
ρ
(
x
)
d
(
x
). Hence, (9a) can also be represented as
(9b)
M
(
t
)
=
F
(
t
)
+
∫
0
t
F
(
t
-
x
)
d
M
(
x
)
=
F
(
t
)
+
∫
0
t
F
(
t
-
x
)
ρ
(
x
)
d
x
.The renewal equation of the type (9b) has been given by many well-known authors such as Xie [20], Murthy et al. [42], Leemis [37], and other notables.In order to obtain the RNF for the Uniform density, we substitute into (9a), for the specific Uniform U
(
0
,
b
) baseline distribution, for which minimum life a
=
0, in order to obtain M
(
t
)
=
t
/
b
+
∫
0
t
M
(
t
-
x
)
d
x
/
b; letting t
-
x
=
τ in this last equation yields
(10)
M
(
t
)
=
t
b
+
1
b
∫
t
0
M
(
τ
)
(
-
d
τ
)
=
t
b
+
1
b
∫
0
t
M
(
τ
)
d
τ
.
The above equation (10) shows that the RNIF is given by ρ
(
t
)
=
d
M
(
t
)
/
d
t
=
1
/
b
+
M
(
t
)
/
b. Solving this last differential equation and applying the boundary condition M
(
t
=
0
)
=
0, we obtain M
(
t
)
=
e
t
/
b
-
1, where time must start at zero (or the last RNL); that is, this last expression is valid only for 0
≤
t
≤
b, b
>
0. Note that Ross [39, Problem 3.7, page 154] gives the same identical M
(
t
) only for the standard U
(
0,1
) 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 F
-
(
x
)
=
1
-
F
(
x
) for the reliability function at lifetime x; 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 theU
(
a
,
b
), we transform the origin from zero to minimum life = a
>
0 by letting τ
=
a
+
(
b
-
a
)
t
/
b
=
a
+
c
t
/
b in the RNF M
(
t
)
=
e
t
/
b
-
1. This yields M
(
τ
)
=
e
(
τ
-
a
)
/
c
-
1, and hence M
(
t
)
=
e
(
t
-
a
)
/
c
-
1, 0
≤
a
≤
t
<
b, and c
=
b
-
a
>
0 is the Uniform-density base. The corresponding RNIF is given by ρ
(
t
)
=
[
e
(
t
-
a
)
/
c
]
/
c, 0
≤
a
≤
t
<
b.Because the Uniform renewal functionM
(
t
)
=
e
(
t
-
a
)
/
c
-
1 is valid only for the interval [
a
,
b
], we will obtain the n-fold convolution of the U
(
a
,
b
)-density which in turn will enable us to obtain M
(
t
) for t
>
b by making use of (1a) that uses the infinite-sum of convolution cdfs, F
(
n
)
(
t
). As stated by Olds [43], the convolutions of Uniform density of equal bases, c, 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 [
-
1
/
2,1
/
2
] were obtained by Maghsoodloo and Hool [44] only for n
=
2-fold, 3, 4, 5, 6, and 8-folds. There are other articles on the Uniform convolutions such as Shiu [45], Killmann and Collani [46] and Kang et al. [47]. We used the procedure in Maghsoodloo and Hool [44] but re-developed each of the n
=
2 through n
=
8 convolutions of U
(
a
,
b
) 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 U
(
a
,
b
) with itself. For example, our 7-fold convolution-density of U
(
a
,
b
) with itself is given below, where c
=
b
-
a
>
0
(11)
f
(
7
)
(
t
)
=
{
(
7
a
-
t
)
6
(
720
c
7
)
,
7
a
≤
t
≤
6
a
+
b
-
[
6
(
7
a
-
t
)
6
+
42
c
(
7
a
-
t
)
5
+
105
c
2
(
7
a
-
t
)
4
+
140
c
3
(
7
a
-
t
)
3
6
a
+
b
≤
t
≤
5
a
+
2
b
+
105
c
4
(
7
a
-
t
)
2
+
42
c
5
(
7
a
-
t
)
+
7
c
6
]
×
(
720
c
7
)
-
1
,
[
15
(
7
a
-
t
)
6
+
210
c
(
7
a
-
t
)
5
+
1155
c
2
(
7
a
-
t
)
4
+
3220
c
3
(
7
a
-
t
)
3
5
a
+
2
b
≤
t
≤
4
a
+
3
b
+
4935
c
4
(
7
a
-
t
)
2
+
3990
c
5
(
7
a
-
t
)
+
1337
c
6
]
×
(
720
c
7
)
-
1
,
-
[
10
(
7
a
-
t
)
6
+
210
c
(
7
a
-
t
)
5
+
1785
c
2
(
7
a
-
t
)
4
+
7840
c
3
(
7
a
-
t
)
3
4
a
+
3
b
≤
t
≤
3
a
+
4
b
+
18795
c
4
(
7
a
-
t
)
2
+
23520
c
5
(
7
a
-
t
)
+
12089
c
6
]
×
(
360
c
7
)
-
1
,
[
15
(
7
a
-
t
)
6
+
420
c
(
7
a
-
t
)
5
+
4830
c
2
(
7
a
-
t
)
4
+
29120
c
3
(
7
a
-
t
)
3
3
a
+
4
b
≤
t
≤
2
a
+
5
b
+
96810
c
4
(
7
a
-
t
)
2
+
168000
c
5
(
7
a
-
t
)
+
119182
c
6
]
×
(
720
c
7
)
-
1
,
-
[
6
(
7
a
-
t
)
6
+
210
c
(
7
a
-
t
)
5
+
3045
c
2
(
7
a
-
t
)
4
+
23380
c
3
(
7
a
-
t
)
3
2
a
+
5
b
≤
t
≤
a
+
6
b
+
100065
c
4
(
7
a
-
t
)
2
+
225750
c
5
(
7
a
-
t
)
+
208943
c
6
]
×
(
720
c
7
)
-
1
,
(
7
b
-
t
)
6
(
720
c
7
)
,
a
+
6
b
≤
t
≤
7
b
.


## 5.2. The Renewal-Intensity Function Approximation When TTF Is Uniform

The Uniform RNIF is given byρ
(
t
)
=
[
e
(
t
-
a
)
/
c
]
/
c only for the interval 0
≤
a
≤
t
<
b. Bartholomew [35] describes ρ
(
t
)
×
Δ
t as 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 at t. However, as described by nearly every author in reliability engineering, the HZF h
(
t
) gives the conditional hazard-rate at time t only amongst survivors of age t; that is, h
(
t
)
×
Δ
t
=
Pr
(
t
≤
T
≤
t
+
Δ
t
)
/
R
(
t
). The hazard function for the U
(
a
,
b
) baseline distribution is given by h
(
t
)
=
1
/
(
b
-
t
), a
≤
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 geometric series for h
(
t
)
=
(
1
/
b
)
/
(
1
-
t
/
b
) and the Maclaurin series for e
t
/
b; that ρ
(
t
)
<
h
(
t
) for all a
<
t
≤
b.In order to compute the RNIFρ
(
t
) for t
>
b, we used two different approximate procedures. First, by directly differentiating (1a) as follows: ρ
(
t
)
=
(
d
/
d
t
)
∑
n
=
1
∞
F
(
n
)
(
t
)
=
∑
n
=
1
∞
f
(
n
)
(
t
), where the exact f
(
1
)
(
t
)
=
f
(
t
), and uniform convolutions f
(
n
)
(
t
), for n
=
2,3
,
…
,
12 can be calculated by our Matlab program. For n
>
12, the program uses the ordinate of N
(
n
μ
,
n
σ
2
) approximation, where μ
=
(
a
+
b
)
/
2 and σ
2
=
(
b
-
a
)
2
/
12. 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 exactn
=
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 for 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 T
12
=
∑
i
=
1
12
X
i, each X
i
~
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. The summary in Table 2 shows the normal approximation to F
(
12
)
(
t
) for intervals of length 0.50
-
σ. Table 2 clearly shows that the worst relative error occurs at σ
/
2 and that the normal approximation improves as Z moves toward the right tail. The accuracy is within 2 decimals up to one-σ and 3 decimals beyond 1.49
σ. Therefore, we conclude that the normal approximation to each of Uniform 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 2
Normal approximation to the 12-fold convolution ofU
(
0,1
).
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.3
E
-
05It should also be noted that the normal approximation to the UniformF
(
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 f
T
n
(
t
)
=
f
(
n
)
(
t
) with those of the corresponding Normal. It can be proven (the proof is available on request) that the skewness of the partial sum T
n
=
∑
i
=
1
n
X
i, X
i’s being iid like X, is given by(12a)
β
3
(
T
n
)
=
μ
3
(
T
n
)
[
μ
2
(
T
n
)
]
3
/
2
=
n
μ
3
(
X
)
(
n
σ
2
)
3
/
2
=
μ
3
(
X
)
(
σ
3
n
)
=
β
3
(
X
)
n
.
Further, the kurtosis of T
n is given by
(12b)
β
4
(
T
n
)
=
α
4
(
X
)
n
+
3
(
n
-
1
)
n
-
3
⟶
β
4
(
T
n
)
=
[
α
4
(
X
)
-
3
]
n
=
β
4
(
X
)
n
,where μ
i
(
i
=
2,3
,
4
) are the universal notation for the i
th central moments. Equation (12b) clearly shows that the kurtosis of the uniform f
(
n
)
(
t
), for n
=
13,14
,
and
15, respectively, is β
4
=
-
1.20
/
13
=
-
0.09231
,
-
1.20
/
14
=
-
0.08571
,
and
-
1.20
/
15
=
-
0.08000, with 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 f
(
n
)
(
t
) to be within 0.01 of Laplace-Gaussian N
(
n
(
a
+
b
)
/
2
,
n
(
b
-
a
)
2
/
12
). 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 f
(
n
)
(
t
) 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β
=
1) does not have a closed-form expression for its n-fold cdf convolution F
(
n
)
(
t
), and hence (1a) cannot directly be used to obtain 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 M
-
(
s
) but it cannot be inverted to yield a closed-form expression for its M
(
t
).It must be highlighted that there have been several articles on approximating the RNF such as Deligonul [49], and also on numerical solutions of M
(
t
) by Tortorella [50], From [24], and other notables. The online supplement by Tortorella provides extensive references on renewal theory and applications. Note that Murthy et al. [42] provide an extensive treatise on Weibull Models, referring to the Weibull with zero minimum life as the standard model. Murthy et al. [42] 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 [1], and others such as Xie [20], with the aid of Matlab to obtain another approximation for M
(
t
).


### 6.2. Discretizing Time in order to Approximate the Renewal Equation

BecauseM
(
t
)
=
F
(
t
)
+
∫
0
t
M
(
t
-
τ
)
f
(
τ
)
d
τ and the underling failure distributions are herein specified, the first term on the right-hand side (RHS) of M
(
t
), F
(
t
), can be easily computed. However, the convolution integral on the RHS, ∫
0
t
M
(
t
-
τ
)
f
(
τ
)
d
τ, except for rare cases, cannot in general be computed and has to be approximated. The discretization method was first applied by Xie [20], where he called his procedure “THE RS-METHOD,” RS for Riemann-Stieltjes. However, Xie [20] 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(
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 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 Δ
t
=
0.50 weeks. As a result, ∫
0
5
M
(
5
-
τ
)
f
(
τ
)
d
τ
≡
∑
i
=
1
10
∫
(
i
-
1
)
/
2
i
/
2
M
(
5
-
τ
)
f
(
τ
)
d
τ, 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 theMean-Value Theorem for Integrals, which states the following: if a function f
(
x
) is continuous over the real closed interval [
a
,
b
], then for certain, there exists a real number x
0 such that ∫
a
b
f
(
x
)
d
x
≡
f
(
x
0
)
×
(
b
-
a
), a
≤
x
0
≤
b, f
(
x
0
) being the ordinate of the integrand at x
0. Because both M
(
t
-
τ
) and the density f
(
t
) are continuous, then, for example, applying the aboveMean-Value Theorem for Integrals to the 4th subinterval, there exists for certain a real number τ
4 such that ∫
3
/
2
2
M
(
5
-
τ
)
f
(
τ
)
d
τ
≡
M
(
5
-
τ
4
)
f
(
τ
4
)
(
2
-
3
/
2
), 3
/
2
≤
τ
4
≤
2. As a result, ∫
0
5
M
(
5
-
τ
)
f
(
τ
)
d
τ
≡
∑
i
=
1
10
M
(
5
-
τ
i
)
f
(
τ
i
)
(
1
/
2
), where 0
≤
τ
1
≤
0.50, 0.5
≤
τ
2
≤
1
,
…, 4.5
≤
τ
10
≤
5. Clearly the exact values of M
(
5
-
τ
i
), i
=
1
,
…
,
10 cannot in general be determined, and because in this example f
(
τ
i
)
(
1
/
2
)
=
∫
(
i
-
1
)
/
2
i
/
2
f
(
τ
)
d
τ
=
Pr
[
(
i
-
1
)
/
2
≤
TTF
≤
i
/
2
], it follows that ∫
0
5
M
(
5
-
τ
)
f
(
τ
)
d
τ
≡
∑
i
=
1
10
M
(
5
-
τ
i
)
∫
(
i
-
1
)
/
2
i
/
2
f
(
τ
)
d
τ. As proposed by Elsayed [1], 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 M
(
5
-
0.50
i
) which results in
(13)
M
(
5
)
≡
F
(
5
)
+
∫
0
5
M
(
5
-
τ
)
f
(
τ
)
d
τ
≅
F
(
5
)
+
∑
i
=
1
10
M
(
5
-
0.50
i
)
∫
(
i
-
1
)
/
2
i
/
2
f
(
τ
)
d
τ
=
∑
i
=
1
10
∫
(
i
-
1
)
/
2
i
/
2
f
(
τ
)
d
τ
+
∑
i
=
1
10
M
(
5
-
0.50
i
)
∫
(
i
-
1
)
/
2
i
/
2
f
(
τ
)
d
τ
=
∑
i
=
1
10
{
[
1
+
M
(
5
-
0.50
i
)
]
×
∫
(
i
-
1
)
/
2
i
/
2
f
(
τ
)
d
τ
}
.
The above equation (13) is similar to that of (7.10) of Elsayed [1], where his subintervals are of length Δ
t
=
1. 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 M
(
0
)
≡
0 at i
=
10 to calculate the last term of (13); further, at i
=
9, (13) yields [
1
+
M
(
5
-
4.5
)
]
∫
8
/
2
(
8
+
1
)
/
2
f
(
τ
)
d
τ
=
[
1
+
M
(
0.50
)
]
∫
4
4.5
f
(
τ
)
d
τ. 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
×
F
(
0.50
)
+
0
×
R
(
0.50
). Hence, at i
=
9, the value of the term before last in (13) reduces approximately to [
1
+
F
(
0.50
)
]
×
∫
4
4.5
f
(
τ
)
d
τ. At i
=
8, the value of (13) is given by [
1
+
M
(
1
)
]
∫
3.5
4
f
(
τ
)
d
τ, where M
(
1
)
=
∑
i
=
1
2
{
[
1
+
M
(
1
-
0.50
i
)
]
×
∫
(
i
-
1
)
/
2
i
/
2
f
(
τ
)
d
τ
}, where M
(
0.50
) has been approximated. Continuing in this manner, we solved (13) backward-recursively in order to approximate M
(
t
). The smaller Δ
t always leads to a better approximation of M
(
t
).


### 6.3. Renewal-Intensity Function Approximation for the Weibull Distribution

Next, after approximating the Weibull RNF, how do we use its approximateM
(
t
) to obtain a fairly accurate value of Weibull RNIF ρ
(
t
)? Because ρ
(
t
)
≡
d
M
(
t
)
/
d
t
≡
Lim
Δ
t
→
0
(
M
(
t
+
Δ
t
)
-
M
(
t
)
)
/
Δ
t, then, for sufficiently small Δ
t
>
0 the approximate ρ
(
t
)
≈
(
M
(
t
+
Δ
t
)
-
M
(
t
)
)
/
Δ
t, which uses the right-hand derivative, and ρ
(
t
)
≈
(
M
(
t
)
-
M
(
t
-
Δ
t
)
)
/
Δ
t, 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 ρ
(
t
) by averaging the two, where t and Δ
t are inputted by the user. It is recommended that the user inputs 0
<
Δ
t
<
0.01
t such that the probability of 2 or more failures during an interval of length Δ
t is almost zero. Further, our program shows that ρ
(
t
)
<
h
(
t
) if and only if the shape β
>
1.


### 6.4. Time-Discretization Accuracy

The accuracy of above discretization method was checked in two different ways.First, we verified that atβ
≅
3.439541, the Weibull mean, median, and mode become almost identical at which the Weibull skewness is β
3
=
0.04052595 and Weibull kurtosis is β
4
=
-
0.288751, with these last 2 standardized-moments being very close to those of Laplace-Gaussian of identically equal to zero. Our Matlab program at β
=
3.439541, δ
=
0, θ
=
α
=
2000, and Δ
t
=
50 yielded the Weibull M
(
4000
)
≅
1.753638831 (with CPU-time = 240.4183 sec), while the corresponding Normal (with MTBF = 1797.84459964 and σ
=
577.9338342) resulted in M
(
4000
)
=
1.77439715 (CPU ≅
48.6). Secondly, we ran our program for the exponential (which is the Weibull with minimum-life δ
=
0, θ
=
α
=
1
/
λ and β
=
1), at λ
=
0.001 and t
=
10000 for varying values of Δ
t. Table 3 depicts the comparisons against the exact M
(
t
)
=
λ
t
=
0.001
×
10000
=
10 expected RNLs.Table 3
Time-discretizing approximation method results.
Δ
t
M
(
t
)
ApproximateM
(
t
)
Relative error
Elapsed time (seconds)
10
10.0000000000
9.9501662508319
−0.49834%
40238.61714
25
10.0000000000
9.8760351886669
−1.23965%
2018.707614
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.066944The above time-discretization-method, using theMean-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 Δ
t
=
0.01
t, the computational time exceeds one minute and the corresponding error may not be acceptable; further, it is not a closed-form approximation.


## 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β
=
1) does not have a closed-form expression for its n-fold cdf convolution F
(
n
)
(
t
), and hence (1a) cannot directly be used to obtain 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 M
-
(
s
) but it cannot be inverted to yield a closed-form expression for its M
(
t
).It must be highlighted that there have been several articles on approximating the RNF such as Deligonul [49], and also on numerical solutions of M
(
t
) by Tortorella [50], From [24], and other notables. The online supplement by Tortorella provides extensive references on renewal theory and applications. Note that Murthy et al. [42] provide an extensive treatise on Weibull Models, referring to the Weibull with zero minimum life as the standard model. Murthy et al. [42] 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 [1], and others such as Xie [20], with the aid of Matlab to obtain another approximation for M
(
t
).


## 6.2. Discretizing Time in order to Approximate the Renewal Equation

BecauseM
(
t
)
=
F
(
t
)
+
∫
0
t
M
(
t
-
τ
)
f
(
τ
)
d
τ and the underling failure distributions are herein specified, the first term on the right-hand side (RHS) of M
(
t
), F
(
t
), can be easily computed. However, the convolution integral on the RHS, ∫
0
t
M
(
t
-
τ
)
f
(
τ
)
d
τ, except for rare cases, cannot in general be computed and has to be approximated. The discretization method was first applied by Xie [20], where he called his procedure “THE RS-METHOD,” RS for Riemann-Stieltjes. However, Xie [20] 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(
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 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 Δ
t
=
0.50 weeks. As a result, ∫
0
5
M
(
5
-
τ
)
f
(
τ
)
d
τ
≡
∑
i
=
1
10
∫
(
i
-
1
)
/
2
i
/
2
M
(
5
-
τ
)
f
(
τ
)
d
τ, 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 theMean-Value Theorem for Integrals, which states the following: if a function f
(
x
) is continuous over the real closed interval [
a
,
b
], then for certain, there exists a real number x
0 such that ∫
a
b
f
(
x
)
d
x
≡
f
(
x
0
)
×
(
b
-
a
), a
≤
x
0
≤
b, f
(
x
0
) being the ordinate of the integrand at x
0. Because both M
(
t
-
τ
) and the density f
(
t
) are continuous, then, for example, applying the aboveMean-Value Theorem for Integrals to the 4th subinterval, there exists for certain a real number τ
4 such that ∫
3
/
2
2
M
(
5
-
τ
)
f
(
τ
)
d
τ
≡
M
(
5
-
τ
4
)
f
(
τ
4
)
(
2
-
3
/
2
), 3
/
2
≤
τ
4
≤
2. As a result, ∫
0
5
M
(
5
-
τ
)
f
(
τ
)
d
τ
≡
∑
i
=
1
10
M
(
5
-
τ
i
)
f
(
τ
i
)
(
1
/
2
), where 0
≤
τ
1
≤
0.50, 0.5
≤
τ
2
≤
1
,
…, 4.5
≤
τ
10
≤
5. Clearly the exact values of M
(
5
-
τ
i
), i
=
1
,
…
,
10 cannot in general be determined, and because in this example f
(
τ
i
)
(
1
/
2
)
=
∫
(
i
-
1
)
/
2
i
/
2
f
(
τ
)
d
τ
=
Pr
[
(
i
-
1
)
/
2
≤
TTF
≤
i
/
2
], it follows that ∫
0
5
M
(
5
-
τ
)
f
(
τ
)
d
τ
≡
∑
i
=
1
10
M
(
5
-
τ
i
)
∫
(
i
-
1
)
/
2
i
/
2
f
(
τ
)
d
τ. As proposed by Elsayed [1], 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 M
(
5
-
0.50
i
) which results in
(13)
M
(
5
)
≡
F
(
5
)
+
∫
0
5
M
(
5
-
τ
)
f
(
τ
)
d
τ
≅
F
(
5
)
+
∑
i
=
1
10
M
(
5
-
0.50
i
)
∫
(
i
-
1
)
/
2
i
/
2
f
(
τ
)
d
τ
=
∑
i
=
1
10
∫
(
i
-
1
)
/
2
i
/
2
f
(
τ
)
d
τ
+
∑
i
=
1
10
M
(
5
-
0.50
i
)
∫
(
i
-
1
)
/
2
i
/
2
f
(
τ
)
d
τ
=
∑
i
=
1
10
{
[
1
+
M
(
5
-
0.50
i
)
]
×
∫
(
i
-
1
)
/
2
i
/
2
f
(
τ
)
d
τ
}
.
The above equation (13) is similar to that of (7.10) of Elsayed [1], where his subintervals are of length Δ
t
=
1. 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 M
(
0
)
≡
0 at i
=
10 to calculate the last term of (13); further, at i
=
9, (13) yields [
1
+
M
(
5
-
4.5
)
]
∫
8
/
2
(
8
+
1
)
/
2
f
(
τ
)
d
τ
=
[
1
+
M
(
0.50
)
]
∫
4
4.5
f
(
τ
)
d
τ. 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
×
F
(
0.50
)
+
0
×
R
(
0.50
). Hence, at i
=
9, the value of the term before last in (13) reduces approximately to [
1
+
F
(
0.50
)
]
×
∫
4
4.5
f
(
τ
)
d
τ. At i
=
8, the value of (13) is given by [
1
+
M
(
1
)
]
∫
3.5
4
f
(
τ
)
d
τ, where M
(
1
)
=
∑
i
=
1
2
{
[
1
+
M
(
1
-
0.50
i
)
]
×
∫
(
i
-
1
)
/
2
i
/
2
f
(
τ
)
d
τ
}, where M
(
0.50
) has been approximated. Continuing in this manner, we solved (13) backward-recursively in order to approximate M
(
t
). The smaller Δ
t always leads to a better approximation of M
(
t
).


## 6.3. Renewal-Intensity Function Approximation for the Weibull Distribution

Next, after approximating the Weibull RNF, how do we use its approximateM
(
t
) to obtain a fairly accurate value of Weibull RNIF ρ
(
t
)? Because ρ
(
t
)
≡
d
M
(
t
)
/
d
t
≡
Lim
Δ
t
→
0
(
M
(
t
+
Δ
t
)
-
M
(
t
)
)
/
Δ
t, then, for sufficiently small Δ
t
>
0 the approximate ρ
(
t
)
≈
(
M
(
t
+
Δ
t
)
-
M
(
t
)
)
/
Δ
t, which uses the right-hand derivative, and ρ
(
t
)
≈
(
M
(
t
)
-
M
(
t
-
Δ
t
)
)
/
Δ
t, 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 ρ
(
t
) by averaging the two, where t and Δ
t are inputted by the user. It is recommended that the user inputs 0
<
Δ
t
<
0.01
t such that the probability of 2 or more failures during an interval of length Δ
t is almost zero. Further, our program shows that ρ
(
t
)
<
h
(
t
) if and only if the shape β
>
1.


## 6.4. Time-Discretization Accuracy

The accuracy of above discretization method was checked in two different ways.First, we verified that atβ
≅
3.439541, the Weibull mean, median, and mode become almost identical at which the Weibull skewness is β
3
=
0.04052595 and Weibull kurtosis is β
4
=
-
0.288751, with these last 2 standardized-moments being very close to those of Laplace-Gaussian of identically equal to zero. Our Matlab program at β
=
3.439541, δ
=
0, θ
=
α
=
2000, and Δ
t
=
50 yielded the Weibull M
(
4000
)
≅
1.753638831 (with CPU-time = 240.4183 sec), while the corresponding Normal (with MTBF = 1797.84459964 and σ
=
577.9338342) resulted in M
(
4000
)
=
1.77439715 (CPU ≅
48.6). Secondly, we ran our program for the exponential (which is the Weibull with minimum-life δ
=
0, θ
=
α
=
1
/
λ and β
=
1), at λ
=
0.001 and t
=
10000 for varying values of Δ
t. Table 3 depicts the comparisons against the exact M
(
t
)
=
λ
t
=
0.001
×
10000
=
10 expected RNLs.Table 3
Time-discretizing approximation method results.
Δ
t
M
(
t
)
ApproximateM
(
t
)
Relative error
Elapsed time (seconds)
10
10.0000000000
9.9501662508319
−0.49834%
40238.61714
25
10.0000000000
9.8760351886669
−1.23965%
2018.707614
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.066944The above time-discretization-method, using theMean-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 Δ
t
=
0.01
t, the computational time exceeds one minute and the corresponding error may not be acceptable; further, it is not a closed-form approximation.


## 7. Conclusions

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ρ
(
t
) and h
(
t
) for t
>
0, except in the only case of CFR for which ρ
(
t
)
≡
h
(
t
)
≡
λ. Further, when h
(
t
) is an IFR, then ρ
(
t
)
<
h
(
t
) so that R
(
t
)
<
e
-
M
(
0
,
t
), t
>
0, for the four baseline failure distributions that we have studied. However, this is not quite consistent with R
(
t
)
=
e
-
M
(
0
,
t
) given by some authors in reliability engineering, such as Ebeling [36], for a NHPP. Further, some authors such as Leemis [37] use the same notation λ
(
t
) for the Weibull RNIF and also use λ
(
t
) 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 analternating renewal process[51, page 350].



---

*Source: 857437-2014-03-19.xml*