# geoestadistica-con-pykrige-variogramas-e-interpolacion-espacial

Geoestadística con PyKrige: Variogramas e Interpolación Espacial

Volver a la portada del curso
Lectura (Markdown)

Validación de Modelos Kriging

Sección 5 · Validación de Modelos Espaciales

Validación de modelos kriging

Un mapa de kriging siempre se ve bien: es suave, continuo y tiene colores agradables. Eso no dice nada sobre si es correcto. La validación responde dos preguntas distintas: ¿qué tan precisas son las estimaciones? y ¿la varianza de kriging describe bien el error real? Un buen modelo cumple ambas.

Validación cruzada dejando uno fuera

Con pocos datos no conviene separar un conjunto de prueba grande. La alternativa clásica es leave-one-out (LOO): se quita un punto, se estima con los demás, se compara con el valor real y se repite para todos. Con cada error $e_i = z_i - \hat{z}_i$ y su varianza de kriging $\sigma_i^2$ se calculan:

  • Error medio (ME): $\frac{1}{n}\sum e_i$. Debe ser cercano a cero; si no, el modelo tiene sesgo.
  • RMSE: $\sqrt{\frac{1}{n}\sum e_i^2}$. Mide la precisión, en las mismas unidades que los datos.
  • Cociente de desviación cuadrática media (MSDR): $\frac{1}{n}\sum e_i^2/\sigma_i^2$. Debe ser cercano a 1. Si es mucho mayor, el modelo subestima su incertidumbre; si es mucho menor, la sobreestima.
import numpy as np
from pykrige.ok import OrdinaryKriging

def loo_kriging(x, y, z, **kw):
    e, s2 = np.empty(len(z)), np.empty(len(z))
    for i in range(len(z)):
        m = np.arange(len(z)) != i
        ok = OrdinaryKriging(x[m], y[m], z[m], **kw)
        est, var = ok.execute("points", x[i:i + 1], y[i:i + 1])
        e[i], s2[i] = z[i] - est[0], var[0]
    return e, s2

e, s2 = loo_kriging(x, y, z, variogram_model="exponential", nlags=12)
print(e.mean(), np.sqrt((e**2).mean()), (e**2 / s2).mean())

En los 150 pozos del curso (profundidad del nivel freático, en metros) este modelo da ME = −0.017 m, RMSE = 1.52 m y MSDR = 0.96. Además, el 96% de los valores reales cae dentro de $\hat{z}_i \pm 1.96\,\sigma_i$, muy cerca del 95% nominal. Observa que execute devuelve la varianza, no la desviación estándar: hay que sacar raíz para construir intervalos.

Comparar contra referencias simples

Un RMSE aislado no significa nada. Compáralo con dos referencias: predecir siempre la media (RMSE igual a la desviación estándar de los datos, 2.22 m aquí) y la distancia inversa ponderada con exponente 2, que en LOO da 1.61 m. El kriging ordinario mejora a ambas; si no lo hiciera, el variograma estaría mal especificado o los datos no tendrían estructura espacial útil.

Las estadísticas Q1, Q2 y cR de PyKrige

Con enable_statistics=True, PyKrige calcula tres diagnósticos propuestos por Kitanidis. Ojo con un detalle: no hacen LOO, sino una validación secuencial. Cada punto $i$ se estima solo con los puntos $0, \dots, i-1$, en el orden en que los entregaste, y se estandariza el error: $\varepsilon_i = e_i/\sigma_i$.

  • Q1 es el valor absoluto de la media de $\varepsilon$. Debe quedar cerca de 0; como guía, por debajo de $2/\sqrt{n-1}$.
  • Q2 es la media de $\varepsilon^2$. Debe quedar cerca de 1, aproximadamente dentro de $1 \pm 2.8/\sqrt{n-1}$.
  • cR combina Q2 con la media geométrica de las varianzas. Entre modelos que pasan Q1 y Q2, conviene el de menor cR.
ok = OrdinaryKriging(x, y, z, variogram_model="exponential", nlags=12,
                     enable_statistics=True)
print(round(ok.Q1, 3), round(ok.Q2, 3), round(ok.cR, 3))  # 0.062 1.026 2.976
ok.print_statistics()

Con $n = 150$, los umbrales son 0.16 para Q1 y [0.77, 1.23] para Q2: el modelo pasa los dos. Como dependen del orden, una permutación de los datos cambia ligeramente Q1. En el notebook de modelos de variograma, un gaussiano mal ajustado da Q2 = 115: esa es la señal de alarma que buscas.

Búsqueda de hiperparámetros con scikit-learn

PyKrige incluye pykrige.rk.Krige, un estimador compatible con scikit-learn. Recibe X con las coordenadas en columnas, así que funciona con GridSearchCV para comparar modelos y valores de nlags con validación cruzada por bloques:

from pykrige.rk import Krige
from sklearn.model_selection import GridSearchCV, KFold

rejilla = {"method": ["ordinary", "universal"],
           "variogram_model": ["linear", "spherical", "exponential", "gaussian"],
           "nlags": [6, 12]}
busqueda = GridSearchCV(Krige(), rejilla, cv=KFold(5, shuffle=True, random_state=0),
                        scoring="neg_root_mean_squared_error")
busqueda.fit(np.c_[x, y], z)
print(busqueda.best_params_, -busqueda.best_score_)

En los pozos gana el kriging ordinario exponencial con nlags=6 (RMSE de 1.55 m en 5 bloques). Un detalle importante: Krige usa por defecto n_closest_points=10 en el kriging ordinario, es decir, una ventana móvil de 10 vecinos. En el universal ese parámetro se ignora y PyKrige lo avisa.

Trampas comunes

  • Validar con los mismos datos que ajustaron el variograma sin reajustarlo en cada pliegue: da métricas optimistas. Para un número honesto, ajusta todo dentro del ciclo.
  • Mirar solo el RMSE: un modelo con buen RMSE y MSDR de 3 entrega intervalos de confianza falsos.
  • Usar exact_values=True con pepita grande: en los puntos medidos el kriging devuelve el dato exacto con varianza 0, aunque haya error de medición. Con exact_values=False, la pepita se trata como ruido y el mapa se suaviza.
  • Datos agrupados: si los pozos se concentran en una zona, el LOO mide sobre todo esa zona. Considera validación por bloques espaciales.

Cierre

Reporta siempre tres números junto a un mapa de kriging: el RMSE de validación cruzada, la comparación contra la referencia más simple y el MSDR (o Q2). Juntos dicen si el mapa es preciso y si su incertidumbre es creíble.

Recursos