{'fraction_non_alphanumeric': 0.034833778619500155, 'fraction_numerical': 0.011295613649430146, 'mean_word_length': 4.498326515433247, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Reliability Analysis of Artificial Intelligence Systems Using Recurrent Events Data from Autonomous Vehicles

Artificial intelligence (AI) systems have become increasingly common and the trend will continue. Examples of AI systems include autonomous vehicles (AV), computer vision, natural language processing, and AI medical experts. To allow for safe and effective deployment of AI systems, the reliability of such systems needs to be assessed. Traditionally, reliability assessment is based on reliability test data and the subsequent statistical modeling and analysis. The availability of reliability data for AI systems, however, is limited because such data are typically sensitive and proprietary. The California Department of Motor Vehicles (DMV) oversees and regulates an AV testing program, in which many AV manufacturers are conducting AV road tests. Manufacturers participating in the program are required to report recurrent disengagement events to California DMV. This information is being made available to the public. In this paper, we use recurrent disengagement events as a representation of the reliability of the AI system in AV, and propose a statistical framework for modeling and analyzing the recurrent events data from AV driving tests. We use traditional parametric models in software reliability and propose a new nonparametric model based on monotonic splines to describe the event process. We develop inference procedures for selecting the best models, quantifying uncertainty, and testing heterogeneity in the event process. We then analyze the recurrent events data from four AV manufacturers, and make inferences on the reliability of the AI systems in AV. We also describe how the proposed analysis can be applied to assess the reliability of other AI systems.


The Problem
With the rapid development of new technology, artificial intelligence (AI) systems are emerging in many areas. Typical applications of AI systems include autonomous vehicles (AV), computer vision, speech recognition, and AI medical experts. The reliability and safety of AI systems need to be assessed before massive deployment in the field. For example, the reliability of AV needs to be demonstrated so that people can use them with confidence.
Traditionally, reliability assessment of products and systems is based on reliability test data collected from the laboratory and the field. Reliability information is then extracted from statistical modeling and analysis of the data.
Commonly used data types for reliability analysis are time-to-event data, degradation data, and recurrent events data. Reliability data collected by manufacturers are highly sensitive and are usually not publicly available. In the area of AV, however, the California Department of Motor Vehicles (DMV) has launched an AV driving program. More details of the study are given in Section 2. Many AV manufacturers are participating in the program and are conducting AV road test in California. Manufacturers are required to report disengagement events and mileage information to the California DMV. The reported events are available for public access. Because of the availability of these recurrent events data, we focus on the reliability analysis of the AI systems in AV units in this paper.
A disengagement event happens when the AI system and/or the backup driver determines that the driver needs to take over the driving. The recurrence rates of disengagement events can be used as a proxy for the reliability of the AI system in AV. A lower occurrence rate (event rate) of disengagement events would indicate a more reliable AI system in the AV. In the reliability literature, parametric forms have been used to describe the event rate through a nonhomogeneous Poisson process (NHPP) model. In practice, some specific questions arise in the analysis of the recurrent disengagement events data, • How to model the event process, and what kind of parametric forms should be used?
• Does the parametric form provide an adequate fit to the data, and are there any other flexible forms for modeling?
• Is there any population heterogeneity in the event processes from multiple test units?
To answer those practical questions, we develop a statistical framework for modeling and analyzing the recurrent events data from AV driving tests.
Specifically, we use the NHPP model to describe the disengagement event process with adjustment for the time-varying mileage data using traditional parametric models in software reliability for the intensity functions. We also propose a new nonparametric model based on monotonic splines to describe the event process. The spline model is flexible enough to fit various patterns in the event process and can also be used to assess if the fully-parametric model provides an adequate fit. We develop inference procedures for selecting the best models and quantifying uncertainty using the fractional-random-weight bootstrap. The parametric models and spline models can be used as complementary tools. In addition, we use the gamma frailty model to quantify and assess heterogeneity in the event process.
From the California driving study, we use data from four manufacturers that have been conducting extensive AV driving tests in California. We apply the developed methods to analyze the recurrent events data from the four AV manufacturers. Based on the modeling, we make inferences on the reliability of the AI systems in AV, and summarize interesting findings on the reliability of the AI systems in AV. Although our analysis focuses on the AI systems in AV, the statistical analysis can also be applied to assess the reliability of other AI systems.

Literature Review
There is only a small amount of literature on reliability analysis of AI systems.  Zhao et al. (2020) proposed a safety framework based on Bayesian inference for critical systems using deep learning models. Alshemali and Kalita (2020) provided a review on improving the reliability of natural language processing. Due to the limited availability of test data, statistical reliability analysis of AI reliability/safety is in an emerging stage.
As an example of the modeling of the reliability and safety of AV, Kalra  In the reliability literature, NHPP models are widely used to analyze recurrent events. Zuo, Meeker, and Wu (2008), and Hong, Li, and Osborn (2015) analyzed recurrent events data with window-observations, which has similar data types with the disengagement events data from the California driving study. Parametric models, such as the Musa-Okumoto model (e.g., Musa and Okumoto 1984) and the Gompertz model (e.g., Huang, Lyu, and Kuo 2003), are commonly used in software reliability applications (e.g., Wood 1996). Ehrlich et al. (1998) used accelerated testing methods to study software reliability. Burke, Jones, and Noufaily (2020) proposed flexible parametric models for time-to-event data analysis. Useful reference books for reliability data analysis for researchers in the AI reliability area include Escobar (1998), andLawless (2003). Overall, parametric models are common in analyzing recurrent events data in the context of reliability study.
We propose to use monotonic splines (Ramsay 1988, andMeyer 2008) as a nonparametric method to model the event process. Although monotonic splines are used in some degradation settings (Xie et al. 2018), the application of monotonic splines to recurrent events in reliability is new. To model population heterogeneity, the gamma frailty model is used. An early use of the gamma frailty model in reliability is found in Lawless (1995). More recently, Shan, Hong, and Meeker (2020) used the gamma frailty model to describe seasonal warrant return data. Duchateau and Janssen (2008) provide a comprehensive review of frailty models.
Due to the complicated structure of the window-observed recurrent events data in our study, we use fractional-random-weight bootstrap as a convenient way to generate bootstrap samples for statistical inference. The idea of fractional-random-weight bootstrap is introduced in Rubin (1981), and some theoretical properties are shown in Jin, Ying, and Wei (2001). Xu et al. (2020) demonstrated the use of fractional-random-weight bootstrap in many complex applications in reliability, survival analysis, and logistic regression. A simultaneous confidence band (SCB) can be used to assess if a parametric model is adequate for fitting the data. Hong, Escobar, and Meeker (2010) showed that an SCB could be obtained from a simultaneous confidence region (SCR) for parameters. However, in the context of bootstrap, it is not straightforward to construct SCR for parameter estimators with multiple dimensions. Hence, we use the idea of the equal-precision band in Nair (1984) and use bootstrap samples to calibrate pointwise confidence intervals to provide SCB to quantify statistical uncertainty.
In summary, we provide a general analytic framework by integrating existing methods and proposing new methods for reliability analysis of data from an AI study. The parametric and nonparametric models, and the statistical interval and testing procedures will be useful tools for practitioners working in the area of AI reliability.

Overview
The rest of the paper is organized as follows. Section 2 describes the California autonomous vehicle driving study and introduces the data. Section 3 describes the parametric model, the spline model, and the gamma frailty model that are used to describe the recurrent disengagement events data. Section 4 describes the parameter estimation procedures and the inference procedures for various models. Section 5 conducts a simulation study to show the statistical performance of the estimation procedures. Section 6 conducts the data analysis and summarizes interesting findings. Section 7 contains some concluding remarks and areas for future research.

Disengagement Events Data
The California DMV requires manufacturers to report when their vehicles disengaged from autonomous mode during tests. A disengagement event occurs when there is an autonomous technology failure, or when situations requiring the test driver to take manual control of the vehicle to operate safely. Disengagement events can be initialized by a warning from the autonomous vehicle system, or by test drivers as the driver thinks it is not safe to continue auto driving. Disengagement reports are provided annually, and include the ID number of vehicles in testing, location and date of disengagement events, description of causing of disengagement, monthly autonomous mileage of each testing vehicle, and annual total autonomous miles of each testing vehicle.
The starting date for the data used for this paper is December 1, 2017. The study period is from December 1, 2017, to November 30, 2019, which is a 24 month or 2 year study period.
The data for the period from December 1, 2018, to November 30, 2019, are available in csv format from the California DMV website, while the data for the period from December 1, 2017, to November 30, 2018, are available in pdf format that need to be manually converted to csv format.
After data cleaning, the event time is computed as the number of days since the starting date. The monthly mileage is divided by the number of days in that month to obtain the daily mileage, under the approximation that the daily mileage is constant over each month. The unit for the mileage is thousands of miles (k-miles). Figure 1 shows a subset of the recurrent events data and mileage data from manufacturer Waymo. Figure  plots the mileage as a function of time for five representative units. Table 1 shows a summary of the recurrent events data and mileage data from the four manufacturers. We can see that Waymo and Cruise have driven more than 1 million miles and the disengagement event rate is around 0.1 events per k-miles during the 24 month testing period. Pony AI and Zoox have smaller driving miles, and the event rates are around 0.2 events per k-miles and 0.6 events per k-miles, respectively.

Notation for Data
The number of AV testing units (vehicles) in the fleet is denoted by n. The total follow-up time is τ = 730 days (i.e., two years). Let t ij , i = 1, . . . , n, j = 1, . . . , n i be event time j for unit i. Here t ij records the number of days since December 1, 2017, and n i is the number of events for unit i. Note that it is possible that n i = 0, indicating that there were no events observed for unit i.
Let x i (t), 0 < t ≤ τ , be the mileage driven for unit i at time (day) t. The unit of x i (t) is k-miles. Because only monthly mileage was reported, the daily average was used for x i (t). Thus, x i (t) is a step function, which can be represented as, Here, n τ = 24 is the number of months in the follow-up period, x il is the daily mileage for unit i during month l, τ l is the ending day since the starting of the study for month l, and Here, λ 0 (t; θ) = λ 0 (t) is the baseline intensity function (BIF) with parameter vector θ. Be- , θ] is the mileage-adjusted event intensity. The BIF can be interpreted as the event rate per k-miles at time t when x i (t) = 1. The baseline cumulative intensity function (BCIF) is, Note that Λ 0 (0; θ) = 0 and Λ 0 (t; θ) is a non-decreasing function of t. The BCIF Λ 0 (t) can be interpreted as the expected number of events from time 0 to t when x(t) = 1 for all t. The cumulative intensity function (CIF) for unit i is, Model In software reliability, the BCIF and the BIF are used as reliability metrics when recurrent events can be collected from repairable systems (e.g., Wood 1996). The trends in the BIF can indicate the evolution of the reliability in the underlying AV system. For example, improvement of the autonomous technology in AV can lead to a decreasing trend in the BIF (i.e., BCIF increasing at a decreasing rate).
Typically in software reliability, one specifies a parametric form for the BCIF Λ 0 (t; θ).
Note that the BIF can be obtained by taking the derivative of BCIF with respect to t.

