1#include "raid/statistics/reductions.hpp"23namespace raid::stats {45Moments moments(std::span<const f64> x) {6 RAID_REQUIRE(x.size() >= 2, insufficient_samples, "moments need n >= 2");7 const f64 n = static_cast<f64>(x.size());8 const f64 m = sum_neumaier(x) / n;9 std::vector<f64> d2(x.size()), d3(x.size()), d4(x.size());10 f64 mn = x[0], mx = x[0];11 for (usize i = 0; i < x.size(); ++i) {12 const f64 d = x[i] - m; d2[i] = d * d; d3[i] = d2[i] * d; d4[i] = d2[i] * d2[i];13 mn = std::min(mn, x[i]); mx = std::max(mx, x[i]);14 }15 const f64 m2 = sum_neumaier(std::span<const f64>(d2)) / n, m3 = sum_neumaier(std::span<const f64>(d3)) / n, m4 = sum_neumaier(std::span<const f64>(d4)) / n;16 const f64 var_unb = sum_neumaier(std::span<const f64>(d2)) / (n - 1.0);17 return {x.size(), m, var_unb, m2 > 0 ? m3 / std::pow(m2, 1.5) : std::numeric_limits<f64>::quiet_NaN(), m2 > 0 ? m4 / (m2 * m2) - 3.0 : std::numeric_limits<f64>::quiet_NaN(), mn, mx};18}1920f64 autocovariance(std::span<const f64> x, usize k, bool demean) {21 RAID_REQUIRE(k < x.size(), insufficient_samples, "autocovariance lag >= n");22 const f64 n = static_cast<f64>(x.size());23 const f64 m = demean ? sum_neumaier(x) / n : 0.0;24 std::vector<f64> p(x.size() - k);25 for (usize i = 0; i + k < x.size(); ++i) p[i] = (x[i] - m) * (x[i + k] - m);26 return sum_neumaier(std::span<const f64>(p)) / n;27}2829std::vector<f64> autocorrelation(std::span<const f64> x, usize max_lag) {30 const f64 c0 = autocovariance(x, 0);31 std::vector<f64> r(max_lag + 1, 1.0);32 for (usize k = 1; k <= max_lag && k < x.size(); ++k) r[k] = c0 > 0 ? autocovariance(x, k) / c0 : std::numeric_limits<f64>::quiet_NaN();33 return r;34}3536Drawdown max_drawdown(std::span<const f64> cum) {37 Drawdown d{0.0, 0, 0};38 if (cum.empty()) return d;39 f64 peak = cum[0]; usize peak_i = 0;40 for (usize i = 0; i < cum.size(); ++i) {41 if (cum[i] > peak) { peak = cum[i]; peak_i = i; }42 const f64 dd = peak != 0.0 ? (peak - cum[i]) / peak : 0.0;43 if (dd > d.max_drawdown) { d.max_drawdown = dd; d.peak_index = peak_i; d.trough_index = i; }44 }45 return d;46}4748} // namespace raid::stats49