In Chapter 19 Lesson 1 and Lesson 2, we defined metrics for localization, SLAM, and navigation robustness. In practice, those metrics are only meaningful if they are measured across a controlled range of operating difficulty. Stress testing systematically perturbs the robot’s sensing, actuation, environment, and compute conditions. Ablation studies then isolate which subsystem is responsible for performance gains (or failures).

For a mission episode, let the measured metric vector be \( \mathbf{m} = [m_1,\dots,m_q]^\top \), for example: success indicator, localization RMSE, minimum collision margin, path efficiency, and time-to-completion. A scalar mission score can be introduced for ranking experiments:

\[ J = w_s m_{\text{succ} } - w_e m_{\text{err} } + w_c m_{\text{clr} } - w_p m_{\text{path} } - w_t m_{\text{time} }, \quad w_i \ge 0 \]

The scalar score \( J \) is convenient for optimization and ablation ranking, while the full metric vector \( \mathbf{m} \) remains necessary for diagnosis.

We define a stress vector \( \boldsymbol{\xi} \) that parameterizes disturbances:

\[ \boldsymbol{\xi} = \begin{bmatrix} \xi_{\text{slip} } & \xi_{\text{drop} } & \xi_{\text{bias} } & \xi_{\text{dyn} } & \xi_{\text{delay} } \end{bmatrix}^{\!\top} \]

A normalized stress level \( \rho \in [0,1] \) is often used for campaign design, where \( \rho=0 \) denotes nominal conditions and \( \rho=1 \) denotes the maximum tested disturbance envelope.

\( \rho \): scalar stress index used to schedule severity sweeps. \( P_\rho(\boldsymbol{\xi}) \): stress-conditioned perturbation distribution. \( \omega \): scenario random seed (map layout, obstacle arrivals, etc.).

Each episode is generated by first sampling \( \omega \) and then sampling perturbations from \( P_\rho \). The performance metric is therefore a random variable \( J(\rho,\omega) \).

\[ J(\rho,\omega) = \mathcal{E}\!\left(\pi,\; \mathcal{M},\; \boldsymbol{\xi}\sim P_\rho,\; \omega\right) \]

where \( \pi \) is the AMR stack (localization, mapping, navigation pipeline) and \( \mathcal{M} \) denotes the environment generator or dataset replay.

flowchart TD
A["Define metrics and thresholds"] --> B["Choose stress variables"]
B --> C["Set severity grid rho = 0..1"]
C --> D["Sample seeds and scenarios"]
D --> E["Run full stack under perturbations"]
E --> F["Collect episode metrics"]
F --> G["Estimate reliability and risk"]
G --> H["Escalate / refine hard regions"]
        

A stress campaign should report more than just mean score. We define reliability, expected performance, and tail-risk measures.

\( R(\rho) \): success probability at stress level \( \rho \):

\[ R(\rho) = \mathbb{P}\!\left(\text{mission success} \mid \rho \right) \]

\( \mu_J(\rho) \): stress-conditioned expected mission score:

\[ \mu_J(\rho) = \mathbb{E}[J(\rho,\omega)] \]

\( L(\rho) \): loss random variable \( L = J_{\max} - J \). Tail risk can be summarized by Conditional Value-at-Risk:

\[ \operatorname{CVaR}_{\alpha}(L \mid \rho) = \mathbb{E}\!\left[L \mid L \ge q_{\alpha}(L)\right], \quad \alpha \in (0,1) \]

For binary success over \( n_\rho \) episodes at stress level \( \rho \), with \( k_\rho \) successful runs, the estimator is \( \hat{R}(\rho)=k_\rho/n_\rho \). A Wilson confidence interval is numerically stable even when success is near 0 or 1:

\[ \hat{R}_W(\rho) = \frac{\hat{R} + \frac{z^2}{2n_\rho} }{1+\frac{z^2}{n_\rho} }, \qquad \Delta_W(\rho) = \frac{z}{1+\frac{z^2}{n_\rho} } \sqrt{\frac{\hat{R}(1-\hat{R})}{n_\rho} + \frac{z^2}{4n_\rho^2} } \]

The reported interval is \( [\hat{R}_W(\rho)-\Delta_W(\rho),\; \hat{R}_W(\rho)+\Delta_W(\rho)] \).

