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>
97 lines
2.9 KiB
C++
97 lines
2.9 KiB
C++
#include "datascience/stats_summary.h"
|
|
|
|
#include <algorithm>
|
|
#include <cmath>
|
|
#include <cstring>
|
|
#include <vector>
|
|
|
|
namespace fn::ds {
|
|
|
|
double stats_sum(const double* data, std::size_t n) {
|
|
if (n == 0 || data == nullptr) return 0.0;
|
|
// Kahan summation — coste despreciable, evita drift en sumas grandes.
|
|
double s = 0.0, c = 0.0;
|
|
for (std::size_t i = 0; i < n; ++i) {
|
|
double y = data[i] - c;
|
|
double t = s + y;
|
|
c = (t - s) - y;
|
|
s = t;
|
|
}
|
|
return s;
|
|
}
|
|
|
|
double stats_mean(const double* data, std::size_t n) {
|
|
if (n == 0) return 0.0;
|
|
return stats_sum(data, n) / static_cast<double>(n);
|
|
}
|
|
|
|
double stats_min(const double* data, std::size_t n) {
|
|
if (n == 0 || data == nullptr) return 0.0;
|
|
double m = data[0];
|
|
for (std::size_t i = 1; i < n; ++i) if (data[i] < m) m = data[i];
|
|
return m;
|
|
}
|
|
|
|
double stats_max(const double* data, std::size_t n) {
|
|
if (n == 0 || data == nullptr) return 0.0;
|
|
double m = data[0];
|
|
for (std::size_t i = 1; i < n; ++i) if (data[i] > m) m = data[i];
|
|
return m;
|
|
}
|
|
|
|
double stats_variance(const double* data, std::size_t n, bool sample) {
|
|
if (n == 0 || data == nullptr) return 0.0;
|
|
if (sample && n < 2) return 0.0;
|
|
// Welford one-pass.
|
|
double mean = 0.0;
|
|
double M2 = 0.0;
|
|
for (std::size_t i = 0; i < n; ++i) {
|
|
double x = data[i];
|
|
double delta = x - mean;
|
|
mean += delta / static_cast<double>(i + 1);
|
|
double delta2 = x - mean;
|
|
M2 += delta * delta2;
|
|
}
|
|
double denom = sample ? static_cast<double>(n - 1) : static_cast<double>(n);
|
|
return M2 / denom;
|
|
}
|
|
|
|
double stats_std(const double* data, std::size_t n, bool sample) {
|
|
return std::sqrt(stats_variance(data, n, sample));
|
|
}
|
|
|
|
void stats_sort(const double* data, std::size_t n, double* out) {
|
|
if (n == 0 || out == nullptr) return;
|
|
if (out != data && data != nullptr) {
|
|
std::memcpy(out, data, n * sizeof(double));
|
|
}
|
|
std::sort(out, out + n);
|
|
}
|
|
|
|
double stats_quantile_sorted(const double* sorted, std::size_t n, double p) {
|
|
if (n == 0 || sorted == nullptr) return 0.0;
|
|
if (p <= 0.0) return sorted[0];
|
|
if (p >= 1.0) return sorted[n - 1];
|
|
// R type-7: h = (n-1) * p; result = sorted[floor(h)] + (h - floor(h)) *
|
|
// (sorted[floor(h)+1] - sorted[floor(h)])
|
|
double h = (static_cast<double>(n) - 1.0) * p;
|
|
std::size_t lo = static_cast<std::size_t>(std::floor(h));
|
|
std::size_t hi = lo + 1;
|
|
if (hi >= n) hi = n - 1;
|
|
double frac = h - static_cast<double>(lo);
|
|
return sorted[lo] + frac * (sorted[hi] - sorted[lo]);
|
|
}
|
|
|
|
double stats_quantile(const double* data, std::size_t n, double p) {
|
|
if (n == 0) return 0.0;
|
|
std::vector<double> tmp(n);
|
|
stats_sort(data, n, tmp.data());
|
|
return stats_quantile_sorted(tmp.data(), n, p);
|
|
}
|
|
|
|
double stats_percentile(const double* data, std::size_t n, double pct) {
|
|
return stats_quantile(data, n, pct * 0.01);
|
|
}
|
|
|
|
} // namespace fn::ds
|