Consider a robotic system executing a task under variable initial states and environments. Each run (or rollout) produces an outcome: success, failure, or partial completion. We model an experiment as a random variable \( Y \in \{0,1\} \), where \( Y=1 \) denotes success and \( Y=0 \) denotes failure under an evaluation distribution \( \mathbb{P} \) over scenarios.

\[ R \;=\; \mathbb{P}(Y=1) \;=\; \mathbb{E}[Y], \quad p_{\mathrm{fail}} \;=\; \mathbb{P}(Y=0) \;=\; 1 - R. \]

A benchmark specifies the scenario distribution and task definition. A failure taxonomy refines the event \( \{Y=0\} \) into interpretable categories; stress testing modifies the scenario distribution to deliberately probe regions of high failure probability or high severity.

flowchart TD
  A["Robot policy (planner + controller)"] --> B["Scenario sampler (benchmark or stress)"]
  B --> C["Execution: simulator or hardware"]
  C --> D["Online monitors (constraints, safety checks)"]
  D --> E["Outcome + logs"]
  E --> F["Failure labeling using taxonomy"]
  F --> G["Metrics (p_fail, severity, CVaR, etc.)"]
  G --> H["Design / update stress tests and benchmarks"]
        

The central idea: metrics live on random variables (outcomes), while taxonomies live on structured labels attached to those outcomes. Stress testing changes how we sample, but the mathematical object of interest remains the same failure indicator.

Let \( Z \) denote a random scenario (task instance, environment, disturbance realizations, etc.), and \( \Omega \) capture any internal randomness (e.g., randomized planner). A policy \( \pi \) induces a random trajectory \( \Gamma(Z,\Omega) \) in state space. A set of task and safety constraints is encoded as a function \( g \) on trajectories such that \( g(\Gamma) \leq 0 \) corresponds to constraint satisfaction.

\[ I_{\mathrm{fail}}(Z,\Omega) \;=\; \mathbf{1}\!\big( g(\Gamma(Z,\Omega)) > 0 \big) \]

is the failure indicator. The induced failure probability under an evaluation distribution \( \mathbb{P} \) over \( (Z,\Omega) \) is

\[ p_{\mathrm{fail}} \;=\; \mathbb{E}_{(Z,\Omega)\sim\mathbb{P}} \left[ I_{\mathrm{fail}}(Z,\Omega) \right] \;=\; \mathbb{P}\big( I_{\mathrm{fail}}(Z,\Omega)=1 \big). \]

If we further define a severity random variable \( C(Z,\Omega) \geq 0 \) (e.g., maximum penetration, minimum distance to collision, cost overruns), then generic risk functionals \( \rho(C) \) such as expectation, variance, or coherent risk measures (Section 4) become the mathematical objects we want to estimate and compare across algorithms.

In practice, \( p_{\mathrm{fail}} \) and associated risk measures are unknown and must be estimated from a finite set of rollouts; Section 5 treats this estimation problem.

A failure taxonomy assigns each failed rollout to a type describing where in the pipeline the failure occurred and what its observable effect was. Mathematically, a taxonomy is a discrete-valued random variable \( F \) taking values in a finite set \( \{1,\dots,K\} \) conditioned on failure:

\[ F \;=\; \phi(Z,\Omega,\Gamma) \quad \text{on} \quad \{ I_{\mathrm{fail}}(Z,\Omega)=1 \}, \]

where \( \phi \) is a labeling rule (possibly implemented by human annotators plus automatic checks).

A common design for robotic systems uses axes such as:

• Stage: perception, planning, control, execution, communication.
• Trigger: environment shift, modeling error, numerical issue, hardware fault.
• Effect: safety violation, task failure without safety issue, performance degradation but task completed.

Let \( F^{\mathrm{stage}} \), \( F^{\mathrm{trigger}} \), \( F^{\mathrm{effect}} \) denote the corresponding categorical random variables on failure runs. Then

\[ p_{s} \;=\; \mathbb{P}\big( F^{\mathrm{stage}} = s \mid I_{\mathrm{fail}} = 1 \big) \]

