Matemáticas y estadística — Ecuaciones de MatLens

Esta sección resume las ecuaciones que MatLens implementa en sus núcleos de estadística inferencial, estadística multivariante, ajuste no lineal, propagación de incertidumbre, inferencia bayesiana, mapas hiperespectrales y series operando, con las convenciones exactas del código.

Corrección por comparaciones múltiples

Código: matlens/core/inference.py::p_adjust. Mismas definiciones que p.adjust de R. Sean $p_{(1)} \le \dots \le p_{(m)}$ los $m$ valores p ordenados.

\text{Bonferroni: } \tilde p_i = \min(1,\, m\,p_i)

\text{Holm: } \tilde p_{(i)} = \min\!\Big(1,\ \max_{j \le i}\,(m-j+1)\,p_{(j)}\Big)

\text{Hochberg: } \tilde p_{(i)} = \min\!\Big(1,\ \min_{j \ge i}\,(m-j+1)\,p_{(j)}\Big)

\text{BH: } \tilde p_{(i)} = \min\!\Big(1,\ \min_{j \ge i}\,\frac{m}{j}\,p_{(j)}\Big), \qquad \text{BY: igual con el factor } q=\sum_{l=1}^{m}\frac{1}{l}

$p_i$ : valor p sin corregir (adimensional, en $[0,1]$ ); $m$ : número de comparaciones; $\tilde p_i$ : valor p ajustado.
Holm y Bonferroni controlan la tasa de error por familia (FWER); BH controla la FDR bajo independencia o dependencia positiva; BY bajo dependencia arbitraria.
Referencias: Dunn 1961, JASA 56:52; Holm 1979, Scand. J. Stat. 6:65; Hochberg 1988, Biometrika 75:800; Benjamini & Hochberg 1995, JRSS-B 57:289; Benjamini & Yekutieli 2001, Ann. Stat. 29:1165.

Interpretación y límites

El ajuste de valores p no cambia los datos ni la prueba: sólo reescala la regla de decisión. Bonferroni y Holm son válidos bajo cualquier dependencia; Hochberg y BH requieren independencia o dependencia positiva (PRDS); BY es válido siempre a costa de un factor $q=\sum 1/l$ ( $\approx \ln m + 0.577$ ).
Algoritmo exacto: se ordenan los $p$ (descendente para Hochberg/BH/BY), se aplica el factor de cada método, se toma el máximo acumulado (Holm) o el mínimo acumulado (los demás) y se recorta a 1; los valores vuelven a su orden original.
Error: ValueError si la lista está vacía, contiene no finitos o valores fuera de $[0,1]$ , o si el método es desconocido.

Ejemplo resuelto

Ocho valores p (p_adjust):

p	Bonferroni	Holm	Hochberg	BH	BY
0.001	0.008	0.008	0.008	0.008	0.0217
0.008	0.064	0.056	0.056	0.032	0.087
0.039	0.312	0.234	0.148	0.0672	0.183
0.041	0.328	0.234	0.148	0.0672	0.183
0.042	0.336	0.234	0.148	0.0672	0.183
0.06	0.48	0.234	0.148	0.08	0.217
0.074	0.592	0.234	0.148	0.0846	0.23
0.205	1	0.234	0.205	0.205	0.557

Pruebas t (Student y Welch)

Código: matlens/core/stats.py::ttest_welch, ttest_one_sample, ttest_paired (vía scipy.stats.ttest_ind, ttest_1samp, ttest_rel).

t = \frac{\bar y_1-\bar y_2}{\sqrt{s_1^2/n_1+s_2^2/n_2}}, \qquad \nu = \frac{\left(s_1^2/n_1+s_2^2/n_2\right)^2}{\dfrac{(s_1^2/n_1)^2}{n_1-1}+\dfrac{(s_2^2/n_2)^2}{n_2-1}}

$\bar y_i$ , $s_i^2$ (con $n_i-1$ en el denominador), $n_i$ : media, varianza y tamaño del grupo $i$ (unidades de los datos y su cuadrado); $\nu$ : grados de libertad de Welch–Satterthwaite.
Por defecto equal_var=False (Welch). Se informa el intervalo de confianza de $\bar y_1-\bar y_2$ al nivel $1-\alpha$ , la $d$ de Cohen con desviación agrupada y el p de Levene (centrado en la mediana, valor por defecto de SciPy) como nota.
Supuestos: observaciones independientes y aproximadamente normales; al menos 2 datos por grupo.
Referencia: implementación propia sobre SciPy (no hay cita en el código).

ANOVA de una vía