Theorem (monotone failure under ordered stress). Suppose the event “failure” is monotone with respect to stress severity and \( P_{\rho_2} \) first-order stochastically dominates \( P_{\rho_1} \) in difficulty for \( \rho_2 \ge \rho_1 \). Then \( R(\rho_2) \le R(\rho_1) \).

Proof sketch. Let \( F(\boldsymbol{\xi}) \in \{0,1\} \) be the failure indicator. By monotonicity, \( F \) is non-decreasing in disturbance severity. Stochastic dominance implies \( \mathbb{E}_{P_{\rho_2} }[F] \ge \mathbb{E}_{P_{\rho_1} }[F] \). Since \( R(\rho)=1-\mathbb{E}[F] \), the result follows.

An ablation removes or degrades a module while holding everything else fixed. For AMR, typical ablations include disabling scan matching, loop closure, IMU fusion, dynamic obstacle prediction, or recovery behaviors.

Let \( A \) denote an ablated system variant and \( \text{Full} \) the complete stack. The estimand is the expected degradation:

\[ \Delta_A(\rho) = \mathbb{E}\!\left[J_{\text{Full} }(\rho,\omega) - J_A(\rho,\omega)\right] \]

To reduce variance, evaluate both systems on the same scenario seed \( \omega_i \) (common random numbers). Define paired differences:

\[ D_i(\rho) = J_{\text{Full} }(\rho,\omega_i) - J_A(\rho,\omega_i), \qquad \hat{\Delta}_A(\rho)=\frac{1}{n}\sum_{i=1}^{n} D_i(\rho) \]

Proposition. \( \hat{\Delta}_A(\rho) \) is unbiased for \( \Delta_A(\rho) \).

Proof. By linearity of expectation,

\[ \mathbb{E}[\hat{\Delta}_A(\rho)] = \frac{1}{n}\sum_{i=1}^{n} \mathbb{E}[D_i(\rho)] = \mathbb{E}[D_1(\rho)] = \Delta_A(\rho) \]

when the pairs are identically distributed.

The variance benefit of pairing is explicit:

\[ \operatorname{Var}(D_i) = \operatorname{Var}(J_{\text{Full} }) + \operatorname{Var}(J_A) - 2\operatorname{Cov}(J_{\text{Full} }, J_A) \]

If the same seed induces positive covariance between the two runs, then pairing reduces variance and narrows confidence intervals.

One-at-a-time ablations are informative but can miss interactions. For three binary module toggles (off/on), use coded variables \( z_j \in \{-1,+1\} \) and fit a factorial model:

\[ Y = \beta_0 + \beta_1 z_1 + \beta_2 z_2 + \beta_3 z_3 + \beta_{12} z_1 z_2 + \beta_{13} z_1 z_3 + \beta_{23} z_2 z_3 + \varepsilon \]

In a balanced \(2^3\) design with coded levels, columns are orthogonal. Therefore, the main effect of module 1 is:

\[ \widehat{\text{Effect} }_1 = 2 \cdot \frac{1}{N}\sum_{i=1}^{N} y_i z_{1,i} \]

and similarly for interactions (replace \( z_{1,i} \) by \( z_{1,i}z_{2,i} \), etc.).

Orthogonality proof idea. In a balanced coded design, each column has equally many \(+1\) and \(-1\) entries, and every pair of columns agrees/disagrees equally often. Hence \( \sum_i z_{a,i}z_{b,i}=0 \) for distinct terms, so the least-squares normal equations decouple term-wise.

Ablation impact usually grows with stress. This is captured by an interaction model across stress levels:

\[ Y_{i a r} = \mu + \alpha_a + \beta_i + (\alpha\beta)_{ia} + u_r + \varepsilon_{iar} \]

where \( \alpha_a \) is ablation effect, \( \beta_i \) is stress-level effect, \( (\alpha\beta)_{ia} \) is the stress-by-ablation interaction, and \( u_r \) is a random seed/scenario effect (paired or repeated campaign blocking).

If the interaction term is significant, the ranking of modules is stress-dependent, which is common in AMR (e.g., IMU fusion may be marginal under nominal conditions but critical under wheel slip).

For multiple ablations, control family-wise error with Holm’s method. Let \( p_{(1)} \le \cdots \le p_{(m)} \) be sorted p-values. Reject hypotheses sequentially while

