1#include <catch2/catch_test_macros.hpp>2#include <catch2/matchers/catch_matchers_floating_point.hpp>3#include "raid/statistics/bootstrap.hpp"4#include "raid/statistics/reductions.hpp"5#include "raid/core/rng.hpp"6#include <cmath>7#include <set>8using namespace raid; using namespace raid::stats; using Catch::Matchers::WithinAbs;910TEST_CASE("Bootstrap indices: valid, deterministic per (seed, replica), block structure", "[bootstrap]") {11 BootstrapOptions o; o.scheme = BootstrapScheme::circular_block; o.mean_block = 5; auto a = bootstrap_indices(100, o, 3), b = bootstrap_indices(100, o, 3), c = bootstrap_indices(100, o, 4);12 CHECK(a == b); CHECK(a != c); for (u32 i : a) CHECK(i < 100);13 usize runs = 0; for (usize i = 1; i < a.size(); ++i) if (a[i] != (a[i - 1] + 1) % 100) ++runs; CHECK(runs <= 21); // ≤ 20 blocks + wrap14 o.scheme = BootstrapScheme::stationary; o.mean_block = 10; auto s = bootstrap_indices(10000, o, 0); usize starts = 1; for (usize i = 1; i < s.size(); ++i) if (s[i] != (s[i - 1] + 1) % 10000) ++starts;15 CHECK_THAT(10000.0 / static_cast<double>(starts), WithinAbs(10.0, 2.5)); // mean block length ≈ 1016}17TEST_CASE("Block bootstrap SE of the mean of an AR(1) exceeds the iid SE", "[bootstrap]") {18 Rng r(77); const usize n = 2000; std::vector<f64> x(n); x[0] = 0; for (usize t = 1; t < n; ++t) x[t] = 0.7 * x[t - 1] + r.normal();19 auto m = [](std::span<const f64> v) { return mean(v); };20 BootstrapOptions bo; bo.replicas = 499; bo.scheme = BootstrapScheme::stationary; bo.mean_block = 20; auto blk = bootstrap(x, m, bo);21 bo.scheme = BootstrapScheme::iid; auto iid = bootstrap(x, m, bo);22 const double true_se = std::sqrt((1.0 / (1 - 0.49)) / (1 - 0.7) / (1 - 0.7) / n); // long-run variance of AR(1): σ²/(1−φ)² with σ²_x = 1/(1−φ²)23 CHECK(blk.se > 1.6 * iid.se); CHECK_THAT(blk.se, WithinAbs(true_se, 0.35 * true_se)); CHECK(blk.failures == 0); CHECK(blk.ci_lo < blk.point); CHECK(blk.ci_hi > blk.point);24}25TEST_CASE("bootstrap_rows resamples (y,X) jointly and is reproducible", "[bootstrap]") {26 BootstrapOptions o; o.replicas = 50; auto st = [](std::span<const u32> idx) { double s = 0; for (u32 i : idx) s += i; return std::vector<f64>{s / static_cast<double>(idx.size()), static_cast<double>(idx[0])}; };27 auto a = bootstrap_rows(500, st, o), b = bootstrap_rows(500, st, o); CHECK(a.statistics == b.statistics); CHECK(a.point[0] == 249.5); CHECK(a.se[0] > 0);28}29TEST_CASE("bootstrap_rv on iid returns brackets the true integrated variance", "[bootstrap]") {30 Rng r(5); std::vector<f64> ret(390); for (auto& v : ret) v = 1e-3 * r.normal(); BootstrapOptions o; o.replicas = 399; o.mean_block = 5; auto b = bootstrap_rv(ret, o);31 const double iv = 390 * 1e-6; CHECK(b.ci_lo < iv); CHECK(b.ci_hi > iv); CHECK_THAT(b.point, WithinAbs(sum_squares(std::span<const f64>(ret)), 1e-18));32}33