Spline Models
Although parametric models can fit certain trends in the event process, they may not be flexible enough to describe the event process for AV testing, as the evolution of the AI technology in an AV system can be complicated, which motivates us to propose the spline model for describing the BCIF. In the spline model, the BCIF is represented as a linear combination of spline bases. That is, which provides a nonparametric method to describe the BCIF. Here θ = (β 1 , . . . , β ns ) is the vector for the spline coefficients, γ l (t)'s are the spline bases, and n s is the number of spline bases. The BIF can be obtained by taking a derivative with respect to t. That is, Because of the constraints that Λ 0 (0; θ) = 0 and that Λ 0 (t; θ) is a non-decreasing function of t, some special considerations are needed in the spline model. We use the I-splines in Ramsay (1988). Figure 2 shows examples of I-spline basis functions (i.e., γ l (t)). We can see that each spline basis takes value zero at t = 0 and is monotonically increasing. By taking non-negative coefficients (i.e., β l ≥ 0), a non-decreasing Λ 0 (t; θ) is obtained.
A brief introduction on the construction of I-spline basis is as follows. The boundary knots are 0 and τ . The b interior knots are denoted by t h+1 , . . . , t h+b for splines of order h. The complete sequence for the knots are denoted by 0 = t 1 = · · · = t h < t h+1 < · · · < t h+b < t h+b+1 = · · · = t 2h+b = τ. The total number of spline bases is n s = h + b. I-splines are obtained by integrating the M-splines; that is, and the M-spline bases of order h are defined recursively. The M-splines of order 1 is for q = 1, · · · , b + 1. The M-splines of order h are obtained as for q = 1, · · · , b + h.

