Template:Grp confidence bounds

Confidence Bounds

In general, in order to obtain the virtual age, the exact occurrence time of each event (failure) should be available (see equations for Type I and Type 2 models). However, the times are unknown until the corresponding events occur. For this reason, there are no closed-form expressions for total failure number and failure intensity, which are functions of failure times and virtual age. Therefore, in Weibull++, a Monte Carlo simulation is used to predict values of virtual time, failure number, MTBF and failure rate. The approximate confidence bounds obtained from simulation are provided. The uncertainty of model parameters is also considered in the bounds.

Bounds on Cumulative Failure (Event) Numbers

The variance of the cumulative failure number N(t) is:

[math]\displaystyle{ Var[N(t)]=Var\left[ E(N(t)|\lambda ,\beta ,q) \right]+E\left[ Var(N(t)|\lambda ,\beta ,q) \right] }[/math]

The first term accounts for the uncertainty of the parameter estimation. The second term considers the uncertainty caused by the renewal process even when model parameters are fixed. However, unless q = 1 , [math]\displaystyle{ Var\left[ E(N(t)|\lambda ,\beta ,q) \right] }[/math] cannot be calculated because E(N(t)) cannot be expressed as a closed-form function of λ,β, and q . In order to consider the uncertainty of the parameter estimation, [math]\displaystyle{ Var\left[ E(N(t)|\lambda ,\beta ,q) \right] }[/math] is approximated by:

[math]\displaystyle{ Var\left[ E(N(t)|\lambda ,\beta ,q) \right]=Var[E(N({{v}_{t}})|\lambda ,\beta )]=Var[\lambda v_{t}^{\beta }] }[/math]

where v_t is the expected virtual age at time t and [math]\displaystyle{ Var[\lambda v_{t}^{\beta }] }[/math] is::

[math]\displaystyle{ \begin{align} & Var[\lambda v_{t}^{\beta }]= & {{\left( \frac{\partial (\lambda v_{t}^{\beta })}{\partial \beta } \right)}^{2}}Var(\hat{\beta })+{{\left( \frac{\partial (\lambda v_{t}^{\beta })}{\partial \lambda } \right)}^{2}}Var(\hat{\lambda }) \\ & +2\frac{\partial (\lambda v_{t}^{\beta })}{\partial \beta }\frac{\partial (\lambda v_{t}^{\beta })}{\partial \lambda }Cov(\hat{\beta },\hat{\lambda }) \end{align} }[/math]

By conducting this approximation, the uncertainty of λ and β are considered. The value of v_t and the value of the second term in the equation for the variance of number of failures are obtained through the Monte Carlo simulation using parameters [math]\displaystyle{ \hat{\lambda },\hat{\beta },\hat{q}, }[/math] which are the ML estimators. The same simulation is used to estimate the cumulative number of failures [math]\displaystyle{ \hat{N}(t)=E(N(t)|\hat{\lambda },\hat{\beta },\hat{q}) }[/math] .

Once the variance and the expected value of N(t) have been obtained, the bounds can be calculated by assuming that N(t) is lognormally distributed as:

[math]\displaystyle{ \frac{\ln N(t)-\ln \hat{N}(t)}{\sqrt{Var(\ln N(t))}}\tilde{\ }N(0,1) }[/math]

The upper and lower bounds for a given confidence level α can be calculated by:

[math]\displaystyle{ N{{(t)}_{U,L}}=\hat{N}(t){{e}^{\pm {{z}_{a}}\sqrt{Var(N(t))}/\hat{N}(t)}} }[/math]

where z_a is the standard normal distribution.

If N(t) is assumed to be normally distributed, the bounds can be calculated by:

[math]\displaystyle{ N{{(t)}_{U}}=\hat{N}(t)+{{z}_{a}}\sqrt{Var(N(t))} }[/math]

[math]\displaystyle{ N{{(t)}_{L}}=\hat{N}(t)-{{z}_{a}}\sqrt{Var(N(t))} }[/math]

In Weibull++, the N(t)_U is the smaller of the upper bounds obtained from lognormal and normal distribution appoximation. The N(t)_L is set to the largest of the lower bounds obtained from lognormal and normal distribution appoximation. This combined method can prevent the out-of-range values of bounds for some small t values.

Bounds of Cumulative Failure Intensity and MTBF

For a given time t , the expected value of cumulative MTBF m_c(t) and cumulative failure intensity λ_c(t) can be calculated using the following equations:

[math]\displaystyle{ {{\hat{\lambda }}_{c}}(t)=\frac{\hat{N}(t)}{t};{{\hat{m}}_{c}}(t)=\frac{t}{\hat{N}(t)} }[/math]