quantifies how often the pipeline fails at stage \( s \) among all failures. Similarly,

\[ p_{\mathrm{fail}} \;=\; \sum_{k=1}^{K} \mathbb{P}(F=k) \]

if \( \{F=k\}_{k=1}^{K} \) forms a partition of the failure event. Ensuring that taxonomy classes are mutually exclusive and collectively exhaustive on \( \{Y=0\} \) is important: otherwise the probabilities \( \mathbb{P}(F=k) \) cannot be interpreted as a decomposition of \( p_{\mathrm{fail}} \).

To incorporate severity, assign each class a cost \( c_k \geq 0 \) (e.g., near-miss vs collision vs catastrophic damage) and define the expected class-weighted loss

\[ \mathbb{E}[C] \;=\; \sum_{k=1}^{K} c_k \, \mathbb{P}(F=k), \]

which couples taxonomy design directly to risk-sensitive evaluation.

Failure probability is only one aspect of performance. Two policies may have similar \( p_{\mathrm{fail}} \) but very different tail behaviour of their severity random variable \( C \). Stress testing is often motivated by tail risk.

A basic tail-risk quantity is the (upper) quantile, or Value-at-Risk:

\[ \mathrm{VaR}_{\alpha}(C) \;=\; \inf\big\{ \ell \in \mathbb{R} : \mathbb{P}(C \leq \ell) \geq \alpha \big\}, \quad \alpha \in (0,1). \]

The Conditional Value-at-Risk (CVaR) at level \( \alpha \) averages the worst \( 1-\alpha \) tail:

\[ \mathrm{CVaR}_{\alpha}(C) \;=\; \inf_{\eta \in \mathbb{R}} \left\{ \eta + \frac{1}{1-\alpha} \mathbb{E}\big[ (C - \eta)_+ \big] \right\}, \quad (x)_+ = \max(x,0). \]

CVaR is a coherent risk measure; it respects properties like monotonicity, positive homogeneity, and subadditivity, which are desirable when aggregating risk across tasks or robot teams.

In stress-oriented evaluation, we consider either:

• Conditional stress: restrict scenarios to a subset \( S_{\mathrm{stress}} \subseteq \mathcal{Z} \) (e.g., low friction, high latency) and study

  \[ p_{\mathrm{fail}}(S_{\mathrm{stress}}) \;=\; \mathbb{P}\big( I_{\mathrm{fail}} = 1 \mid Z \in S_{\mathrm{stress}} \big), \]

• Distributional stress: replace the nominal scenario distribution \( \mathbb{P}_0 \) by a new distribution \( \mathbb{Q} \) that over-samples near failure regions, and evaluate \( p_{\mathrm{fail}}^{\mathbb{Q}} \) and tail-risk under \( \mathbb{Q} \).

The latter connects to distributionally robust evaluation, where one may define a set of plausible distributions \( \mathcal{Q} \) around \( \mathbb{P}_0 \) and study a worst-case risk:

\[ \rho_{\mathcal{Q}} \;=\; \sup_{\mathbb{Q} \in \mathcal{Q}} \mathbb{E}_{\mathbb{Q}}[C]. \]

While full distributionally robust optimization is beyond the scope of this course, the basic idea underlies many practical stress-testing heuristics: explore regions of scenario space where \( C \) is large or where taxonomy labels correspond to severe failure modes.

Suppose we execute \( N \) independent rollouts under some (possibly stress-biased) scenario distribution, and record failure indicators \( I_{\mathrm{fail}}^{(1)},\dots,I_{\mathrm{fail}}^{(N)} \). The empirical failure probability is

\[ \hat{p}_{\mathrm{fail}} \;=\; \frac{1}{N} \sum_{i=1}^{N} I_{\mathrm{fail}}^{(i)}. \]

Under independence and identical distribution assumptions, Hoeffding's inequality for bounded random variables yields a distribution-free confidence bound:

