Asignatura: Estadística Aplicada con Python y R Programa: Ingeniería Agrícola — Universidad de Sucre
Un ingeniero agrícola se enfrenta constantemente a la misma pregunta bajo dos disfraces diferentes:
En ambos casos el procedimiento es el mismo método científico, aplicado a dos objetos matemáticos distintos: una distribución de probabilidad \(f(x;\boldsymbol\theta)\) en el primer caso, y una función determinista \(y = g(t;\boldsymbol\theta)\) en el segundo. Este documento desarrolla la teoría común a ambos y, al final, muestra por qué —en más de un caso real— son literalmente la misma cosa vista desde dos ángulos.
El método tiene seis pasos, y ninguno es opcional:
Se tiene una muestra \(x_1, x_2, \dots, x_n\) de una variable aleatoria continua \(X\), y se supone que proviene de una familia paramétrica \(f(x;\boldsymbol\theta)\) (Normal, Gamma, Weibull, Gumbel…), donde \(\boldsymbol\theta\) es un vector de parámetros desconocidos (por ejemplo, \(\boldsymbol\theta = (\mu,\sigma)\) en la Normal). Ajustar es encontrar el valor de \(\boldsymbol\theta\) que hace más compatible al modelo con la muestra observada.
Elegir la familia no es un paso estadístico sino físico: la variable, ¿es siempre positiva? ¿simétrica o sesgada? ¿es un máximo de bloque (justifica una familia de valores extremos)? ¿es un tiempo hasta un evento (justifica una familia de supervivencia)? La estadística estima los parámetros de la familia elegida; no elige la familia por el estudiante.
Método de los momentos. Se igualan los momentos muestrales (media, varianza…) con las expresiones teóricas de esos momentos en función de \(\boldsymbol\theta\), y se despeja. Es sencillo y da un punto de partida razonable, pero no es eficiente estadísticamente (desperdicia información de la muestra) y a veces produce estimaciones fuera del espacio admisible de parámetros.
Máxima verosimilitud (MLE). Es el método estándar. La función de verosimilitud mide qué tan plausible es haber observado la muestra si los parámetros fueran \(\boldsymbol\theta\):
\[L(\boldsymbol\theta) = \prod_{i=1}^{n} f(x_i;\boldsymbol\theta), \qquad \ell(\boldsymbol\theta) = \ln L(\boldsymbol\theta) = \sum_{i=1}^{n}\ln f(x_i;\boldsymbol\theta)\]
El estimador de máxima verosimilitud es el que maximiza esa log-verosimilitud:
\[\hat{\boldsymbol\theta}_{\text{MLE}} = \arg\max_{\boldsymbol\theta} \; \ell(\boldsymbol\theta)\]
Para algunas familias (la Normal) existe solución cerrada; para la mayoría (Gamma, Weibull, Gumbel, GEV) no la hay, y el computador maximiza \(\ell\) numéricamente mediante algoritmos de optimización iterativos. Bajo condiciones de regularidad razonables, el MLE es consistente (converge al valor verdadero cuando \(n\to\infty\)), asintóticamente eficiente (tiene la menor varianza posible entre estimadores insesgados, en el límite) y asintóticamente normal, lo que permite construir errores estándar e intervalos de confianza para los parámetros a partir de la curvatura de \(\ell\) en su máximo (la matriz de información de Fisher).
L-momentos (mención). En hidrología y ciencias ambientales es común un tercer método, basado en combinaciones lineales de los datos ordenados, más robusto que los momentos clásicos frente a valores atípicos y frecuentemente usado para distribuciones de valores extremos (Gumbel, GEV).
El histograma engaña. Su forma cambia con el número de barras (bins) elegido arbitrariamente; el mismo conjunto de datos puede parecer uni o bimodal, simétrico o sesgado, solo por esa elección. Es un buen punto de partida exploratorio, nunca la prueba final de un ajuste.
Función de distribución empírica (ECDF).
\[F_n(x) = \frac{\#\{x_i \le x\}}{n}\]
Es una función escalonada que no depende de ninguna decisión arbitraria de agrupación, y que —si el modelo es correcto— debe acompañar de cerca a la CDF teórica \(F(x;\hat{\boldsymbol\theta})\).
Gráfico cuantil-cuantil (Q-Q). Es el diagnóstico más informativo. Se ordenan los datos \(x_{(1)} \le \dots \le x_{(n)}\), se les asigna una posición de graficación (una probabilidad acumulada estimada para cada dato, por ejemplo \(p_i = \frac{i-0.44}{n+0.12}\) o \(p_i=\frac{i}{n+1}\)), se calcula el cuantil teórico \(F^{-1}(p_i;\hat{\boldsymbol\theta})\) y se grafica el par \(\big(F^{-1}(p_i),\,x_{(i)}\big)\). Si el modelo es correcto, los puntos caen sobre la recta de 45°. La forma de la desviación es diagnóstica:
| Patrón observado | Interpretación |
|---|---|
| Puntos sobre la recta en todo el rango | Buen ajuste |
| Cola derecha por encima de la recta | Los datos tienen valores extremos más grandes de lo que predice el modelo (el modelo subestima el riesgo) |
| Cola derecha por debajo de la recta | El modelo sobrestima los extremos |
| Curvatura sistemática (forma de “S”) | Asimetría o curtosis no capturada por la familia elegida |
| Recta paralela mas no coincidente | Parámetro de ubicación o escala mal estimado |
El gráfico P-P (probabilidad-probabilidad, \(F_n(x)\) contra \(F(x;\hat{\boldsymbol\theta})\)) es un complemento: es más sensible al centro de la distribución, mientras que el Q-Q es más sensible a las colas — la parte que casi siempre le interesa al ingeniero (diseño para eventos extremos, tiempos de proceso largos).
Todas se enmarcan en una prueba de hipótesis:
\[H_0:\; \text{la muestra proviene de } F(x;\hat{\boldsymbol\theta}) \qquad \text{vs.} \qquad H_1:\; \text{no proviene de esa distribución}\]
con la regla usual: si el valor \(p\) es menor que el nivel de significancia \(\alpha\) (típicamente 0.05), se rechaza \(H_0\).
Kolmogorov-Smirnov (KS). Se basa en la distancia máxima vertical entre la ECDF y la CDF teórica:
\[D_n = \sup_{x} \big| F_n(x) - F(x;\hat{\boldsymbol\theta}) \big|\]
Es sensible principalmente al centro de la distribución.
Anderson-Darling (AD). Una versión ponderada de la misma idea, que da más peso a las colas:
\[A^2 = -n - \frac{1}{n}\sum_{i=1}^{n}(2i-1)\Big[\ln F(x_{(i)}) + \ln\big(1-F(x_{(n+1-i)})\big)\Big]\]
Suele ser la prueba de elección en ingeniería, porque el interés casi siempre está en los extremos (crecidas, tiempos máximos de proceso).
Chi-cuadrado de Pearson. Se agrupan los datos en clases, se comparan las frecuencias observadas \(O_j\) con las esperadas bajo el modelo \(E_j\):
\[\chi^2 = \sum_{j} \frac{(O_j - E_j)^2}{E_j}\]
que sigue aproximadamente una distribución \(\chi^2\) con \(g - 1 - k\) grados de libertad (\(g\) = número de clases, \(k\) = número de parámetros estimados). Hereda el problema del histograma: su resultado depende de cómo se agrupen los datos.
Advertencia estadística importante. Cuando los parámetros \(\hat{\boldsymbol\theta}\) se estiman con la misma muestra que luego se usa para la prueba, el valor \(p\) “clásico” de KS es demasiado optimista (problema conocido desde Lilliefors, 1967): el modelo se ajustó a la medida de esos datos, así que es natural que parezca ajustar bien. La solución moderna es re-calcular el valor \(p\) por simulación Monte Carlo: se generan muchas muestras sintéticas del modelo ajustado, se re-estiman sus parámetros y se recalcula el estadístico en cada una, construyendo así la distribución nula correcta del estadístico bajo esa contaminación.
Límites de las pruebas formales. (i) No rechazar \(H_0\) no prueba que el modelo sea verdadero, solo que no hay evidencia suficiente en contra; con muestras pequeñas la prueba tiene poca potencia. (ii) Con muestras muy grandes casi cualquier modelo se rechaza, aunque sea útil en la práctica (“todos los modelos son incorrectos, algunos son útiles” — George Box). (iii) Evaluar muchas familias candidatas a la vez aumenta la probabilidad de que alguna “pase” por puro azar (data dredging); restringir las candidatas por argumento físico es, por tanto, parte del método, no un atajo.
Las pruebas de hipótesis responden “¿rechazo este modelo?”; los criterios de información responden una pregunta distinta y más útil cuando hay varias candidatas razonables: “¿cuál equilibra mejor el ajuste con la simplicidad?”
\[\text{AIC} = 2k - 2\hat\ell \qquad \text{BIC} = k\ln n - 2\hat\ell \qquad \text{AICc} = \text{AIC} + \frac{2k(k+1)}{n-k-1}\]
donde \(k\) es el número de parámetros estimados y \(\hat\ell\) el log-verosimilitud maximizado. El término \(2k\) (o \(k\ln n\) en BIC) penaliza la complejidad: agregar parámetros siempre mejora \(\hat\ell\), pero no siempre mejora el modelo — es la formalización estadística de la navaja de Ockham. Menor AIC/BIC es mejor. Con muestras pequeñas (\(n/k < 40\)) se prefiere el AICc, que corrige el sesgo de AIC en ese régimen.
Para comparar modelos se usa la diferencia \(\Delta_i = \text{AIC}_i - \text{AIC}_{\min}\): valores de \(\Delta_i < 2\) indican soporte sustancial similar al mejor modelo; \(4\)–\(7\), soporte considerablemente menor; \(>10\), prácticamente sin soporte. Estos criterios no dicen si un modelo es “correcto”; solo ordenan candidatas entre sí.
El objetivo final no es el ajuste en sí, sino un número que sirva para decidir: el cuantil de diseño \(x_p = F^{-1}(p;\hat{\boldsymbol\theta})\). En variables de valores extremos (precipitación máxima, caudal pico) se expresa en función del período de retorno \(T\) (años):
\[x_T = F^{-1}\!\left(1-\frac{1}{T}\right)\]
y el riesgo de que ese valor sea excedido al menos una vez en una vida útil de \(N\) años es \(R = 1-(1-1/T)^N\) — una obra diseñada para \(T=50\) años con vida útil de 50 años tiene, en realidad, un riesgo de falla cercano al 64 %, no del 2 %.
Como los parámetros \(\hat{\boldsymbol\theta}\) se estimaron a partir de una muestra finita, \(x_T\) también tiene incertidumbre de muestreo, que se cuantifica típicamente con un bootstrap paramétrico: se simulan muchas muestras del tamaño original desde el modelo ajustado, se re-estima \(\boldsymbol\theta\) y se recalcula \(x_T\) en cada una, obteniendo así un intervalo de confianza. Reportar \(x_T\) sin su incertidumbre es, en ingeniería, una afirmación incompleta.
| Tipo de modelo | Origen | Ventaja | Limitación |
|---|---|---|---|
| Teórico | Se deriva resolviendo las ecuaciones físicas del fenómeno (p. ej. la segunda ley de Fick para difusión de humedad) | Explica el mecanismo; extrapola razonablemente | Requiere geometría, propiedades termofísicas y suele tener solución en serie infinita, costosa de evaluar |
| Semiempírico (semiteórico) | Se obtiene simplificando un modelo teórico (truncando una serie, o linealizando) | Conserva parte del significado físico de sus parámetros; más simple de ajustar | Válido solo dentro del rango de condiciones (temperatura, humedad relativa, velocidad de aire…) con que se derivó |
| Empírico | Relación matemática conveniente entre la variable de interés y el tiempo (o la variable independiente), sin partir de ningún principio físico | Muy flexible, fácil de ajustar, buen ajuste dentro del rango de datos | No informa sobre el mecanismo; puede comportarse de forma físicamente absurda fuera del rango de datos (extrapolación peligrosa) |
Un ejemplo trabajado en el curso: la difusión de humedad en una placa (ley de Fick) da una solución en serie infinita; truncarla en su primer término da el modelo semiempírico de Henderson y Pabis; elevar el tiempo a una potencia libre (\(t^n\)) para ganar flexibilidad de ajuste, sin que ese exponente salga de ninguna derivación física, da el modelo empírico de Page.
Se tiene ahora una variable respuesta \(y_i\) medida en puntos \(t_i\) (por ejemplo, humedad relativa \(MR\) medida cada hora), y se propone un modelo determinista \(y = g(t;\boldsymbol\theta)\) con parámetros \(\boldsymbol\theta\) desconocidos. Se estiman minimizando la suma de cuadrados de los residuos:
\[\hat{\boldsymbol\theta}_{\text{MCO}} = \arg\min_{\boldsymbol\theta} \; \text{SSE}(\boldsymbol\theta) = \sum_{i=1}^{n} \big[y_i - g(t_i;\boldsymbol\theta)\big]^2\]
Este es, de nuevo, un problema de máxima verosimilitud. Si se asume que los errores de medición son independientes y se distribuyen \(\varepsilon_i \sim N(0,\sigma^2)\), la log-verosimilitud es una función decreciente de SSE, así que minimizar SSE es exactamente maximizar la verosimilitud bajo ese supuesto. Los dos hilos de este documento —ajuste de distribuciones y ajuste de curvas— son, en el fondo, la misma optimización.
A diferencia de la regresión lineal, \(g\) no es lineal en \(\boldsymbol\theta\) (por ejemplo, \(t^n\) con \(n\) desconocido), así que no hay solución cerrada: se resuelve con algoritmos iterativos (Gauss-Newton, Levenberg-Marquardt, o variantes con restricciones como Trust Region Reflective), que parten de valores iniciales razonables y convergen mejorando la solución paso a paso. Una consecuencia práctica importante: los valores iniciales importan. Un mal punto de partida puede hacer que el algoritmo no converja o converja a un óptimo local en vez del global, especialmente con modelos de varios parámetros correlacionados entre sí (como el modelo de dos exponenciales).
Coeficiente de determinación (\(R^2\)).
\[R^2 = 1 - \frac{\text{SSE}}{\sum_i (y_i-\bar y)^2}\]
Mide qué fracción de la variabilidad de \(y\) explica el modelo. Su limitación central: nunca disminuye al agregar parámetros, así que no sirve para comparar modelos con distinto número de parámetros — un modelo con más parámetros casi siempre tendrá \(R^2\) mayor o igual, aunque sea peor en el sentido de parsimonia (y de comportamiento fuera del rango de datos).
Raíz del error cuadrático medio (RMSE). \(\sqrt{\text{SSE}/n}\), en las mismas unidades que \(y\); más interpretable que \(R^2\) como “error típico de predicción”.
Análisis de residuos. El diagnóstico más importante y el más olvidado. Se grafican los residuos \(y_i - \hat y_i\) contra \(t_i\) (o contra \(\hat y_i\)): deben verse aleatorios, sin patrón, con varianza aproximadamente constante (homocedasticidad). Un patrón sistemático —una curva, ondas, un abanico que se abre— indica que el modelo omite algo, aunque el \(R^2\) sea altísimo. Un gráfico Q-Q de los residuos contra la Normal (y una prueba como Shapiro-Wilk) evalúa si el supuesto de errores normales —el que justifica que mínimos cuadrados equivalga a máxima verosimilitud— es razonable.
Con varios modelos candidatos que describen la misma curva, se aplica el mismo criterio de la sección 1.5, adaptado a la log-verosimilitud gaussiana:
\[\text{AIC} = n\Big[\ln(2\pi) + 1 + \ln\big(\text{SSE}/n\big)\Big] + 2(p+1), \qquad \text{AICc} = \text{AIC} + \frac{2k(k+1)}{n-k-1},\; k=p+1\]
donde \(p\) es el número de parámetros del modelo (el \(+1\) cuenta la varianza del error, que también se estima). Con pocos datos por curva —algo frecuente en ensayos de secado o de crecimiento, con una decena de mediciones— el AICc es casi obligatorio: penaliza con más fuerza los modelos sobreparametrizados quer un \(R^2\) engañosamente alto no delataría.
Un modelo empírico se ajusta dentro del rango de los datos observados; fuera de ese rango, su comportamiento no tiene por qué tener sentido físico. Ejemplos típicos: un polinomio ajustado a una curva de secado puede, para tiempos largos, predecir que el producto se humedece de nuevo; una curva sigmoide ajustada a pocos puntos de crecimiento puede no capturar una fase de senescencia que aún no se ha observado. Un modelo con base física (teórico o semiempírico) suele comportarse mejor fuera del rango de calibración, porque su forma no es arbitraria sino heredada del mecanismo. Por eso, todo ajuste empírico debe reportarse con su rango de validez explícito, y cualquier uso fuera de ese rango debe señalarse como extrapolación, no como predicción.
Los dos hilos de este documento convergen en una identidad útil. Sea \(T\) el tiempo aleatorio en que ocurre cierto evento (un componente falla, una partícula de agua sale de un grano). Su función de supervivencia es la probabilidad de que el evento aún no haya ocurrido en el instante \(t\):
\[S(t) = P(T > t) = 1 - F(t)\]
Si \(T\) tiene distribución Weibull con forma \(n\) y escala \(\lambda\), esa supervivencia es:
\[S(t) = \exp\!\left[-\left(\frac{t}{\lambda}\right)^{n}\right]\]
que, reescrita con \(k = \lambda^{-n}\), es exactamente la ecuación de Page usada para describir curvas de secado: \(MR(t) = e^{-k t^n}\). La razón de humedad de una curva de secado es la función de supervivencia de un “tiempo de salida del agua”; la velocidad de secado, \(-dMR/dt\), es la densidad de probabilidad de ese tiempo. Con \(n=1\), la Weibull se reduce a la Exponencial, y Page se reduce a Lewis.
Esta equivalencia no es una curiosidad matemática: significa que todo lo desarrollado en la sección 1 (estimación por MLE, diagnóstico Q-Q, pruebas de bondad de ajuste, AIC) aplica igual de bien si se reformula el problema de la curva como un problema de distribución (por ejemplo, mirando la variabilidad de tiempos de proceso entre lotes en vez de la curva promedio), y viceversa: todo lo desarrollado en la sección 2 (mínimos cuadrados, residuos, rango de validez) es MLE bajo un supuesto de error particular. El método es uno solo; los nombres cambian según el objeto que se ajusta.
| Paso del método | En una distribución | En una curva empírica/semiempírica |
|---|---|---|
| 1. Explorar | Histograma, ECDF, asimetría | Gráfico \(y\) vs. \(t\), transformaciones (p. ej. \(\ln y\) vs. \(t\)) |
| 2. Proponer candidatas | Familias con argumento físico (positiva, extremos, tiempo hasta evento…) | Modelos teóricos, semiempíricos y empíricos, en ese orden de preferencia |
| 3. Estimar | Máxima verosimilitud (dist.fit()) |
Mínimos cuadrados no lineales (curve_fit()) |
| 4. Diagnosticar | Q-Q, P-P | Curva ajustada sobre los datos, residuos vs. \(t\), Q-Q de residuos |
| 5. Contrastar | KS, Anderson-Darling, \(\chi^2\); AIC/BIC entre candidatas | \(R^2\)/RMSE (con cautela), AIC/AICc entre candidatas |
| 6. Usar | Cuantiles, período de retorno, con intervalo de confianza | Predicción dentro del rango de validez, con intervalo de confianza de los parámetros |