\[ p_{(k)} \le \frac{\alpha}{m-k+1} \]

and stop at the first failure. This is stricter than naive pairwise testing and avoids inflated false discoveries.

A reproducible campaign should pre-register metrics, thresholds, seeds, stress levels, and ablations. Use the same map seeds across all ablations (paired design), and report both aggregate metrics and per-scenario failure cases.

• Fix evaluation metrics and mission success criteria before running experiments.
• Use a stress grid (e.g., \( \rho=0,0.2,\dots,0.8 \)) and justify the upper bound.
• Use common random seeds for all variants to reduce estimator variance.
• Log raw episode outputs (not only means) for CI estimation and postmortem analysis.
• Report stress-by-ablation interactions, not only average ablation impact.

flowchart TD
A["Full stack baseline"] --> B["Choose ablations (1-factor or factorial)"]
B --> C["Use same seeds for all variants"]
C --> D["Run across stress grid"]
D --> E["Compute paired deltas and CI"]
E --> F["Test stress-by-ablation interactions"]
F --> G["Rank critical modules"]
G --> H["Design mitigations and retest"]
        

The Python script below implements a compact synthetic AMR stress harness with paired ablations and bootstrap confidence intervals for \( \Delta_A(\rho) \).

Code: Chapter19_Lesson4.py

# Chapter19_Lesson4.py
# AMR Stress Testing + Paired Ablation (compact reference)

from dataclasses import dataclass
import math
import numpy as np

@dataclass
class Episode:
    stress: float
    success: int
    rmse: float
    margin: float
    path_ratio: float
    ttc: float

WEIGHTS = {"scan": 0.9, "imu": 0.8, "predict": 0.7}

def simulate(seed: int, stress: float, mods: dict) -> Episode:
    rng = np.random.default_rng(seed)
    dyn = rng.poisson(2 + 6 * stress)
    dropout = rng.binomial(1, min(0.85, 0.05 + 0.75 * stress))
    slip = abs(rng.normal(0.0, 0.02 + 0.12 * stress))
    bias = abs(rng.normal(0.0, 0.01 + 0.08 * stress))
    missing = sum(WEIGHTS[k] * (1 - mods[k]) for k in WEIGHTS)

    rmse = 0.05 + 0.09 * stress + 0.04 * dyn + 1.4 * slip + 1.6 * bias + 0.10 * dropout + 0.12 * missing + rng.normal(0, 0.02)
    rmse = max(0.01, rmse)

    margin = 0.55 - 0.18 * stress - 0.01 * dyn - 0.07 * dropout - 0.05 * missing + 0.02 * mods["predict"] + rng.normal(0, 0.03)
    path_ratio = 1.03 + 0.14 * stress + 0.02 * dyn + 0.08 * dropout + 0.05 * missing - 0.03 * mods["scan"] + rng.normal(0, 0.02)
    path_ratio = max(1.0, path_ratio)
    ttc = 34 + 7 * path_ratio + 2.0 * dyn + 8.0 * stress + 3.5 * missing + rng.normal(0, 1.5)

    logit = 3.2 - 2.6 * stress - 2.2 * rmse + 1.0 * margin - 0.5 * (path_ratio - 1) - 0.7 * missing
    p = 1.0 / (1.0 + math.exp(-logit))
    success = int(rng.uniform() < p)
    return Episode(stress, success, rmse, margin, path_ratio, max(5.0, ttc))

def score(ep: Episode) -> float:
    y = 100.0
    y -= 35.0 * (1 - ep.success)
    y -= 16.0 * max(0.0, ep.rmse - 0.10)
    y -= 14.0 * max(0.0, 0.25 - ep.margin)
    y -= 5.0 * max(0.0, ep.path_ratio - 1.10)
    y -= 0.20 * max(0.0, ep.ttc - 45.0)
    return max(0.0, y)

def bootstrap_ci(diffs: np.ndarray, n_boot: int = 2000, alpha: float = 0.05):
    rng = np.random.default_rng(7)
    means = []
    n = len(diffs)
    for _ in range(n_boot):
        idx = rng.integers(0, n, size=n)
        means.append(float(np.mean(diffs[idx])))
    lo = np.quantile(means, alpha / 2)
    hi = np.quantile(means, 1 - alpha / 2)
    return float(lo), float(hi)