\[ \mathbb{P}\Big( \big| \hat{p}_{\mathrm{fail}} - p_{\mathrm{fail}} \big| > \varepsilon \Big) \;\leq\; 2 \exp\!\big( -2 N \varepsilon^2 \big), \quad \varepsilon > 0. \]

Solving for \( N \) to guarantee accuracy \( \varepsilon \) with probability at least \( 1-\delta \) gives

\[ N \;\geq\; \frac{1}{2\varepsilon^2} \log\!\left( \frac{2}{\delta} \right). \]

This provides concrete guidance when planning stress tests: for example, to estimate \( p_{\mathrm{fail}} \) within \( \varepsilon = 0.01 \) with confidence \( 1-\delta = 0.95 \), we require on the order of \( N \approx 18{,}500 \) runs.

flowchart TD
  S["Specify scenario distribution and stress subset"] --> N["Choose N, epsilon, delta via Hoeffding bound"]
  N --> G["Generate N scenarios (possibly biased toward stress region)"]
  G --> R["Run rollouts in simulator/hardware"]
  R --> M["Compute p_hat, CVaR, class frequencies"]
  M --> D["Diagnose weaknesses via taxonomy & adjust design"]
        

For rare but critical failures (e.g., \( p_{\mathrm{fail}} \ll 10^{-3} \)), naive Monte Carlo becomes inefficient; variance-reduction and rare-event methods (importance sampling, cross-entropy) are then needed, but the basic empirical estimator and confidence machinery remain the core.

We now implement a simple Monte Carlo stress test in Python. The code focuses on the evaluation logic; connecting it to a robotics simulator such as pybullet or MoveIt is a matter of replacing a stub function with actual simulation calls.

import numpy as np

# Example: stress scenarios over friction, payload mass, and control latency.
# In practice, replace simulate_trial(...) with calls into PyBullet, MoveIt, etc.

rng = np.random.default_rng(seed=42)

def sample_stress_scenarios(N):
    """Draw N stress scenarios.
    friction: biased toward low values using a Beta distribution.
    mass: log-normal (kg).
    latency: exponential tail (s).
    """
    friction = rng.beta(a=0.5, b=2.0, size=N)        # more mass near 0
    mass = rng.lognormal(mean=0.0, sigma=0.25, size=N)
    latency = rng.exponential(scale=0.08, size=N)    # mean 80 ms
    return friction, mass, latency

def simulate_trial(friction, mass, latency):
    """Return (failed, taxonomy_label).

    failed: bool
    taxonomy_label: string representing failure class, e.g.
      "control:overshoot", "planning:infeasible", "sensing:occlusion", ...
    Here we implement a toy algebraic model instead of a real simulator.
    """
    # Constraint 1: stopping distance constraint violated if
    # effective friction * mass is too small.
    stopping_violation = (friction * mass) < 0.6

    # Constraint 2: closed-loop stability degrades if latency is too large.
    latency_violation = latency > 0.15  # seconds

    if stopping_violation and latency_violation:
        return True, "control+environment:multi-factor"
    if stopping_violation:
        return True, "control:insufficient_braking"
    if latency_violation:
        return True, "control:latency_instability"

    # Example of a non-safety failure: task time-out.
    timeout_violation = (mass > 2.5) and (latency > 0.10)
    if timeout_violation:
        return True, "performance:timeout"

    return False, "success"

def run_stress_test(N, delta=0.05):
    friction, mass, latency = sample_stress_scenarios(N)

    failures = 0
    labels = []
    for f, m, L in zip(friction, mass, latency):
        failed, label = simulate_trial(f, m, L)
        if failed:
            failures += 1
            labels.append(label)

    p_hat = failures / N
    # Hoeffding half-width for the two-sided confidence interval
    eps = np.sqrt(np.log(2.0 / delta) / (2.0 * N))
    ci_low = max(0.0, p_hat - eps)
    ci_high = min(1.0, p_hat + eps)

    # Empirical taxonomy distribution on failures
    from collections import Counter
    counts = Counter(labels)
    tax_probs = {k: v / failures for k, v in counts.items()} if failures > 0 else {}

    return {
        "N": N,
        "failures": failures,
        "p_hat": p_hat,
        "ci": (ci_low, ci_high),
        "taxonomy_distribution": tax_probs,
    }

