---
title: "Gaussian-Process Regression from First Principles"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Gaussian-Process Regression from First Principles}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 7,
  fig.height = 4
)
library(gaussianprocesses)
```

# Model

Let

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

where \(m(x)\) is a mean function and \(k(x,x')\) is a positive
semidefinite covariance kernel.

For training inputs

\[
X = (x_1,\ldots,x_n)^\top
\]

and observations

\[
y_i = f(x_i) + \varepsilon_i,
\qquad
\varepsilon_i \sim \mathcal N(0,\sigma_i^2),
\]

define

\[
K_{ff} = K(X,X),
\qquad
\Sigma_\varepsilon =
\operatorname{diag}(\sigma_1^2,\ldots,\sigma_n^2).
\]

Then

\[
y \sim \mathcal N\left(
m(X),
K_{ff}+\Sigma_\varepsilon
\right).
\]

The package represents the same objects explicitly:

| Mathematics | Package object |
|---|---|
| \(m\) | `zero_mean()`, `constant_mean()` |
| \(k\) | kernel specification such as `rbf_kernel()` |
| \(K(X,X)\) | `evaluate_kernel(kernel, X)` |
| \(\Sigma_\varepsilon\) | `noise_variance` in `fit_gp()` |
| fitted conditional model | `fit_gp()` |

A scalar `noise_variance` gives the usual homoscedastic model. A vector gives
one known observation variance per training point.

# Conditioning

At prediction inputs \(X_*\), define

\[
K_{f*} = K(X,X_*),
\qquad
K_{**} = K(X_*,X_*).
\]

Write

\[
C = K_{ff}+\Sigma_\varepsilon.
\]

The posterior latent mean is

\[
\mu_* =
m(X_*) +
K_{*f}C^{-1}\left(y-m(X)\right),
\]

and the latent covariance is

\[
\Sigma_* =
K_{**} -
K_{*f}C^{-1}K_{f*}.
\]

The implementation never forms \(C^{-1}\). It computes a Cholesky
factorization

\[
C = R^\top R
\]

and solves triangular systems.

The stored vector

\[
\alpha =
C^{-1}\left(y-m(X)\right)
\]

corresponds directly to `model$alpha`.

# Minimal reproducible fit

```{r}
x <- c(-1.5, -0.5, 0.5, 1.5)
y <- c(0.9, -0.1, 0.2, 1.0)

kernel <- rbf_kernel(
  variance = 1.2,
  length_scale = 0.7
)

model <- fit_gp(
  x,
  y,
  kernel = kernel,
  noise_variance = 0.05
)

model
```

The fitted object stores the Cholesky factor \(R\) and the vector
\(\alpha\). The covariance matrices are not stored, because the kernel
recomputes them exactly, and \(R^\top R\) reproduces \(C\) plus any
numerical jitter (here none):

```{r}
K <- evaluate_kernel(model$kernel, model$x)
C <- K + diag(model$training_noise_variance)
C

model$cholesky
all.equal(crossprod(model$cholesky), C + diag(model$jitter, nrow(C)))

model$alpha
```

`?gaussianprocesses_model` describes every field of a fitted model, how
models are versioned, and how to save them.

# Posterior prediction

```{r}
x_new <- seq(-2, 2, length.out = 80)

prediction <- predict_gp(
  model,
  x_new,
  interval_level = 0.95
)

