1#include "raid/statistics/rolling.hpp"2#include "raid/statistics/reductions.hpp"3#include "raid/core/error.hpp"4#include <cmath>56namespace raid::stats {78namespace {9struct RunningSum {10 f64 s{0}, c{0};11 void add(f64 v) { const f64 t = s + v; if (std::abs(s) >= std::abs(v)) c += (s - t) + v; else c += (v - t) + s; s = t; }12 f64 value() const { return s + c; }13 void reset(f64 exact) { s = exact; c = 0; }14};15inline usize minp(const RollingOptions& o) { return o.min_periods == 0 ? o.window : o.min_periods; }16} // namespace1718std::vector<f64> rolling_sum(std::span<const f64> x, const RollingOptions& o) {19 RAID_REQUIRE(o.window >= 1, invalid_argument, "window must be >= 1");20 const usize n = x.size(), w = o.window, mp = minp(o);21 std::vector<f64> out(n, kNaN);22 RunningSum rs; usize since_resum = 0;23 for (usize i = 0; i < n; ++i) {24 rs.add(x[i]);25 if (i >= w) rs.add(-x[i - w]);26 if (++since_resum >= o.resum_every) { const usize b = i + 1 >= w ? i + 1 - w : 0; rs.reset(sum_neumaier(x.subspan(b, i + 1 - b))); since_resum = 0; }27 const usize count = i + 1 < w ? i + 1 : w;28 if (count >= mp) out[i] = rs.value();29 }30 return out;31}3233std::vector<f64> rolling_mean(std::span<const f64> x, const RollingOptions& o) {34 auto s = rolling_sum(x, o);35 for (usize i = 0; i < s.size(); ++i) { const usize count = i + 1 < o.window ? i + 1 : o.window; s[i] /= static_cast<f64>(count); }36 return s;37}3839std::vector<f64> rolling_variance(std::span<const f64> x, const RollingOptions& o, int ddof) {40 // Two accumulators (sum, sum of squares about a per-window shift) would still cancel for large means;41 // we recompute each window by Welford from scratch when the window is small (<= 64), otherwise use the42 // shifted-sums formula with periodic exact re-summation. Both are FP64 references; speed is not the goal here.43 const usize n = x.size(), w = o.window, mp = minp(o);44 std::vector<f64> out(n, kNaN);45 if (w <= 64) {46 for (usize i = 0; i < n; ++i) {47 const usize b = i + 1 >= w ? i + 1 - w : 0; const usize count = i + 1 - b;48 if (count < mp || count <= static_cast<usize>(ddof)) continue;49 Welford wf; for (usize j = b; j <= i; ++j) wf.add(x[j]);50 out[i] = wf.variance(ddof);51 }52 return out;53 }54 const f64 shift_c = x.empty() ? 0.0 : x[0];55 RunningSum s1, s2; usize since = 0;56 for (usize i = 0; i < n; ++i) {57 const f64 d = x[i] - shift_c; s1.add(d); s2.add(d * d);58 if (i >= w) { const f64 e = x[i - w] - shift_c; s1.add(-e); s2.add(-e * e); }59 if (++since >= o.resum_every) {60 const usize b = i + 1 >= w ? i + 1 - w : 0; f64 a = 0, q = 0;61 std::vector<f64> t1(i + 1 - b), t2(i + 1 - b);62 for (usize j = b; j <= i; ++j) { t1[j - b] = x[j] - shift_c; t2[j - b] = t1[j - b] * t1[j - b]; }63 a = sum_neumaier(std::span<const f64>(t1)); q = sum_neumaier(std::span<const f64>(t2)); s1.reset(a); s2.reset(q); since = 0;64 }65 const usize count = i + 1 < w ? i + 1 : w;66 if (count >= mp && count > static_cast<usize>(ddof)) {67 const f64 cnt = static_cast<f64>(count), m = s1.value() / cnt;68 out[i] = std::max(0.0, (s2.value() - cnt * m * m) / (cnt - static_cast<f64>(ddof)));69 }70 }71 return out;72}7374std::vector<f64> rolling_covariance(std::span<const f64> x, std::span<const f64> y, const RollingOptions& o, int ddof) {75 RAID_REQUIRE(x.size() == y.size(), invalid_argument, "rolling_covariance: size mismatch");76 const usize n = x.size(), w = o.window, mp = minp(o);77 std::vector<f64> out(n, kNaN);78 for (usize i = 0; i < n; ++i) {79 const usize b = i + 1 >= w ? i + 1 - w : 0; const usize count = i + 1 - b;80 if (count < mp || count <= static_cast<usize>(ddof)) continue;81 f64 mx = sum_neumaier(x.subspan(b, count)) / static_cast<f64>(count), my = sum_neumaier(y.subspan(b, count)) / static_cast<f64>(count);82 f64 s = 0, c = 0;83 for (usize j = b; j <= i; ++j) { const f64 v = (x[j] - mx) * (y[j] - my), t = s + v; if (std::abs(s) >= std::abs(v)) c += (s - t) + v; else c += (v - t) + s; s = t; }84 out[i] = (s + c) / static_cast<f64>(count - static_cast<usize>(ddof));85 }86 return out;87}8889std::vector<f64> rolling_correlation(std::span<const f64> x, std::span<const f64> y, const RollingOptions& o) {90 auto cxy = rolling_covariance(x, y, o, 1); auto vx = rolling_variance(x, o, 1); auto vy = rolling_variance(y, o, 1);91 for (usize i = 0; i < cxy.size(); ++i) cxy[i] = (vx[i] > 0 && vy[i] > 0) ? cxy[i] / std::sqrt(vx[i] * vy[i]) : kNaN;92 return cxy;93}9495std::vector<f64> expanding_mean(std::span<const f64> x, usize min_periods) {96 std::vector<f64> out(x.size(), kNaN); Welford w;97 for (usize i = 0; i < x.size(); ++i) { w.add(x[i]); if (w.n >= min_periods) out[i] = w.mean; }98 return out;99}100std::vector<f64> expanding_variance(std::span<const f64> x, usize min_periods, int ddof) {101 std::vector<f64> out(x.size(), kNaN); Welford w;102 for (usize i = 0; i < x.size(); ++i) { w.add(x[i]); if (w.n >= min_periods) out[i] = w.variance(ddof); }103 return out;104}105std::vector<f64> shift(std::span<const f64> x, long k) {106 std::vector<f64> out(x.size(), kNaN);107 for (usize i = 0; i < x.size(); ++i) { const long j = static_cast<long>(i) - k; if (j >= 0 && j < static_cast<long>(x.size())) out[i] = x[static_cast<usize>(j)]; }108 return out;109}110111} // namespace raid::stats112