if __name__ == "__main__":
    result = run_stress_test(N=5000, delta=0.05)
    print("N:", result["N"])
    print("Estimated failure probability:", result["p_hat"])
    print("95% Hoeffding CI:", result["ci"])
    print("Failure taxonomy distribution:")
    for k, v in result["taxonomy_distribution"].items():
        print(f"  {k}: {v:.3f}")
      

To connect this to a robotic system, one would replace simulate_trial by code that:

• sets environment parameters in a simulator (e.g., PyBullet),
• executes a motion plan from MoveIt or OMPL,
• runs safety monitors to compute failed, and
• assigns a taxonomy label based on logged data.

In C++, stress tests often integrate with planning libraries such as OMPL and MoveIt. Below is a minimal skeleton that samples stress scenarios and records empirical failure probability and a simple text-based taxonomy. The core loop can be wrapped around calls to an OMPL planner and a dynamics simulator.

#include <random>
#include <iostream>
#include <string>

struct TrialResult {
    bool failed;
    std::string label;
};

// Toy stress model as in the Python example.
TrialResult simulate_trial(double friction, double mass, double latency) {
    bool stopping_violation = (friction * mass) < 0.6;
    bool latency_violation = latency > 0.15;
    bool timeout_violation = (mass > 2.5) && (latency > 0.10);

    if (stopping_violation && latency_violation) {
        return {true, "control+environment:multi-factor"};
    }
    if (stopping_violation) {
        return {true, "control:insufficient_braking"};
    }
    if (latency_violation) {
        return {true, "control:latency_instability"};
    }
    if (timeout_violation) {
        return {true, "performance:timeout"};
    }
    return {false, "success"};
}

int main() {
    const int N = 5000;
    const double delta = 0.05;

    std::mt19937 rng(42);
    std::gamma_distribution<double> gamma_f(0.5, 1.0);   // for friction
    std::lognormal_distribution<double> logn_m(0.0, 0.25);
    std::exponential_distribution<double> exp_lat(1.0 / 0.08);

    int failures = 0;

    for (int i = 0; i < N; ++i) {
        double friction = std::min(1.0, gamma_f(rng));  // truncated
        double mass = logn_m(rng);
        double latency = exp_lat(rng);

        TrialResult tr = simulate_trial(friction, mass, latency);
        if (tr.failed) {
            ++failures;
            // Here we could increment a map<std::string,int> to build taxonomy stats.
        }
    }

    double p_hat = static_cast<double>(failures) / N;
    double eps = std::sqrt(std::log(2.0 / delta) / (2.0 * N));
    double ci_low = std::max(0.0, p_hat - eps);
    double ci_high = std::min(1.0, p_hat + eps);

    std::cout << "N: " << N << "\n";
    std::cout << "Estimated failure probability: " << p_hat << "\n";
    std::cout << "Hoeffding 95% CI: [" << ci_low << ", " << ci_high << "]\n";

    return 0;
}
      

Within a MoveIt-based benchmarking tool, the body of the loop would typically:

• sample world parameters (obstacle poses, friction),
• call an OMPL planner via MoveIt to compute a trajectory,
• simulate execution in Gazebo or a dynamics engine, and
• label failures according to a shared taxonomy definition.

Java is less central to manipulation planning but is used in some robotics stacks (e.g., ROSJava, custom warehouse robots). The example below mirrors the Monte Carlo pattern, ready to be integrated with a simulator accessed via Java APIs.

import java.util.Random;

public class StressTest {

    private static final Random rng = new Random(42L);

    // Sample stress scenarios (simplified distributions).
    static double[] sampleScenario() {
        double friction = Math.max(0.0, Math.min(1.0, rng.nextGaussian() * 0.2 + 0.3));
        double mass = Math.exp(rng.nextGaussian() * 0.25);  // log-normal-ish
        double latency = -0.08 * Math.log(1.0 - rng.nextDouble());
        return new double[]{friction, mass, latency};
    }

