Examples

Three examples illustrate advantages of the suggested regression method for recurrent event prediction.

Example 1.
The input data for this example is generated by Monte Carlo method for Kijima model 1 with restoration factor q=1 and Weibull underlying function with shape parameter beta=2.0 and scale parameter eta=1.0. There is the exact solution in this case for CIF function described as x2. In total we had 20 components. The first 5 components were suspended at time=1.5, the second set of 5 components were suspended at time=2 and the rest of the set (10 items) were suspended at time=3.0. Therefore the effect of censoring played a significant role in this exercise. The input file is here. The expected error of non-parametric estimation (as input for regression analyis) was between 0.05 and 0.1. This error is small if CIF>3.
Output:
Calculation #1. Mon Dec 25 10:34:03 EST 2023

Entered number of coefficients of Pade function is: 3
Required error of approximation is: 0.01
Probability function parameter beta=2.5087905 eta=1.0574883
Calculated correlation coefficient is 0.9983209551863421
Coefficient of determination is 0.9906786098460901
Maximum error is 0.12417020305677487 at x=1.89133269927662
Standard deviation (residual error) is 0.1766787910431362
Number of iterations without significant changing of Standard Deviation is 1

Approximation formula:
(x/1.0574883)2.5087905)(1.5663906-0.3242146X+0.029424906X2)/(1.0),
where X=x1.55

Red solid curve in figure corresponds to the exact solution, red signs represent input data generated by simulation, black solid curve is the regression solution which is compared with the result of maximum likelihood estimation (MLE) by Kijima model with constant restoration factor (dash-dotted line). One can see that regression method yields good approximation of input data. The standard deviation of this approximation from exact solution is 0.1953 while standard deviation of input data from exact solution is 0.2341. Increasing number of coefficients of Pade function to 4 we obtained the approximation with maximum error 0.1207 at x=1.89133269927662 and the standard deviation (residual error) 0.1709. Note that in this example the input data were generated using the same Kijima model, therefore the MLE procedure based on this model produced good result.
The output file is here

Example 2.
The input data for this example is generated by Monte Carlo method for Kijima model 1 with restoration factor q=0 and Weibull underlying function with shape parameter beta=2.0 and scale parameter eta=1.0. In total we had 20 components. All of them were suspended at time=4. The input file is here
Output:
Calculation #1. Thu Dec 14 14:35:37 EST 2023

Entered number of coefficients of Pade function is: 3
Required error of approximation is: 0.01
Probability function parameter beta=2.0727682 eta=1.0019438
Calculated correlation coefficient is 0.9990280583608917
Coefficient of determination is 0.9855101997954171
Maximum error is 0.07114915429082248 at x=3.91619793862935
Standard deviation (residual error) is 0.054700980486738555
Number of iterations without significant changing of Standard Deviation is 1

Approximation formula:
(x/1.0019438)2.0727682(1.2918677+0.06436612X)/(1.0+0.5501749X),
where X=x2

Red solid curve in figure corresponds to the exact solution, red signs represents input data generated by simulation, black solid curve is the regression solution. The standard deviation of obtained approximation from exact solution is also less than the standard deviation of input data from the exact solution.
The output file is here

Example 3.
This example from paper [1] shows the possibility of data extrapolation using regression analysis. We consider the Kijima model with increasing restoration factor. We have simulated a statistical sample of 100 trajectories (components) from a general renewal model with the underlying Weibull distribution (scale parameter eta=1 and shape parameter beta= 2) and with Kijima I restoration factor depending on failure number, i, in the following functional form: qi = qi-1 +0.5(1- qi-1). Specifically, with q1 = 0.2, q2 = 0.6, q3 = 0.8, ... The estimation interval was truncated at 2 units of time, corresponding to ~3 failures of the system. The input file is here

Output:
Calculation #1. Thu Dec 14 14:35:37 EST 2023

Entered number of coefficients of Pade function is: 5
Required error of approximation is: 0.01
Probability function parameter beta=2.0404139 eta=1.0431365
Calculated correlation coefficient is 0.9998009792440604
Coefficient of determination is 0.9996658781294799
Maximum error is 0.03442081447507269 at x=1.45346888328358
Standard deviation (residual error) is 0.016067707661711172
Number of iterations without significant changing of Standard Deviation is 0

Approximation formula:

(x/1.0431365)2.0404139)(0.088154934+3.2770913X-3.9587612X2+1.7845159X3-0.2693846X4)/(1),
where X=x0.75

Red curve in the figure corresponds to the exact solution. Black solid line represents obtained approximation; dash-dotted curve is the solution obtained by Kijima model with constant restoration factor. The example shows that by increasing the number of parameters of approximation function we can obtain better data extrapolation. However, be advised that increasing this number can lead to overfitting phenomena.
The output file is here

Example 4.
This example shows the possibility of Unavailability Q calculation as Failure Intensity (FI) function multiplied by mean time to repair. We consider the exponential function for failure time distribution. We simulated a statistical sample of 400 trajectories (components) from a general renewal model with the underlying Weibull distribution (scale parameter eta=1 and shape parameter beta= 1). The estimation interval was truncated at 3 units of time. The input file is here

Output:
Calculation #1. Thu Dec 14 14:35:37 EST 2023

Entered number of coefficients of Pade function is: 5
Required error of approximation is: 0.01
Probability function parameter beta=1.1864227 eta=1.15588
Calculated correlation coefficient is 0.9999434503387546
Coefficient of determination is 0.9992598657108345
Maximum error is 0.02204465 at x=2.99686113713495
Standard deviation (residual error) is 0.008039685093776072
Number of iterations without significant changing of Standard Deviation is 0

The output file is here
The result is shown by dash-dotted line. The red curve in the figure corresponds to the exact solution. The example shows that the model can be used for Unavailability calculation. Note that in some cases for approximation of the derivative of the function it is required more parameters (for example if Failure Intensity function has maximum).

References
1. A. Yevkin, V. Krivtsov, A generalized model for recurrent failures prediction, Reliability Engineering and System Safety 204 (2020).
2. A. Yevkin, V. Krivtsov, Modeling Recurrent Failure Processes using Padé Approximants, RAMS 2025. Proceedings of Reliability and Maintainability Symposium, RAMS 2025, Florida, US.