def run():
    stress_grid = [0.0, 0.2, 0.4, 0.6, 0.8]
    seeds = list(range(100, 150))

    full = {"scan": 1, "imu": 1, "predict": 1}
    ablations = {
        "Full": full,
        "NoScan": {"scan": 0, "imu": 1, "predict": 1},
        "NoIMU": {"scan": 1, "imu": 0, "predict": 1},
        "NoPredict": {"scan": 1, "imu": 1, "predict": 0},
    }

    records = {}
    for s in stress_grid:
        for name, mods in ablations.items():
            vals = []
            for seed in seeds:
                ep = simulate(seed, s, mods)
                vals.append((seed, ep, score(ep)))
            records[(s, name)] = vals

    print("stress, ablation, success_rate, mean_score, mean_rmse")
    for s in stress_grid:
        for name in ablations:
            vals = records[(s, name)]
            succ = np.mean([v[1].success for v in vals])
            mean_score = np.mean([v[2] for v in vals])
            mean_rmse = np.mean([v[1].rmse for v in vals])
            print(f"{s:.1f}, {name:9s}, {succ:.3f}, {mean_score:.2f}, {mean_rmse:.3f}")

    print("\nPaired deltas (Full - Ablated) with bootstrap CI")
    for s in stress_grid:
        full_by_seed = {seed: sc for seed, _, sc in records[(s, "Full")]}
        for name in ["NoScan", "NoIMU", "NoPredict"]:
            diffs = np.array([full_by_seed[seed] - sc for seed, _, sc in records[(s, name)]], dtype=float)
            lo, hi = bootstrap_ci(diffs)
            print(f"stress={s:.1f}, {name:9s}, delta={diffs.mean():6.2f}, CI95=({lo:6.2f},{hi:6.2f})")

if __name__ == "__main__":
    run()
      

This C++17 version provides a lightweight stress-test harness suitable for embedding into simulator regression tests or ROS-based benchmarks.

Code: Chapter19_Lesson4.cpp

// Chapter19_Lesson4.cpp
// Compact AMR stress-test harness (C++17)

#include <cmath>
#include <iomanip>
#include <iostream>
#include <map>
#include <random>
#include <string>
#include <vector>

struct Episode {
    int success;
    double rmse;
    double margin;
    double score;
};

Episode simulate(int seed, double stress, const std::map<std::string, int>& mods) {
    std::mt19937 gen(seed);
    std::normal_distribution<double> n01(0.0, 1.0);
    std::uniform_real_distribution<double> u01(0.0, 1.0);
    auto norm = [&](double mu, double sigma) { return mu + sigma * n01(gen); };

    std::poisson_distribution<int> pd(2 + static_cast<int>(6 * stress));
    int dyn = pd(gen);
    int dropout = (u01(gen) < std::min(0.85, 0.05 + 0.75 * stress)) ? 1 : 0;

    double slip = std::fabs(norm(0.0, 0.02 + 0.12 * stress));
    double bias = std::fabs(norm(0.0, 0.01 + 0.08 * stress));
    double missing = 0.9 * (1 - mods.at("scan")) + 0.8 * (1 - mods.at("imu")) + 0.7 * (1 - mods.at("predict"));

    double rmse = 0.05 + 0.09 * stress + 0.04 * dyn + 1.4 * slip + 1.6 * bias + 0.10 * dropout + 0.12 * missing + norm(0.0, 0.02);
    if (rmse < 0.01) rmse = 0.01;

    double margin = 0.55 - 0.18 * stress - 0.01 * dyn - 0.07 * dropout - 0.05 * missing + 0.02 * mods.at("predict") + norm(0.0, 0.03);
    double pathRatio = 1.03 + 0.14 * stress + 0.02 * dyn + 0.08 * dropout + 0.05 * missing - 0.03 * mods.at("scan") + norm(0.0, 0.02);
    if (pathRatio < 1.0) pathRatio = 1.0;
    double ttc = 34 + 7 * pathRatio + 2.0 * dyn + 8.0 * stress + 3.5 * missing + norm(0.0, 1.5);
    if (ttc < 5.0) ttc = 5.0;

    double logit = 3.2 - 2.6 * stress - 2.2 * rmse + 1.0 * margin - 0.5 * (pathRatio - 1.0) - 0.7 * missing;
    double p = 1.0 / (1.0 + std::exp(-logit));
    int success = (u01(gen) < p) ? 1 : 0;

    double score = 100.0;
    score -= 35.0 * (1 - success);
    score -= 16.0 * std::max(0.0, rmse - 0.10);
    score -= 14.0 * std::max(0.0, 0.25 - margin);
    score -= 5.0 * std::max(0.0, pathRatio - 1.10);
    score -= 0.20 * std::max(0.0, ttc - 45.0);
    if (score < 0.0) score = 0.0;

    return {success, rmse, margin, score};
}