    static class TrialResult {
        boolean failed;
        String label;
        TrialResult(boolean f, String l) { failed = f; label = l; }
    }

    static TrialResult simulateTrial(double friction, double mass, double latency) {
        boolean stoppingViolation = (friction * mass) < 0.6;
        boolean latencyViolation = latency > 0.15;
        boolean timeoutViolation = (mass > 2.5) && (latency > 0.10);

        if (stoppingViolation && latencyViolation) {
            return new TrialResult(true, "control+environment:multi-factor");
        }
        if (stoppingViolation) {
            return new TrialResult(true, "control:insufficient_braking");
        }
        if (latencyViolation) {
            return new TrialResult(true, "control:latency_instability");
        }
        if (timeoutViolation) {
            return new TrialResult(true, "performance:timeout");
        }
        return new TrialResult(false, "success");
    }

    public static void main(String[] args) {
        int N = 5000;
        double delta = 0.05;
        int failures = 0;

        for (int i = 0; i < N; ++i) {
            double[] s = sampleScenario();
            TrialResult tr = simulateTrial(s[0], s[1], s[2]);
            if (tr.failed) {
                failures++;
                // In practice, increment a Map<String,Integer> of taxonomy counts.
            }
        }

        double pHat = (double) failures / N;
        double eps = Math.sqrt(Math.log(2.0 / delta) / (2.0 * N));
        double ciLow = Math.max(0.0, pHat - eps);
        double ciHigh = Math.min(1.0, pHat + eps);

        System.out.println("N = " + N);
        System.out.println("Estimated failure probability: " + pHat);
        System.out.println("Hoeffding 95% CI: [" + ciLow + ", " + ciHigh + "]");
    }
}
      

To use a Java-based simulator, replace simulateTrial with calls into the robot API (e.g., via ROSJava services) and evaluate constraint violations and taxonomy labels using its logs.

MATLAB and Simulink are widely used in control and robotics. The script below illustrates a Monte Carlo parameter sweep. In a real setup, simulate_trial would call sim on a Simulink model that includes a robot arm, controller, and environment blocks, possibly using Robotics System Toolbox.

function stress_test_matlab()
    N = 5000;
    delta = 0.05;

    failures = 0;

    for i = 1:N
        [friction, mass, latency] = sample_scenario();
        [failed, label] = simulate_trial(friction, mass, latency); %#ok<NASGU>
        if failed
            failures = failures + 1;
            % Use containers.Map to maintain taxonomy counts if desired.
        end
    end

    p_hat = failures / N;
    eps = sqrt(log(2.0 / delta) / (2.0 * N));
    ci_low = max(0.0, p_hat - eps);
    ci_high = min(1.0, p_hat + eps);

    fprintf('N = %d\n', N);
    fprintf('Estimated failure probability: %.4f\n', p_hat);
    fprintf('Hoeffding 95%% CI: [%.4f, %.4f]\n', ci_low, ci_high);
end

function [friction, mass, latency] = sample_scenario()
    friction = betarnd(0.5, 2.0);   % requires Statistics and Machine Learning Toolbox
    mass = lognrnd(0.0, 0.25);
    latency = exprnd(0.08);
end

function [failed, label] = simulate_trial(friction, mass, latency)
    stopping_violation = (friction * mass) < 0.6;
    latency_violation = latency > 0.15;
    timeout_violation = (mass > 2.5) && (latency > 0.10);

    if stopping_violation && latency_violation
        failed = true;
        label = 'control+environment:multi-factor';
    elseif stopping_violation
        failed = true;
        label = 'control:insufficient_braking';
    elseif latency_violation
        failed = true;
        label = 'control:latency_instability';
    elseif timeout_violation
        failed = true;
        label = 'performance:timeout';
    else
        failed = false;
        label = 'success';
    end

    % In a Simulink workflow, this function would:
    %  - configure model parameters (friction, mass, latency),
    %  - run sim('model') to obtain trajectories,
    %  - compute constraint violations and taxonomy labels from logged signals.
end
      

