Presentan
Georgina Zoroastro Ramírez García
Daniela Márquez Morales
José Carlos Iñiguez Álvarez

Preliminares

Modelos de regresión en espacios de características

Los modelos de regresión lineal nos dan el primer ejemplo de lo que definiremos como procesos Gaussianos; en particular nos proveen la motivación para conceptos como los kerneles y sus efectos en los distintos modelos. Para comenzar, consideremos una colección de puntos \((\mathbf{x}_i, y_i)\), \(i = 1,\dots,n\), donde \(\mathbf{x}_i \in \mathbb{R}^D\) es un vector de entrada (input) y \(y_i \in \mathbb{R}\) es su correspondiente valor escalar de salida (output). Si \(X \in \mathbb{R}^{D\times n}\) es la matriz formada por los \(n\) inputs y \(\mathbf{y}\) es el vector de outputs, definimos el conjunto de entrenamiento como

\[\mathcal{D} = (X, \mathbf{y}) = \left\{(\mathbf{x}_i, y_i)\;\vert\; i = 1,\dots,n\right\}.\]

El objetivo es realizar inferencia sobre la relación entre inputs y outputs utilizando la distribución condicional \(p(\mathbf{y} \mid X)\); la inferencia no se realiza directamente sobre los inputs.

Un breve recordatorio a la inferencia Bayesiana

El enfoque que utilizaremos para nuestros modelos de regresión lineal (y en particular para la regresión en procesos Gaussianos) es el de la inferencia Bayesiana. Para utilizarla debemos considerar la verosimilitud del modelo y suponer que nuestros parámetros se distribuyen de cierta manera; esta distribución se denomina a priori y nuestro objetivo es encontrar la distribución a posteriori.

Si nuestro vector de parámetros es \(\theta\) y nuestra verosimilitud la denotamos como \(p(X \mid \theta)\), nuestra distribución a priori es \(p(\theta)\). La distribución a posteriori se obtiene mediante el teorema de Bayes:

\[p(\theta \mid X) = \frac{p(X \mid \theta)\; p(\theta)}{p(X)},\]

donde \(p(X) = \int p(X \mid \theta)\, p(\theta)\, d\theta\) es la constante de normalización. Podemos prescindir de este término durante la inferencia, lo cual nos permite escribir:

\[p(\theta \mid X) \propto p(X \mid \theta)\; p(\theta).\]

Una de las principales características de la inferencia Bayesiana es su capacidad para modelar la incertidumbre del modelo. Podemos utilizar la distribución a posteriori para cuantificar el comportamiento de datos que no están en el conjunto de entrenamiento mediante la distribución predictiva.

Distribución predictiva a priori

La distribución predictiva a priori marginaliza sobre los parámetros sin considerar las observaciones:

\[p(x_\text{new}) = \int p(x_\text{new} \mid \theta)\; p(\theta)\; d\theta.\]

Esta integral se interpreta como el promedio ponderado de la verosimilitud \(p(x_\text{new} \mid \theta)\) bajo la a priori \(p(\theta)\): refleja lo que esperamos observar antes de ver cualquier dato.

Distribución predictiva a posteriori

La distribución predictiva a posteriori incorpora la información de las observaciones \(\mathbf{x}\) actualizando la a priori por la a posteriori:

\[p(x_\text{new} \mid \mathbf{x}) = \int p(x_\text{new} \mid \theta)\; p(\theta \mid \mathbf{x})\; d\theta.\]

Para ver por qué esta integral es la elección natural, notemos que podemos escribirla como una marginalización conjunta:

\[p(x_\text{new} \mid \mathbf{x}) = \int p(x_\text{new}, \theta \mid \mathbf{x})\; d\theta = \int p(x_\text{new} \mid \theta,\mathbf{x})\; p(\theta \mid \mathbf{x})\; d\theta.\]

Si asumimos que \(x_\text{new}\) es condicionalmente independiente de \(\mathbf{x}\) dado \(\theta\) (es decir, toda la información relevante sobre \(x_\text{new}\) está capturada por \(\theta\)), entonces \(p(x_\text{new} \mid \theta, \mathbf{x}) = p(x_\text{new} \mid \theta)\) y recuperamos la expresión anterior.

Caso gaussiano

Cuando tanto la verosimilitud como la a priori son gaussianas, la integral anterior tiene solución analítica. Supongamos:

\[\theta \sim \mathcal{N}(\mu_0, \Sigma_0), \qquad x_\text{new} \mid \theta \sim \mathcal{N}(A\theta, \Sigma_\varepsilon).\]

Entonces la distribución predictiva a priori es también gaussiana:

\[p(x_\text{new}) = \mathcal{N}(x_\text{new} \mid A\mu_0,\; A\Sigma_0 A^T + \Sigma_\varepsilon).\]

De manera análoga, una vez observados los datos \(\mathbf{x}\) y obtenida la a posteriori \(p(\theta \mid \mathbf{x}) = \mathcal{N}(\theta \mid \mu_n, \Sigma_n)\), la distribución predictiva a posteriori es:

\[p(x_\text{new} \mid \mathbf{x}) = \mathcal{N}(x_\text{new} \mid A\mu_n,\; A\Sigma_n A^T + \Sigma_\varepsilon).\]

La estructura es idéntica en ambos casos: la media predictiva es la transformación lineal de la media de \(\theta\), y la varianza predictiva suma la incertidumbre sobre \(\theta\) (propagada por \(A\)) más el ruido irreducible \(\Sigma_\varepsilon\). Este resultado es el que usaremos directamente en el modelo de regresión lineal Bayesiano, donde \(A = \phi(x_{test})^T\) y \(\theta = w\). #### Intervalos de la distribución predictiva

Dado que la distribución predictiva a posteriori es gaussiana, podemos construir intervalos de credibilidad de manera directa. Para un punto nuevo \(x_\text{new}\), la media y varianza predictivas son:

\[\mu_\text{pred} = A\mu_n, \qquad \sigma^2_\text{pred} = A\Sigma_n A^T + \Sigma_\varepsilon.\]

El intervalo de credibilidad al \(100(1-\alpha)\%\) es:

\[\mu_\text{pred} \;\pm\; z_{\alpha/2}\;\sigma_\text{pred},\]

donde \(z_{\alpha/2}\) es el cuantil correspondiente de la distribución normal estándar. Para \(\alpha = 0.05\) se tiene \(z_{0.025} = 1.96\), lo que da el intervalo de credibilidad al 95%:

\[A\mu_n \;\pm\; 1.96\;\sqrt{A\Sigma_n A^T + \Sigma_\varepsilon}.\]

Es importante distinguir este intervalo del intervalo de confianza frecuentista:

  • El intervalo de credibilidad afirma que, dado los datos observados, \(x_\text{new}\) cae en el intervalo con probabilidad \(1 - \alpha\). Es una declaración probabilística sobre \(x_\text{new}\).

  • El intervalo de confianza frecuentista afirma que, bajo muestreo repetido, el \(100(1-\alpha)\%\) de los intervalos construidos de esta forma contendrían el valor verdadero. No es una declaración probabilística sobre \(x_\text{new}\) dado este experimento específico.

La varianza predictiva \(\sigma^2_\text{pred} = A\Sigma_n A^T + \Sigma_\varepsilon\) tiene además una descomposición interpretable en dos fuentes de incertidumbre:

  • \(A\Sigma_n A^T\): incertidumbre epistémica, debida al desconocimiento de \(\theta\). Se reduce conforme se agregan más datos, ya que \(\Sigma_n \to 0\) cuando \(n \to \infty\).

  • \(\Sigma_\varepsilon\): incertidumbre aleatoria (ruido irreducible), inherente al proceso generador de datos. No se reduce con más observaciones.

Esta distinción será relevante en el modelo de regresión lineal Bayesiano, donde \(\sigma^2_n\) juega el papel de \(\Sigma_\varepsilon\) y la incertidumbre epistémica queda capturada por \(\phi(x_{test})^T A^{-1} \phi(x_{test})\). ### Función de características y modelo lineal