Modeling Heterogeneity
It is maybe possible that some AV units are more likely to generate more events even after accounting for the mileage driven, resulting in heterogeneity in the event process. Similar to the approach in Chapter 3 of Cook and Lawless (2007), a frailty term can be added to the parametric intensity function to model this extra heterogeneity. The gamma frailty model is, Here, the frailty term u i is a random variable that follows a gamma distribution with mean one and variance φ. The probability density function (pdf) of u i is, The gamma frailty model is popular in reliability and survival applications because there is a closed-form expression for the marginal likelihood for recurrent events data based on the model in (7).

Model Estimation and Inference
This section presents parameter estimation procedures for the parametric and spline models, and the corresponding statistical inference procedures. This section also contains parameter estimation and hypothesis testing for the gamma frailty model.

Parameter Estimation
We use the maximum likelihood (ML) methods for parameter estimation. The likelihood function is, with the convention that 0 j=1 (·) = 1. Here, the intensity function and the CIF are defined in (2) and (4), respectively. The log-likelihood function is obtained by taking the logarithm of the L(θ) in (8). That is, with the convention that 0 j=1 (·) = 0. The last step in (9) is obtained by substituting the expression for x i (t) in (1).
For parametric models, the functional forms of Λ 0 (t) and λ 0 (t) in Table 2 can be substituted into (9) to evaluate the log-likelihood function. The ML estimate of θ, denoted by θ, can be obtained by finding the value of θ that maximizes l(θ). For parametric models, the dimension of θ is typically 2 or 3. We used the R optim() function with the "Nelder-Mead" option to do the optimization. For the spline model, the functional forms of Λ 0 (t) and λ 0 (t) in (5) and (6), respectively, can be substituted into (9) to evaluate the log-likelihood function. To estimate the parameter θ for the spline model, one needs to specify the locations of the knots and the number of knots, and address the non-negativity constraints on θ as indicated in (5).
We use I-splines of degree 3, which is generally smooth enough to fit the BCIF. The To select the best number of knots, we use the Akaike information criterion (AIC), which is computed as Here, df is the number of degrees of freedom for the model. For the spline model, df is the number of non-zero coefficients. To address the non-negativity constraints on θ, we use the "L-BFGS-B" option in the R optim() function. The "L-BFGS-B" algorithm allows users to specify an interval for the variable to be optimized. We set the interval to be [0, ∞) for the spline coefficients. Using the "L-BFGS-B" algorithm, some of the elements of the estimates θ could be set to zero.