Simulink's parameter sweep mechanisms (e.g., Simulink.SimulationInput) can automate large stress campaigns, while MATLAB scripts aggregate the results into failure probabilities and taxonomy statistics.

Mathematica is useful for both symbolic manipulations (e.g., deriving sample-complexity bounds) and numeric Monte Carlo experiments. The following code mirrors the previous examples in a functional style.

(* Toy stress scenario distributions *)
sampleScenario[] := Module[
  {friction, mass, latency},
  friction = RandomVariate[BetaDistribution[0.5, 2.0]];
  mass = RandomVariate[LogNormalDistribution[0.0, 0.25]];
  latency = RandomVariate[ExponentialDistribution[1/0.08]];
  {friction, mass, latency}
];

(* Toy failure model *)
simulateTrial[{friction_, mass_, latency_}] := Module[
  {stoppingViolation, latencyViolation, timeoutViolation},
  stoppingViolation = (friction*mass) < 0.6;
  latencyViolation = latency > 0.15;
  timeoutViolation = (mass > 2.5) && (latency > 0.10);

  Which[
    stoppingViolation && latencyViolation, {True, "control+environment:multi-factor"},
    stoppingViolation, {True, "control:insufficient_braking"},
    latencyViolation, {True, "control:latency_instability"},
    timeoutViolation, {True, "performance:timeout"},
    True, {False, "success"}
  ]
];

stressTest[N_Integer?Positive, delta_?Positive] := Module[
  {scenarios, trials, failures, pHat, eps, ciLow, ciHigh},
  scenarios = Table[sampleScenario[], {N}];
  trials = simulateTrial /@ scenarios;

  failures = Count[trials[[All, 1]], True];
  pHat = N[failures / N];

  eps = Sqrt[Log[2.0/delta]/(2.0*N)];
  ciLow = Max[0.0, pHat - eps];
  ciHigh = Min[1.0, pHat + eps];

  << "Return association with stats and taxonomy distribution" >>
  << "In a full implementation, build Counts[trials[[All, 2]]] / failures" >>

  << "For brevity, we return only core statistics here." >>
  <|
    "N" -> N,
    "Failures" -> failures,
    "pHat" -> pHat,
    "CI" -> {ciLow, ciHigh}
  |>
];

(* Example call *)
stressTest[5000, 0.05]
      

Mathematica is particularly convenient for analytically studying bounds such as Hoeffding's inequality, or for prototyping new risk measures before porting them to a real-time robotics stack.

Problem 1 (Taxonomy as a Partition): Let \( I_{\mathrm{fail}} \) be the failure indicator and \( F \in \{1,\dots,K\} \) a taxonomy label that is only defined when \( I_{\mathrm{fail}}=1 \). Suppose the events \( E_k = \{F = k\} \) are pairwise disjoint and their union equals \( \{I_{\mathrm{fail}}=1\} \). Show that \( p_{\mathrm{fail}} = \sum_{k=1}^{K} \mathbb{P}(F=k) \).

Solution: By assumption, the sets \( E_k \) form a partition of the failure event:

\[ \{ I_{\mathrm{fail}} = 1 \} \;=\; \biguplus_{k=1}^{K} E_k, \]

where \( \biguplus \) denotes disjoint union. Using finite additivity of probability measures on disjoint sets,

\[ p_{\mathrm{fail}} \;=\; \mathbb{P}(I_{\mathrm{fail}} = 1) \;=\; \mathbb{P}\!\left( \bigcup_{k=1}^{K} E_k \right) \;=\; \sum_{k=1}^{K} \mathbb{P}(E_k) \;=\; \sum_{k=1}^{K} \mathbb{P}(F=k). \]

Thus the taxonomy frequencies \( \mathbb{P}(F=k) \) decompose the total failure probability.

Problem 2 (Sample Complexity via Hoeffding): Derive the bound \( N \geq \frac{1}{2\varepsilon^2} \log\!\left(\frac{2}{\delta}\right) \) needed to guarantee \( \mathbb{P}(|\hat{p}_{\mathrm{fail}} - p_{\mathrm{fail}}| \leq \varepsilon) \geq 1-\delta \) for a Bernoulli failure indicator.