En la regresión lineal estándar consideramos el modelo \(y = w^T x + \varepsilon\), donde \(w\) es el vector de pesos que buscamos estimar. Sin embargo, podemos generalizar este modelo trabajando en espacios más expresivos. En lugar de usar directamente \(\mathbf{x}\), consideramos una función de características \(\phi\colon \mathbb{R}^D \to \mathbb{R}^N\) (con \(N \geq D\), inclusive \(N \gg D\)), que transforma cada input a un espacio de mayor dimensión. Algunos ejemplos son:

  • \(\phi(\mathbf{x}) = \mathbf{x}\) (identidad, recupera la regresión lineal estándar),
  • \(\phi(x) = (1, x, \dots, x^N)^T\) (características polinomiales),
  • \(\phi(x) = \frac{1}{\sqrt{N}}\bigl(1, \cos(\omega x), \sin(\omega x), \dots, \cos(N\omega x), \sin(N\omega x)\bigr)^T\) (bases de Fourier).

El modelo de regresión lineal en el espacio de características queda definido como

\[y = w^T\phi(\mathbf{x}) + \varepsilon, \qquad \varepsilon \sim \mathcal{N}(0, \sigma_n^2),\]

donde \(w \in \mathbb{R}^N\) es el vector de pesos del modelo. La matriz de diseño se define como \(\Phi = \Phi(X) \in \mathbb{R}^{n \times N}\), cuya \(i\)-ésima fila es \(\phi(\mathbf{x}_i)^T\).

Verosimilitud

Dado \(w\), las observaciones son independientes y la verosimilitud del vector de outputs es

\[p(\mathbf{y} \mid X, w) = \prod_{i=1}^n p(y_i \mid \mathbf{x}_i, w) = \mathcal{N}(\mathbf{y} \mid \Phi w,\, \sigma_n^2 I).\]

Enfoque Bayesiano

Distribución a priori y a posteriori

A diferencia del enfoque frecuentista, que produce una estimación puntual \(\hat{w}\), el enfoque Bayesiano trata \(w\) como una variable aleatoria y le asigna una distribución a priori. Elegimos una a priori gaussiana

\[p(w) = \mathcal{N}(w \mid 0, \Sigma_N),\]

cuya elección no es arbitraria: la familia gaussiana es conjugada con la verosimilitud gaussiana, lo que garantiza que la distribución a posteriori también sea gaussiana y tenga forma analítica. Aplicando el teorema de Bayes,

\[p(w \mid \mathbf{y}, X) \propto p(\mathbf{y} \mid X, w)\, p(w),\]

Derivación analítica de la distribución a posteriori

Para demostrar formalmente que la distribución resultante es una normal multivariada, trabajamos con el logaritmo de la densidad y completamos el cuadrado. Las densidades del modelo son:

\[p(w) = \frac{1}{(2\pi)^{N/2} |\Sigma_N|^{1/2}} \exp \left( -\frac{1}{2} w^T \Sigma_N^{-1} w \right),\]

\[p(\mathbf{y} \mid X, w) = \frac{1}{(2\pi \sigma_n^2)^{n/2}} \exp \left( -\frac{1}{2\sigma_n^2} (\mathbf{y} - \Phi w)^T (\mathbf{y} - \Phi w) \right).\]

Tomando el logaritmo y omitiendo constantes que no dependen de \(w\):

\[\ln p(w \mid \mathbf{y}, X) = -\frac{1}{2} \left[ \frac{1}{\sigma_n^2} (\mathbf{y} - \Phi w)^T (\mathbf{y} - \Phi w) + w^T \Sigma_N^{-1} w \right] + \text{cte.}\]

Expandiendo el término cuadrático de la verosimilitud y usando que \(\mathbf{y}^T \Phi w = w^T \Phi^T \mathbf{y}\) (escalar igual a su transpuesta):

\[\ln p(w \mid \mathbf{y}, X) = -\frac{1}{2} \left[ w^T \left( \sigma_n^{-2} \Phi^T \Phi + \Sigma_N^{-1} \right) w - 2 w^T \left( \sigma_n^{-2} \Phi^T \mathbf{y} \right) \right] + \text{cte.}\]

Igualando con la estructura de la log-densidad de una Gaussiana \(\mathcal{N}(\mu_w, A^{-1})\), es decir \(-\frac{1}{2}\bigl[w^T A w - 2w^T A\mu_w\bigr] + \text{cte.}\), deducimos por comparación de coeficientes:

  1. Matriz de precisión: \(A = \sigma_n^{-2}\Phi^T\Phi + \Sigma_N^{-1}\),
  2. Vector de medias: \(\mu_w = \sigma_n^{-2} A^{-1}\Phi^T\mathbf{y}\).

Por lo tanto la distribución a posteriori es:

\[p(w \mid \mathbf{y}, X) = \mathcal{N}(w \mid \mu_w,\, A^{-1}).\]

La diferencia conceptual con el enfoque frecuentista es ilustrativa: mientras OLS produce una única recta ajustada, la distribución a posteriori Bayesiana representa una distribución sobre funciones, cada una compatible con los datos observados en distinta medida. La siguiente figura ilustra esta diferencia.

Comparación entre el enfoque frecuentista (OLS) y el Bayesiano. A la izquierda, una única estimación puntual con banda de confianza al 95\%. A la derecha, 300 funciones muestreadas de la distribución a posteriori.

Comparación entre el enfoque frecuentista (OLS) y el Bayesiano. A la izquierda, una única estimación puntual con banda de confianza al 95%. A la derecha, 300 funciones muestreadas de la distribución a posteriori.

Comparación de bandas de incertidumbre

Una pregunta natural es si las bandas de incertidumbre de ambos enfoques son comparables. La respuesta es que se pueden visualizar juntas, pero miden cosas conceptualmente distintas:

  • Banda frecuentista (intervalo de confianza): mide la incertidumbre sobre \(\hat{w}\) bajo muestreo repetido. Si repitieras el experimento muchas veces, el 95% de los intervalos construidos contendrían el valor verdadero de \(f(x)\). No es una declaración probabilística sobre \(f(x)\) dado este conjunto de datos.

  • Banda Bayesiana (intervalo de credibilidad): mide la incertidumbre a posteriori sobre \(f(x_{test})\) dados los datos observados. Es una declaración probabilística directa: \(f(x_{test})\) cae en este intervalo con probabilidad 0.95 dado lo observado.

Con pocos datos, la banda Bayesiana suele ser más ancha porque incorpora la incertidumbre de la a priori \(\Sigma_N\). Con muchos datos, ambas bandas convergen al mismo resultado, pues la verosimilitud domina a la a priori.

Comparación de bandas de incertidumbre al 95\%. Verde: intervalo de confianza frecuentista. Morado: intervalo de credibilidad Bayesiano.

Comparación de bandas de incertidumbre al 95%. Verde: intervalo de confianza frecuentista. Morado: intervalo de credibilidad Bayesiano.

Distribución predictiva y el puente hacia los Procesos Gaussianos

Predicción en un punto nuevo

Dado un punto de prueba \(x_{test}\), queremos inferir \(f_{test} = w^T\phi(x_{test})\) marginalizando sobre la incertidumbre en \(w\):

\[p(f_{test} \mid x_{test}, X, \mathbf{y}) = \int p(f_{test} \mid x_{test}, w)\; p(w \mid X, \mathbf{y})\; dw.\]

Como \(f_{test}\) es una transformación lineal de \(w\) y la distribución a posteriori es gaussiana, \(f_{test}\) también se distribuye de forma normal:

\[f_{test} \mid x_{test}, X, \mathbf{y} \;\sim\; \mathcal{N}\!\left(\mu_w^T\phi(x_{test}),\; \phi(x_{test})^T A^{-1} \phi(x_{test})\right).\]

El intervalo de credibilidad al 95% para \(f_{test}\) es entonces:

\[\mu_w^T\phi(x_{test}) \;\pm\; 1.96\;\sqrt{\phi(x_{test})^T A^{-1} \phi(x_{test})}.\]