int main() {
    std::vector<double> stressGrid = {0.0, 0.4, 0.8};
    std::map<std::string, int> full = { {"scan", 1}, {"imu", 1}, {"predict", 1} };
    std::map<std::string, int> noScan = { {"scan", 0}, {"imu", 1}, {"predict", 1} };

    std::cout << std::fixed << std::setprecision(3);
    for (double s : stressGrid) {
        double succFull = 0.0, scoreFull = 0.0, rmseFull = 0.0;
        double deltaSum = 0.0;
        int n = 0;
        for (int seed = 100; seed < 160; ++seed) {
            Episode a = simulate(seed, s, full);
            Episode b = simulate(seed, s, noScan); // same seed -> paired
            succFull += a.success;
            scoreFull += a.score;
            rmseFull += a.rmse;
            deltaSum += (a.score - b.score);
            n++;
        }
        std::cout << "stress=" << s
                  << " fullSucc=" << succFull / n
                  << " fullScore=" << scoreFull / n
                  << " fullRMSE=" << rmseFull / n
                  << " pairedDelta(Full-NoScan)=" << deltaSum / n
                  << "\n";
    }
    return 0;
}
      

The Java code estimates main and interaction effects in a balanced \(2^3\) factorial ablation design. This is useful for backend evaluation services or educational JVM toolchains.

Code: Chapter19_Lesson4.java

// Chapter19_Lesson4.java
// 2^3 factorial ablation effect estimation (Java)

import java.util.ArrayList;
import java.util.List;
import java.util.Locale;
import java.util.Random;

public class Chapter19_Lesson4 {

    static class Run {
        int z1, z2, z3; // coded levels: -1 off, +1 on
        double y;
        Run(int z1, int z2, int z3, double y) { this.z1 = z1; this.z2 = z2; this.z3 = z3; this.y = y; }
    }

    static double simulateScore(int z1, int z2, int z3, double stress, int seed) {
        Random rng = new Random(seed);
        double noise = rng.nextGaussian();

        // Main effects + interactions + stress interaction
        double y = 80.0
                + 4.5 * z1 + 3.8 * z2 + 2.9 * z3
                + 1.5 * z1 * z2 + 1.0 * z1 * z3 + 0.6 * z2 * z3
                - 11.0 * stress - 1.8 * stress * z1 - 1.2 * stress * z3
                + 1.2 * noise;

        if (y < 0.0) y = 0.0;
        if (y > 100.0) y = 100.0;
        return y;
    }

    static List<Run> designRuns(double stress, int reps) {
        int[] levels = {-1, 1};
        List<Run> out = new ArrayList<>();
        int seed = 1000;
        for (int z1 : levels) {
            for (int z2 : levels) {
                for (int z3 : levels) {
                    for (int r = 0; r < reps; r++) {
                        out.add(new Run(z1, z2, z3, simulateScore(z1, z2, z3, stress, seed++)));
                    }
                }
            }
        }
        return out;
    }

    static double effect(List<Run> runs, String term) {
        double acc = 0.0;
        for (Run run : runs) {
            int x;
            switch (term) {
                case "z1": x = run.z1; break;
                case "z2": x = run.z2; break;
                case "z3": x = run.z3; break;
                case "z1z2": x = run.z1 * run.z2; break;
                case "z1z3": x = run.z1 * run.z3; break;
                case "z2z3": x = run.z2 * run.z3; break;
                default: throw new IllegalArgumentException("unknown term");
            }
            acc += run.y * x;
        }
        // For balanced coded 2^k design: effect = 2 * mean(y*x)
        return 2.0 * acc / runs.size();
    }

