{'fraction_non_alphanumeric': 0.03176562688184993, 'fraction_numerical': 0.007497290136095387, 'mean_word_length': 5.127859778597786, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Closed-Loop Statistical Verification of Stochastic Nonlinear Systems Subject to Parametric Uncertainties

This paper proposes a statistical verification framework using Gaussian processes (GPs) for simulation-based verification of stochastic nonlinear systems with parametric uncertainties. Given a small number of stochastic simulations, the proposed framework constructs a GP regression model and predicts the system's performance over the entire set of possible uncertainties. Included in the framework is a new metric to estimate the confidence in those predictions based on the variance of the GP's cumulative distribution function. This variance-based metric forms the basis of active sampling algorithms that aim to minimize prediction error through careful selection of simulations. In three case studies, the new active sampling algorithms demonstrate up to a 35% improvement in prediction error over other approaches and are able to correctly identify regions with low prediction confidence through the variance metric.


I. INTRODUCTION
Ever increasing demands for high performance necessitates the adoption of advanced, nonlinear control techniques in a wide variety of engineering systems. The complexity of these controllers, coupled with the increasing complexity of the open-loop systems themselves, complicates verification that the closed-loop system actually satisfies the performance requirements at off-nominal conditions, called parametric uncertainties. Control system verification analyzes the closedloop system and attempts to identify at which parametric uncertainties the resulting trajectory fails to meet the requirements.
Traditional analytical and numerical verification techniques [1], [2] provably verify the system's performance at various uncertainties, but require modeling assumptions that restrict their suitability to complex systems with highly nonlinear dynamics and hybrid controllers. Recent work [3]- [5] addressed some of these concerns, but are limited by the size of the state and parameter spaces. Stochastic dynamics further challenge these approaches and restrict their applicability. For comparison, statistical verification techniques [6], [7] replace provable guarantees with looser probabilistic bounds, but are capable of handling a much broader range of problems with fewer restrictions. Statistical verification methods rely upon extensive simulations of the system to test the likelihood of requirement satisfaction at various conditions and are well-suited to stochastic systems.
Recent work developed data-driven statistical verification techniques [8], [9] to efficiently verify deterministic closedloop systems without relying upon exhaustive simulations. These data-driven methods construct machine learning models from small sets of simulation trajectories and predict the response at unobserved parameter settings. The latest work [10] introduced Gaussian process (GP) prediction models and novel active sampling procedures to improve the accuracy of predictions given limited amounts of simulation data. The deterministic GP-based approach also quantifies the prediction confidence online without the need for external validation datasets and carefully selects simulations to maximize the prediction confidence.
This paper extends data-driven methods to stochastic verification problems and provides three contributions: a GP-based statistical verification framework for stochastic systems, a metric to estimate the accuracy of the framework's predictions, and active sampling algorithms specifically developed for closed-loop statistical verification. Given a small set of stochastic simulations, the first two contributions will predict the probability of requirement satisfaction at all possible parametric uncertainties and compute the confidence in those predictions. The third contribution exploits this prediction confidence to select subsequent simulations and minimize prediction error. Three case studies are provided to demonstrate the effectiveness of the approach to maximize accuracy given a constraint on the number of simulations.

II. PROBLEM DESCRIPTION
Consider the following nonlinear systeṁ where x(t) ∈ R n is the state vector, u(t) ∈ R m is the control input vector, θ ∈ R p represents the parametric uncertainties, and w(t) ∈ R n is the stochastic noise. A candidate controller generates control inputs according to some predetermined function u(t) = g(x(t)). Given this controller, the resulting closed-loop dynamics reduce tȯ and the trajectory Φ(x(t)|x 0 , θ) is a function of nominal initial state x 0 , parametric uncertainties θ, and stochastic noise w(t). We assume the nominal initial state x 0 is known and fixed. The parametric uncertainties θ can arise from a variety of sources. For instance, these parameters may include system properties such as vehicle mass and inertia that will typically vary across different operational scenarios. The main impact of the parametric uncertainties is as variable initial conditions that determine the evolution of the subsequent state trajectory Φ(x(t)|x 0 , θ). These parametric uncertainties may not be known during real-world execution of the system, but are assumed to fall within some bounded, compact set Θ that is known in advance.

Assumption 1
The set of all possible parametric uncertainties θ ∈ Θ is a known, compact set Θ ∈ R p . Assumption 1 ensures the system has known feasible bounds on the uncertain conditions for the verification procedure to examine. Most physical systems will naturally have bounds on the set of allowable operating conditions. For instance, aircraft will typically fly within a maximum take-off weight set by structural limitations and a minimum (empty) weight.
The major difference between (2) and the deterministic systems found in earlier problems [8]- [10] is the addition of stochastic noise w(t) in the dynamics. This noise term encapsulates any source of randomness, such as process and/or measurement noise and ultimately breaks the previous works' fundamental assumption of deterministic (noise-free) dynamics. The randomness will cause multiple trajectories starting from the same x 0 and θ to obtain different results. The noise could have a large effect on the satisfaction of the performance requirements at a given x 0 and θ.

A. Measurements of Trajectory Performance
A trajectory must meet certain performance requirements in order for it to be considered "satisfactory." As before [10], these requirements are supplied by relevant certification authorities and may include a wide range of potential specifications with varying complexity. Commonly, the requirements can be written as signal temporal logic (STL) [11] specifications, which also offers a convenient method to quantify the satisfaction of the requirements. Central to STL is the computation of a scalar robustness degree ρ ϕ ∈ R that quantifies the minimum robustness of the trajectory with respect to the requirement specification ϕ. Therefore, the satisfaction of the requirements along the entire trajectory is captured with a single scalar variable. Building upon that concept, this work assumes the robustness of a trajectory Φ(x(t)|x 0 , θ) to a given performance requirement is measured with a single measurement variable y(θ), where y(θ) = ρ ϕ in problems with STL specifications. Other applications will specify the performance requirements and quantify their satisfaction through different methods, as will be seen in the task allocation example from Section V-B.

Assumption 2
The simulation model provides a single, scalar output y(θ) ∈ R that quantifies the minimum robustness of trajectory Φ(x(t)|x 0 , θ) to a specified performance requirement. The sign of y(θ) indicates satisfaction of the requirements, where y > 0 signifies "satisfactory" performance and y ≤ 0 signifies "unsatisfactory" performance.
Since the closed-loop trajectory is a function of nominal (fixed) x 0 and θ, we write the trajectory robustness measurements as an explicit function of θ to emphasize the effect of the parametric uncertainties upon the robustness of the trajectory. As θ changes, so will the satisfaction of the performance requirement.

B. Probabilistic Satisfaction of Requirements
Although a closed-loop trajectory and its corresponding robustness measurement are functions of parametric uncertainties, stochasticity will cause multiple trajectories at the same θ setting to obtain different y(θ). If a large number of simulations are performed at this setting, then the corresponding y(θ) measurements will produce a continuous distribution of possible robustness values. The underlying true distribution of y(θ) may take a variety of different forms, but this work addresses the case with a Gaussian distribution.
While this assumption does restrict the class of problems, it simplifies the implementation details to improve reader clarity and better highlight the contributions of the approach. This paper is intended to introduce the statistical verification framework and discuss the challenges and contributions on clearer, and still relevant, verification problems. Upcoming extensions will relax Assumption 3 to consider spatiallyvarying y [12] and non-Gaussian distributions [13], but also require significantly more complex and sensitive statistical inference techniques to model those distributions. The implementation details will change in these extensions, but the fundamental approach remains the same.
From the verification perspective, where each trajectory either satisfies the requirement or does not, the Gaussian distribution of y(θ) creates a Bernoulli distribution for the likelihood an arbitrary trajectory at θ will satisfy the requirement.

Definition 1
The satisfaction probability function p sat (θ) ∈ [0, 1] defines the probability an arbitrary trajectory initialized with θ will satisfy the performance requirement.
This probability of satisfaction is the cumulative distribution of positive y(θ) measurements, While p sat (θ) defines the probabilistic satisfaction of the requirements, the challenge for stochastic verification is p sat (θ) will generally be unknown in advance. Instead, the statistical verification goal is to compute an estimated p sat (θ) for all possible conditions in Θ and minimize the difference between the predictions and p sat (θ).
Simultaneously, statistical verification must also contend with computational costs. In most applications, external constraints restrict the computational budget allocated to the statistical verification process. For instance, verification typically relies upon high-fidelity models of complex systems and subsystems that drive up the computational time for a single trajectory. It is infeasible to simply saturate Θ with multiple simulations at every θ because each simulation trajectory requires a certain amount of time and computational resources. Regardless of the source, each application has some feasible limit on the amount of time, resources, or money allocated to the statistical verification process. This work models the feasible limit as an upper bound N total on the number of simulation trajectories.

III. STATISTICAL VERIFICATION FRAMEWORK
Given the verification objective, there are multiple possible approaches to compute p sat (θ) from a finite set of simulation-based observations. The most basic form is to treat each trajectory as a Bernoulli trial with binary evaluations {"satisfactory", "unsatisfactory"} from the sign of y(θ). While this approach can address all types of systems, including those that do not meet Assumptions 2 and 3, a binary-based approach requires multiple simulations at every θ to construct a binomial distribution. Each repetition counts against the sampling budget |L| ≤ N total and means one few simulation at another location. Although the cost of each simulation is less of a problem for applications with simple simulation models, the cost of training a prediction model given a large number of simulations at each location is non-negligible. Therefore, the large number of simulations required for practical predictions with binomial distributions restricts the tractability of those approaches.
Instead, this work entirely avoids the need for expensive binomial distributions in the relevant class of systems through the direct use of y(θ) and its Gaussian probability density function (PDF). The main idea is to exploit a single stochastic trajectory at each training location and its noisy measurement y(θ) to infer the underlying PDF, which defines the cumulative distribution p sat (θ). Although various computational methods are possible [14], [15], this problem is particularly well-suited to Gaussian process regression [16] due to the Gaussian distribution of y(θ).

A. GP-based Prediction Model
The Gaussian process regression model follows a similar format to the earlier work in verification of deterministic systems [10]. A finite collection of N total simulation trajectories forms a training dataset L = {D, y} consisting of their parameter settings D = {θ 1 , θ 2 , . . . , θ N } and corresponding robustness values y = [y(θ 1 ), . . . , y(θ N )] T . Unlike the deterministic approach, these stochastic measurements require the introduction of a Gaussian likelihood model, where the GP does not model y(θ) directly, but infers the latent mean y(θ).
The training process constructs the GP regression model from the information provided by the training dataset L.
More details are found in [16], but fundamentally the training procedure uses Bayesian inference to place a posterior probability distribution on latent meanȳ given L. Assuming a zero-mean prior N (ȳ|0, K) and likelihood model N (y|ȳ, 2 y I), the posterior predictive distribution at an arbitrary location θ * follows a Gaussian distribution P(ȳ(θ * )|L, θ * ) = N (µ(θ * ), Σ(θ * )). The posterior predictive mean µ(θ * ) and covariance Σ(θ * ) are given by where Different choices are possible, but this work uses the common squared exponential kernel with automatic relevance determination (SE-ARD) [16] for the kernel function.
The choice in GP hyperparameters may drastically change the predictive distribution, even with the same training set L, since these hyperparameters control the kernel function κ and the likelihood model. Unfortunately, the optimal choice of hyperparameters that perfectly replicatesȳ(θ) will not be known in advance and must be estimated online. This work uses maximum likelihood estimation (MLE) [16] to optimize the kernel hyperparameters given only the information provided by L. The use of the SE-ARD kernel also enables the MLE procedure to adjust the sensitivity of the kernel to each element of θ asȳ(θ) may have a higher sensitivity to certain elements than others. Additionally, the likelihood model contains a hyperparameter to estimate standard deviation y as this term may also be completely unknown. If y is unknown, then an estimate y may be computed in the same MLE optimization process as the kernel hyperparameters or separately by repeated sampling at a training locations in D.

B. Prediction Confidence
While (5) computes the expected probability of requirement satisfaction, p sat (θ) will likely fail to perfectly model p sat (θ). Therefore, the statistical verification framework must not only provide p sat (θ), but also indicate where its confidence in the accuracy of these predictions is low. Unsurprisingly, the true prediction error p sat (θ) = p sat (θ)− p sat (θ) is unknown, but offline validation methods using external validation datasets [14] can estimate p sat (θ) by comparing the predictions against the validation set's known, true values. However, these external validation methods are wasteful as they generally require the removal of valuable training data for the independent validation set.
For comparison, probabilistic inequalities provide theoretically-justified bounds on p sat (θ) without the use of validation sets. In particular, Chebyshev's inequality bounds the prediction error by the variance of the CDF, where a > 0. Lower CDF variance will translate into lower probabilistic bounds on p sat (θ) and thus higher confidence in the accuracy of p sat (θ). The primary challenge with (6) is the variance lacks an analytical closed-form solution.
Fortunately, the variance can be approximated using a 1st or 2nd order Taylor series expansion of the nonlinear random function [17]. The 1st order approximation is which also happens to conservatively upper bound the 2nd order approximation. For simplicity, we subsequently refer to the approximate CDF variance in (7) as V (θ|L).
Despite the convenience of (7), the accuracy of an approximation of a nonlinear random function is limited. This inaccuracy makes it inadvisable to blindly substitute (7) into (6) without careful consideration. However, the approximation in (7) does still provide a perfect metric to identify regions of Θ where the confidence in p sat (θ) is low and provide intuition on the sensitivity of p sat (θ). The power of (7) is its value as a computationally-efficient, online validation tool for signifying prediction confidence. For instance, approximate variance is highest in regions with small |µ(θ)| and large Σ(θ), meaning the confidence in the predictions should be quite low. The next section will exploit the approximate CDF variance to derive a closed-loop verification process that seeks to minimize prediction errors.
IV. CLOSED-LOOP STATISTICAL VERIFICATION As discussed in Section II, statistical verification will have some feasible limit on the amount of time or computational resources allocated to the overall process. This limit is approximated as a cap on the number of simulation trajectories, |L| ≤ N total . Given this restriction on the size of training dataset L, the ideal scenario would only perform informative simulations and carefully select all θ ∈ D to minimize the prediction error over all θ ∈ Θ. However, the informationmaximizing training dataset is not apparent until after all the trajectories have been obtained. To address that issue, this work applies active learning [18] to iteratively select informative θ settings for future simulations and minimize the expected prediction error. To emphasize the iterative, feedback-based nature of the active learning procedure, we label the process closed-loop statistical verification.
Active learning describes a wide variety of different procedures [18]- [23], each with their own definition for the "best" sample to run next. Most of these procedures focus on a particular aspect of the Gaussian PDF. For instance, procedures focused on the PDF mean [9], [19] favor points withȳ(θ) near zero, meaning the best location is θ = argmin|µ(θ)|. Although originally intended for binary classification with support vector machines, such an approach does correctly label points with low |µ(θ)| as informative since the CDF variance (7) is high there. Likewise, PDF variance-based approaches and extensions [20]- [23] are significantly more common for GP methods and aim to reduce the PDF variance Σ(θ). These approaches favor points with high variance, θ = argmax Σ(θ), which also correctly emphasizes points with comparatively high CDF variance since Σ(θ) will be large. Although they both correctly emphasize certain aspects, neither of those two approaches explicitly minimizes the posterior CDF variance.

A. Reduction in CDF Variance
Earlier work in deterministic closed-loop verification [10] specifically developed new sample selection metrics to maximize prediction confidence. Even though the implementation details have changed, this work utilizes the same underlying motivation and attempts to minimize the approximate CDF variance (7) in order to maximize the confidence in p sat (θ). The ideal selection criterion would minimize the maximum posterior CDF variance V (θ|L + ) after the new sample data at θ was added, but this requires the posterior training set L + = L∪{θ, y(θ)} to be known apriori. Prior knowledge of L + is an impossible proposition since y(θ) cannot be known before a simulation has actually been performed there. Even if the expected posterior set L + = L ∪ {θ, µ(θ)} replaces infeasible L, the high computational cost of retraining the GP at every prospective sample location, nominally an O(N 3 ) operation and O(N 2 ) at best, renders the approach computationally intractable.
A more computationally tractable approach maximizes the local improvement in posterior CDF variance, effectively selecting the θ setting which will experience the greatest reduction in CDF variance if a simulation is performed at that location. The local change in CDF variance is labeled by V (θ|L) = V (θ|L)−V (θ| L + ). Although the change V (θ|L) requires expected posterior information L + , it does not require the GP model to be retrained, which was the source of the previous computational intractability. The local posterior CDF variance at location θ * after a measurement there can be written purely in terms of the current information. After the Woodbury matrix identity, the expected posterior covariance Σ(θ * ) + reduces to while mean µ(θ * ) + = µ(θ * ). Ultimately, the local change in CDF variance is given by Note that the CDF variance reduction criterion in (9) favors points with low |µ(θ)| and high Σ(θ), just as in (7). Perform simulation at θ, obtain measurement y(θ) 6: Add {θ, y(θ)} to training set L, remove θ from U 7: Retrain GP model with updated L 8: end for 9: Return: expected p sat (θ) and variance V θ|L

B. Sampling Algorithms
The section criterion from (9) forms the basis of the closed-loop stochastic verification framework. Since the actual set Θ is infinite, the framework constructs a high resolution lattice Θ d to replicate Θ and then selects all training locations from Θ d rather than Θ directly (D ⊂ Θ d ). Set U = Θ d \ D contains all remaining points not in training set L available for future simulations. This paper presents two versions of the closed-loop verification framework.
1) Sequential Sampling: The most straightforward implementation of closed-loop verification is sequential sampling, described in Algorithm 1. The process starts with an initial training set L of passively-selected samples, usually obtained through random sampling or another open-loop design of experiments approach. This initial training set of size |L| = N 0 produces the initial GP so active sampling can be performed. The sequential procedure then selects one sample from U according to (9) before it performs one simulation at the selected θ vector, obtains y(θ), and updates the GP prediction model. This iterative procedure repeats until the number of remaining samples T = N total − N 0 has been reached.
2) Batch Sampling: In comparison to Algorithm 1, batch sampling [18] offers further computational efficiency by selecting multiple samples between retraining steps. Assuming the same limit N total , batch selection reduces the number of GP retraining processes and can also fully exploit any parallel computing capabilities of the simulation environment. Batch approaches select M samples at once and perform their simulations in parallel. While this provides computational savings, it also introduces the possibility of redundant samples if the batch sample set S does not possess adequate diversity.
First proposed for deterministic closed-loop verification in [10], determinantal point processes (DPPs) [24] present efficient probabilistic methods for encouraging diversity within S without severe computational overhead. The batch framework first converts V (θ|L) into a probability distribution P V (θ) = V (θ|L)/Z V , where Z V is the normalization constant. A relatively small number of M T samples are obtained from Θ d according to P V (θ) and used to construct a DPP. This paper uses M T = 1000, but these are simply Transform V (θ|L) into distribution P V (θ) 6: Form k-DPP from M T random samples of P V (θ) 7: Generate M random samples from DPP, add to S 8: Run simulation ∀θ ∈ S, obtain measurements y S 9: Add {S, y S } to training set L, remove S from U 10: Retrain GP model with updated L 11: end for 12: Return: expected p sat (θ) and variance V θ|L samples of θ locations from the distribution and not actual simulations are performed. The DPP randomly selects M locations from U for S based upon a modified version of P V that penalizes similarities in S and subsequently spreads the datapoints out across high-valued regions in Θ d with significantly less redundancy. Algorithm 2 details the batch closed-loop verification framework. Rather than T number of additional simulations, Algorithm 2 operates in T batches of M simulations, assuming M · T ≤ N total − N 0 .