La media \(\mu_w^T\phi(x_{test})\) es la predicción puntual óptima bajo pérdida cuadrática, mientras que la varianza \(\phi(x_{test})^T A^{-1} \phi(x_{test})\) cuantifica la incertidumbre a posteriori en \(x_{test}\): es grande cuando \(x_{test}\) está lejos de los datos de entrenamiento y pequeña cuando hay muchas observaciones cercanas.

La identidad de Woodbury y la derivación del kernel

El problema con la expresión anterior es que \(A \in \mathbb{R}^{N\times N}\), cuya inversión cuesta \(\mathcal{O}(N^3)\). Cuando \(N \to \infty\) esto es completamente intratable. La identidad de Woodbury permite reescribir:

\[A^{-1} = \Sigma_N - \Sigma_N \Phi^T \!\left(\Phi\Sigma_N\Phi^T + \sigma_n^2 I\right)^{-1}\!\Phi\Sigma_N.\]

Observemos que \(A^{-1}\) sigue viviendo en \(\mathbb{R}^{N\times N}\): Woodbury no elimina la dependencia en \(N\) por sí sola. Lo que sí hace es reorganizar los términos de forma que, al sustituir en la distribución predictiva, todo se expresa en términos de productos de la forma \(\phi(x)^T\Sigma_N\phi(x')\). Veamos esto explícitamente.

Derivación de la media predictiva

Sustituyendo \(A^{-1}\) en la media \(\mu_w^T\phi(x_{test})\), con \(\mu_w = \sigma_n^{-2}A^{-1}\Phi^T\mathbf{y}\):

\[\mathbb{E}[f_{test}] = \phi(x_{test})^T \mu_w = \sigma_n^{-2}\phi(x_{test})^T A^{-1} \Phi^T \mathbf{y}.\]

Aplicando Woodbury:

\[\mathbb{E}[f_{test}] = \sigma_n^{-2}\phi(x_{test})^T \left[\Sigma_N - \Sigma_N\Phi^T\!\left(\Phi\Sigma_N\Phi^T + \sigma_n^2 I\right)^{-1}\!\Phi\Sigma_N\right] \Phi^T \mathbf{y}.\]

Distribuyendo y factorizando, los términos se agrupan como:

\[\mathbb{E}[f_{test}] = \underbrace{\phi(x_{test})^T\Sigma_N\Phi^T}_{= \mathbf{k}_*^T} \left[\sigma_n^{-2}I - \sigma_n^{-2}\left(\Phi\Sigma_N\Phi^T + \sigma_n^2 I\right)^{-1}\Phi\Sigma_N\Phi^T\sigma_n^{-2}\cdots\right]\mathbf{y},\]

que tras simplificación algebraica (usando la identidad de Sherman-Morrison sobre el término entre corchetes) resulta en:

\[\mathbb{E}[f_{test}] = \mathbf{k}_*^T \left(K + \sigma_n^2 I\right)^{-1} \mathbf{y},\]

donde hemos definido:

\[(\mathbf{k}_*)_i = \phi(x_i)^T \Sigma_N\, \phi(x_{test}) = k(x_i, x_{test}), \qquad K_{ij} = \phi(x_i)^T\Sigma_N\,\phi(x_j) = k(x_i, x_j).\]

Derivación de la varianza predictiva

Para la varianza, sustituimos \(A^{-1}\) directamente en \(\phi(x_{test})^T A^{-1} \phi(x_{test})\):

\[\mathbb{V}[f_{test}] = \phi(x_{test})^T \left[\Sigma_N - \Sigma_N\Phi^T\!\left(\Phi\Sigma_N\Phi^T + \sigma_n^2 I\right)^{-1}\!\Phi\Sigma_N\right] \phi(x_{test}).\]

Distribuyendo los dos términos:

\[\mathbb{V}[f_{test}] = \underbrace{\phi(x_{test})^T\Sigma_N\phi(x_{test})}_{= k(x_{test}, x_{test})\;=\;k_{**}} - \underbrace{\phi(x_{test})^T\Sigma_N\Phi^T}_{= \mathbf{k}_*^T} \left(\Phi\Sigma_N\Phi^T + \sigma_n^2 I\right)^{-1} \underbrace{\Phi\Sigma_N\phi(x_{test})}_{= \mathbf{k}_*}.\]

Es decir:

\[\mathbb{V}[f_{test}] = k_{**} - \mathbf{k}_*^T\left(K + \sigma_n^2 I\right)^{-1}\mathbf{k}_*.\]

Resultado: \(N\) desaparece

En ambas expresiones, \(\phi\) y \(\Sigma_N\) han quedado completamente absorbidos dentro del kernel \(k(x, x') = \phi(x)^T\Sigma_N\phi(x')\). La dimensión \(N\) ya no aparece: \(K \in \mathbb{R}^{n\times n}\) y \(\mathbf{k}_* \in \mathbb{R}^n\), por lo que la inversión cuesta \(\mathcal{O}(n^3)\) independientemente de \(N\).

Distribución predictiva final e intervalos

La distribución predictiva completa en términos del kernel es:

\[p(f_{test} \mid x_{test}, X, \mathbf{y}) = \mathcal{N}\!\left(\mathbf{k}_*^T(K + \sigma_n^2 I)^{-1}\mathbf{y},\;\; k_{**} - \mathbf{k}_*^T(K + \sigma_n^2 I)^{-1}\mathbf{k}_*\right),\]

donde \(k_{**} = k(x_{test}, x_{test})\). El intervalo de credibilidad al 95% en términos del kernel es entonces:

\[\mathbf{k}_*^T(K + \sigma_n^2 I)^{-1}\mathbf{y} \;\pm\; 1.96\;\sqrt{k_{**} - \mathbf{k}_*^T(K + \sigma_n^2 I)^{-1}\mathbf{k}_*}.\]

La varianza predictiva \(k_{**} - \mathbf{k}_*^T(K+\sigma_n^2 I)^{-1}\mathbf{k}_*\) tiene una interpretación intuitiva: empieza con la incertidumbre a priori \(k_{**}\) sobre \(x_{test}\) y la reduce en función de cuánta información aportan los datos de entrenamiento sobre ese punto a través de \(\mathbf{k}_*^T(K+\sigma_n^2 I)^{-1}\mathbf{k}_*\).

El puente

El resultado anterior tiene una consecuencia fundamental: dado que \(\phi\) solo aparece a través de productos de la forma \(k(x, x') = \phi(x)^T\Sigma_N\phi(x')\), podemos reemplazar \(\phi\) por cualquier función kernel válida \(k\colon\mathbb{R}^D\times\mathbb{R}^D\to\mathbb{R}\) sin necesidad de conocerla explícitamente. Cuando \(\dim(\phi) \to \infty\), trabajar con \(w\) es intratable, pero el kernel sigue siendo computable y las ecuaciones predictivas anteriores siguen siendo exactas. Este es precisamente el argumento que motiva los Procesos Gaussianos: en lugar de poner una a priori sobre un vector finito de parámetros \(w\), se pone directamente una a priori sobre funciones, con función de covarianza \(k(x, x')\):

\[\underbrace{w^T\phi(\mathbf{x})}_{\dim(\phi) < \infty} \;\xrightarrow{\;\dim(\phi)\to\infty\;}\; \underbrace{\text{Proceso Gaussiano}}_{\text{\textit{a priori} sobre funciones}.}\]

Regresión de Procesos Gaussianos

Ya analizamos la regresión del modelo lineal bajo el enfoque bayesiano, pero realmento lo que nos interesa es liberarnos de las limitaciones del modelo lineal y, en general, de cualquier modelo paramétrico, ya que la calidad de sus predicciones depende de la elección del modelo elegido. Para ello, el conocimiento previo se codifica no en la distribución de los parámetro de algún modelo paramétrico, sino directamente sobre la distribución del espacio de funciones posibles, es decir, en la distribución de $f(x)$ para cada $x$; lo cual nos conduce a la regresión de procesos estocásticos. Aquí consideramos únicamente a los procesos gaussianos (GP).

Un proceso gaussiano se define como una colección de variables aleatorias tales que cualquier subconjunto finito de ellas tiene una distribución gaussiana conjunta. Tal proceso queda completamente determinado por una función media:

\[ m(\mathbf{x}) = \mathbb{E}[f(\mathbf{x})], \]

y una función de covarianza:\[ k(\mathbf{x},\mathbf{x}') = \mathbb{E} \left[ (f(\mathbf{x})-m(\mathbf{x})) (f(\mathbf{x}')-m(\mathbf{x}')) \right]. \]

Si $f(x)$ es un proceso gaussiano lo escribimos como\[ f(x)\sim \mathcal{GP}(m(x),k(x,x')) \]

Notemos que un proceso gaussiano satisface la propiedad fundamental de la consistencia marginal, que dice que examinar un conjunto de variables no debe cambiar la distribución marginal de un subconjunto menor. Esta propiedad nos será de mucha utilidad a la hora de hacer la regresión de GP.

Entonces, queremos obtener una distribución a posteriori de la función latente \(f(x)\). Para ello tenemos un conjunto de entrenamiento con $n$ observaciones:

$$

\mathcal{D}=(X_{tr},y_{tr})=\{(x_{tr,i},y_{tr,i})\mid i=1,\ldots,n\}.

$$

Aquí $y$ es la variable de respuesta, y $x$ es el vector de variables predictoras de dimensión $D$. La matriz de diseño se denota por $X$.

Consideremos primero el caso en el que los outputs están libres de ruido:\[ y_i=f(x_i) \]

Usamos el hecho de que todo subconjunto finito de valores de un proceso gaussiano tiene distribución normal multivariada y, por lo tanto, la priori conjunta del conjunto de entrenamiento y el conjunto de prueba es\[ \begin{bmatrix}f_{tr}\\f_\text{ts}\end{bmatrix}\sim\mathcal{N}\left(0,\begin{bmatrix}K(X_{tr},X_{tr}) & K(X_{tr},X_\text{ts})\\K(X_\text{ts},X_{tr}) & K(X_\text{ts},X_\text{ts})\end{bmatrix}\right), \]

donde \[ f_{tr} =\begin{bmatrix}f(x_{tr,1})\\\vdots\\f(x_{tr,n})\end{bmatrix},\qquad f_\text{ts} =\begin{bmatrix}f(x_{\text{ts}1})\\\vdots\\f(x_{\text{ts}n})\end{bmatrix}, \]

y donde \(K\) es el kernel que representa la matriz de covarianza de sus entradas.

Usamos el resultado del condicionamiento gaussiano: si dos vectores, \(\alpha\) y \(\beta\), siguen una distribución normal multivariada tal que

\[ \begin{bmatrix}\mathbf{\alpha} \\ \mathbf{\beta} \end{bmatrix} \sim\mathcal{N} \left( \begin{bmatrix} \boldsymbol{\mu}_{\alpha} \\ \boldsymbol{\mu}_{\beta} \end{bmatrix}, \begin{bmatrix} A & C \\C^\top & B\end{bmatrix} \right), \]

la distribución condicional de \(\alpha\) dada \(\beta\) es

\[ \alpha\mid \beta\sim\mathcal{N}\left(\boldsymbol{\mu}_{\alpha}+CB^{-1}(\beta-\boldsymbol{\mu}_{\beta}),\,A-CB^{-1}C^\top\right) \]

Aplicándolo a nuestro problema, tenemos que

\[ f_\text{ts} \mid X_\text{ts},X_{tr},f_{tr}\sim\mathcal{N}\left(\bar f_\text{ts},\operatorname{cov}(f_\text{ts})\right), \]

donde la media predictiva es

\[ \bar {f_\text{ts}} = \mathbb{E}[f_\text{ts} \mid X_\text{ts},X_{tr},f_{tr}] = K(X_\text{ts},X_{tr})K(X_{tr},X_{tr})^{-1}f \]

y la covarianza predictiva es

\[ \operatorname{cov}(f_\text{ts}) = K(X_\text{ts},X_\text{ts}) - K(X_\text{ts},X_{tr})K(X_{tr},X_{tr})^{-1}K(X_{tr},X_\text{ts}). \]

Para el caso en el que las observaciones sí contienen ruido gaussiano \(\epsilon_x\sim\mathcal{N\left(0,\sigma_n^2\right)}\), las observaciones se expresan como

\[ y = f(x)+\varepsilon_x \]

y por lo tanto, la covarianza prior de las observaciones es

\[ \operatorname{cov}(y_{tr})=K(X_{tr},X_{tr})+\sigma_n^2 I \]

Se procede de manera similar al caso libre de ruido, i.e., usando la distribución conjunta de las observaciones de entrenamiento y los valores latentes de prueba:

\[ \begin{bmatrix}y_{tr}\\f_\text{ts}\end{bmatrix}\sim\mathcal{N}\left(0,\begin{bmatrix}K(X_{tr},X_{tr})+\sigma_n^2I & K(X_{tr},X_\text{ts})\\K(X_\text{ts},X_{tr}) & K(X_\text{ts},X_\text{ts})\end{bmatrix}\right) \]

\[\implies f_\text{ts} \mid X_{tr},y_{tr},X_\text{ts}\sim\mathcal{N}\left(\bar f_\text{ts},\operatorname{cov}(f_\text{ts})\right), \]

donde

\[ \bar {f_\text{ts}}=\mathbb{E}[f_\text{ts} \mid X_{tr},y_{tr},X_\text{ts}]=K(X_\text{ts},X_{tr})\left[K(X_{tr},X_{tr})+\sigma_n^2I\right]^{-1} \]

y

\[ \operatorname{cov}(f_\text{ts})=K(X_\text{ts},X_\text{ts})-K(X_\text{ts},X_{tr})\left[K(X_{tr},X_{tr})+\sigma_n^2I\right]^{-1}K(X_{tr},X_\text{ts}) \]

Estas expresiones constituyen la solución al problema de regresión bajo el modelo de procesos gaussianos (RMPG).

Revisitando al modelo lineal Bayesiano, notemos que este es un caso particular de los proceso Gaussianos, donde

\[ f(x)=\phi(x)^\top \mathbb{w} \]

y, dado que especificamos la priori $\mathbb{w}$ como \(\mathcal{N}(\mathbb{0},\Sigma_p)\), su matriz de covarianza es

\[ \mathbb{E}\big[(f(x)\cdot f(x')\big] = \phi(x)^\top \Sigma_p\ \phi(x') \]

Este es precisamente un kernel.

Esto dice que para cualquier conjunto de funciones base de algún espacio de características se puede construir una función de covarianza (kernel), trasladando así el modelo paramétrico al modelo no-paramétrico de un proceso gaussiano, donde en vez de trabajar con parámetros, y por lo tanto con una matriz de covarianza finita, se trabaja con una función de covarianza.

Aún más, más adelante veremos que el teorema de Mercer establece la inversa de la afirmación anterior: toda función de covarianza admite una expansión (posiblemente infinita) en términos de funciones base.

Funciones de Covarianza

Ya se estableció que, en regresión de procesos gaussianos, la definición del kernel es lo que arranca el proceso de aprendizaje, ya que es a través de la función de covarianza que se codifica la información inicial que se tiene sobre la función latente de la que provienen los datos que investigamos. Esta información especifica de alguna manera la noción de cercanía entre cada par de puntos de la función. Por ejemplo, una suposición básica podría ser que a entradas \(x\) cercanas, valores de $f(x)$, y por lo tanto de \(y\), cercanos también (o de alguna forma similares). En otras palabras, la función de covarianza define la suavidad de las funciones del proceso. Es así que el kernel expresa nuestras hipótesis a priori acerca de la estructura de la función que se desea modelar. A continuación formalizaremos esto.

Sea \(\mathcal{X}\) el espacio de entrada. Una función \(k:\mathcal{X}^2\to \mathbb{R}\) se denomina si satisface dos propiedades fundamentales:

  • Simetría: \(k(x,x')=k(x',x)\) para cualquier par de puntos \(x,x'\in \mathcal{X}\).

  • Dado un conjunto de $n$ puntos \(\{x_i|\ i=1,\dots,n\}\), la matriz de covarianzas Gram \(K_G\) correspondiente (\(K_G \ \ni \ K_{_Gij}=k(x_i,x_j)\)) es positiva semi-definida, es decir, para cualquier conjunto finito de puntos \({x_1,…,x_n}⊂X{x_1,…,x_n}⊂X\) y cualquier vector de coeficientes reales \(c_1,…,c_n\) se cumple que \(\sum_{i=1}^n\sum_{j=1}^nc_ic_jk(x_i,x_j)\ge0\).

La última propiedad es equivalente a decir que el kernel es positivo semi-definido:

\[ \int k(x,x')f(x)f(x')\,d\mu(x)\,d\mu(x')\geq 0. \]

En este contexto, cada entrada \(K_{ij}\) mide la similitud entre los puntos \(x_i\) y \(x_j\) de acuerdo con el kernel elegido.

Para que un kernel califique como función de covarianzas, este debe ser un kernel de Mercer.

Intuitivamente, los kernels de Mercer son una generalización de las matrices de covarianza en espacios de dimensión finita.

Espacios de Hilbert con kernel reproductor

Un (RKHS) \(\mathcal{H}_k\) es un espacio de Hilbert de funciones $f:$ asociado a un kernel \(k\): para todo \(x\in\mathcal{X}\) \[ k(x,\cdot)\in \mathcal{H}_k, \]

y satisface la : \[ f(x)=\langle f,k(x,\cdot)\rangle_{\mathcal{H}_k}. \]

Esta última propiedad implica que el valor de la función \(f\) en el punto \(x\) puede obtenerse mediante un producto interno entre \(f\) y el kernel centrado en \(x\). En este sentido, el kernel define una geometría sobre el espacio de funciones.

Las funciones de una RKHS pueden expresarse como combinaciones lineales de kernels: \[ f(\cdot)=\sum_{i=1}^{n}\alpha_i k(x_i,\cdot). \] Esto demuestra que las bases del espacio \(\mathcal{H}_k\) son kernels.

De lo anterior, si \[ f(\cdot)=\sum_{i=1}^{n}\alpha_i k(x_i,\cdot), \qquad g(\cdot)=\sum_{j=1}^{n}\beta_j k(y_j,\cdot), \] entonces el producto interno en la RKHS está dado por \[ \langle f,g\rangle_{\mathcal{H}_k} = \sum_{i=1}^{n}\sum_{j=1}^{n} \alpha_i\beta_j k(x_i,y_j). \] La norma inducida por este producto interno es \[ \|f\|_{\mathcal{H}_k} = \sqrt{\langle f,f\rangle_{\mathcal{H}_k}}. \]

Cabe remarcar que cada kernel genera un único RKHS, y viceversa.

El Teorema de Mercer

De aquí en adelante veremos cómo varias nociones y resultados fundamentales de matrices finitas se pueden extender a una noción de matrices infinitas y varios resultados análogos. Empezamos definiendo el análogo a los eigenvectores para espacios de funciones.

A una función \(\phi(\cdot)\) que satisface\[ \int k(\mathbf{x},\mathbf{x}') \psi(\mathbf{x})\,d\mu(\mathbf{x}) = \lambda \psi(\mathbf{x}') \] se le llama eigenfunción del kernel \(k\) con eigenvalor \(\lambda\) con respecto a la medida \(\mu\).

\begin{theorem}{El Teorema de Mercer} Sea \((\mathcal{X},\mu)\) un espacio de medida finito; sea \(k\in L_\infty(\mathcal{X}^2,\mu^2)\) un kernel que es una función de covarianzas. Entonces existe un conjunto ortonormal de eigenfunciones \(\{\psi_i\}_{i=1}^{\infty}\) del operador $T_k:L_2(\mathcal{X},\mu)\to L_2(\mathcal{X},\mu)$:

\[ (T_k f)(x):= \int_{\mathcal{X}} k(x,u)f(u)\,du, \]

con eigenvalores no negativos \(\{\lambda_i\}_{i=1}^{\infty}\) tales que \[ k(x,y) = \sum_{i=1}^{\infty} \lambda_i \psi_i(x)\psi_i(y), \] donde la convergencia es absoluta y uniforme casi seguramente en \(\mu^2\).

\end{theorem}

Esta descomposición es el análogo de dimensión infinita de la diagonalización de una matriz hermitiana.

En general, un kernel \(k\) tiene un número infinito de eigenfunciones (aunque también puede tener sólo un número finito de eigenfunciones diferentes a cero, en cuyo caso se dice que el kernel es degenerado y de rango finito).

El teorema de Mercer permite definir un mapa de características \[ \phi:\mathcal{X}\to\mathcal{H}, \] dado por \[ \phi(x) = \bigl[ \sqrt{\lambda_1}\psi_1(x), \sqrt{\lambda_2}\psi_2(x), \ldots \bigr]^\top. \] Con esta definición, el kernel puede escribirse como un producto interno: \[ k(x,y) = \langle \phi(x),\phi(y)\rangle_{\mathcal{H}} = \phi(x)^\top\phi(y). \] Lo de arriba ilustra cómo el kernel mide la similitud entre puntos después de ser transformados al espacio de características. Enfatizamos que dicho espacio puede ser de dimensión muy alta o incluso infinita, pero el permite calcular productos internos sin construir explícitamente \(\phi(x)\).

La expansión de Mercer suele ser infinita. Para hacerla operativa, se aproxima el kernel conservando solo los primeros \(m\) términos espectrales, usualmente aquellos asociados con los autovalores más grandes:

\[ k_m(x,y)= \sum_{i=1}^{m}\lambda_i\psi_i(x)\psi_i(y). \]

El término que se descarta es el residuo espectral:

\[ r_m(x,y)=k(x,y)-k_m(x,y) = \sum_{i=m+1}^{\infty}\lambda_i\psi_i(x)\psi_i(y). \]

Bajo las hipótesis del Teorema de Mercer, la expansión converge de forma absoluta y uniforme; por tanto, al aumentar \(m\), la aproximación \(k_m\) se acerca progresivamente al kernel original. Además, como los autovalores son no negativos, los primeros términos suelen concentrar la mayor parte de la variabilidad inducida por el kernel.

Este truncamiento equivale a reemplazar el mapa de características infinito

\[ \phi(x)= \left[\sqrt{\lambda_1}\psi_1(x), \sqrt{\lambda_2}\psi_2(x),\ldots\right]^\top \]

por un mapa finito de dimensión \(m\):

\[ \phi_m(x)= \left[\sqrt{\lambda_1}\psi_1(x), \ldots, \sqrt{\lambda_m}\psi_m(x)\right]^\top. \]

Así,

\[ k(x,y)\approx k_m(x,y)=\langle \phi_m(x),\phi_m(y)\rangle. \]

En una muestra finita, si la matriz kernel admite la descomposición espectral

\[ K=U\Lambda U^\top, \]

entonces una aproximación de rango \(m\) se obtiene como

\[ K_m=U_m\Lambda_m U_m^\top, \]

donde \(U_m\) contiene los \(m\) autovectores principales y \(\Lambda_m\) los \(m\) autovalores más grandes. Un criterio práctico para elegir \(m\) es conservar una proporción dada de la energía espectral:

\[ \rho_m= \frac{\sum_{i=1}^{m}\lambda_i}{\sum_{i=1}^{n}\lambda_i} \geq 1-\alpha, \]

En regresión con procesos gaussianos, esto permite aproximar la matriz de covarianza como

\[ K(X,X)\approx K_m(X,X), \]

lo cual reduce el costo computacional y conecta la teoría de Mercer con aproximaciones de bajo rango.

Ejemplos de Funciones de Covarianza

Existe una amplia variedad de kernels que cumplen con ser funciones de covarianza. Cada uno caracteriza de manera diferente el espacio de funciones compatibles con el proceso gaussiano, pero muchos de ellos comparten ciertas características que determinan propiedades importantes del espacio de funciones. Por ejemplo, cuando la función de covarianza depende únicamente de

  • la diferencia entre entradas, \(\tau=x-x'\) se dice que es estacionaria, pues es invariante respecto a traslaciones en el espacio de entradas.

  • de \(r=|x-x'|\), se dice que es isotrópica, pues es invariante respecto a cualquier movimiento rígido en el espacio de entradas. En este caso, se dice que el kernel \(k\) es una (RBF).

  • del producto punto \(x\cdot x'\), se denomina función de covarianza de producto punto. Estas funciones son invariantes ante rotaciones alrededor del origen, pero no ante traslaciones.

Suavidad de Procesos Estocásticos y la Función de Covarianza

Una medida de suavidad del proceso gaussiano en cuestión es considerar la continuidad y la diferenciabilidad en media cuadrática del proceso.

Definición: Sea \(x_1, x_2, \ldots\) una sucesión de puntos en el espacio $^D$, y sea \(x_\text{ts}\) un punto fijo en $^D$ tal que \(|x_k - x_\text{ts}| \to 0\) cuando \(k \to \infty\). Entonces un proceso \(f(x)\) es en \(x_\text{ts}\) si \(\mathbb{E}[|f(x_k) - f(x_\text{ts})|^2] \to 0\) cuando \(k \to \infty\).

Un proceso estocástico es continuo en media cuadrática en \(x_\text{ts}\) si y solo si su función de covarianza \(k(x, x')\) es continua en el punto \(x = x' = x_\text{ts}\). Es fácil verificar esta afirmación mediante la siguiente identidad

\[ k(x,x)−k(x_∗,x_∗)=\mathbb{E}\big[|f(x)−f(x_∗)|^2\big]−2\big(k(x_∗,x_∗)−k(x,x_∗)\big), \]

la cual se obtiene de \[ \mathbb{E}\big[|f(x)−f(x_∗)|^2\big]=k(x,x)+k(x_\text{ts},x_\text{ts})−2k(x,x_\text{ts}). \]

Para funciones de covarianza estacionarias esto se reduce a verificar la continuidad en \(k(0)\).

Definición: La derivada en media cuadrática de \(f(x)\) en la dirección \(i\)-ésima se define como \[ \frac{\partial f(x)}{\partial x_i} = \operatorname{l.i.m.}_{h \to 0} \frac{f(x + h e_i) - f(x)}{h}, \] donde \(\operatorname{l.i.m.}\) denota el límite en media cuadrática y \(e_i\) es el vector unitario en la dirección \(i\)-ésima.

Se puede demostrar que la función de covarianza de \(\partial f(x)/\partial x_i\) está dada por \(\partial^2 k(x, x')/\partial x_i \partial x_j'\). Estas definiciones pueden extenderse a derivadas de orden superior.

Para procesos estacionarios, si la derivada parcial de orden $2k$ \(\partial^{2k} f(x)/\partial^2 x_i \ldots \partial^2 x_k\) existe y es finita en \(x = 0\), entonces la derivada parcial de orden \(k\) \(\partial^k f(x)/\partial x_i \ldots x_k\) existe para todo \(x \in \mathbb{R}^D\) como un límite en media cuadrática.

Ejemplos de Funciones de Covarianza

A continuación mencionamos algunas de las funciones de covarianza más comunes.

  • La función de covarianza exponencial cuadrática (SE) tiene la forma \[ k_{\mathrm{SE}}(r)=\exp\!\left(-\frac{r^2}{2\ell^2}\right), \] donde \(\ell\) define la característica de escala longitud.

    En una dimensión, el número medio de cruces por cero para un proceso SE es proporcional a \((2\pi \ell)^{-1}\). Esta función de covarianza es infinitamente diferenciable, por lo que el proceso gaussiano asociado tiene derivadas de media cuadrática de todos los órdenes y es por lo tanto muy suave.

    \[ \Downarrow \]

  • \(k(r)^t\) es un kernel válido para todo \(t>0\), cuando $k(r)$ es una SE.

  • La clase de covarianzas de Matérn. Un kernel de esta familia es de la forma \[ k_{\mathrm{Matern}}(r)= \frac{2^{1-\nu}}{\Gamma(\nu)} \left(\frac{\sqrt{2\nu}r}{\ell}\right)^\nu K_\nu\!\left(\frac{\sqrt{2\nu}r}{\ell}\right), \] con parámetros positivos \(\nu\) y \(\ell\), y \(K_\nu\) es una función de Bessel modificada.

    \[ \Downarrow \]

  • El caso especial \(\nu=1/2\) dentro de la clase de Matérn produce la función de covarianza exponencial: \[ k(r)=\exp(-r/\ell). \]

    \[ \Downarrow \]

  • En dimensión \(D=1\), la covarianza exponencial corresponde al proceso de Ornstein–Uhlenbeck (OU).

Ejemplos

La “receta”

Finalmente, al momento de querer realizar una regresión con procesos gausianos, se siguen de manera general los siguientes pasos:

  1. Definir los conjuntos de datos
  • Datos de entrenamiento: Un conjunto \(\mathcal{D} = (X_{tr}, y) = \{(x_i, y_i) \mid i = 1, \dots, n\}\).
  • Datos de prueba: Los puntos de entrada \(X_{ts}\) donde deseas estimar la salida \(f_{ts}\).
  1. Elegir una función kernel

3.Construir la matriz de kernel \(K\in\mathbb{R}^{n\times n}\), el vector de kernel \(\textbf{k}_{*}\) y calcular \(k_{* *}\).
Utilizando la función de kernel elegida, se construyen:

  • \(K \in \mathbb{R}^{n \times n}\): La matriz de covarianza entre todos los puntos de entrenamiento, donde \(K_{ij} = k(x_i, x_j)\).

  • \(\textbf{k}_* \in \mathbb{R}^n\): El vector de covarianza entre los puntos de entrenamiento y el punto de prueba, donde \((k_*)_i = k(x_i, x_{ts})\).

  • \(k_{* *} = k(x_{ts}, x_{ts})\).

  1. Calcular la media y varianza predictivas:
  • Media predictiva:

\[\mathbb{E}[f_{ts}] = k_*^T (K + \sigma_n^2 I)^{-1} y\]

  • Varianza predictiva

\[\text{Var}[f_{ts}] = k(x_{ts}, x_{ts}) - k_*^T (K + \sigma_n^2 I)^{-1} k_*\]

Ejemplo 1: millas por galón vs caballos de fuerza

Auto MPG es un conjunto de datos que contiene el consumo de combustible, el rendimiento y las características técnicas de varios modelos de automóviles de los años 70 y principios de los 80.

Consideremos:

  • Variable Predictora (\(x\)): hp (Horsepower / Caballos de fuerza). Representa la potencia del motor.
  • Variable Respuesta (\(y\)): mpg (Miles per Gallon / Millas por galón). Mide la eficiencia de combustible del vehículo.

Dado que suponemos que el comportamiento de los autos es continuo, de tal forma que si un auto con \(100\text{ hp}\) rinde \(25\text{ mpg}\), uno con \(101\text{ hp}\) tendrá un rendimiento muy similar, para este ejemplo consideraremos el Kernel RBF / Exponencial Cuadrado / SE (Squared Exponencial), el cual es el kernel por excelencia para funciones suaves: \[K(x, x') = \exp(-\gamma ||x - x'||^2),\]

Parte importante a destacar de esta elección es que captura la relación de que los puntos cercanos en el espacio de entrada (\(x\)) deben tener salidas muy parecidas en el espacio de destino (\(y\)).

# 1. Cargar la librería necesaria para graficar
if(!require(ggplot2)) install.packages("ggplot2")
library(ggplot2)

# 2. Descargar y preparar los datos desde el repositorio de la UCI
url <- "https://archive.ics.uci.edu/ml/machine-learning-databases/auto-mpg/auto-mpg.data"
dataFrame <- read.table(url, header = FALSE, na.strings = "?")
names(dataFrame) <- c("mpg", "cylinders", "displacement", "hp", "weight", 
                      "acceleration", "model_year", "origin", "car_name")

# Limpieza básica: Eliminar filas con valores faltantes (NA) en 'hp'
dataFrame <- na.omit(dataFrame)

# Para efectos prácticos del cálculo de matrices, ordenamos por 'hp'
dataFrame <- dataFrame[order(dataFrame$hp), ]

# 3. Definimos las variables de entrenamiento
X_tr <- dataFrame$hp
y    <- dataFrame$mpg

# Hiperparámetros del modelo (ajustados para la escala de 'hp')
sigma_n2 <- 3.5      # Varianza del ruido de las observaciones (ancho mínimo en los puntos)
l        <- 15       # Longitud de escala del kernel (qué tan suave es la curva)
sigma_f2 <- 150      # Varianza de la función (amplitud vertical de la incertidumbre)

# Función para calcular el Kernel Gaussiano (Squared Exponential)
kernel_gaussiano <- function(x1, x2, l, sigma_f2) {
  dist_cuadrado <- outer(x1, x2, "-")^2
  return(sigma_f2 * exp(-dist_cuadrado / (2 * l^2)))
}

# Construcción de las matrices esenciales
K <- kernel_gaussiano(X_tr, X_tr, l, sigma_f2)
I <- diag(length(X_tr))
K_inv <- solve(K + sigma_n2 * I) # Inversión de la matriz de entrenamiento

# 4. Puntos de prueba para generar las curvas e incertidumbres
X_ts <- seq(min(X_tr), max(X_tr), length.out = 300)

# Inicializar vectores para almacenar la media y la varianza predictiva
pred_media <- numeric(length(X_ts))
pred_var   <- numeric(length(X_ts))

# Paso a paso para cada punto de prueba
for (i in 1:length(X_ts)) {
  k_star <- kernel_gaussiano(X_tr, X_ts[i], l, sigma_f2)
  k_starstar <- kernel_gaussiano(X_ts[i], X_ts[i], l, sigma_f2)
  
  # Media predictiva (Ecuación de la receta)
  pred_media[i] <- t(k_star) %*% K_inv %*% y
  
  # Varianza predictiva (Incertidumbre latente)
  pred_var[i] <- k_starstar - t(k_star) %*% K_inv %*% k_star
}

# Crear el data frame con los resultados predictivos
puntos_prueba <- data.frame(
  hp = X_ts,
  Pred_GPR = pred_media,
  se = sqrt(pmax(0, pred_var)) # Desviación estándar predictiva
)

# Calcular los límites reales de la banda (l_n y u_n)
puntos_prueba$l_n <- puntos_prueba$Pred_GPR - (1.96 * puntos_prueba$se)
puntos_prueba$u_n <- puntos_prueba$Pred_GPR + (1.96 * puntos_prueba$se)

# 5. Ajustar también la regresión lineal para mantener la comparación
modelo_lineal <- lm(mpg ~ hp, data = dataFrame)
puntos_prueba$Pred_Lineal <- predict(modelo_lineal, newdata = puntos_prueba)

# 6. Graficar el resultado con ggplot2
ggplot(data = dataFrame, aes(x = hp, y = mpg)) +
  # Banda de incertidumbre
  geom_ribbon(data = puntos_prueba, aes(y = Pred_GPR, ymin = l_n, ymax = u_n), 
              fill = "royalblue", alpha = 0.40) +
  
  # Líneas discontinuas para marcar los límites l_n y u_n de la banda
  geom_line(data = puntos_prueba, aes(y = l_n), color = "royalblue", linetype = "dashed", alpha = 0.6) +
  geom_line(data = puntos_prueba, aes(y = u_n), color = "royalblue", linetype = "dashed", alpha = 0.6) +
  
  # Puntos observados de entrenamiento
  geom_point(alpha = 0.7, color = "black", size = 2) +
  
  # Recta del modelo lineal clásico
  geom_line(data = puntos_prueba, aes(y = Pred_Lineal, color = "Regresión Lineal"), 
            linewidth = 1.2) +
  
  # Curva del Proceso Gaussiano (Media predictiva)
  geom_line(data = puntos_prueba, aes(y = Pred_GPR, color = "Proceso Gaussiano (GPR)"), 
            linewidth = 1.2) +
  
  # Colores y etiquetas
  scale_color_manual(name = "Modelos", 
                     values = c("Regresión Lineal" = "firebrick", 
                                "Proceso Gaussiano (GPR)" = "royalblue")) +
  labs(title = "Proceso Gaussiano con Bandas de Incertidumbre Reales",
       subtitle = "La varianza predictiva se reduce drásticamente cerca de los datos observados",
       x = "Caballos de Fuerza (hp)",
       y = "Millas por Galón (mpg)") +
  theme_minimal() +
  theme(legend.position = "bottom",
        plot.title = element_text(face = "bold", size = 14),
        plot.subtitle = element_text(size = 11))

Ejemplo 2: extracción de cafeína

En Effect of time and water temperature on caffeine extraction from coffee (Nhan, P. P. & Phu, N. T., 2012) se realizaron diferentes experimentos para determinar una manera “óptima” de extracción de cafeína.

En particular consideramos experimento con los siguientes registros:

Time (minutes) C. arabica (Caffeine %) C. robusta (Caffeine %)
5 1.16 1.94
10 1.17 2.00
15 1.20 2.07
30 1.21 2.07 
# Instalar los paquetes necesarios
# install.packages("DiceKriging")
# install.packages("ggplot2")

library(DiceKriging)
library(ggplot2)

# ==============================================================================
# DATASET: Tabla 2 - Efecto del Tiempo en la Extracción de Cafeína a 100°C
# ==============================================================================

# Crear el dataframe basado en los datos reales de la Tabla 2 del paper
datos_tabla2 <- data.frame(
  Tiempo = c(5, 10, 15, 30),
  Arabica = c(1.16, 1.17, 1.20, 1.21),
  Robusta = c(1.94, 2.00, 2.07, 2.07)
)

print("--- MODELOS PARA LA TABLA 2 (Tiempo vs Cafeína) ---")
## [1] "--- MODELOS PARA LA TABLA 2 (Tiempo vs Cafeína) ---"
# ==============================================================================
# 1. MODELADO CON PROCESOS GAUSSIANOS (Kriging)
# ==============================================================================

# Ajustar GP para C. arabica con kernel Gaussiano
modelo_gp_arabica <- km(formula = ~1, 
                        design = data.frame(x = datos_tabla2$Tiempo), 
                        response = data.frame(y = datos_tabla2$Arabica),
                        covtype = "gauss", nugget = 1e-6)
## 
## optimisation start
## ------------------
## * estimation method   : MLE 
## * optimisation method : BFGS 
## * analytical gradient : used
## * trend model : ~1
## * covariance model : 
##   - type :  gauss 
##   - nugget : 1e-06 
##   - parameters lower bounds :  1e-10 
##   - parameters upper bounds :  50 
##   - variance bounds :  5.656667e-05 0.00699 
##   - best initial criterion value(s) :  10.37254 
## 
## N = 2, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate     0  f=      -10.373  |proj g|=      0.18755
## At iterate     1  f =      -10.405  |proj g|=    0.00069372
## At iterate     2  f =      -10.437  |proj g|=      0.038277
## At iterate     3  f =      -10.442  |proj g|=      0.096882
## At iterate     4  f =      -10.446  |proj g|=      0.069568
## At iterate     5  f =      -10.446  |proj g|=      0.073467
## At iterate     6  f =      -10.446  |proj g|=      0.073974
## At iterate     7  f =      -10.446  |proj g|=       0.07409
## At iterate     8  f =      -10.446  |proj g|=      0.074753
## At iterate     9  f =      -10.446  |proj g|=      0.075503
## At iterate    10  f =      -10.446  |proj g|=      0.076824
## At iterate    11  f =      -10.446  |proj g|=      0.078656
## At iterate    12  f =      -10.446  |proj g|=      0.081121
## At iterate    13  f =      -10.447  |proj g|=      0.083539
## At iterate    14  f =       -10.45  |proj g|=      0.083347
## At iterate    15  f =      -10.455  |proj g|=      0.073426
## At iterate    16  f =      -10.465  |proj g|=       0.04333
## At iterate    17  f =      -10.465  |proj g|=      0.016732
## At iterate    18  f =      -10.471  |proj g|=    0.00056798
## At iterate    19  f =      -10.471  |proj g|=      0.006455
## At iterate    20  f =      -10.471  |proj g|=     0.0064538
## At iterate    21  f =      -10.471  |proj g|=    0.00047955
## At iterate    22  f =      -10.471  |proj g|=    0.00047954
## 
## iterations 22
## function evaluations 25
## segments explored during Cauchy searches 23
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.000479537
## final function value -10.4712
## 
## F = -10.4712
## final  value -10.471224 
## converged
# Ajustar GP para C. robusta con kernel Gaussiano
modelo_gp_robusta <- km(formula = ~1, 
                        design = data.frame(x = datos_tabla2$Tiempo), 
                        response = data.frame(y = datos_tabla2$Robusta),
                        covtype = "gauss", nugget = 1e-6)
## 
## optimisation start
## ------------------
## * estimation method   : MLE 
## * optimisation method : BFGS 
## * analytical gradient : used
## * trend model : ~1
## * covariance model : 
##   - type :  gauss 
##   - nugget : 1e-06 
##   - parameters lower bounds :  1e-10 
##   - parameters upper bounds :  50 
##   - variance bounds :  0.0003632333 0.03932333 
##   - best initial criterion value(s) :  6.676638 
## 
## N = 2, M = 5 machine precision = 2.22045e-16
## At X0, 0 variables are exactly at the bounds
## At iterate     0  f=      -6.6766  |proj g|=      0.27009
## At iterate     1  f =      -6.7976  |proj g|=      0.078844
## At iterate     2  f =      -6.7993  |proj g|=      0.089109
## At iterate     3  f =      -6.7999  |proj g|=       0.10048
## At iterate     4  f =      -6.7999  |proj g|=       0.09953
## At iterate     5  f =         -6.8  |proj g|=      0.097486
## At iterate     6  f =      -6.8002  |proj g|=      0.094147
## At iterate     7  f =      -6.8008  |proj g|=      0.088455
## At iterate     8  f =      -6.8022  |proj g|=      0.078911
## At iterate     9  f =      -6.8056  |proj g|=      0.063049
## At iterate    10  f =      -6.8136  |proj g|=      0.037943
## At iterate    11  f =       -6.831  |proj g|=     0.0072606
## At iterate    12  f =      -6.8616  |proj g|=      0.014224
## At iterate    13  f =      -6.8634  |proj g|=      0.033678
## At iterate    14  f =      -6.8644  |proj g|=     0.0055571
## At iterate    15  f =      -6.8644  |proj g|=     0.0055364
## At iterate    16  f =      -6.8644  |proj g|=      0.012546
## At iterate    17  f =      -6.8644  |proj g|=     0.0021285
## 
## iterations 17
## function evaluations 20
## segments explored during Cauchy searches 18
## BFGS updates skipped 0
## active bounds at final generalized Cauchy point 0
## norm of the final projected gradient 0.0021285
## final function value -6.86444
## 
## F = -6.86444
## final  value -6.864441 
## converged
# ==============================================================================
# 2. PREDICCIÓN Y REGIONES DE INCERTIDUMBRE
# ==============================================================================

# Grid de tiempo para evaluar las curvas (de 2 a 33 minutos)
pred_grid_t2 <- data.frame(x = seq(2, 33, length.out = 300))

# Predicciones de Media y Desviación Estándar para ambas especies
preds_arabica <- predict(modelo_gp_arabica, newdata = pred_grid_t2, type = "UK")
preds_robusta <- predict(modelo_gp_robusta, newdata = pred_grid_t2, type = "UK")

# Estructurar los resultados en el dataframe de predicciones para C. arabica
pred_grid_t2$media_arabica    <- preds_arabica$mean
pred_grid_t2$inferior_arabica <- preds_arabica$mean - 2 * preds_arabica$sd
pred_grid_t2$superior_arabica <- preds_arabica$mean + 2 * preds_arabica$sd

# Estructurar los resultados en el dataframe de predicciones para C. robusta
pred_grid_t2$media_robusta    <- preds_robusta$mean
pred_grid_t2$inferior_robusta <- preds_robusta$mean - 2 * preds_robusta$sd
pred_grid_t2$superior_robusta <- preds_robusta$mean + 2 * preds_robusta$sd

# ==============================================================================
# 3. REGRESIÓN LINEAL TRADICIONAL PARA COMPARACIÓN
# ==============================================================================
lm_arabica <- lm(Arabica ~ Tiempo, data = datos_tabla2)
lm_robusta <- lm(Robusta ~ Tiempo, data = datos_tabla2)

pred_grid_t2$lm_arabica <- predict(lm_arabica, newdata = data.frame(Tiempo = pred_grid_t2$x))
pred_grid_t2$lm_robusta <- predict(lm_robusta, newdata = data.frame(Tiempo = pred_grid_t2$x))

# ==============================================================================
# 4. GRÁFICO COMPARATIVO CON ggplot2
# ==============================================================================

ggplot() +
  # --- BANDAS DE INCERTIDUMBRE (INTERVALOS DE CONFIANZA DEL 95%) ---
  # C. robusta (Franja Azul Transparente)
  geom_ribbon(data = pred_grid_t2, aes(x = x, ymin = inferior_robusta, ymax = superior_robusta), 
              fill = "royalblue", alpha = 0.15) +
  # C. arabica (Franja Roja Transparente)
  geom_ribbon(data = pred_grid_t2, aes(x = x, ymin = inferior_arabica, ymax = superior_arabica), 
              fill = "red2", alpha = 0.15) +
  
  # --- MEDIAS PREDICTIVAS DE PROCESOS GAUSSIANOS ---
  geom_line(data = pred_grid_t2, aes(x = x, y = media_robusta, color = "Media GP (C. robusta)"), size = 1.2) +
  geom_line(data = pred_grid_t2, aes(x = x, y = media_arabica, color = "Media GP (C. arabica)"), size = 1.2) +
  
  # --- LÍNEAS DE REGRESIÓN LINEAL TRADICIONAL ---
  geom_line(data = pred_grid_t2, aes(x = x, y = lm_robusta, color = "Regresión Lineal (C. robusta)"), size = 0.8, linetype = "dashed") +
  geom_line(data = pred_grid_t2, aes(x = x, y = lm_arabica, color = "Regresión Lineal (C. arabica)"), size = 0.8, linetype = "dashed") +
  
  # --- DATOS REALES DE ENTRADA ---
  geom_point(data = datos_tabla2, aes(x = Tiempo, y = Robusta, color = "Observaciones (C. robusta)"), size = 3.5, stroke = 1.2) +
  geom_point(data = datos_tabla2, aes(x = Tiempo, y = Arabica, color = "Observaciones (C. arabica)"), size = 3.5, stroke = 1.2) +
  
  labs(title = "Análisis de la Tabla 2: Tiempo vs Extracción de Cafeína a 100°C",
       subtitle = "La varianza predictiva disminuye a cero en los puntos de entrenamiento (5, 10, 15 y 30 min)",
       x = "Tiempo de Extracción (minutos)",
       y = "Contenido de Cafeína (%)",
       color = "Modelos y Datos por Especie") +
  
  scale_color_manual(values = c(
    "Media GP (C. robusta)"         = "blue",
    "Regresión Lineal (C. robusta)" = "cyan3",
    "Observaciones (C. robusta)"    = "darkblue",
    "Media GP (C. arabica)"         = "red",
    "Regresión Lineal (C. arabica)" = "orange",
    "Observaciones (C. arabica)"    = "darkred"
  )) +
  
  theme_minimal() +
  theme(legend.position = "right")

Referencias