Statistical Inference
For inference based on parametric models, the normal approximation based on large sample theory can be used. Because the inference for parametric models is relatively straightforward, the focus of this section on the inference for the spline model. For the spline model, the normal approximation is inappropriate because the parameter estimates can occur at the boundary of the parameter space (i.e., β l can be zero). Thus, we use the fractional-random-weight bootstrap (e.g., Xu et al. 2020). Different from other bootstrap methods, the fractionalrandom-weight bootstrap provides a convenient way to quantify uncertainty for data with complex structure. For our case, we need to handle recurrent events with window observations and time-varying mileage information, which complicate the structure.
The implementation of the fractional-random-weight bootstrap is straightforward. The log-likelihood function in (9) is re-weighted as where the random weights are independently generated from an exponential distribution with mean one. The bootstrap algorithm for generating bootstrap versions of the estimate of BCIF Λ * 0 (t) is summarized as follows.
1. Generate random weights w i , i = 1, . . . , n, independently from an exponential distribution with mean one.
3. Use AIC to select the best number of knots based on the log-likelihood function in (10).
Based on the bootstrap estimates Λ * b 0 (t), b = 1, . . . , B, one can construct approximate 100(1 − α p )% pointwise confidence interval (PCI) for Λ 0 (t) for a given t, Here, Λ * (b) 0 (t) is the ordered version of Λ * b 0 (t), and [ · ] is the rounding function. Based on the spline model, one can construct a simultaneous confidence band (SCB) for the BCIF, which can be used to assess the adequacy of the parametric models in fitting the data. An approximate 100 where t L and t U are boundaries to be specified. If an estimated parametric model is contained in the SCB in (11) the parametric model is statistically consistent with the data (i.e., there is no statistical evidence to reject the parametric model). Thus, the SCB constructed by the spline model provides a tool to check if a parametric model is adequate.
The 100(1 − α p )% PCIs provide a structure for computing SCBs for Λ 0 (t) for all in the t ∈ [t L , t U ] period. We use bootstrap samples to calibrate the PCIs so that it can approximately provide the nominal 100(1 − α)% CP, similar to the idea of the equal-precision SCB for a cdf in Nair (1984). The CP can be estimated as By setting CP(α p ) = 1 − α, one can find the solution to be α c . Thus, the SCB in (11) can be computed as which is time-efficient because the bootstrap samples have already been obtained.