    public static void main(String[] args) {
        double[] stressGrid = {0.0, 0.4, 0.8};
        int reps = 16;

        for (double s : stressGrid) {
            List<Run> runs = designRuns(s, reps);
            double mean = 0.0;
            for (Run r : runs) mean += r.y;
            mean /= runs.size();

            System.out.println("=== stress = " + s + " ===");
            System.out.printf(Locale.US, "mean score: %.3f%n", mean);
            System.out.printf(Locale.US, "effect(z1=ScanMatch): %.3f%n", effect(runs, "z1"));
            System.out.printf(Locale.US, "effect(z2=IMUFusion): %.3f%n", effect(runs, "z2"));
            System.out.printf(Locale.US, "effect(z3=Predictor): %.3f%n", effect(runs, "z3"));
            System.out.printf(Locale.US, "effect(z1z2): %.3f%n", effect(runs, "z1z2"));
            System.out.printf(Locale.US, "effect(z1z3): %.3f%n", effect(runs, "z1z3"));
            System.out.printf(Locale.US, "effect(z2z3): %.3f%n", effect(runs, "z2z3"));
            System.out.println();
        }
    }
}
      

The MATLAB script performs paired ablation sweeps and plots a reliability curve. The same logic can be wrapped in a Simulink test harness (e.g., scenario blocks + logged outputs) for hardware-in-the-loop experiments.

Code: Chapter19_Lesson4.m

% Chapter19_Lesson4.m
% AMR stress testing + paired ablation (MATLAB / Simulink support script)

clear; clc; rng(42);

stressGrid = [0.0 0.2 0.4 0.6 0.8];
seeds = 100:149;

modsFull = [1 1 1];      % [scan imu predict]
modsNoScan = [0 1 1];
modsNoIMU = [1 0 1];

fprintf('stress   pSucc(full)   score(full)   pairedDelta(NoScan)   pairedDelta(NoIMU)\n');

for s = stressGrid
    fullScores = zeros(size(seeds));
    noScanScores = zeros(size(seeds));
    noIMUScores = zeros(size(seeds));
    fullSucc = zeros(size(seeds));

    for k = 1:numel(seeds)
        sd = seeds(k);
        eFull = simulateEpisode(sd, s, modsFull);
        eNoScan = simulateEpisode(sd, s, modsNoScan); % paired seed
        eNoIMU = simulateEpisode(sd, s, modsNoIMU);   % paired seed

        fullScores(k) = missionScore(eFull);
        noScanScores(k) = missionScore(eNoScan);
        noIMUScores(k) = missionScore(eNoIMU);
        fullSucc(k) = eFull.success;
    end

    fprintf('%4.1f      %7.3f       %8.2f          %8.2f            %8.2f\n', ...
        s, mean(fullSucc), mean(fullScores), mean(fullScores - noScanScores), mean(fullScores - noIMUScores));
end

% Reliability curve plot
pSucc = zeros(size(stressGrid));
for i = 1:numel(stressGrid)
    s = stressGrid(i);
    tmp = zeros(size(seeds));
    for k = 1:numel(seeds)
        tmp(k) = simulateEpisode(seeds(k), s, modsFull).success;
    end
    pSucc(i) = mean(tmp);
end
figure; plot(stressGrid, pSucc, '-o', 'LineWidth', 1.5); grid on;
xlabel('Stress level'); ylabel('Success probability');
title('AMR reliability curve under stress');

function ep = simulateEpisode(seed, stress, mods)
    rng(seed, 'twister');
    dyn = poissrnd(2 + 6*stress);
    dropout = rand < min(0.85, 0.05 + 0.75*stress);
    slip = abs(normrnd(0, 0.02 + 0.12*stress));
    bias = abs(normrnd(0, 0.01 + 0.08*stress));
    missing = 0.9*(1-mods(1)) + 0.8*(1-mods(2)) + 0.7*(1-mods(3));

    rmse = 0.05 + 0.09*stress + 0.04*dyn + 1.4*slip + 1.6*bias + 0.10*dropout + 0.12*missing + normrnd(0,0.02);
    rmse = max(0.01, rmse);
    margin = 0.55 - 0.18*stress - 0.01*dyn - 0.07*dropout - 0.05*missing + 0.02*mods(3) + normrnd(0,0.03);
    pathRatio = 1.03 + 0.14*stress + 0.02*dyn + 0.08*dropout + 0.05*missing - 0.03*mods(1) + normrnd(0,0.02);
    pathRatio = max(1.0, pathRatio);
    ttc = max(5.0, 34 + 7*pathRatio + 2.0*dyn + 8.0*stress + 3.5*missing + normrnd(0,1.5));

    logitVal = 3.2 - 2.6*stress - 2.2*rmse + 1.0*margin - 0.5*(pathRatio - 1) - 0.7*missing;
    p = 1/(1 + exp(-logitVal));
    ep = struct('success', rand < p, 'rmse', rmse, 'margin', margin, 'pathRatio', pathRatio, 'ttc', ttc);