V. EXAMPLES
This section examines statistical verification applied to three stochastic systems and demonstrates the improved performance of Algorithms 1 and 2 over existing techniques.

A. Example 1: Model Reference Adaptive Control System
The first example considers a stochastic version of the concurrent learning model reference adaptive control (CL-MRAC) system from earlier deterministic work [9], [10]. The CL-MRAC example examines a second order linear system with two uncertain parameters θ = [θ 1 , θ 2 ] T . The adaptive control system estimates these parameters online and attempts to track a desirable reference trajectory. The nonlinearities associated with adaptation result in a highly nonlinear closed-loop system.
The verification goal is to predict the probability the actual trajectory will remain within 1 unit of the reference trajectory at all times between 0 ≤ t ≤ 40 seconds, given in STL format as where tracking error e 1 (t) is the difference between actual position x 1 (t) and reference position x m1 (t). The STL operator  N (0, I). The underlying true distribution was obtained by repeated sampling at each of the points in Θ d . To avoid any potential issues, Gaussian distributions were fit to the raw data and the variance was averaged across Θ d to return y = 0.0372. Figure 1 compares the performance of Algorithm 2 against similar procedures using the existing selection metrics discussed in Section IV as well as an open-loop, random sampling procedure. These procedures all start with an initial training set of 50 simulations and select batches of M = 10 points until a total of 450 simulations has been reached. Neither the kernel hyperparameters nor y are assumed to be known so the procedures estimate these online using maximum likelihood estimation. In order to fairly compare each procedure, the algorithms all start from the same 100 randomly-chosen initial training sets with the same random seed. At the conclusion of the process, Algorithm 2 demonstrates a 29%, 31%, and 35% improvement in average mean absolute error (MAE) over the PDF mean, PDF variance, and random sampling approaches. Additionally, the performance of the algorithms in each of the 100 test cases can be directly compared since they all start with the same L and random seed. At this level, Figure 2 illustrates that Algorithm 2 will either eventually match or outperform the existing sampling strategies nearly 100% of the time. Although the exact numbers will change for different distributions, these results highlight the value of the selection criteria from (9) to further reduce prediction error given a fixed number of samples. Figure 3 demonstrates the utility of CDF variance to identify regions of large prediction errors, regardless of the particular sampling strategy. As discussed in Section III-B, the true CDF variance is unavailable, but the approximate CDF variance (7) still provides a meaningful metric to compare prediction confidence across Θ d . One of the best uses of (7) is to rank points in Θ d according to their CDF variance in order to identify which predictions to trust the least. Figure 3 displays the reduction in prediction error when the points with the top 5% of CDF variance are removed. The 12-22% improvement in MAE proves that the CDF variance did indeed correctly identify and remove regions with large p sat (θ). When the top 10% is removed, the improvement jumps to 40%.

