Imagine product or part lifetime as the service time in a self-service system. Suppose you could observe periodic input and output counts of that self-service system, without identifying individual service times. How would you estimate the service time distribution without following individual products or parts from input to output? The maximum likelihood reliability estimator for an M/G/Infinity self-service-time distribution function from “ships” (inputs) and “returns” (outputs) counts also works for nonstationary arrival process M(t)/G/Infinity self-service systems, under a condition! A constant or increasing arrival (inputs or ships) rate satisfies the condition. If you identify returns by failure mode then you could estimate make nonparametric estimates of field reliability by failure mode and quantify reliability growth, without life data, if the condition is satisfied.
Failure?
Recently a friend sent me his annual sales and returns.
Table 1. Nevada table of ships (sales) and returns (sums in bottom row)
| Year | Sales | 2020 | 2021 | 2022 | 2023 | 2024 | 2025 |
| 2020 | 393 | ||||||
| 2021 | 502 | ||||||
| 2022 | 592 | ||||||
| 2023 | 141 | ||||||
| 2024 | 0 | ||||||
| 2025 | 0 | ||||||
| Returns | 14 | 9 | 11 | 12 | 3 | 5 |
I computed the nonparametric reliability estimate using my 30+ year old spreadsheet, https://docs.google.com/spreadsheets/d/1fb_pjwygp_2CUeayNEV4PAX6wyR8sdWqPr76qeJ_n4I/edit?gid=0#gid=0/.
Table 2. Nonparametric maximum likelihood reliability estimators from successively earlier cohorts
| Age, years | All cohorts | -2025 | -‘24&’25 | Etc. | Etc. | 2020 |
| 1 | 0.9771 | 0.9771 | 0.9771 | 0.9771 | 0.9743 | 0.9644 |
| 2 | 0.9771 | 0.9771 | 0.9771 | 0.9771 | 0.9743 | |
| 3 | 0.9771 | 0.9771 | 0.9771 | 0.9771 | ||
| 4 | 0.9149 | 0.9149 | 0.9149 | |||
| 5 | 0.8936 | 0.8936 | ||||
| 6 | 0.8582 |
Because the sales were zero in 2024 and 2025, I became suspicious that the PAV (Penalize Adjacent Values) algorithm for nonparametric maximum likelihood estimator (npmle) of a nondecreasing cumulative distribution function G(s) might fail [Marshall and Proschan]. So I also computed the least-squares reliability estimator (nplse). (Thanks to Dick Mensing for suggesting this in 1991, https://accendoreliability.com/uncertainty-in-population-estimates/#more-503100/.) Notice how the reliability npmle goes nuts in ages 4, 5, and 6? That’s because a condition for the PAV algorithm is violated.
The npmle is valid if E[λ(t-S)]=∫λ(t-s)dG(s), integrated over s=1,2,…,t, is non-decreasing [George 1994]. Table 3 uses the nplse probability density function estimate, not that from the npmle, because the PAV algorithm produces non-decreasing E[λ(t-S)], whether the reliability estimate is valid or not.
Table 3. Condition for validity of PAV npmle algorithm: entries λ(t-s) are the ships, E[λ(t-S)] is approximated by ∑p(s)*λ(t-s), s=1,2,..,t where p(s)=G(s+1)-G(s).
| t | λ(t-s) | Etc. | Etc. | Etc. | Etc. | Etc. | p(s) | E[λ(t-S)] | NonDecreasing? |
| 1 | 393 | 0.02 | 8.68 | Eh? | |||||
| 2 | 502 | 393 | 0.00 | 9.03 | TRUE | ||||
| 3 | 592 | 502 | 393 | 0 | 13.5 | TRUE | |||
| 4 | 141 | 592 | 502 | 393 | 0.01 | 7.8 | FALSE | ||
| 5 | 0 | 141 | 592 | 502 | 393 | 0 | 5.5 | FALSE | |
| 6 | 0 | 0 | 141 | 592 | 502 | 393 | 0 | 6.33 | TRUE |
Reliability Estimation Without Life Data?
[The following article dates from 1994, unpublished.]
Except for high-cost or repairable products or parts, there is little incentive to record when returns or failures were shipped, so there is no time-to-return data from which to estimate reliability. Since 21 CFR Part 821, “Medical Device Tracking Requirements”, pacemakers have been tracked by serial number from implant to failure. Medtronic used to publish pacemaker sales and failure counts and their reliability estimates from their lifetime data. I made nonparametric reliability estimates from their pacemaker sales and failure counts. Our estimates agreed, and I sent them to Medtronic. Medtronic quit publishing sales, failure counts, and reliability estimates from their pacemaker lifetime data!
Product lives may be regarded as service times of an M/G/Infinity self-service system. (“M” means inputs are a stationary Poisson process with service time distribution G(t).) System inputs are “ships”: product sales, installations, or shipments, by accounting interval. System outputs are “returns”: complaints, failures, repairs, replacements, spares ships or returns.
Larry Baxter derived a reliability function estimator from ships and renewals data, when all the ships occurred at the start of reporting. Unfortunately ships occur every month, and returns come back without a clue about how old they are.
Ships and returns counts provide statistically sufficient data to estimate products’ and parts’ reliability functions. A reliability function is one minus the service-time distribution function G(t) of an M/G/Infinity self-service system, because ships “serve” for a random time [George and Agrawal].
The output from an M/G/Infinity service system is a Poisson process with mean rate λG(t) where λ is the Poisson mean input rate [Mirasol]. The Poisson output mean rate λG(t) is analogous to the failure rate function of a time-to-failure random variable. Ulf Grenander and later Albert Marshall and Frank Proschan derived the Pool Adjacent Violators (PAV) nonparametric maximum likelihood estimator (npmle) of an nondecreasing failure rate function. The npmle of λG(t) gives the npmle of G(t) by dividing λG(t) by the input rate λ, because, under reasonable conditions, a function of an npmle is also an npmle. Reality may violate conditions to represent returns as output from a M/G/Infinity system.
What if Ships Process is Not Stationary?
Nonparametric maximum likelihood reliability function estimates, from ships and returns counts population data, derivable from data required by generally accepted accounting principles, are desirable for their properties and require no unwarranted distribution assumptions. They could be valid even if the input Poisson process rate λ(t) is nonstationary.
Theorem 1 of the article by Eick, Massey and Whitt says, “The departure process (of a nonstationary M(t)/G/Infinity service system) is a Poisson process with time dependent output rate δ(t)=E[λ(t-S)].” S is the service-time random variable, the expected value E[λ(t-S)] is ∫λ(t-s)dG(s), integrated from s=0 to t for periods t=1,2,…, and λ(t-S) is the ships input rate at time t-S. G(s) is the cumulative distribution function of service time S. If E[λ(t-S)] is nondecreasing, then the PAV nonparametric maximum likelihood estimator may be used to estimate reliability 1-G(s).
For example, suppose the Poisson ships rate is a constant λ. Then, E[λ(t-S)]=∫λdG(s)=λG(t)) and is nondecreasing because G(t) is a distribution function, nondecreasing in age t. For another example, suppose the Poisson rate λ(t) is a linearly increasing function, λo+λt, for t ≥ 0. Then, E[λ(t-S)] = ∫(λo+λ(t-s))dG(s) = λ(λo+λG(t)) and is also nondecreasing.
The npmle of E[λ(t-S)] puts mass at a subset of M(t)/G/Infinity departure times (return times in the reliability context). Those masses may be solved for the discrete npmle of the probability density function g(t). The npmle of E[λ(t-S)] with linearly increasing λ(t) is the npmle of λo +λ[G(t) – G(0) – ∑ug(u)], where the sum from u>0 to t runs over the set of return times, u, at which the npmle of E[λ(t-S)] puts mass. This discrete function can be solved for G(t) at jumps. The npmle of G(t) is constant elsewhere.
Suppose the first return occurs at time t1>0 and the npmle puts mass g(t1) at t1. Then the npmle of G(0) is 0 and the npmle of g(t1) is (g(t1)-λo)/(λ(λ-t1)) if nonnegative. Remaining values of g(u) may be solved recursively for values of u where the npmle puts mass.
An Example of Grouped Data
Table 4 shows the numbers of ships and returns by quarter for a product. Figure 2 shows the npmle and nonparametric least squares estimates (nplse) [George 2008] of the cumulative distribution function G(t) of ages-at-returns. There is not enough data to test whether the ships process is nonstationary Poisson. The condition for using the npmle with a nonstationary Poisson input is satisfied for all observed times. Table 5 lists increasing values of E[λ(t-S)].
Nonparametric least squares estimators (nplse) are computed under two assumptions: returns are from dead forever products or products that undergo renewal. Both estimators agree on G(t) after the third quarter (figure 2). The nplse puts all failure probability in the first quarter, and the npmle divided it between the third and fourth quarter.
Suppose the input rate λ(t) doesn’t keep up with the output rate, so ships < returns in a period. This violates the conditions for the maximum likelihood estimator [George 1993]. The returns may be a renewal or other kind of replacement process and need a different model of ships and returns such as a renewal process shown in figure 2 [Tortorella, George 2008].
Table 4. Ships and returns data. E[λ(t-S)] should be nondecreasing for a valid npmle.
| Quarter | Ships | Returns | E[λ(t-S)] |
| 1 | 892 | 0 | 0 |
| 2 | 64537 | 52 | 0 |
| 3 | 352133 | 937 | 0.7187 |
| 4 | 851677 | 2826 | 52.58 |
| 5 | 1452959 | 3913 | 325.6 |
| 6 | 2469820 | 4175 | 914.6 |
| 7 | 3815927 | 7975 | 1723 |
| 8 | 5288164 | 9326 | 2932 |
| 9 | 6837147 | 6962 | 4676 |
| 10 | 8268514 | 9436 | 6735 |
| 11 | 6001803 | 5838 | 8938 |
What Happens in Steady State?
Information from the M/G/Infinity output dwindles as it approaches steady state. Eventually G(t)->1 and output rate λG(t)-> λ, so steady-state outputs have no information about the service time distribution. Furthermore, a nonstationary Poisson mean input rate function λ(t) could invalidate the condition for the PAV (Pool Adjacent Violators) algorithm for the maximum likelihood estimator of λ(t)G(t). Ships and returns counts can be derived from sales, warranty, repair cost, or spares sales data required by generally accepted accounting principles. But, I’ve seen data where period ships were negative for unknown reasons: product returns, data entry errors, other?
What if There are Multiple Failure Modes, Product Changes?
[Added after COVID-19.]
If the Poisson outputs from an M(t)/G/Infinity self-service system are categorized by failure mode or COVID-19 cases are categorized by variant or by recovery or death, then it is possible to estimate reliability or survival functions conditional on mode, variant, or outcome. It’s tricky, because observed failures, recoveries, or deaths become mixtures of variants.
If inputs and outputs were also categorized by failure mode or COVID-19 recovery, death, or variant, then the PAV algorithm could be used separately on each category. If not, it is possible to estimate “survival” functions for both COVID-19 time to death and time to recovery from case, recovery, and death counts, because categorized outputs have Poisson distributions conditional on category. Maximum likelihood and least squares both work to estimate conditional reliability or survival functions, including Paretos or Case Fatality Ratios in the optimization (figure 3). In the reliability context, it is possible to estimate reliability functions, from ships and returns counts, the latter categorized by failure mode.
COVID-19 variants seem to have become less lethal, but survival functions estimates for successive variants have to account for deaths and recoveries from earlier variants. The same methods work for quantifying product reliability growth due to product changes, not MTBF growth as in [Duane, AMSAA Crow]. Maximize the log likelihood function or minimize the sum of squared errors between observed failure mode counts and actuarial hindcasts ∑a(s|mode)n(t-s), s=1,2,…,t, of past failure mode counts, simultaneously, as functions of the conditional distribution functions G(t|category) (likelihood) or conditional actuarial rates a(s|mode)=(G(t|mode)-G(t-1|mode))/G(t|mode).
If you would like examples of multivariate survival function or reliability function estimation, without life data, ask or send data.
References
L. A. Baxter, “Estimation from Quasi Life Tables,” Biometrika, Vol. 81, No. 3, pp. 567-577, Aug., 1994
L. H. Crow, “Methods for Assessing Reliability Growth Potential,” IEEE Proceedings Annual Reliability and Maintainability Symposium, pp. 484-489, 1984
J. T. Duane, “Learning Curve Approach To Reliability Monitoring,” IEEE Transactions on Aerospace, Vol. 2, pp. 563-566, 1964
S. G. Eick, W. A. Massey, and W. Whitt, “The Physics of the Mt/G/Infinity Queue,” Ops. Res. , vol. 41, pp. 731-742, 1993
L. L. George, “Estimate Reliability Functions Without Life Data,” ASQ Reliability Review, vol. 13, pp. 21-25, 1993
L. L. George, “A Note on Estimation of a Hidden Service Distribution of an M/G/Infinity Service System,” Random Samples, ASQC Santa Clara Valley, June 1994
L. L. George and A. Agrawal, “Estimation of a Hidden Service Distribution of an M/G/Infinity Service System,” Naval Research Logistics Quarterly, vol. 20, pp. 549-555, 1973
L. L. George, “Renewal Estimation Without Renewal Counts,” ASQ Tech Briefs, Vol. 2, 2008
U. Grenander, “On the Theory of Mortality Measurement,” Skandinavisk Aktuarietidskrift, 39, parts I and II, pp. 70-96 and 125-153, 1956
Albert W. Marshall and Frank Proschan, “Maximum Likelihood Estimation for Distributions with Monotone Failure Rate,” Annals of Math. Statist., vol. 36, pp. 69-77, 1965
Noel Mirasol, “The Output of an M/G/Infinity Queuing System is Poisson,” Operations Research, Vol. 11, No. 2, pp. 282-284, Mar.-Apr., 1963
Mike Tortorella, “Life Estimation from Pooled Discrete Renewal Counts,” published in: Jewell, N.P., et al., (eds) Lifetime Data: Models in Reliability and Survival Analysis, Springer, Boston, MA, 1995
Leave a Reply