Regresión por Procesos Gaussianos (GPR)

La Regresión por Procesos Gaussianos es un método no paramétrico y bayesiando. A diferencia de la regresión lineal, donde buscamos los parámetros de una función, en GPR buscamos una distribución sobre funciones.

1. Definición Formal

Un Proceso Gaussiano (GP) es una colección de variables aleatorias, tal que cualquier subconjunto finito de ellas tiene una distribución normal conjunta.

Se define completamente por su función de media \(m(x)\) y su función de covarianza (kernel) \(k(x, x')\):

\[f(x) \sim \mathcal{GP}(m(x), k(x, x'))\]

En la práctica, solemos asumir por simplicidad que \(m(x) = 0\).

2. El Kernel (Función de Covarianza)

El kernel define la “forma” de las funciones (su suavidad, periodicidad, etc.). El más común es el Kernel de Exponencial Cuadrática (RBF):

\[k(x, x') = \sigma_f^2 \exp \left( -\frac{\|x - x'\|^2}{2\ell^2} \right)\]

\(\sigma_f^2\): Varianza de la señal (amplitud).
\(\ell\): Escala de longitud (qué tan rápido cambia la función).

3. Demostración: Predicciones con GPR

Queremos predecir el valor \(f_*\) para un punto nuevo \(x_*\), dado un conjunto de datos de entrenamiento \(\mathbf{X}\) y sus etiquetas \(\mathbf{y}\).

A. La Distribución Conjunta

Asumiendo que las observaciones tienen ruido gaussiano \(\epsilon \sim \mathcal{N}(0, \sigma_n^2)\), los datos observados son \(y = f(x) + \epsilon\). La distribución conjunta entre los datos de entrenamiento \(\mathbf{y}\) y el valor crítico \(f_*\) es:

\[\begin{bmatrix} \mathbf{y} \\ f_* \end{bmatrix} \sim \mathcal{N} \left( \mathbf{0}, \begin{bmatrix} K(X, X) + \sigma_n^2 I & K(X, x_*) \\ K(x_*, X) & k(x_*, x_*) \end{bmatrix} \right)\]

Para simplificar la notación, definamos \(K = K(X, X) + \sigma_n^2 I\).

B. Condicionamiento (Inferencia)

Para obtener la predicción, calculamos la probabilidad condicional \(P(f_* \mid X, \mathbf{y}, x_*)\). Aplicando el teorema de la distribución condicional de una normal multivariada, obtenemos:

Media Predictiva (\(\mu_*\)):

\[\bar{f}_* = K(x_*, X) [K(X, X) + \sigma_n^2 I]^{-1} \mathbf{y}\]

Varianza Predictiva (\(\text{var}(f_*)\)):

\[\text{var}(f_*) = k(x_*, x_*) - K(x_*, X) [K(X, X) + \sigma_n^2 I]^{-1} K(X, x_*)\]

4. Optimización: Máxima Verosimilitud Marginal

¿Cómo elegimos los hiperparámetros \(\theta\) (como \(\ell\) y \(\sigma_f\))? Minimizamos el logaritmo de la verosimilitud marginal negativa (LML):

\[\log P(\mathbf{y} \mid X, \theta) = -\frac{1}{2} \mathbf{y}^\top K_y^{-1} \mathbf{y} - \frac{1}{2} \log |K_y| - \frac{n}{2} \log 2\pi\]

Donde \(K_y = K(X, X) + \sigma_n^2 I\).

Término de ajuste: \(-\frac{1}{2} \mathbf{y}^\top K_y^{-1} \mathbf{y}\) (qué tan bien explica el modelo los datos).
Penalización por complejidad: \(-\frac{1}{2} \log |K_y|\) (evita el sobreajuste).

5. Ventajas y Desventajas

Ventajas	Desventajas
Incertidumbre: Proporciona un intervalo de confianza para cada predicción.	Escalabilidad: Requiere invertir una matriz \(n \times n\), con costo \(O(n^3)\).
Flexibilidad: Puedes modelar casi cualquier comportamiento eligiendo el kernel adecuado.	Memoria: Almacenar la matriz de covarianza es \(O(n^2)\).
Bayesiano: Funciona bien con pocos datos de entrenamiento.	Selección de Kernel: Requiere conocimiento previo del dominio.

Algoritmo

Definir un kernel y sus hiperparámetros iniciales.
Calcular la matriz de covarianza \(K\) entre todos los puntos de entrenamiento.
Invertir \(K\) (normalmente usando la descomposición de Cholesky para estabilidad numérica).
Calcular la media y varianza para los nuevos puntos de prueba usando las fórmulas del Paso 3.
Optimizar los hiperparámetros si es necesario mediante gradiente descendente sobre la LML.

1. La Intuición: Del Espacio de Pesos al Espacio de Funciones

Imagina una regresión lineal estándar: \(f(x) = \mathbf{w}^\top \mathbf{x}\). Aquí, la incertidumbre está en los pesos \(\mathbf{w}\). En GPR, movemos esa incertidumbre directamente a los valores de la función.

La “Vista de Pesos” (Interpretación Dual)

Si tenemos un modelo \(f(x) = \phi(x)^\top \mathbf{w}\) con un prior \(\mathbf{w} \sim \mathcal{N}(0, \Sigma_p)\), la media y la covarianza de la función son:

Media: \(E[f(x)] = \phi(x)^\top E[\mathbf{w}] = 0\)
Covarianza: \(E[f(x)f(x')] = \phi(x)^\top E[\mathbf{w}\mathbf{w}^\top] \phi(x') = \phi(x)^\top \Sigma_p \phi(x')\)

Aquí ocurre la magia: definimos \(k(x, x') = \phi(x)^\top \Sigma_p \phi(x')\). No necesitamos conocer las características \(\phi(x)\) explícitamente, solo cómo se relacionan los puntos a través del kernel \(k\). Esto es el Kernel Trick.

2. Demostración: La Distribución Condicional Posterior

Esta es la parte más importante: ¿por qué la predicción tiene esa forma?

El Teorema de la Normal Condicional

Si tenemos un vector gaussiano conjunto:

\[\begin{bmatrix} \mathbf{a} \\ \mathbf{b} \end{bmatrix} \sim \mathcal{N} \left( \begin{bmatrix} \mu_a \\ \mu_b \end{bmatrix}, \begin{bmatrix} \Sigma_{aa} & \Sigma_{ab} \\ \Sigma_{ba} & \Sigma_{bb} \end{bmatrix} \right)\]

La distribución de \(\mathbf{a}\) dado \(\mathbf{b}\) es:

\[P(\mathbf{a} \mid \mathbf{b}) = \mathcal{N}(\mu_a + \Sigma_{ab}\Sigma_{bb}^{-1}(\mathbf{b} - \mu_b), \Sigma_{aa} - \Sigma_{ab}\Sigma_{bb}^{-1}\Sigma_{ba})\]

Aplicación a GPR

En nuestro caso, queremos \(P(f_* \mid \mathbf{y})\). Sustituimos:

\(\mathbf{a} \to f_*\) (lo que queremos predecir)
\(\mathbf{b} \to \mathbf{y}\) (lo que observamos)
\(\Sigma_{bb} \to K(X, X) + \sigma_n^2 I\) (covarianza de los datos con ruido)
\(\Sigma_{ab} \to K(x_*, X)\) (correlación entre test y entrenamiento)

Al sustituir en la fórmula del teorema, obtenemos la Media Predictiva:

\[\bar{f}_* = \underbrace{K(x_*, X) K_y^{-1}}_{\text{Pesos de similitud}} \mathbf{y}\]

Intuición: La predicción es una combinación lineal de las etiquetas de entrenamiento \(\mathbf{y}\). Si el nuevo punto \(x_*\) es muy parecido a un punto de entrenamiento \(x_i\), el valor \(k(x_*, x_i)\) será alto, dándole más peso a esa \(y_i\) en el resultado final.

3. ¿Qué significa realmente la Varianza?

La varianza predictiva es:

\[\text{var}(f_*) = k(x_*, x_*) - K(x_*, X) K_y^{-1} K(X, x_*)\]

\(k(x_*, x_*)\): Es la incertidumbre a priori (máxima).
\(K(x_*, X) K_y^{-1} K(X, x_*)\): Es la información ganada al observar los datos.
Conclusión: Cuanto más cerca esté \(x_*\) de los datos observados \(X\), más grande será el sustraendo y menor será la incertidumbre. Si te alejas de los datos, la varianza vuelve a su valor máximo original.

4. Selección de Kernel y Propiedades Históricas

El kernel no es solo una fórmula, codifica nuestras suposiciones sobre el mundo:

Estacionariedad: Si \(k(x, x') = k(x - x')\), el proceso es invariante a traslaciones (la relación solo depende de la distancia).
Diferenciabilidad: El kernel RBF (Exponencial Cuadrática) es infinitamente diferenciable, lo que genera funciones extremadamente suaves. Si los datos son rugosos (como precios de acciones), el RBF suele fallar y se prefiere el kernel Matérn.

El Kernel Matérn (Generalización)

\[k_{Matern}(r) = \frac{2^{1-\nu}}{\Gamma(\nu)} \left( \sqrt{2\nu} \frac{r}{\ell} \right)^\nu K_\nu \left( \sqrt{2\nu} \frac{r}{\ell} \right)\]

Donde \(\nu\) controla la suavidad. Si \(\nu = 1/2\), obtenemos el proceso de Ornstein-Uhlenbeck (muy rugoso, no diferenciable).

5. Implementación Práctica: Estabilidad Numérica

En la teoría usamos \(K^{-1}\), pero en computación nunca se invierte una matriz directamente. Se utiliza la Descomposición de Cholesky:

Hallar \(L\) tal que \(L L^\top = K_y\).
Resolver \(L \boldsymbol{\alpha} = \mathbf{y}\) por sustitución hacia adelante.
La media se vuelve \(\bar{f}_* = K(x_*, X) (L^\top)^{-1} \boldsymbol{\alpha}\).

Esto es mucho más estable y rápido que calcular la inversa directamente, que es propensa a errores de redondeo.

6. “Filosofía GP”

Un GP es como un “filtro de funciones”. Empezamos con una nube infinita de funciones posibles (el prior). Al observar puntos específicos, “aplastamos” esa nube de funciones para que solo pasen las que cruzan (o están cerca) de nuestros datos observados. Lo que queda es el posterior.

Para entender la naturaleza Bayesiana de los Procesos Gaussianos (GP), debemos pasar de la idea de “ajustar una curva” a la idea de “actualizar creencias sobre funciones”.

En la estadística bayesiana tradicional, estimamos parámetros \(\theta\). En GP, realizamos inferencia en espacios de funciones. Aquí te explico la mecánica bayesiana detrás del modelo.

1. El Marco Bayesiano: Prior, Likelihood y Posterior

Aplicamos el Teorema de Bayes directamente sobre la función \(f\):

\[P(f \mid \mathbf{y}, X) = \frac{P(\mathbf{y} \mid f, X) P(f)}{P(\mathbf{y} \mid X)}\]

A. El Prior (\(P(f)\))

Antes de ver cualquier dato, definimos qué tipo de funciones esperamos. Esto se hace mediante el Proceso Gaussiano. Decimos que, para cualquier conjunto de puntos \(X\), la distribución de sus valores de función es:

\[f(X) \sim \mathcal{N}(\mathbf{0}, K(X, X))\]

Aquí, el Kernel \(K\) codifica nuestras creencias previas: ¿Es la función suave? ¿Es periódica? ¿Es ruidosa?

B. El Likelihood (\(P(\mathbf{y} \mid f, X)\))

Es la probabilidad de observar los datos \(\mathbf{y}\) dado que la función real es \(f\). Si asumimos ruido gaussiano blanco \(\epsilon \sim \mathcal{N}(0, \sigma_n^2)\):

\[\mathbf{y} \mid f \sim \mathcal{N}(f(X), \sigma_n^2 I)\]

C. El Posterior (\(P(f \mid \mathbf{y}, X)\))

Es la actualización de nuestra creencia. Al combinar el prior (todas las funciones posibles) con el likelihood (qué tan bien encajan con los datos), el posterior restringe el espacio de funciones a solo aquellas que pasan cerca de los puntos observados.

2. La Verosimilitud Marginal (Evidencia)

En el denominador de Bayes tenemos \(P(\mathbf{y} \mid X)\). En GP, esto se llama Verosimilitud Marginal y es fundamental para el “Bayesianismo” del modelo, ya que permite la selección de modelos sin validación cruzada.

La integral para obtenerla es:

\[P(\mathbf{y} \mid X) = \int P(\mathbf{y} \mid f, X) P(f \mid X) df\]

Como ambas distribuciones son gaussianas, la integral tiene una solución analítica cerrada:

\[\log P(\mathbf{y} \mid X) = \underbrace{-\frac{1}{2} \mathbf{y}^\top K_y^{-1} \mathbf{y}}_{\text{Ajuste de datos}} \underbrace{-\frac{1}{2} \log |K_y|}_{\text{Complejidad (Occam's Razor)}} - \frac{n}{2} \log 2\pi\]

¿Por qué es Bayesiana esta parte?

El segundo término (\(-\frac{1}{2} \log |K_y|\)) penaliza la complejidad. Un kernel muy flexible intentará explicar cada punto de ruido, pero este término aumentará, reduciendo la probabilidad total. El GP elige automáticamente el modelo más simple que explique bien los datos.

3. Inferencia de Jerarquía (Empirical Bayes)

En un flujo de trabajo bayesiano puro, deberíamos poner priors sobre los hiperparámetros del kernel (\(\ell, \sigma_f\)) y usar MCMC. Sin embargo, en GPR solemos usar Bayes Empírico:

Nivel 1: Inferencia sobre la función \(f\) (Exacta mediante matrices).
Nivel 2: Inferencia sobre hiperparámetros \(\theta\) (Maximizando la Verosimilitud Marginal).

4. Comparativa: Bayesiana vs. Frecuentista

Concepto	Regresión Tradicional (Frecuentista)	Regresión GP (Bayesiana)
Resultado	Un solo valor óptimo \(\hat{y}\).	Una distribución completa \(P(f_* \mid \mathbf{y})\).
Incertidumbre	Basada en intervalos de confianza (bootstrapping).	Natural: la varianza crece donde no hay datos.
Parámetros	Los pesos \(\mathbf{w}\) son fijos pero desconocidos.	La función \(f\) es una variable aleatoria.
Sobreajuste	Se controla con regularización (\(L_1, L_2\)).	Se controla con la Verosimilitud Marginal (Navaja de Occam).

Si no hay datos cerca de un punto \(x_*\), el likelihood no tiene fuerza para corregir el prior. Por lo tanto, la predicción “colapsa” hacia la media del prior (usualmente 0) y la incertidumbre vuelve al máximo.

A medida que agregamos datos, la “mancha” de incertidumbre se estrecha alrededor de las observaciones. Esto es la esencia del aprendizaje bayesiano: La transición de la ignorancia (prior) al conocimiento (posterior).

Para elevar el nivel de la teoría, vamos a entrar en la visión funcional profunda. Hasta ahora hemos visto el GP como una “extensión de la normal multivariada”, pero su verdadera potencia reside en el Teorema de Representación de Mercer y en la Dualidad de Pesos.

Teoría Avanzada de Procesos Gaussianos

1. El Espacio de Hilbert Reproducible (RKHS)

La Regresión por Procesos Gaussianos no ocurre en un vacío; vive en un espacio matemático llamado Reproducing Kernel Hilbert Space (RKHS).

La Propiedad de Reproducción

Para cada kernel \(k(\cdot, \cdot)\), existe un espacio único de funciones \(\mathcal{H}\) tal que:

Para todo \(x\), la función \(k(\cdot, x)\) pertenece a \(\mathcal{H}\).
Propiedad de Reproducción: Para cualquier función \(f \in \mathcal{H}\), el valor de la función en un punto \(x\) se puede obtener mediante el producto interno:

\[f(x) = \langle f, k(\cdot, x) \rangle_{\mathcal{H}}\]

Intuición: El kernel actúa como un “evaluador”. La similitud entre dos funciones en este espacio es equivalente a qué tan correlacionadas están según el kernel.

2. Demostración: La Verosimilitud Marginal y la Navaja de Occam

Mencionamos que el GP elige el modelo más simple. Vamos a demostrarlo analíticamente desde la Evidencia del Modelo.

Queremos maximizar \(p(\mathbf{y} \mid X, \theta)\). Al integrar sobre todas las posibles funciones \(f\):

\[p(\mathbf{y} \mid X, \theta) = \int p(\mathbf{y} \mid \mathbf{f}) p(\mathbf{f} \mid X, \theta) d\mathbf{f}\]

Siendo ambas Gaussianas:

\(p(\mathbf{f} \mid X, \theta) = \mathcal{N}(\mathbf{0}, K)\)
\(p(\mathbf{y} \mid \mathbf{f}) = \mathcal{N}(\mathbf{f}, \sigma_n^2 I)\)

La log-verosimilitud marginal resultante es:

\[\mathcal{L}(\theta) = \underbrace{-\frac{1}{2} \mathbf{y}^\top K_y^{-1} \mathbf{y}}_{\text{Ajuste de datos}} - \underbrace{\frac{1}{2} \log |K_y|}_{\text{Complejidad}} - \frac{n}{2} \log 2\pi\]

Análisis de los términos:

Ajuste de datos (\(-\mathbf{y}^\top K_y^{-1} \mathbf{y}\)): Es una medida de qué tan bien la estructura de covarianza explica la varianza observada en \(\mathbf{y}\). Si el kernel es muy rígido y los datos varían mucho, este término será muy negativo.
Penalización de Complejidad (\(\log |K_y|\)): El determinante de la matriz de covarianza mide el “volumen” de funciones que el modelo puede generar.

Si el kernel es muy flexible (muchos grados de libertad), \(|K_y|\) es grande, lo que resta valor a la verosimilitud.
Este término es el que evita automáticamente el sobreajuste sin necesidad de un set de validación.

3. El Concepto de “No Paramétrico”

A menudo se dice que el GP es no paramétrico. Esto es una sutileza técnica:

Modelos Paramétricos: Tienen un número fijo de parámetros \(\mathbf{w}\) (ej. Regresión Lineal). La complejidad es constante independientemente del tamaño de los datos.
Modelos No Paramétricos (GP): El número de “parámetros” crece con el número de datos. En GPR, cada punto de entrenamiento se convierte en un parámetro (un centro de masa para la predicción).

4. Demostración: El Límite de las Redes Neuronales

Uno de los resultados más fascinantes de la teoría moderna (Neal, 1996) es la relación entre GP y Deep Learning:

Teorema: Una red neuronal de una sola capa oculta, con un número infinito de neuronas (\(N \to \infty\)) y cuyos pesos se distribuyen de forma independiente según una Gaussiana, converge exactamente a un Proceso Gaussiano.

¿Por qué importa esto?

Significa que los GPs son el límite teórico de las redes neuronales. Cuando usamos un GP, estamos usando una red neuronal de ancho infinito. Esto ha dado lugar al estudio del Neural Tangent Kernel (NTK) en la IA actual.

5. El Problema de la Inversión (Estabilidad Numérica)

La teoría dice \(K^{-1}\), pero en la práctica usamos la Descomposición de Cholesky.

Demostración de eficiencia: Calcular \(K^{-1}\) directamente es numéricamente inestable si la matriz es casi singular (puntos muy cercanos).

Factorizamos \(K + \sigma_n^2 I = LL^\top\), donde \(L\) es triangular inferior.
La media predictiva \(\bar{f}_* = k_*^\top K_y^{-1} \mathbf{y}\) se resuelve en dos pasos:

Paso 1: Resolver \(L\mathbf{\alpha} = \mathbf{y}\) (Forward substitution).
Paso 2: Resolver \(L^\top \mathbf{\beta} = \mathbf{\alpha}\) (Backward substitution).

\(\bar{f}_* = k_*^\top \mathbf{\beta}\).

Este método reduce el error de punto flotante y es la base de todas las librerías como Scikit-Learn o GPyTorch.

6. Modelado de Ruido No-Gaussiano

¿Qué pasa si mis datos no son Gaussianos (ej. clasificación)? Aquí la teoría se complica porque la integral del posterior ya no tiene solución analítica. Se usan dos métodos:

Aproximación de Laplace: Se hace una expansión de Taylor de segundo orden alrededor del máximo del posterior.
Inferencia Variacional (VI): Se busca la distribución Gaussiana \(q(f)\) que sea más similar (en términos de divergencia KL) al posterior verdadero no-gaussiano.

Para cerrar el ciclo teórico, vamos a profundizar en la arquitectura de kernels y la Inferencia Variacional, que es lo que permite que los Procesos Gaussianos funcionen en el mundo real con millones de datos.

1. Álgebra de Kernels: Construyendo Estructuras Complejas

La teoría de los GPs permite combinar kernels simples para modelar comportamientos complejos. Esto se basa en que la suma y el producto de funciones de covarianza válidas también son kernels válidos.

Suma (\(k_1 + k_2\)): Modela la suma de procesos independientes. Por ejemplo, una tendencia a largo plazo (Kernel Lineal) más variaciones locales (Kernel RBF).
Producto (\(k_1 \times k_2\)): Modela interacciones. Por ejemplo, un kernel RBF multiplicado por uno Periódico resulta en una función que es periódica, pero cuya forma cambia lentamente con el tiempo (decaimiento de la periodicidad).

2. Inferencia Variacional Esparcida (Sparse GP)

El mayor problema teórico del GP es su costo computacional de \(O(n^3)\) debido a la inversión de la matriz. Para datos masivos, usamos Inducing Points (Puntos de Inducción).

La Idea Intuitiva

En lugar de calcular la covarianza entre los \(n\) puntos (donde \(n=1,000,000\)), seleccionamos un conjunto pequeño de \(m\) puntos “maestros” (\(m \ll n\)), llamados puntos de inducción \(\mathbf{u}\). Suponemos que estos puntos son suficientes para resumir toda la forma de la función.

Demostración Matemática (Lower Bound)

Buscamos maximizar el ELBO (Evidence Lower Bound). En lugar de la verosimilitud marginal exacta, optimizamos una aproximación:

\[\mathcal{L}_{ELBO} = \mathbb{E}_{q(\mathbf{f})} [\log p(\mathbf{y} \mid \mathbf{f})] - KL(q(\mathbf{u}) \| p(\mathbf{u}))\]

\(\mathbb{E}[\log p(\mathbf{y} \mid \mathbf{f})]\): Qué tan bien explican los puntos de inducción a los datos observados.
\(KL(q(\mathbf{u}) \| p(\mathbf{u}))\): Una penalización (Divergencia de Kullback-Leibler) que asegura que nuestra aproximación no se aleje demasiado de nuestras creencias previas (el prior).

3. Optimización de Hiperparámetros (Gradientes)

Para hallar la mejor escala de longitud \(\ell\) o la varianza \(\sigma_f\), necesitamos las derivadas de la Log-Verosimilitud Marginal (\(\mathcal{L}\)) respecto a los hiperparámetros \(\theta_j\):

\[\frac{\partial \mathcal{L}}{\partial \theta_j} = \frac{1}{2} \mathbf{y}^\top K^{-1} \frac{\partial K}{\partial \theta_j} K^{-1} \mathbf{y} - \frac{1}{2} \text{tr} \left( K^{-1} \frac{\partial K}{\partial \theta_j} \right)\]

El primer término empuja al modelo a reducir el error de ajuste.
El segundo término (la traza) actúa como un regulador de la complejidad, controlando cuánto “volumen” de funciones estamos permitiendo.

4. Clasificación con Procesos Gaussianos (GPC)

En regresión, la salida es continua. En clasificación, la salida debe estar en \([0, 1]\). Aquí la teoría Bayesiana se vuelve “no analítica”.

El proceso:

Se define una función latente \(f(x)\) mediante un GP.
Se pasa esa función por una función sigmoide o probit: \(\pi(x) = \sigma(f(x))\).
Como la distribución posterior \(p(f \mid \mathbf{y})\) ya no es Gaussiana (la sigmoide rompe la normalidad), usamos la Aproximación de Laplace:

Encontramos el máximo a posteriori (MAP) de \(f\).
Hacemos una expansión de Taylor de segundo orden alrededor de ese máximo para “fingir” que es una Gaussiana y poder operar.

5. Resumen Teórico Final: La Jerarquía de la Inferencia

Para ser un experto en GPR, debes visualizar estos tres niveles de abstracción:

Nivel de la Función: Inferencia de \(f\) (el valor de la curva). Se resuelve con Álgebra Lineal.
Nivel de los Hiperparámetros: Optimización de \(\ell, \sigma\) (la forma de la curva). Se resuelve con Gradiente Descendente sobre la Verosimilitud Marginal.
Nivel del Modelo: Elección del Kernel (¿RBF, Matérn, Lineal?). Se resuelve comparando la Evidencia del Modelo (Bayes Factor).

Dato Curioso: Los Procesos Gaussianos son considerados “Gold Standard” en la optimización de hiperparámetros de otros modelos (como Redes Neuronales Profundas) mediante una técnica llamada Optimización Bayesiana.

# 1. Configuración de datos
set.seed(42)
x_train <- matrix(seq(-5, 5, length.out = 10), ncol = 1)
y_train <- sin(x_train) + rnorm(10, sd = 0.1)
x_test  <- matrix(seq(-6, 6, length.out = 100), ncol = 1)

# 2. Definición del Kernel RBF (Exponencial Cuadrática)
rbf_kernel <- function(x1, x2, l = 1, sigma_f = 1) {
  dist_sq <- outer(as.vector(x1)^2, as.vector(x2)^2, "+") - 2 * x1 %*% t(x2)
  return(sigma_f^2 * exp(-0.5 * dist_sq / l^2))
}

# 3. Parámetros
l <- 1          # Escala de longitud
sigma_f <- 1    # Amplitud
sigma_n <- 0.1  # Ruido observado

# 4. Cálculo de Matrices de Covarianza
K   <- rbf_kernel(x_train, x_train, l, sigma_f) + sigma_n^2 * diag(nrow(x_train))
Ks  <- rbf_kernel(x_train, x_test, l, sigma_f)
Kss <- rbf_kernel(x_test, x_test, l, sigma_f)

# 5. Inferencia Bayesiana (Media y Varianza Posterior)
# Invertimos K (usando solve para este ejemplo didáctico)
K_inv <- solve(K)
mu_s <- t(Ks) %*% K_inv %*% y_train
cov_s <- Kss - t(Ks) %*% K_inv %*% Ks

# 6. Preparación de intervalos de confianza (95%)
std_dev <- sqrt(diag(cov_s))
upper <- mu_s + 1.96 * std_dev
lower <- mu_s - 1.96 * std_dev

# 7. Visualización
plot(x_train, y_train, pch = 19, col = "red", xlim = c(-6, 6), ylim = c(-2, 2),
     main = "GPR desde Cero en R", xlab = "x", ylab = "f(x)")
lines(x_test, mu_s, col = "blue", lwd = 2)
polygon(c(x_test, rev(x_test)), c(upper, rev(lower)), 
        col = rgb(0, 0, 1, 0.2), border = NA)
grid()

# Instalar si es necesario: install.packages("kernlab")
library(kernlab)

# 1. Ajuste del modelo
# gausspr utiliza optimización interna para los parámetros del kernel
modelo_gp <- gausspr(x_train, y_train, kernel = "rbfdot", kpar = list(sigma = 0.5), variance.model = TRUE)

# 2. Predicción
predicciones <- predict(modelo_gp, x_test)

# 3. Obtener la varianza (incertidumbre)
# kernlab permite obtener el intervalo de predicción
incertidumbre <- predict(modelo_gp, x_test, type = "probabilities") 

# Visualización rápida con ggplot2
library(ggplot2)

## 
## Adjuntando el paquete: 'ggplot2'

## The following object is masked from 'package:kernlab':
## 
##     alpha

df_plot <- data.frame(
  x = as.vector(x_test),
  y = as.vector(predicciones),
  lower = as.vector(predicciones - 1.96 * 0.2), # Aproximación visual
  upper = as.vector(predicciones + 1.96 * 0.2)
)

ggplot(df_plot, aes(x = x, y = y)) +
  geom_ribbon(aes(ymin = lower, ymax = upper), fill = "steelblue", alpha = 0.3) +
  geom_line(color = "darkblue", size = 1) +
  geom_point(data = data.frame(x=x_train, y=y_train), aes(x=x, y=y), color = "red", size = 3) +
  theme_minimal() +
  labs(title = "Regresión por Proceso Gaussiano (kernlab)",
       subtitle = "El área sombreada representa la incertidumbre epistémica")

## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once per session.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

# Función para calcular la Log-Verosimilitud Marginal Negativa
nlogmag_lik <- function(params, x, y) {
  l <- exp(params[1])       # Usamos exp para asegurar valores positivos
  sigma_f <- exp(params[2])
  sigma_n <- 0.1
  
  K <- rbf_kernel(x, x, l, sigma_f) + sigma_n^2 * diag(nrow(x))
  
  # Log-verosimilitud (fórmula teórica vista anteriormente)
  # -1/2 * y' * K^-1 * y - 1/2 * log|K| - n/2 * log(2pi)
  term1 <- 0.5 * t(y) %*% solve(K) %*% y
  term2 <- 0.5 * determinant(K, logarithm = TRUE)$modulus
  term3 <- (nrow(x) / 2) * log(2 * pi)
  
  return(as.numeric(term1 + term2 + term3))
}

# Optimización usando optim()
opt <- optim(par = c(0, 0), fn = nlogmag_lik, x = x_train, y = y_train)
best_params <- exp(opt$par)

print(paste("Mejor escala de longitud (l):", round(best_params[1], 4)))

## [1] "Mejor escala de longitud (l): 1.9169"

print(paste("Mejor amplitud (sigma_f):", round(best_params[2], 4)))

## [1] "Mejor amplitud (sigma_f): 0.9775"

# install.packages("DiceKriging")
library(DiceKriging)

## Warning: package 'DiceKriging' was built under R version 4.5.3

# Creamos un kernel: Exponencial Cuadrático + Lineal
# Esto permite capturar la curvatura y la tendencia general
modelo_compuesto <- km(design = data.frame(x = x_train), 
                       response = y_train, 
                       covtype = "matern5_2", # Un kernel más robusto que el RBF
                       iso = TRUE)

## 
## optimisation start
## ------------------
## * estimation method   : MLE 
## * optimisation method : BFGS 
## * analytical gradient : used
## * trend model : ~1
## * covariance model : 
##   - type :  matern5_2 
##   - nugget : NO
##   - parameters lower bounds :  1e-10 
##   - parameters upper bounds :  20 
##   - best initial criterion value(s) :  -6.753224 
## 
## N = 1, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate     0  f=       6.7532  |proj g|=      0.44523
## At iterate     1  f =       6.6943  |proj g|=       0.22952
## At iterate     2  f =       6.6799  |proj g|=       0.03329
## At iterate     3  f =       6.6795  |proj g|=     0.0019416
## At iterate     4  f =       6.6795  |proj g|=    1.8439e-05
## At iterate     5  f =       6.6795  |proj g|=    1.0033e-08
## 
## iterations 5
## function evaluations 6
## segments explored during Cauchy searches 5
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 1.00325e-08
## final function value 6.67952
## 
## F = 6.67952
## final  value 6.679521 
## converged

# Predicción con el modelo compuesto
p <- predict(modelo_compuesto, newdata = data.frame(x = x_test), type = "UK")

# Graficar resultados
plot(x_test, p$mean, type = "l", col = "blue", lwd = 2, ylim = c(-2, 2))
lines(x_test, p$lower95, col = "gray", lty = 2)
lines(x_test, p$upper95, col = "gray", lty = 2)
points(x_train, y_train, col = "red", pch = 19)

# Función de Adquisición: Expected Improvement (EI)
calc_ei <- function(x_new, model, y_best) {
  p <- predict(model, newdata = data.frame(x = x_new), type = "UK")
  mu <- p$mean
  s <- p$sd
  
  z <- (mu - y_best) / s
  ei <- (mu - y_best) * pnorm(z) + s * dnorm(z)
  return(ei)
}

# 1. Encontrar el mejor valor observado hasta ahora
y_best <- max(y_train)

# 2. Calcular EI para todos los puntos potenciales (x_test)
ei_values <- sapply(x_test, calc_ei, model = modelo_compuesto, y_best = y_best)

# 3. El siguiente punto a experimentar es el que maximiza el EI
next_point <- x_test[which.max(ei_values)]

# Visualización del EI (Donde el modelo sugiere explorar)
plot(x_test, ei_values, type = "l", col = "darkgreen", main = "Función de Adquisición (EI)")
abline(v = next_point, col = "orange", lty = 3)

Regresión por Procesos Gaussianos (GPR)

Jhonier Rangel

2026-05-08

Regresión por Procesos Gaussianos (GPR)

1. Definición Formal

2. El Kernel (Función de Covarianza)

3. Demostración: Predicciones con GPR

A. La Distribución Conjunta

B. Condicionamiento (Inferencia)

4. Optimización: Máxima Verosimilitud Marginal

5. Ventajas y Desventajas

Algoritmo

1. La Intuición: Del Espacio de Pesos al Espacio de Funciones

La “Vista de Pesos” (Interpretación Dual)

2. Demostración: La Distribución Condicional Posterior

El Teorema de la Normal Condicional

Aplicación a GPR

3. ¿Qué significa realmente la Varianza?

4. Selección de Kernel y Propiedades Históricas

El Kernel Matérn (Generalización)

5. Implementación Práctica: Estabilidad Numérica

6. “Filosofía GP”

1. El Marco Bayesiano: Prior, Likelihood y Posterior

A. El Prior (\(P(f)\))

B. El Likelihood (\(P(\mathbf{y} \mid f, X)\))

C. El Posterior (\(P(f \mid \mathbf{y}, X)\))

2. La Verosimilitud Marginal (Evidencia)

¿Por qué es Bayesiana esta parte?

3. Inferencia de Jerarquía (Empirical Bayes)

4. Comparativa: Bayesiana vs. Frecuentista

Teoría Avanzada de Procesos Gaussianos

1. El Espacio de Hilbert Reproducible (RKHS)

La Propiedad de Reproducción

2. Demostración: La Verosimilitud Marginal y la Navaja de Occam

Análisis de los términos:

3. El Concepto de “No Paramétrico”

4. Demostración: El Límite de las Redes Neuronales

¿Por qué importa esto?

5. El Problema de la Inversión (Estabilidad Numérica)

6. Modelado de Ruido No-Gaussiano

1. Álgebra de Kernels: Construyendo Estructuras Complejas

2. Inferencia Variacional Esparcida (Sparse GP)

La Idea Intuitiva

Demostración Matemática (Lower Bound)

3. Optimización de Hiperparámetros (Gradientes)

4. Clasificación con Procesos Gaussianos (GPC)

El proceso:

5. Resumen Teórico Final: La Jerarquía de la Inferencia