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>
66 lines
2.8 KiB
Markdown
66 lines
2.8 KiB
Markdown
---
|
|
name: rhat_ess
|
|
kind: function
|
|
lang: cpp
|
|
domain: datascience
|
|
version: "1.0.0"
|
|
purity: pure
|
|
signature: "double rhat(const double* chains, size_t m, size_t n); double rhat_split(const double* chains, size_t m, size_t n); double ess_basic(const double* chains, size_t m, size_t n, size_t max_lag, double cutoff)"
|
|
description: "Diagnosticos multi-chain MCMC: Gelman-Rubin R-hat (clasico y split), y Effective Sample Size basico. Cadenas en layout row-major chains[j*n + i]. Convergencia tipica R_hat < 1.01."
|
|
tags: [mcmc, rhat, ess, gelman_rubin, convergence, datascience]
|
|
uses_functions: ["autocorr_cpp_datascience"]
|
|
uses_types: []
|
|
returns: []
|
|
returns_optional: false
|
|
error_type: ""
|
|
imports: [cstddef, cmath, vector]
|
|
tested: false
|
|
tests: []
|
|
test_file_path: ""
|
|
file_path: "cpp/functions/datascience/rhat_ess.cpp"
|
|
params:
|
|
- name: chains
|
|
desc: "Buffer row-major chains[j * n + i] con m cadenas de n samples cada una. Layout estable de cualquier MCMC sampler que produzca arrays apilados."
|
|
- name: m
|
|
desc: "Numero de cadenas (>=2 para rhat clasico)."
|
|
- name: n
|
|
desc: "Samples por cadena (>=2 clasico, >=4 split)."
|
|
- name: max_lag
|
|
desc: "(ess_basic) lag maximo de la ACF para tau_int. Default 200."
|
|
- name: cutoff
|
|
desc: "(ess_basic) umbral |r(k)| para truncar la suma. Default 0.05."
|
|
output: "rhat / rhat_split: escalar >= 1; valores < 1.01 indican convergencia razonable. ess_basic: suma de ESS por cadena (no corregido por between-chain variance)."
|
|
---
|
|
|
|
# rhat_ess
|
|
|
|
Diagnosticos estandar de convergencia para MCMC, esenciales en `mcmc-lab` (que precisamente computa R-hat sobre multi-chain).
|
|
|
|
## R-hat clasico vs split
|
|
|
|
`rhat` es la formula original de Gelman-Rubin (1992). `rhat_split` (Stan / pymc moderno) parte cada cadena en dos mitades antes del calculo — detecta cadenas que parecen estables pero estan stuck en modos distintos al inicio vs fin. **Recomendado: usar rhat_split por defecto.**
|
|
|
|
## Patron tipico
|
|
|
|
```cpp
|
|
constexpr int M = 8, N = 10000;
|
|
std::vector<double> chains(M * N);
|
|
// ... rellenar con tu sampler ...
|
|
|
|
double r = fn::ds::rhat_split(chains.data(), M, N);
|
|
double ess = fn::ds::ess_basic(chains.data(), M, N);
|
|
|
|
if (r > 1.01) {
|
|
// No convergido — correr mas iteraciones o reseed.
|
|
}
|
|
```
|
|
|
|
## Layout
|
|
|
|
`chains[j * n + i]` = sample `i` de la cadena `j`. Asi cada cadena es un slice contiguo, accesible con `chains + j * n`. Esto es lo que produce naturalmente un sampler que persiste un SSBO de samples (cada chain ocupa un bloque contiguo).
|
|
|
|
## Notas
|
|
|
|
- ESS basico no usa la formula multi-chain de Stan (que combina B y W). Es conservador (subestima ESS cuando las cadenas estan bien mezcladas) pero suficiente como primer indicador.
|
|
- Para inferencia formal en produccion, considerar pymc/arviz que tienen las versiones bias-corrected y rank-normalized R-hat. Aqui buscamos diagnostico interactivo, no certificacion estadistica.
|