B. Example 2: Robust Multi-Agent Task Allocation
The second example addresses the robust multi-agent task allocation problem from [25] with added stochasticity. The task allocation problem attempts to assign a team of four UAVs to complete fire surveillance tasks that will take longer or shorter depending on the wind speed θ 1 and direction θ 2 . Task durations are also corrupted by zero-mean Gaussian noise. Unforeseen time delays will compound and may potentially lead the UAVs to miss the completion of tasks within their assigned window, thus lowering the overall realized mission score. The verification goal is to determine whether the team of UAV agents will sufficiently complete the ordered tasking and achieve a minimum mission score at different wind settings. Sampling grid Θ d spans the set    Figure 6 demonstrates a similar reduction in MAE as witnessed in Figure 3. Once the points with the top 5% of CDF variance are removed, the MAE of the predictions in the remaining 95% of the data reduces by up to 12-16%. If the top 10% are removed, the average reduction reaches at least 28% for all the sampling strategies. This supports the conclusion that CDF variance identifies points with low prediction confidence where p sat (θ) may be large.

C. Example 3: Lateral-Directional Autopilot
Lastly, the third example considers the altitude-hold requirement of a lateral-directional autopilot [26]. In this problem, a Dryden wind turbulence model [27] augments the original nonlinear aircraft dynamics model to introduce stochasticity into the closed-loop dynamics. The performance requirement expects the autopilot to maintain altitude x(t) within a set threshold of the initial altitude when the autopilot is engaged. This work explores a relaxed version of the original requirement [26] in which the aircraft must remain within a 55 foot window around the initial altitude x(0), The performance measurement is the STL robustness degree ρ ϕ . This example tests the satisfaction of the requirement against initial Euler angles for roll  Figures 3 and 6 is skipped but displays the same exact behavior. This example begins with 100 simulations in the initial training set L and compares the performance over 120 random initializations. As before, Algorithm 2 outperforms the PDF variance-focused and random sampling procedures with a 27% reduction in average MAE. Unlike the previous two examples, Algorithm 2 only slightly edges out the competing mean-focused strategy. This new result is due to the comparatively low y for noisy y(θ) from the low-altitude turbulence model. The resulting tighter distribution only experiences changes in p sat (θ) near y(θ) = 0, which favors the mean-focused selection criteria. If the turbulence is increased then the distribution widens dramatically and the performance of the PDF mean-focused metric degrades significantly like the previous two examples. These case studies all highlight the improved performance of the novel CDF variance algorithms compared to the existing sampling strategies. Just as important, Algorithms 1 and 2 will remain the best-performing sampling strategy even as the distribution changes, whereas the performance of the other strategies will vary according to changes in y .

VI. CONCLUSION
This work introduced machine learning methods for simulation-based statistical verification of stochastic nonlinear systems. In particular, the paper presented a GP-based statistical verification framework to efficiently predict the probability of requirement satisfaction over the full space of possible parametric uncertainties given a limited amount of simulation data. Additionally, Section III-B developed new criteria based upon the variance of the cumulative distribution function to qualify confidence in the predictions. This metric provides a simple validation tool for online identification of regions where the prediction confidence is low. Section IV builds upon this metric and introduces sequential and batch sampling algorithms for efficient closedloop verification. These new verification procedures demonstrate up to a 35% improvement in prediction error over competing approaches in the three examples in Section V. The examples also serve to highlight the utility of the CDF variance to correctly identify low-confidence regions in the parameter space without external validation datasets. While the paper only considers Gaussian distributions for the performance measurements, this work lays the foundation for more complex statistical verification frameworks capable of handing a wider range of distributions. Upcoming work will adapt the CDF variance metric and closed-loop verification procedures to recent developments in modeling of spatially-varying standard deviations [12] and non-Gaussian distributions [13]. Additionally, this work can be extended to high-dimensional GP representations [28]. While the implementation considers Gaussian distributions, the fundamental concepts have broader utility.