head(
  data.frame(
    x = x_new,
    mean = prediction$mean,
    latent_sd = prediction$latent_sd,
    observation_sd = prediction$observation_sd
  )
)
```

The package reports two uncertainty layers.

The latent variance is

\[
\operatorname{Var}\left[f(x_*)\mid y\right].
\]

The observation variance is

\[
\operatorname{Var}\left[y_*\mid y\right]
=
\operatorname{Var}\left[f(x_*)\mid y\right]
+
\sigma_*^2.
\]

They are deliberately not conflated.

```{r, fig.alt = "Posterior mean of the fitted Gaussian process for x from -2 to 2, with a hatched 95% prediction band and the four training observations as points."}
plot(
  x_new,
  prediction$mean,
  type = "l",
  xlab = "x",
  ylab = "response",
  ylim = range(prediction$prediction_interval)
)
polygon(
  c(x_new, rev(x_new)),
  c(
    prediction$prediction_interval[, "lower"],
    rev(prediction$prediction_interval[, "upper"])
  ),
  border = NA,
  density = 15
)
lines(x_new, prediction$mean)
points(x, y)
```

# Function draws

A mean and an interval summarize each input separately. Draws show whole
functions, including how values at nearby inputs move together.

A draw from \(\mathcal N(\mu, \Sigma)\) is

\[
f = \mu + R^\top z,
\qquad
\Sigma = R^\top R,
\qquad
z \sim \mathcal N(0, I),
\]

where \(R\) is the Cholesky factor of the covariance. No inverse is needed.
For the prior, \(\mu = m(X_*)\) and \(\Sigma = K_{**}\). For the latent
posterior, \(\mu\) and \(\Sigma\) are the posterior mean and covariance from
`predict_gp()`. Closely spaced inputs make \(\Sigma\) numerically singular,
so the factorization uses the package's jitter policy and reports the jitter.

```{r, fig.alt = "Two panels. Left: five smooth prior function draws spread around zero across x from -2 to 2. Right: five posterior function draws that pass close to the four training observations and fan out between and beyond them.", fig.width = 7, fig.height = 3.5}
prior_draws <- sample_gp_prior(x_new, kernel, n_draws = 5, seed = 1)
posterior_draws <- sample_gp_posterior(model, x_new, n_draws = 5, seed = 1)

old_par <- par(mfrow = c(1, 2), mar = c(4, 4, 2, 1))
matplot(
  x_new, prior_draws$latent,
  type = "l", lty = 1, col = "grey40",
  xlab = "x", ylab = "f(x)", main = "Prior draws"
)
matplot(
  x_new, posterior_draws$latent,
  type = "l", lty = 1, col = "grey40",
  xlab = "x", ylab = "f(x)", main = "Posterior draws"
)
points(x, y, pch = 19)
par(old_par)
```

The posterior draws agree near the observations and spread apart away from
them, which is what the pointwise intervals summarize.

The 80 closely spaced inputs make both covariance matrices numerically
singular, so both factorizations needed a small jitter:

```{r}
c(prior = prior_draws$jitter, posterior = posterior_draws$jitter)
```

Latent draws are values of the function \(f\). A future observation also
includes noise, \(y_* = f_* + \varepsilon_*\). With `type = "observation"`,
each noisy draw is generated from a latent draw plus independent noise, and
both are returned:

```{r, fig.alt = "One latent posterior function draw as a line, with the corresponding noisy observation draws as small points scattered around it, and the training observations as larger points."}
noisy <- sample_gp_posterior(
  model,
  x_new,
  n_draws = 1,
  type = "observation",
  seed = 2
)

plot(
  x_new, noisy$observation[, 1],
  pch = 20, cex = 0.6, col = "grey50",
  xlab = "x", ylab = "response"
)
lines(x_new, noisy$latent[, 1])
points(x, y, pch = 19)
```

Over many draws, the latent draws reproduce the posterior covariance and the
observation draws reproduce that covariance plus the noise variance:

```{r}
many <- sample_gp_posterior(
  model,
  c(-1, 0.5),
  n_draws = 5000,
  type = "observation",
  seed = 3
)

round(
  rbind(
    latent_sample = apply(many$latent, 1, var),
    latent_theory = diag(many$latent_covariance),
    observation_sample = apply(many$observation, 1, var),
    observation_theory = diag(many$latent_covariance) +
      many$observation_noise_variance
  ),
  4
)
```

# Known heteroscedastic observation variance

Exact conditioning also supports a known diagonal observation covariance:

```{r}
known_noise <- c(0.02, 0.04, 0.08, 0.16)

heterogeneous_fit <- fit_gp(
  x,
  y,
  kernel = kernel,
  noise_variance = known_noise
)

heterogeneous_fit$training_noise_variance
```

For unknown input-dependent noise, see the package's heteroscedastic GP
approximation. The exact model above is still the underlying conditioning
mechanism once a diagonal noise estimate has been supplied.