The Frailty Model
To estimate the gamma frailty in (7), one needs to calculate the marginal likelihood function.
The marginal likelihood function is, The ML estimates of θ and φ are obtained by finding the values that maximize the log-likelihood function, log[L(θ, φ)]. We again use the R optim() function with the "Nelder-Mead" option to do the optimization.
To check if there is population heterogeneity in the event process, one can use the likelihood ratio test. The test statistic is constructed as, which has a χ 2 1 distribution under the null hypothesis that φ = 0. Here, l( θ) is obtained by evaluating (9) at θ. The null hypothesis is rejected if the test statistic is larger than χ 2 1, (1−α) , indicating that there is population heterogeneity.

Simulation Study
The purpose of the simulation study is to show the properties of the ML estimator for the spline model, to check the CP of the SCB procedure based on the spline model, and see if the SCB can correctly accept or reject a particular parametric model.

Setting
In the simulation study, the spline model is considered as the true underlying model. The spline bases are shown in Figure 2. We consider three scenarios. Figure 3 shows the true BCIFs used in the three simulation scenarios. To simplify the setting, we only consider one parametric model, which is the Gompertz model.
• Scenario 1: we choose θ = (6, 16, 23, 11, 4) as the coefficients. Using the spline bases in Figure 2, the true BCIF is shown as the black/solid line in Figure 3. The Gompertz model fits this BCIF perfectly.
• Scenario 2: we choose θ = (8, 12, 28, 0, 12) as the coefficients. The corresponding true BCIF is shown as the red/dash line in Figure 3. The Gompertz model fits this BCIF well except for the late stage. We evaluate the relative root mean squared errors (RMSE) for the BCIF estimator, the CP for the SCB, and the acceptance probability for the parametric model. In particular, the relative RMSE (RelRMSE) is computed as where Λ 0l (t) is the estimated BCIF using the spline model based on the lth simulated dataset, and Λ 0 (t) is the true BCIF. We use the relative RMSE to remove the scale effect of Λ 0 (t). That is, RMSE tends to be large if Λ 0 (t) is large at a particular t. The CP is estimated by the proportion that the SCB constructed by using the spline model captures the true BCIF.
The acceptance probability is estimated by the proportion of times that the spline-based SCB captures the estimated BCIF from the Gompertz model. In summary, the ML estimator for the spline model works as expected. When a parametric BCIF is adequate, it tends to have higher statistical efficiency. When the parametric model is not adequate, there could be large RelRMSE due to bias in the estimation. The results show that it is important to assess the adequacy of parametric models. The spline model, however, is flexible at the price of losing some statistical efficiency. Figure 5 shows the plot of CP and acceptance probability as a function of sample size under the three scenarios. From Figure 5(a), the CP for the SCB procedure based on the spline model and bootstrap is quite similar for the three scenarios. The CP improves when the sample size increases. The CP is less than nominal when the sample size is small and is getting closer to the nominal CP when the sample is larger than 200. From the plot of the acceptance probability in Figure 5(b), the parametric model is generally accepted when the sample size is small and there is little or no departure from the true model. When there is a large departure from the true model, as in Scenario 3, the SCB does not capture the estimated Gompertz model with high probability when the sample size is large.