The bounds can be easily obtained from the corresponding bounds of N(t).

[math]\displaystyle{ \begin{align} & {{{\hat{\lambda }}}_{c}}{{(t)}_{L}}= & \frac{\hat{N}{{(t)}_{L}}}{t};\text{ }{{{\hat{\lambda }}}_{c}}{{(t)}_{L}}=\frac{\hat{N}{{(t)}_{L}}}{t};\text{ } \\ & {{{\hat{m}}}_{c}}{{(t)}_{L}}= & \frac{t}{\hat{N}{{(t)}_{U}}};\text{ }{{{\hat{m}}}_{c}}{{(t)}_{U}}=\frac{t}{\hat{N}{{(t)}_{L}}} \end{align} }[/math]

Bounds on Instantaneous Failure Intensity and MTBF

The instantaneous failure intensity is given by:

[math]\displaystyle{ {{\lambda }_{i}}(t)=\lambda \beta v_{t}^{\beta -1} }[/math]

where v_t is the virtual age at time t. When [math]\displaystyle{ q\ne 1, }[/math] it is obtained from simulation. When q = 1 , v_t = t from model Type I and Type 2.

The variance of instantaneous failure intensity can be calculated by:

[math]\displaystyle{ \begin{align} & Var({{\lambda }_{i}}(t))= {{\left( \frac{\partial {{\lambda }_{i}}(t)}{\partial \beta } \right)}^{2}}Var(\hat{\beta })+{{\left( \frac{\partial {{\lambda }_{i}}(t)}{\partial \lambda } \right)}^{2}}Var(\hat{\lambda }) \\ & +2\frac{\partial {{\lambda }_{i}}(t)}{\partial \beta }\frac{\partial {{\lambda }_{i}}(t)}{\partial \lambda }Cov(\hat{\beta },\hat{\lambda })+{{\left( \frac{\partial {{\lambda }_{i}}(t)}{\partial v(t)} \right)}^{2}}Var({{{\hat{v}}}_{t}}) \end{align} }[/math]

The expected value and variance of v_t are obtained from the Monte Carlo simulation with parameters [math]\displaystyle{ \hat{\lambda },\hat{\beta },\hat{q}. }[/math] Because of the simulation accuracy and the convergence problem in calculation of [math]\displaystyle{ Var(\hat{\beta }),Var(\hat{\lambda }) }[/math] and [math]\displaystyle{ Cov(\hat{\beta },\hat{\lambda }), }[/math] V'a'r(λ_i(t)) can be a negative value at some time points. When this case happens, the bounds of instantaneous failure intensity are not provided.

Once the variance and the expected value of λ_i(t) are obtained, the bounds can be calculated by assuming that  λ_i(t) is lognormally distributed as:

[math]\displaystyle{ \frac{\ln {{\lambda }_{i}}(t)-\ln {{{\hat{\lambda }}}_{i}}(t)}{\sqrt{Var(\ln {{\lambda }_{i}}(t))}}\tilde{\ }N(0,1) }[/math]

The upper and lower bounds for a given confidence level α can be calculated by:

[math]\displaystyle{ {{\lambda }_{i}}(t)={{\hat{\lambda }}_{i}}(t){{e}^{\pm {{z}_{a}}\sqrt{Var({{\lambda }_{i}}(t))}/{{{\hat{\lambda }}}_{i}}(t)}} }[/math]

where z_a is the standard normal distribution.

If λ_i(t) is assumed to be normally distributed, the bounds can be calculated by:

[math]\displaystyle{ {{\lambda }_{i}}{{(t)}_{U}}={{\hat{\lambda }}_{i}}(t)+{{z}_{a}}\sqrt{Var(N(t))} }[/math]

[math]\displaystyle{ {{\lambda }_{i}}{{(t)}_{L}}={{\hat{\lambda }}_{i}}(t)-{{z}_{a}}\sqrt{Var(N(t))} }[/math]

In Weibull++, λ_i(t)_U is set to the smaller of the two upper bounds obtained from the above lognormal and normal distribution appoximation. λ_i(t)_L is set to the largest of the two lower bounds obtained from the above lognormal and normal distribution appoximation. This combination method can prevent the out of range values of bounds when t values are small.

For a given time t , the expected value of cumulative MTBF m_i(t) is:

[math]\displaystyle{ {{\hat{m}}_{i}}(t)=\frac{1}{{{{\hat{\lambda }}}_{i}}(t)}\text{ } }[/math]

