1// OLS / WLS / ridge with stable linear algebra (Householder QR by default, SVD for rank-deficient designs).2#pragma once3#include <vector>4#include <span>5#include <optional>6#include <Eigen/Dense>7#include "raid/core/types.hpp"8#include "raid/core/estimate.hpp"9#include "raid/statistics/covariance.hpp"1011namespace raid::est {1213enum class Solver { qr, svd, normal_cholesky /* reference only: ill-conditioned designs will be flagged */ };1415struct OlsOptions {16 stats::CovOptions cov{};17 Solver solver{Solver::qr};18 f64 ridge_lambda{0.0}; // L2 penalty on all coefficients except those in `unpenalized` (e.g. the intercept)19 std::vector<usize> unpenalized{}; // column indices not penalised20 std::span<const f64> weights{}; // WLS weights w_i (variance ∝ 1/w_i); empty = OLS21 f64 rcond_warn{1e-10}; // warn when the design's reciprocal condition number is below this22 std::vector<std::string> names{};23};2425struct OlsResult {26 Estimate estimate; // coefficients, inference, loss = RSS27 Eigen::VectorXd residuals;28 Eigen::VectorXd fitted;29 Eigen::MatrixXd xtx_inv; // (X'WX + λP)^{-1}30 f64 sigma2{0}; // RSS/(n−k)31 f64 r2{0}, r2_adj{0};32 Eigen::Index rank{0};33 f64 rcond{0};34 Eigen::VectorXd leverage; // h_ii35};3637OlsResult ols(const Eigen::Ref<const Eigen::MatrixXd>& X, const Eigen::Ref<const Eigen::VectorXd>& y, const OlsOptions& opt = {});3839/// Sufficient statistics of a (possibly partial) regression block: X'X, X'y, y'y, n. Batched/GPU accumulation40/// produces these and the solve happens in FP64 on the CPU.41struct SufficientStats {42 Eigen::MatrixXd xtx; Eigen::VectorXd xty; f64 yty{0}; usize n{0};43 void merge(const SufficientStats& o) { if (n == 0) { *this = o; return; } xtx += o.xtx; xty += o.xty; yty += o.yty; n += o.n; }44};45SufficientStats accumulate(const Eigen::Ref<const Eigen::MatrixXd>& X, const Eigen::Ref<const Eigen::VectorXd>& y);46/// Solve from sufficient statistics (Cholesky with symmetric pivoting fallback to eigen-decomposition); classical SE only.47OlsResult ols_from_stats(const SufficientStats& s, const OlsOptions& opt = {});4849/// Rolling OLS: windows [i−w+1, i] for i ≥ w−1; returns k × n_windows coefficients (column j = window ending at index w−1+j).50struct RollingOlsResult { Eigen::MatrixXd coefficients; std::vector<f64> rss; std::vector<Eigen::Index> ranks; usize window; };51RollingOlsResult rolling_ols(const Eigen::Ref<const Eigen::MatrixXd>& X, const Eigen::Ref<const Eigen::VectorXd>& y, usize window, usize step = 1);5253} // namespace raid::est5455namespace raid::est {56/// Two-sided p-value of a t statistic (df ≤ 0 → standard normal). FP64, validated against scipy in tests/reference.57f64 t_pvalue_two_sided(f64 t, f64 df);58}59