end

function y = missionScore(ep)
    y = 100 - 35*(1-ep.success) ...
          - 16*max(0, ep.rmse - 0.10) ...
          - 14*max(0, 0.25 - ep.margin) ...
          - 5*max(0, ep.pathRatio - 1.10) ...
          - 0.20*max(0, ep.ttc - 45.0);
    y = max(0, y);
end
      

The Mathematica notebook script below symbolically verifies paired estimator moments and performs a compact synthetic reliability/ablation analysis.

Code: Chapter19_Lesson4.nb

(* Chapter19_Lesson4.nb *)
(* Wolfram Mathematica: paired ablation estimator + symbolic moments *)

ClearAll["Global`*"];

(* Symbolic result for paired estimator mean/variance *)
n = Symbol["n"];
mu = Symbol["mu"];
sigma = Symbol["sigma"];

meanExpr = FullSimplify[
  Expectation[(1/n) Sum[Subscript[D, i], {i, 1, n}],
    Table[Subscript[D, i] \[Distributed] NormalDistribution[mu, sigma], {i, 1, n}]
  ],
  Assumptions -> {n > 0, Element[n, Integers], sigma > 0}
];

varExpr = FullSimplify[
  Variance[(1/n) Sum[Subscript[D, i], {i, 1, n}],
    Table[Subscript[D, i] \[Distributed] NormalDistribution[mu, sigma], {i, 1, n}]
  ],
  Assumptions -> {n > 0, Element[n, Integers], sigma > 0}
];

Print["E[mean(D)] = ", meanExpr];
Print["Var[mean(D)] = ", varExpr];

(* Synthetic reliability / ablation study *)
SeedRandom[2026];
stressGrid = {0.0, 0.2, 0.4, 0.6, 0.8};
seeds = Range[100, 149];

simulateScore[seed_, stress_, missing_] := Module[
  {dyn, dropout, slip, bias, rmse, margin, pathRatio, ttc, logit, p, success, y},
  SeedRandom[seed];
  dyn = RandomVariate[PoissonDistribution[2 + 6 stress]];
  dropout = Boole[RandomReal[] < Min[0.85, 0.05 + 0.75 stress]];
  slip = Abs[RandomVariate[NormalDistribution[0, 0.02 + 0.12 stress]]];
  bias = Abs[RandomVariate[NormalDistribution[0, 0.01 + 0.08 stress]]];

  rmse = Max[0.01, 0.05 + 0.09 stress + 0.04 dyn + 1.4 slip + 1.6 bias + 0.10 dropout + 0.12 missing + RandomVariate[NormalDistribution[0, 0.02]]];
  margin = 0.55 - 0.18 stress - 0.01 dyn - 0.07 dropout - 0.05 missing + RandomVariate[NormalDistribution[0, 0.03]];
  pathRatio = Max[1.0, 1.03 + 0.14 stress + 0.02 dyn + 0.08 dropout + 0.05 missing + RandomVariate[NormalDistribution[0, 0.02]]];
  ttc = Max[5.0, 34 + 7 pathRatio + 2 dyn + 8 stress + 3.5 missing + RandomVariate[NormalDistribution[0, 1.5]]];

  logit = 3.2 - 2.6 stress - 2.2 rmse + 1.0 margin - 0.5 (pathRatio - 1) - 0.7 missing;
  p = 1/(1 + Exp[-logit]);
  success = Boole[RandomReal[] < p];

  y = 100 - 35 (1 - success) - 16 Max[0, rmse - 0.10] - 14 Max[0, 0.25 - margin] - 5 Max[0, pathRatio - 1.10] - 0.20 Max[0, ttc - 45.0];
  Max[0, y]
];