Código: matlens/core/stats.py::anova_oneway.

SS_b=\sum_{i=1}^{k} n_i(\bar y_i-\bar y)^2,\quad SS_w=\sum_i\sum_j (y_{ij}-\bar y_i)^2,\quad F=\frac{SS_b/(k-1)}{SS_w/(N-k)}

$k$ : número de grupos; $N$ : total de observaciones; $p = P(F_{k-1,N-k} > F)$ .
Si $k\ge 3$ y cada grupo tiene al menos 2 datos se añade Tukey HSD. Se informan Levene (sobre los grupos) y Shapiro–Wilk sobre los residuos $y_{ij}-\bar y_i$ .
Referencia: implementación propia.

Supuestos, algoritmo y límites

Modelo: $y_{ij}=\mu_i+\varepsilon_{ij}$ , $\varepsilon_{ij}\sim N(0,\sigma^2)$ independientes y con varianza común. Bajo $H_0$ ( $\mu_1=\dots=\mu_k$ ), $F\sim F_{k-1,N-k}$ . Con varianzas heterogéneas y tamaños desiguales la tasa de error tipo I se distorsiona; use Welch o Games–Howell.
Algoritmo: se eliminan los no finitos de cada grupo; error si hay menos de 2 grupos o si $N-k<1$ (ningún grupo con réplicas). Si $MS_w=0$ , $F=\infty$ . Tukey HSD se calcula con scipy.stats.tukey_hsd (IC al nivel $1-\alpha$ ).
La prueba t de Welch (ttest_welch) usa los grados de libertad de Welch–Satterthwaite $\nu=\dfrac{(s_1^2/n_1+s_2^2/n_2)^2}{\frac{(s_1^2/n_1)^2}{n_1-1}+\frac{(s_2^2/n_2)^2}{n_2-1}}$ e informa además Levene (centro en la mediana, valor por defecto de SciPy) y $d$ de Cohen.

Ejemplo resuelto

Tres grupos de 8 datos normales sintéticos (medias 10, 11, 12.5; DE 1, 1, 1.5; semilla fija).

Fuente	SS	gl	MS	F	p
Entre grupos	17.34	2	8.671	4.486	0.0239
Dentro de grupos	40.6	21	1.933
Total	57.94	23

Welch A vs C (ttest_welch): $t$ = -2.815, $\nu$ = 13.61, $p$ = 0.0141, diferencia = -2.081 [-3.67, -0.4914] (IC 95 %), $d$ de Cohen = -1.408.

Tukey HSD, Games–Howell y Dunnett

Código: matlens/core/inference.py::tukey_hsd, games_howell, dunnett.

Tukey HSD (Tukey–Kramer con $n$ desiguales) se calcula con scipy.stats.tukey_hsd; se informa diferencia, IC al nivel $1-\alpha$ y p. Una advertencia aparece si Brown–Forsythe (Levene con mediana) da $p<0.05$ .

Games–Howell:

q_{ij}=\frac{|\bar y_i-\bar y_j|}{\sqrt{\tfrac12\left(s_i^2/n_i+s_j^2/n_j\right)}},\qquad p = P\big(Q_{k,\nu_{ij}} > q_{ij}\big),\qquad \bar y_i-\bar y_j \pm q_{1-\alpha;k,\nu_{ij}}\sqrt{\tfrac12\left(s_i^2/n_i+s_j^2/n_j\right)}

$Q_{k,\nu}$ : distribución del rango studentizado con $k$ grupos y $\nu_{ij}$ de Welch–Satterthwaite (fórmula anterior con los grupos $i,j$ ).
Advertencia si algún grupo tiene menos de 6 datos.

Dunnett (muchos contra un control) usa scipy.stats.dunnett con integración cuasi-Monte Carlo de semilla fija (rng=0).

Referencias: Tukey 1949; Kramer 1956, Biometrics 12:307; Games & Howell 1976, J. Educ. Stat. 1:113; Dunnett 1955, JASA 50:1096.

ANOVA de medidas repetidas de una vía

Código: matlens/core/inference.py::rm_anova_oneway. Tabla de $n$ sujetos × $k$ condiciones; se descartan filas incompletas.

SS_c=n\sum_j(\bar y_{\cdot j}-\bar y)^2,\ SS_s=k\sum_i(\bar y_{i\cdot}-\bar y)^2,\ SS_e=SS_t-SS_c-SS_s,\ F=\frac{SS_c/(k-1)}{SS_e/[(k-1)(n-1)]}