Solution: Hoeffding's inequality states that for i.i.d. Bernoulli variables \( X_i \in [0,1] \) with mean \( \mu \),

\[ \mathbb{P}\big( |\bar{X}_N - \mu| > \varepsilon \big) \;\leq\; 2 \exp(-2 N \varepsilon^2), \quad \varepsilon > 0. \]

Here \( X_i = I_{\mathrm{fail}}^{(i)} \), \( \bar{X}_N = \hat{p}_{\mathrm{fail}} \), and \( \mu = p_{\mathrm{fail}} \). We require

\[ 2 \exp(-2 N \varepsilon^2) \leq \delta. \]

Taking logarithms and rearranging,

\[ -2 N \varepsilon^2 \leq \log(\delta/2) \quad\Rightarrow\quad N \geq \frac{1}{2\varepsilon^2} \log\!\left( \frac{2}{\delta} \right). \]

Therefore this choice of \( N \) suffices to guarantee the desired confidence.

Problem 3 (Class-Weighted Loss): Consider a taxonomy with three failure classes: collision (class 1), near-miss (class 2), and timeout (class 3). Let costs be \( c_1 = 10 \), \( c_2 = 3 \), \( c_3 = 1 \). You observe empirical class probabilities \( \hat{p}_1 = 0.01 \), \( \hat{p}_2 = 0.03 \), \( \hat{p}_3 = 0.06 \). Compute the estimated class-weighted loss \( \widehat{\mathbb{E}[C]} \) and interpret.

Solution: The estimator of the expected cost is

\[ \widehat{\mathbb{E}[C]} \;=\; \sum_{k=1}^{3} c_k \hat{p}_k \;=\; 10 \cdot 0.01 + 3 \cdot 0.03 + 1 \cdot 0.06. \]

Numerically,

\[ \widehat{\mathbb{E}[C]} \;=\; 0.10 + 0.09 + 0.06 \;=\; 0.25. \]

Thus, on average, the policy incurs an expected normalized loss of 0.25 per rollout under the benchmark distribution, dominated by rare but highly costly collisions (class 1).

Problem 4 (Empirical CVaR for Discrete Outcomes): Let \( C_1,\dots,C_N \) be i.i.d. rollout costs. To estimate \( \mathrm{CVaR}_{\alpha}(C) \) empirically, a standard approach is to: (i) sort the samples, (ii) take the empirical \( \mathrm{VaR}_{\alpha} \), and (iii) average the worst \( (1-\alpha) \)-fraction. Derive the formula for \( \widehat{\mathrm{CVaR}}_{\alpha} \) when \( N(1-\alpha) \) is an integer and provide an explicit expression.

Solution: Let \( C_{(1)} \leq \dots \leq C_{(N)} \) denote the order statistics. If \( m = N(1-\alpha) \) is an integer, the empirical \( \mathrm{VaR}_{\alpha} \) is \( C_{(N-m)} \). The empirical CVaR averages the top \( m \) values:

\[ \widehat{\mathrm{CVaR}}_{\alpha}(C) \;=\; \frac{1}{m} \sum_{i=N-m+1}^{N} C_{(i)}. \]

This estimator is consistent under mild conditions and is widely used in risk-sensitive evaluation of robotic policies under stochastic disturbances.

In this lesson we formalized robotic failures as random variables driven by scenario distributions and policy-induced trajectories. We constructed failure taxonomies as discrete labels that partition the failure event, enabling class-wise analysis and class-weighted loss metrics. We introduced tail-risk measures such as VaR and CVaR and linked them to stress testing, both as conditional evaluation on stress subsets and as distributional shifts toward rare but severe scenarios.

Finally, we derived finite-sample accuracy guarantees for empirical failure probabilities using Hoeffding's inequality and implemented taxonomy-aware stress tests in Python, C++, Java, MATLAB/Simulink, and Mathematica. These tools provide a rigorous foundation for the benchmarking and reproducibility practices developed in the rest of this chapter.