The upper and lower bounds can be easily obtained from the corresponding bounds of λ_i(t):

[math]\displaystyle{ {{\hat{m}}_{i}}{{(t)}_{U}}=\frac{1}{{{{\hat{\lambda }}}_{i}}{{(t)}_{L}}} }[/math]

[math]\displaystyle{ {{\hat{m}}_{i}}{{(t)}_{L}}=\frac{1}{{{{\hat{\lambda }}}_{i}}{{(t)}_{U}}} }[/math]

Bounds on Conditional Reliability

Given mission start time t₀ and mission time T , the conditional reliability can be calculated by:

[math]\displaystyle{ R(T|{{t}_{0}})=\frac{R(T+{{v}_{0}})}{R({{v}_{0}})}={{e}^{-\lambda [{{({{v}_{0}}+T)}^{\beta }}-{{v}_{0}}]}} }[/math]

v₀ is the virtual age corresponding to time t₀ . The expected value and the variance of v₀ are obtained from Monte Carlo simulation. The variance of the conditional reliability R(T | t₀) is:

[math]\displaystyle{ \begin{align} & Var(R)= {{\left( \frac{\partial R}{\partial \beta } \right)}^{2}}Var(\hat{\beta })+{{\left( \frac{\partial R}{\partial \lambda } \right)}^{2}}Var(\hat{\lambda }) \\ & +2\frac{\partial R}{\partial \beta }\frac{\partial R}{\partial \lambda }Cov(\hat{\beta },\hat{\lambda })+{{\left( \frac{\partial R}{\partial {{v}_{0}}} \right)}^{2}}Var({{{\hat{v}}}_{0}}) \end{align} }[/math]

Because of the simulation accuracy and the convergence problem in calculation of [math]\displaystyle{ Var(\hat{\beta }),Var(\hat{\lambda }) }[/math] and [math]\displaystyle{ Cov(\hat{\beta },\hat{\lambda }), }[/math] V'a'r(R) can be a negative value at some time points. When this case happens, the bounds are not provided.

The bounds are based on:

[math]\displaystyle{ \log \text{it}(\hat{R}(T))\tilde{\ }N(0,1) }[/math]

[math]\displaystyle{ \log \text{it}(\hat{R}(T))=\ln \left\{ \frac{\hat{R}(T)}{1-\hat{R}(T)} \right\} }[/math]

The confidence bounds on reliability are given by:

[math]\displaystyle{ R=\frac{{\hat{R}}}{\hat{R}+(1-\hat{R}){{e}^{\pm \sqrt{Var(R)}/[\hat{R}(1-\hat{R})]}}} }[/math]

It will be compared with the bounds obtained from:

[math]\displaystyle{ R=\hat{R}{{e}^{\pm {{z}_{a}}\sqrt{Var(R)}/\hat{R}}} }[/math]

The smaller of the two upper bounds will be the final upper bound and the larger of the two lower bounds will be the final lower bound.

Example 4: The following table gives the failure times for the air conditioning unit of an aircraft. The observation ended by the time the last failure occurred, as discussed in Cox [3].

[math]\displaystyle{ \begin{matrix} \text{50} & \text{329} & \text{811} & \text{991} & \text{1489} \\ \text{94} & \text{332} & \text{899} & \text{1013} & \text{1512} \\ \text{196} & \text{347} & \text{945} & \text{1152} & \text{1525} \\ \text{268} & \text{544} & \text{950} & \text{1362} & \text{1539} \\ \text{290} & \text{732} & \text{955} & \text{1459} & {} \\ \end{matrix}\,\! }[/math]

1. Estimate the GRP model parameters using the Type I virtual age option.

2. Plot the failure number and instantaneous failure intensity vs. time with 90% two-sided confidence bounds.

3. Plot the conditional reliability vs. time with 90% two-sided confidence bounds. The mission start time is 40 and mission time is varying.

4. Using the QCP, calculate the expected failure number and expected instantaneous failure intensity by time 1800.

Solution

Enter the data into a parametric RDA folio in Weibull++. On the control panel, select the 3 parameters option and the Type I setting. Keep the default simulation settings. Click Calculate.

1. The estimated parameters are [math]\displaystyle{ \hat{\beta }=1.1976\,\! }[/math], [math]\displaystyle{ \hat{\lambda }=4.94E-03\,\! }[/math], [math]\displaystyle{ \hat{q}=0.1344\,\! }[/math].

2. The following plots show the cumulative number of failures and instantaneous failure intensity, respectively.

3. The following plot shows the conditional reliability.

4. Using the QCP, the failure number and instantaneous failure intensity are: