
Unbiasing Beta for the CrowAMSAA (NHPP) Model
The CrowAMSAA (NHPP) model is probably the most popular model in reliability growth and repairable system analysis. With this model, maximum likelihood estimation (MLE) is used to estimate the parameters. However, it is well known that the MLE solution for the shape parameter in any model is biased (i.e., the estimated beta value is too high), which affects the accuracy of the results. This is a particular concern when working with small sample sizes. This article will demonstrate how the bias is corrected for CrowAMSAA in RGA. This is also applicable to the Crow Extended, Crow Extended  Continuous Evaluation and power law models.
CrowAMSAA (NHPP) Model
The instantaneous failure intensity of a system is a power function in the CrowAMSAA model. It is:
(1) 
where λ and β are the two model parameters. Sometimes this model is also called the power law NHPP model. The function for the expected number of failures is:
(2) 
The CrowAMSAA model can be used for both time and failure terminated data. For example, if a test ends at 300 hours even though the last failure was at 286.1 hours, then it is a time terminated test. If the test ended at the time of the last failure (i.e., 286.1 hours), then it is a failure terminated test.
The data in the next table are from a time terminated test. The F events are failures, and the End event is the end of the test.
Table 1: Data from a time terminated test
Event 
Time to Event

Event 
Time to Event

F 
2.60 
F 
98.10 
F 
16.50 
F 
101.10 
F 
16.50 
F 
132.00 
F 
17.00 
F 
142.20 
F 
21.40 
F 
147.70 
F 
29.10 
F 
149.00 
F 
33.30 
F 
167.20 
F 
56.50 
F 
190.70 
F 
63.10 
F 
193.00 
F 
70.60 
F 
198.70 
F 
73.00 
F 
251.90 
F 
77.70 
F 
282.50 
F 
93.90 
F 
286.10 
F 
95.50 
End 
300.00 
The following data are from a failure terminated test. In this case, the last failure time is the treated as the end of the test.
Table 2: Data from a failure terminated test
Event 
Time to Event 
Event 
Time to Event

F 
0.7 
F 
785.9 
F 
3.7 
F 
887 
F 
13.2 
F 
1010.7 
F 
17.6 
F 
1029.1 
F 
54.5 
F 
1034.4 
F 
99.2 
F 
1136.1 
F 
112.2 
F 
1178.9 
F 
120.9 
F 
1259.7 
F 
151 
F 
1297.9 
F 
163 
F 
1419.7 
F 
174.5 
F 
1571.7 
F 
191.6 
F 
1629.8 
F 
282.8 
F 
1702.3 
F 
355.2 
F 
1928.9 
F 
486.3 
F 
2072.3 
F 
490.5 
F 
2525.2 
F 
513.3 
F 
2928.5 
F 
558.4 
F 
3016.4 
F 
678.1 
F 
3181 
F 
688 
F 
3256.3 
In order to use MLE, we need to get the likelihood function for a set of data. The probability density function (pdf) of the ith failure given that the (i1)th failure occurred at t_{i1} is:
(3) 
For failure terminated data, the likelihood function is the product of the pdf of each failure. It is:
(4) 
where n is the total number of failures.
For time terminated data, the system is working from the last failure time to the test end time T. Therefore, the following conditional reliability equation is used to model the likelihood that no failure occurs by time T.
(5) 
Multiplying Eqn. (5) with the pdf of each failure, we get the likelihood function:
(6) 
Let T* be either T (time terminated) or t_{n }(failure terminated). Eqns. (4) and (6) can be written as:
(7) 
Taking the natural log on both sides:
(8) 
Differentiate Eqn. (8) with respect to:
(9) 
Set the above equation equal to zero and solve for β to get the MLE solution:
(10) 
Bias Correction
ln(T*/t_{i})^{β}, i = 1,…, n are ordered statistics from an exponential distribution with a parameter of 1. This can be explained from the simulation point of view.
The reliability function for an exponential distribution with a parameter of 1 is P_{i} = exp(x_{i}). For time terminated data, it can be seen that P_{i} is a number between 0 and 1. So x_{i} (i = 1,…,n ) is a random number from an exponential distribution with a parameter of 1.
For failure terminated data, x_{n} can be removed and i ranges from 1 to n  1. This is because:
The sum of k random exponential variables with a parameter of 1 is a Gamma random variable with parameters of k and 1. Therefore, for failure terminated data:
(11) 
One variable is eliminated because ln(t_{n}/t_{n}) = 0, so there are only n  1 random variables in the sum in Eqn. (11). Since Y follows a gamma distribution Gamma(n  1, 1), 1/Y is an inverse gamma distribution, 1/Y ~InverseGamma(n  1, 1). The expected value of 1/Y is:
(12) 
From Eqn. (10), we know:
Therefore:
(13) 
Eqn. (13) shows is biased because its expected value is not β. Therefore, the factor (n  2)/n is used to remove the bias.
For time terminated data, since:
(14) 
We can also get:
(15) 
Example:
The Calculate Unbiased Beta option on the Calculations page in the Application Setup window is used to select whether the software will correct the bias, when applicable.
If this option is cleared when you calculate the time terminated data in Table 1, the (biased) MLE beta value will be calculated, as shown next in the Results area.
For an MLE beta of 0.7163 and with 27 data points (n = 27), we can calculate the unbiased beta as follows:
If you select the Calculate Unbiased Beta option in the Application Setup window and recalculate the folio, the unbiased value will appear, as shown next and indicated with the "Beta (UnB)" label:
Conclusion
In this article, we explained why the MLE beta of the
CrowAMSAA model is biased and how to correct it. When the sample size
is large, there may be only a small difference between the biased and
unbiased beta,
as illustrated in the example. When the sample size is small (e.g., less
than 5), the difference is relatively large, so the unbiased beta is
preferred.