egutierrez
d115d8e830
Nuevo dominio cpp/functions/datascience con primitivas puras CPU para post- proceso de samples Monte Carlo y diagnostico de cadenas MCMC. Diseñadas como gemelas CPU de los kernels GPU (rng pareja con gpu_rng_glsl, MH 1D/ND con mc_metropolis_hastings_gpu) para validar numericamente y para datasets pequeños donde el dispatch GPU no compensa. - rng: xoshiro256++ con uniform / normal (Box-Muller) / below (Lemire) / categorical. Determinista bit-exacto dado seed. - stats_summary: sum (Kahan), mean, var/std (Welford one-pass), min, max, quantile / percentile (R type-7). - autocorr: r(k), ACF, tau_int (Sokal) — diagnostico ACF y ESS. - rhat_ess: Gelman-Rubin clasico y split + ESS basico (multi-chain). - beta_dist: lgamma (Lanczos), beta_pdf, beta_cdf (continued fraction), beta_quantile, mean/var/std — para inferencia Beta-Binomial. - drawdown: max_dd absoluto/pct + underwater series para sesiones simuladas y backtests. - samples_to_grid_2d: binning 2D CPU para alimentar heatmap_cpp_viz / contour_cpp_viz desde samples (x[], y[]). - metropolis_hastings: MH 1D y ND con target log-pdf como std::function (no normalizada). Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com>
89 lines
2.5 KiB
C++
89 lines
2.5 KiB
C++
#include "datascience/autocorr.h"
|
|
|
|
#include <cmath>
|
|
|
|
namespace fn::ds {
|
|
|
|
static double mean_inline(const double* x, std::size_t n) {
|
|
if (n == 0) return 0.0;
|
|
double s = 0.0;
|
|
for (std::size_t i = 0; i < n; ++i) s += x[i];
|
|
return s / static_cast<double>(n);
|
|
}
|
|
|
|
static double var_inline(const double* x, std::size_t n, double mu) {
|
|
if (n == 0) return 0.0;
|
|
double s = 0.0;
|
|
for (std::size_t i = 0; i < n; ++i) {
|
|
double d = x[i] - mu;
|
|
s += d * d;
|
|
}
|
|
return s / static_cast<double>(n);
|
|
}
|
|
|
|
double autocorr_lag(const double* x, std::size_t n, std::size_t k) {
|
|
if (x == nullptr || n <= 1 || k >= n) return 0.0;
|
|
double mu = mean_inline(x, n);
|
|
double var = var_inline(x, n, mu);
|
|
if (var <= 0.0) return (k == 0 ? 0.0 : 0.0);
|
|
|
|
double cov = 0.0;
|
|
std::size_t m = n - k;
|
|
for (std::size_t i = 0; i < m; ++i) {
|
|
cov += (x[i] - mu) * (x[i + k] - mu);
|
|
}
|
|
cov /= static_cast<double>(m);
|
|
return cov / var;
|
|
}
|
|
|
|
void autocorr_acf(const double* x, std::size_t n,
|
|
std::size_t max_lag, double* out) {
|
|
if (out == nullptr || max_lag == 0) return;
|
|
if (x == nullptr || n <= 1) {
|
|
for (std::size_t k = 0; k < max_lag; ++k) out[k] = 0.0;
|
|
return;
|
|
}
|
|
double mu = mean_inline(x, n);
|
|
double var = var_inline(x, n, mu);
|
|
if (var <= 0.0) {
|
|
for (std::size_t k = 0; k < max_lag; ++k) out[k] = 0.0;
|
|
return;
|
|
}
|
|
for (std::size_t k = 0; k < max_lag; ++k) {
|
|
if (k >= n) { out[k] = 0.0; continue; }
|
|
double cov = 0.0;
|
|
std::size_t m = n - k;
|
|
for (std::size_t i = 0; i < m; ++i) {
|
|
cov += (x[i] - mu) * (x[i + k] - mu);
|
|
}
|
|
cov /= static_cast<double>(m);
|
|
out[k] = cov / var;
|
|
}
|
|
}
|
|
|
|
double autocorr_tau(const double* x, std::size_t n,
|
|
std::size_t max_lag, double cutoff) {
|
|
if (x == nullptr || n <= 1) return 1.0;
|
|
if (max_lag == 0) max_lag = 1;
|
|
double mu = mean_inline(x, n);
|
|
double var = var_inline(x, n, mu);
|
|
if (var <= 0.0) return 1.0;
|
|
|
|
double sum = 0.0;
|
|
std::size_t kmax = (max_lag < n ? max_lag : n - 1);
|
|
for (std::size_t k = 1; k < kmax; ++k) {
|
|
double cov = 0.0;
|
|
std::size_t m = n - k;
|
|
for (std::size_t i = 0; i < m; ++i) {
|
|
cov += (x[i] - mu) * (x[i + k] - mu);
|
|
}
|
|
cov /= static_cast<double>(m);
|
|
double r = cov / var;
|
|
if (std::fabs(r) < cutoff) break;
|
|
sum += r;
|
|
}
|
|
return 1.0 + 2.0 * sum;
|
|
}
|
|
|
|
} // namespace fn::ds
|