Data Analysis
In this section, we present the data analysis for the recurrent disengagement events data.

Model Fitting
We start by fitting the parametric models in Table 2    To further check how well the models fit the data, Figure 7 shows  Table 4 shows the parameter estimates, standard errors, and approximate 95% confidence intervals (CIs) for the best parametric models for the four manufacturers. The CIs for parameters are based on a normal (Wald) approximation, and some transformations (e.g., a logarithm transformation for positive parameters) are used to improve the performance. Note that the estimate for θ 1 for the Zoox/Musa-Okumoto is set to zero because the model is degenerate at θ 1 = 0, indicating a constant rate situation.

Results Interpretation
We also tested population heterogeneity using the procedure in Section 4.3. Table 5 shows the summary of the gamma frailty models for the four manufacturers. All of the p-values are  Table 5 also lists the names for the best parametric models for each manufacturer.
To further visualize the results, Figure 8 plots the BIFs from the best parametric models for each of the four manufacturers. From the plot, we see the trends for different manufacturers.
A decreasing trend means an improvement in AI technology and an increase in AV reliability. for different manufacturers, we can see that the overall AV reliability is improving over the two-year period.
As a comparison, Figure 9 shows the estimated BIFs based on the spline model and the parametric models, and the 95% PCIs based on the spline model. We can see that the spline model shows more variation but the general trends are the same as the parametric models.
The estimates for Zoox have a large amount of variability due to the limited sample size, as discussed previously.

Concluding Remarks
This paper focuses on the reliability analysis of AV technology using the recurrent disengagement events from the California driving study. We propose a statistical framework for modeling and analyzing the recurrent events data from AV driving tests. We use both parametric models and a nonparametric model to describe the event processes. Based on the spline model, we can select the best model, quantify uncertainty, and test heterogeneity in the event process. We want to point out that the parametric models and spline models are proposed as complementary tools for modeling and inference.
Through the simulation and data analysis, we show that the proposed spline model is flexible for describing the recurrent events data from four AV manufacturers, and parametric models are adequate for data from most manufacturers. It is worth noting that the best parametric models can be different for different manufacturers. The population heterogeneity   in the event process is also low. From the data analysis, we found that the overall AV reliability is improving over the two-year study period.
The currently available data do not include covariates, such as the driving speed when the event occurred, the test environment (e.g., busy street versus freeway), and vehicle models.
In the future, it would be interesting and useful to collect more covariate information. Our proposed modeling framework can be extended to analyze recurrent disengagement events data with covariates. The CA DMV recently started a driverless program where cars on the road do not have to have a driver. It would be interesting to analyze the driverless study data in the future when enough event data are available.
Because reliability is a property that evolves over time, all kinds of AI systems need to be tested over time to quantify their reliability. Although our analysis focuses on AV reliability, our proposed data analytic framework can also be applied to assess the reliability of other AI systems where departures from desired behavior provide recurrent events data. With an appropriate definition of time scale and events, the parametric and spline models discussed in this paper can be extended to analyze data from the reliability testing of other AI systems.
Computing hardware reliability is also an important aspect of AI reliability. For example, GPU's are widely used in AI model computing. Ostrouchov et al. (2020) considered the lifetimes of GPUs used supercomputers. It would be interesting to study GPU reliability, or more broadly computing hardware reliability, and its relationship to AI reliability.