Esfericidad sobre la covarianza $\mathbf S$ de los contrastes ortonormales ( $k\times(k-1)$ , por QR), con autovalores $\lambda_l$ y $p'=k-1$ :

\hat\varepsilon_{GG}=\frac{(\sum\lambda_l)^2}{p'\sum\lambda_l^2},\qquad \hat\varepsilon_{HF}=\min\!\left(1,\ \frac{n p'\hat\varepsilon_{GG}-2}{p'\,(n-1-p'\hat\varepsilon_{GG})}\right)

W=\frac{\det\mathbf S}{(\operatorname{tr}\mathbf S/p')^{p'}},\quad \chi^2=-(n-1)\left(1-\frac{2p'^2+p'+2}{6p'(n-1)}\right)\ln W,\quad \nu=\tfrac{p'(p'+1)}{2}-1

Los valores p corregidos usan $F$ con grados de libertad multiplicados por $\varepsilon$ . Mauchly solo se calcula si $k>2$ y $n>k-1$ ; si rechaza (p < α) se recomienda GG o HF. Se informa $\eta_p^2 = SS_c/(SS_c+SS_e)$ .
Referencias: Mauchly 1940, Ann. Math. Stat. 11:204; Greenhouse & Geisser 1959, Psychometrika 24:95; Huynh & Feldt 1976, J. Educ. Stat. 1:69.

Modelo lineal mixto

Código: matlens/core/inference.py::mixed_lm. $y=\mathbf X\boldsymbol\beta+u_g+\varepsilon$ , intercepto aleatorio por grupo $u_g\sim N(0,\sigma_u^2)$ , $\varepsilon\sim N(0,\sigma^2)$ , ajustado por REML con statsmodels.MixedLM (dependencia opcional). Advertencia si no converge. Referencia: Lindstrom & Bates 1988, JASA 83:1014.

Tablas de contingencia

Código: matlens/core/inference.py::chi2_independence, cramers_v, fisher_exact, mcnemar.

\chi^2=\sum_{ij}\frac{(O_{ij}-E_{ij})^2}{E_{ij}},\qquad V=\sqrt{\frac{\chi^2}{n\,(\min(r,c)-1)}},\qquad \chi^2_{\text{McNemar}}=\frac{(|b-c|-\delta)^2}{b+c}

$O,E$ : frecuencias observadas y esperadas; corrección de Yates en 2×2 por defecto (correction=True); $V$ se calcula siempre sin corrección. Regla de Cochran: advertencia si algún $E<1$ o más del 20 % de $E<5$ .
McNemar: $b,c$ pares discordantes; $\delta=1$ solo si correction=True (por defecto 0); si $b+c<25$ se usa la prueba binomial exacta con estadístico $\min(b,c)$ .
Fisher exacto: solo 2×2, devuelve razón de probabilidades y p.
Referencias: Pearson 1900; Cramér 1946, §21.9; Cochran 1954, Biometrics 10:417; Fisher 1935; McNemar 1947, Psychometrika 12:153.

Pruebas no paramétricas

Código: matlens/core/inference.py::kruskal_dunn, dunn_posthoc, friedman, wilcoxon_signed_rank.

z_{ij}=\frac{\bar R_i-\bar R_j}{\sqrt{\left(\dfrac{N(N+1)}{12}-\dfrac{\sum (t^3-t)}{12(N-1)}\right)\left(\dfrac1{n_i}+\dfrac1{n_j}\right)}},\qquad p=2\,[1-\Phi(|z|)]

\varepsilon^2_H=\frac{H(N+1)}{N^2-1},\qquad W_{\text{Kendall}}=\frac{\chi^2_F}{n(k-1)},\qquad r_{rb}=\frac{\sum R^+-\sum R^-}{\sum R}

$\bar R_i$ : rango medio del grupo $i$ ; $t$ : tamaños de los empates; los p de Dunn se ajustan con p_adjust (Holm por defecto). Kruskal–Wallis advierte con grupos de menos de 5 datos; Wilcoxon descarta diferencias nulas y advierte con menos de 6 pares.
Referencias: Kruskal & Wallis 1952, JASA 47:583; Dunn 1964, Technometrics 6:241; Friedman 1937, JASA 32:675; Wilcoxon 1945, Biometrics Bull. 1:80; King & Minium 2003.

Supervivencia: Kaplan–Meier y log-rank

Código: matlens/core/inference.py::kaplan_meier, logrank. event = 1 evento, 0 censurado.

\hat S(t)=\prod_{t_i\le t}\left(1-\frac{d_i}{n_i}\right),\qquad \widehat{\mathrm{Var}}\,\hat S(t)=\hat S(t)^2\sum_{t_i\le t}\frac{d_i}{n_i(n_i-d_i)}

\text{IC log(−log): } \exp\!\left[-\exp\!\left(\ln(-\ln\hat S)\pm z\,\frac{\sqrt{\sum d_i/[n_i(n_i-d_i)]}}{|\ln\hat S|}\right)\right]

$t_i$ : tiempos de evento distintos; $n_i$ : en riesgo ( $t\ge t_i$ ); $d_i$ : eventos en $t_i$ ; $z=\Phi^{-1}(0.5+\text{conf}/2)$ , conf = 0.95. Límites no finitos se sustituyen por 0 y 1. Mediana: primer $t_i$ con $\hat S\le 0.5$ .

Log-rank con $k$ grupos: en cada $t_i$ , $E_g \mathrel{+}= d\,n_g/n$ y

V_{gh} \mathrel{+}= \frac{d(n-d)}{n^2(n-1)}\left(n\,n_g\,\delta_{gh}-n_g n_h\right),\qquad \chi^2=(\mathbf O-\mathbf E)^\top_{1..k-1}\,\mathbf V^{+}_{1..k-1}\,(\mathbf O-\mathbf E)_{1..k-1},\ \nu=k-1

$\mathbf V^+$ : pseudoinversa. Advertencia si algún $E_g<5$ .
Referencias: Kaplan & Meier 1958, JASA 53:457; Greenwood 1926; Kalbfleisch & Prentice 2002, §1.4; Mantel 1966, Cancer Chemother. Rep. 50:163; Peto & Peto 1972.

Análisis de potencia

Código: matlens/core/inference.py::power_ttest, power_anova, power_proportions.

\text{t: } 1-\beta = P(t'_{\nu,\delta}>t_c)+P(t'_{\nu,\delta}<-t_c),\quad \delta=d\sqrt{n/2},\ \nu=2n-2\ \text{(dos muestras)};\ \delta=d\sqrt n,\ \nu=n-1\ \text{(una/pareada)}

\text{ANOVA: } 1-\beta=P\big(F'_{k-1,k(n-1),\lambda}>F_{1-\alpha}\big),\ \lambda=f^2kn;\qquad \text{proporciones: } h=2\arcsin\sqrt{p_1}-2\arcsin\sqrt{p_2},\ n=\left\lceil 2\left(\frac{z_{1-\alpha/s}+z_{1-\beta}}{h}\right)^2\right\rceil

$d$ , $f$ , $h$ : tamaños de efecto de Cohen (adimensionales); $n$ : tamaño por grupo; $s=2$ bilateral, 1 unilateral; $t_c=t_{1-\alpha/s,\nu}$ .
Para t y ANOVA, $n$ se obtiene duplicando y luego por bisección entera hasta alcanzar la potencia pedida ( $n\ge 2$ ); para proporciones con la fórmula cerrada normal.
Referencia: Cohen 1988, Statistical Power Analysis, cap. 2, 6 y 8.

Tamaños de efecto

Código: matlens/core/inference.py::cohen_d, hedges_g, eta_omega_squared.

d=\frac{\bar y_1-\bar y_2}{s_p},\ \ \mathrm{SE}(d)=\sqrt{\frac{n_1+n_2}{n_1n_2}+\frac{d^2}{2(n_1+n_2)}};\qquad g=J\,d,\ J=1-\frac{3}{4\,\nu-1}

\eta^2=\frac{SS_b}{SS_t},\qquad \omega^2=\frac{SS_b-(k-1)MS_w}{SS_t+MS_w},\qquad \eta^2_{\text{IC}}=\frac{\lambda}{\lambda+N}

Una muestra/pareada: $d=(\bar y-\mu)/s$ , $\mathrm{SE}=\sqrt{1/n+d^2/(2n)}$ , $\nu=n-1$ . IC de $d$ normal: $d\pm z\,\mathrm{SE}$ ; el de $g$ se escala por $J$ .
IC de $\eta^2$ : se invierte la F no central (Brent) para obtener $\lambda$ ; $\omega^2$ comparte ese mismo IC.
Nota (diferencia con la literatura): en $J$ el código usa $\nu=n_1+n_2-2$ (o $n-1$ ), forma aproximada de Hedges y Olkin, no la función gamma exacta.
Referencias: Cohen 1988, §2.2; Hedges & Olkin 1985, ec. 5.10 y 5.14; Hedges 1981, J. Educ. Stat. 6:107; Hays 1963; Olejnik & Algina 2003, Psychol. Methods 8:434; Steiger 2004, Psychol. Methods 9:164.

Comprobación de supuestos

Código: matlens/core/inference.py::shapiro_wilk, levene, brown_forsythe. Shapiro–Wilk (algoritmo AS R94 vía SciPy, $n\ge3$ ); Levene con center='mean'; Brown–Forsythe = Levene con la mediana. Advertencia si $p<\alpha$ ( $\alpha=0.05$ ). Referencias: Shapiro & Wilk 1965, Biometrika 52:591; Royston 1995; Levene 1960; Brown & Forsythe 1974, JASA 69:364.

PCA con límites T² y Q

Código: matlens/core/multivariate.py::pca, scale_matrix.

\mathbf Z=\mathbf U\boldsymbol\Sigma\mathbf V^\top,\quad \mathbf T=\mathbf U_A\boldsymbol\Sigma_A,\quad \lambda_a=\frac{\sigma_a^2}{n-1}

T^2_i=\sum_{a=1}^{A}\frac{t_{ia}^2}{\lambda_a},\qquad T^2_{\lim}=\frac{A(n-1)}{n-A}\,F_{1-\alpha}(A,\,n-A)

Q_i=\lVert\mathbf z_i-\mathbf t_i\mathbf P_A^\top\rVert^2,\quad Q_{\lim}=\theta_1\left[\frac{z_{1-\alpha}\sqrt{2\theta_2h_0^2}}{\theta_1}+1+\frac{\theta_2h_0(h_0-1)}{\theta_1^2}\right]^{1/h_0},\quad h_0=1-\frac{2\theta_1\theta_3}{3\theta_2^2}

$\mathbf Z$ : datos preprocesados (por defecto autoscale; también none, center, pareto = división por $\sqrt{s}$ , con $s$ de $n-1$ ); $A$ : componentes retenidos (por defecto $\min(n-1,p)$ ); $\theta_j=\sum_{a>A}\lambda_a^j$ ; $\alpha=0.05$ .
Advertencia si $p\ge n$ (límites T² aproximados).
Referencias: Jolliffe 2002; Jackson 1991, ec. 1.7.5; Jackson & Mudholkar 1979, Technometrics 21:341; van den Berg et al. 2006, BMC Genomics 7:142.

Regresión PLS (NIPALS)

Código: matlens/core/multivariate.py::pls_nipals, pls_regression.

Para cada variable latente, partiendo de $\mathbf u$ = columna de $\mathbf Y$ con mayor varianza, se itera hasta $\lVert\mathbf u_{\text{new}}-\mathbf u\rVert<10^{-12}\max(1,\lVert\mathbf u\rVert)$ (máx. 500):

\mathbf w=\frac{\mathbf X^\top\mathbf u}{\lVert\mathbf X^\top\mathbf u\rVert},\ \mathbf t=\mathbf X\mathbf w,\ \mathbf c=\frac{\mathbf Y^\top\mathbf t}{\mathbf t^\top\mathbf t},\ \mathbf u=\frac{\mathbf Y\mathbf c}{\mathbf c^\top\mathbf c};\quad \mathbf p=\frac{\mathbf X^\top\mathbf t}{\mathbf t^\top\mathbf t},\ \mathbf X\leftarrow\mathbf X-\mathbf t\mathbf p^\top,\ \mathbf Y\leftarrow\mathbf Y-\mathbf t\mathbf c^\top

\mathbf B=\mathbf W(\mathbf P^\top\mathbf W)^{+}\mathbf Q^\top,\qquad Q^2_A=1-\frac{\mathrm{PRESS}_A}{SS_{\text{tot}}},\quad \mathrm{RMSECV}_A=\sqrt{\frac{\mathrm{PRESS}_A}{nm}}

El número de variables latentes es el de mínimo PRESS en validación cruzada de $\min(7,n)$ bloques contiguos (sin barajar); máximo $\min(10,p,n-2)$ . Se usa scikit-learn PLSRegression si está instalado; si no, NIPALS con autoescalado de $\mathbf X$ y centrado de $\mathbf Y$ .
Advertencia con $n<10$ . Referencia: Wold et al. 2001, Chemom. Intell. Lab. Syst. 58:109.

Agrupamiento

Código: matlens/core/multivariate.py::kmeans, hierarchical, silhouette. Silueta $s(i)=(b-a)/\max(a,b)$ con distancia euclídea ( $s=0$ en grupos unitarios). k-means con siembra k-means++, n_init=10, semilla 0. Jerárquico por scipy.cluster.hierarchy.linkage (Ward por defecto), advertencia si la correlación cofenética es menor que 0.7. Referencias: Rousseeuw 1987, J. Comput. Appl. Math. 20:53; Lloyd 1982, IEEE Trans. Inf. Theory 28:129; Arthur & Vassilvitskii 2007; Ward 1963, JASA 58:236; Sokal & Rohlf 1962.

Ajuste no lineal de curvas

Código: matlens/core/fitting.py::fit (lmfit, Levenberg–Marquardt por defecto).

\min_{\boldsymbol\theta}\ \chi^2=\sum_{i=1}^{N} w_i^2\,[y_i-f(x_i;\boldsymbol\theta)]^2,\qquad w_i=\frac1{\sigma_i}\ \text{(si se da yerr; si no, } w_i=1)

\chi^2_\nu=\frac{\chi^2}{N-k},\quad R^2=1-\frac{RSS}{TSS},\quad R^2_{\text{aj}}=1-(1-R^2)\frac{N-1}{N-k},\quad \mathbf C=\begin{cases}(\mathbf J^\top\mathbf W^2\mathbf J)^{-1} & \sigma_i \text{ absolutas}\\ \chi^2_\nu\,(\mathbf J^\top\mathbf W^2\mathbf J)^{-1} & \text{sin } \sigma_i \text{ o pesos relativos}\end{cases}

$N$ : puntos válidos (finitos, dentro de x_range, con $\sigma_i>0$ ); $k$ : parámetros libres; $RSS$ y $R^2$ se calculan siempre sin pesos.
Valores iniciales: rutina automática por función (pico, decaimiento, sigmoide…), sobrescribible por initial o por el árbol de parámetros (que también aporta límites y parámetros fijos).
Incertidumbre: errores estándar $k=1$ (raíz de la diagonal de $\mathbf C$ ). La banda de confianza usa eval_uncertainty con sigma=2 por defecto (≈95 %).
Covarianza: con yerr (incertidumbres conocidas) la opción absolute_sigma vale True por defecto y la covarianza no se escala por $\chi^2_\nu$ (Bevington y Robinson, Data Reduction and Error Analysis, 3.ª ed.; Press et al., Numerical Recipes, §15.4); con absolute_sigma=False las $\sigma_i$ se tratan como pesos relativos y la covarianza se escala por $\chi^2_\nu$ , como sin pesos. En el diálogo de ajuste la casilla "σ absolutas" controla esta opción.
Referencia: implementación propia sobre lmfit; nombres y fórmulas según la convención NLFit de Origin.

Ejemplo resuelto

El ejemplo compara la covarianza escalada por $\chi^2_\nu$ (sin yerr) con la de σ absolutas. Cuando las σ dadas coinciden con el ruido real, $\chi^2_\nu\approx1$ y ambas incertidumbres son casi iguales; si las σ declaradas fueran demasiado pequeñas, las incertidumbres con σ absolutas quedarían subestimadas. Las incertidumbres no incluyen el sesgo de modelo (una función inadecuada da errores estándar pequeños y parámetros erróneos) ni la correlación entre residuos.

101 puntos de $y = 0.2 + 3\,e^{-x/2.5}$ con ruido gaussiano σ = 0.05 (semilla fija), ajustados con fit("ExpDec1", ...) (errores estándar $k=1$ ):

Ponderación	y0	A1	t1	χ²_ν	R²_aj
sin σ (covarianza × χ²_ν)	0.19527 ± 0.014	2.9781 ± 0.021	2.5434 ± 0.047	0.002994	0.995065
σ = 0.05 absolutas	0.19527 ± 0.013	2.9781 ± 0.02	2.5434 ± 0.043	1.198	0.995065

Propagación de incertidumbre (GUM)

Código: matlens/core/uncertainty.py::gum_linear, monte_carlo, combine.

u_c^2(y)=\sum_i\left(\frac{\partial f}{\partial x_i}\,u_i\right)^2,\quad \frac{\partial f}{\partial x_i}\approx\frac{f(x_i+h)-f(x_i-h)}{2h},\ h=10^{-6}\max(|x_i|,1),\qquad U=k\,u_c

Combinación lineal $y=\sum c_ix_i$ : $u_c^2=\sum(c_iu_i)^2$ . Se informa la contribución dominante.
Monte Carlo (JCGM 101): $n=20000$ muestras normales independientes (semilla 0); valor = media, $u_c$ = desviación típica ( $n-1$ ), $U$ = semiancho del intervalo central del 95,45 %.
Factor de cobertura por defecto $k=2$ .
Supuestos: entradas no correlacionadas (no se incluyen covarianzas); GUM lineal válida si $f$ es aproximadamente lineal en $\pm u_i$ .
Nota (diferencia con la literatura): en Monte Carlo las "contribuciones" no se separan por entrada; cada una guarda la varianza total y la "dominante" es simplemente la primera entrada.
Referencias: JCGM 100:2008; JCGM 101:2008.

Validez de la aproximación lineal

La fórmula de primer orden es exacta para $f$ lineal; para $f$ no lineal el error relativo crece con $u_i/x_i$ y con la curvatura (en un cociente $m/V$ , el término de segundo orden desplaza la media en $\approx \rho\,(u_V/V)^2$ ). Monte Carlo no tiene esa limitación pero supone, igual que la GUM lineal implementada, entradas normales e independientes. Ninguno de los dos métodos incluye errores sistemáticos no declarados como entradas.

Ejemplo resuelto

Densidad $\rho = m/V$ con $m = 2.500 \pm 0.002$ g, $V = 0.950 \pm 0.010$ cm³ (incertidumbres típicas). GUM (gum_linear): $\rho$ = 2.63158 g/cm³, $u_c$ = 0.0278, $U(k=2)$ = 0.0556, contribución dominante: V. Monte Carlo (monte_carlo, $2\times10^4$ muestras, semilla 0): media = 2.63189, $u$ = 0.028, semiancho del intervalo central 95.45 % = 0.0563.

Ajuste bayesiano por MCMC

Código: matlens/core/bayes.py::curve_fit_bayes, ensemble_sample, split_rhat.

Posterior (verosimilitud gaussiana, priors uniformes dentro de los límites):

\ln p(\boldsymbol\theta\mid\mathbf y)=-\frac12\sum_i\frac{r_i^2}{\sigma_i^2}+\text{cte},\qquad \ln p(\boldsymbol\theta,\ln\sigma\mid\mathbf y)=-\frac{1}{2}e^{-2\ln\sigma}\sum_i r_i^2-N\ln\sigma+\text{cte}

$r_i=y_i-f(x_i;\boldsymbol\theta)$ . La segunda forma se usa cuando no se da $\sigma$ : $\ln\sigma$ es un parámetro más con prior uniforme en $(-30,30)$ . Fuera de límites o con residuos no finitos, $\ln p=-\infty$ .

Movimiento de estiramiento (stretch move), con caminantes divididos en dos mitades que se actualizan por turnos:

\mathbf Y=\mathbf X_j+z(\mathbf X_k-\mathbf X_j),\quad g(z)\propto\frac1{\sqrt z}\ \text{en}\ [1/a,a],\quad z=\frac{[(a-1)u+1]^2}{a},\ u\sim U(0,1)

P_{\text{acept}}=\min\!\left(1,\ z^{D-1}\,\frac{p(\mathbf Y)}{p(\mathbf X_k)}\right)

$a=2$ ; $D$ : dimensión; número de caminantes $W=\max(4D,16)$ redondeado a par (se exige $W\ge 2D$ ); inicio en una bola gaussiana de radio relativo $10^{-3}$ alrededor de p0 (y de $\ln s_0$ , con $s_0$ la desviación de los residuos iniciales); 3000 pasos, 30 % de quemado, semilla 0.

Split- $\hat R$ (cada cadena se parte en dos mitades de longitud $n$ ; $m=2W$ cadenas):

B=n\,\mathrm{Var}_j(\bar\theta_j),\quad W_v=\overline{s_j^2},\quad \hat V=\frac{n-1}{n}W_v+\frac{B}{n},\quad \hat R=\sqrt{\hat V/W_v}

Resultados: mediana y percentiles 2,5/16/84/97,5 (intervalos creíbles del 68 % y 95 %, no "valor ± error"). Advertencias si la fracción de aceptación es menor que 0,15 o si $\hat R>1{,}05$ .
Nota (diferencia con la literatura): el $\hat R$ se calcula con cada caminante como cadena y no incluye la normalización por rangos de Vehtari et al. (2021).
Referencias: Goodman & Weare, Commun. Appl. Math. Comput. Sci. 5, 65 (2010); Foreman-Mackey et al., PASP 125, 306 (2013); Gelman et al., Bayesian Data Analysis, 3.ª ed., 2013, §11.4.

Mapas hiperespectrales: MCR-ALS

Código: matlens/core/hyperspectral.py::mcr_als, purest_spectra, pca_map.

\mathbf D=\mathbf C\mathbf S^\top+\mathbf E

$\mathbf D$ : cubo $(n_y,n_x,n_c)$ desplegado a píxeles × canales (intensidad en unidades del instrumento); $\mathbf C$ : concentraciones (píxeles × $k$ ); $\mathbf S$ : espectros puros (canales × $k$ ); $\mathbf E$ : residuo.

Mínimos cuadrados alternados con no negatividad (NNLS por columnas):

\mathbf C=\arg\min_{\mathbf C\ge0}\lVert\mathbf D-\mathbf C\mathbf S^\top\rVert_F^2,\qquad \mathbf S=\arg\min_{\mathbf S\ge0}\lVert\mathbf D-\mathbf C\mathbf S^\top\rVert_F^2

Restricciones: no negatividad en $\mathbf C$ y $\mathbf S$ ; cierre opcional ( $\sum_j C_{ij}=1$ , por normalización de filas). Sin cierre, cada espectro se normaliza a norma unidad y $\mathbf C$ absorbe la escala.
Inicio: espectros de los píxeles más puros (tipo SIMPLISMA): pureza $s_i/(\mu_i+0{,}05\max\mu)$ , y los siguientes maximizan pureza × $\det(\mathbf M\mathbf M^\top)$ de los ya elegidos; o espectros dados por el usuario.
Convergencia: $|r_{t-1}-r_t|<10^{-7}\max(r_t,10^{-12})$ , máximo 200 iteraciones, con $r=\sum E^2/\sum D^2$ .

\mathrm{LOF}\,(\%)=100\sqrt{\frac{\sum E_{ij}^2}{\sum D_{ij}^2}},\qquad R^2\,(\%)=100\left(1-\frac{\sum E_{ij}^2}{\sum D_{ij}^2}\right)

Siempre se advierte de la ambigüedad rotacional: la solución es factible, no única. No se calcula incertidumbre de $\mathbf C$ ni $\mathbf S$ .
pca_map: PCA de $\mathbf D$ centrado, varianza explicada $\sigma_a^2/\sum\sigma^2$ , para elegir $k$ (por defecto 6 componentes).
Referencias: Tauler, Chemom. Intell. Lab. Syst. 30, 133 (1995); de Juan & Tauler, Anal. Chim. Acta 500, 195 (2003); Windig & Guilment, Anal. Chem. 63, 1425 (1991).

Series operando: ajuste secuencial

Código: matlens/core/operando.py::sequential_fit.

\hat{\boldsymbol\theta}_t=\arg\min_{\boldsymbol\theta}\sum_i\big[f(x_i;\boldsymbol\theta)-y_{t,i}\big]^2,\qquad \boldsymbol\theta^{(0)}_t=\hat{\boldsymbol\theta}_{t-1}\ \text{(si el cuadro anterior convergió)}

s^2=\frac{\sum r_i^2}{\max(N-k,1)},\qquad \mathbf C=s^2(\mathbf J^\top\mathbf J)^{+},\qquad \mathrm{rms}=s

Cada cuadro (espectro o patrón frente a tiempo, temperatura o potencial) se ajusta con scipy.optimize.least_squares (sin pesos), con límites opcionales; el primero parte de p0. reverse=True recorre la serie desde el último cuadro. Se ignoran puntos no finitos.
Errores estándar con $k=1$ (1σ).
Referencia: implementación propia.

Series operando: puntos de cambio (PELT)

Código: matlens/core/operando.py::change_points, _seg_cost.

F(t)=\min_{s\in\mathcal R_t}\big[F(s)+\mathcal C(y_{s:t})+\beta\big],\quad F(0)=-\beta,\qquad \mathcal C(y_{a:b})=(b-a)\ln\hat\sigma^2_{a:b}

$\hat\sigma^2_{a:b}$ : varianza (máxima verosimilitud) del segmento con su propia media, con suelo $10^{-12}$ ; coste = log-verosimilitud gaussiana negativa salvo constantes (media y varianza cambian).
Penalización por defecto $\beta=3\ln n$ (tipo BIC, dos parámetros por segmento); longitud mínima de segmento 3. Poda: se descartan candidatos $s$ con $F(s)+\mathcal C(y_{s:t})>F(t)$ .
Devuelve los índices donde empieza cada nuevo segmento. Sin incertidumbre en las posiciones.
Nota (diferencia con la literatura): el coste habitual de PELT para media y varianza es $n\ln\hat\sigma^2$ más una constante; aquí se usa sin el término constante y con suelo de varianza, y la penalización $3\ln n$ es una elección del código.
Referencia: Killick, Fearnhead & Eckley, J. Am. Stat. Assoc. 107, 1590 (2012).