pairedTable = Table[
  Module[{full, noScan, d},
    full = Table[simulateScore[sd, s, 0.0], {sd, seeds}];
    noScan = Table[simulateScore[sd, s, 0.9], {sd, seeds}];
    d = full - noScan;
    {s, Mean[d], StandardDeviation[d], Mean[UnitStep[full - 55]]}
  ],
  {s, stressGrid}
];

Print["{stress, meanDelta, sdDelta, successRate}:"]; Print /@ pairedTable;

ListLinePlot[
  pairedTable[[All, {1, 4}]],
  Frame -> True,
  PlotMarkers -> Automatic,
  FrameLabel -> {"Stress", "Success Probability"},
  PlotLabel -> "Reliability Curve (Full Stack)"
]
      

Problem 1 (Unbiased paired ablation estimator): Let \( D_i = J_{\text{Full} }(\rho,\omega_i) - J_A(\rho,\omega_i) \) for paired seeds \( \omega_i \). Prove that \( \hat{\Delta}_A(\rho)=\frac{1}{n}\sum_i D_i \) is unbiased.

Solution: Assuming the pairs are i.i.d. under the fixed stress level \( \rho \),

\[ \mathbb{E}[\hat{\Delta}_A(\rho)] = \frac{1}{n}\sum_{i=1}^{n}\mathbb{E}[D_i] = \mathbb{E}[D_1] = \Delta_A(\rho) \]

Hence the sample mean of paired differences is an unbiased estimator of the expected ablation degradation.

Problem 2 (Variance reduction by pairing): Show why using common seeds reduces variance compared with independent runs of Full and Ablated stacks.

Solution: For a paired difference,

\[ \operatorname{Var}(D) = \operatorname{Var}(J_{\text{Full} }) + \operatorname{Var}(J_A) - 2\operatorname{Cov}(J_{\text{Full} },J_A) \]

Common seeds induce positive covariance because both systems see the same map complexity and obstacle randomness. Therefore the covariance term is positive and subtracts from total variance, making the estimator more precise.

Problem 3 (Factorial effect in coded design): In a balanced \(2^3\) coded design with \( z_1,z_2,z_3 \in \{-1,+1\} \), derive the estimator for the main effect of factor 1.

Solution: Orthogonality implies the least-squares coefficient satisfies

\[ \hat{\beta}_1 = \frac{1}{N}\sum_{i=1}^{N} y_i z_{1,i} \]

and the usual “effect” scale (difference between average high and low levels) is

\[ \widehat{\text{Effect} }_1 = 2\hat{\beta}_1 = 2\cdot \frac{1}{N}\sum_{i=1}^{N} y_i z_{1,i} \]

Problem 4 (Wilson interval for reliability): At stress level \( \rho=0.6 \), suppose \( n=50 \) and \( k=32 \) successful episodes. Write the Wilson center and radius for a 95% interval.

Solution: First compute \( \hat{R}=32/50=0.64 \) and use \( z=1.96 \).

\[ \hat{R}_W = \frac{0.64 + \frac{1.96^2}{2\cdot 50} }{1+\frac{1.96^2}{50} }, \qquad \Delta_W = \frac{1.96}{1+\frac{1.96^2}{50} } \sqrt{\frac{0.64(1-0.64)}{50} + \frac{1.96^2}{4\cdot 50^2} } \]

The confidence interval is \( [\hat{R}_W-\Delta_W,\; \hat{R}_W+\Delta_W] \).

Problem 5 (Interpreting stress-by-ablation interaction): In the model \( Y_{iar}=\mu+\alpha_a+\beta_i+(\alpha\beta)_{ia}+u_r+\varepsilon_{iar} \), what does a large positive value of \( (\alpha\beta)_{ia} \) mean when the response is a loss metric (larger is worse)?

Solution: It means the ablation \( a \) becomes disproportionately harmful at stress level \( i \). In other words, the module removed in that ablation is especially important under that stress regime, and the interaction reveals a condition-specific vulnerability.

Stress testing quantifies how AMR performance degrades as sensing, dynamics, or environmental conditions worsen; ablation studies quantify which subsystem causes the degradation. The key methodological ideas are: paired seeds (variance reduction), reliability curves with confidence intervals, factorial ablations for interactions, and stress-by-ablation interaction models for condition-specific failure analysis. These tools directly support the full-system evaluation report in the next lesson.