# Implementing Gaussian Process in Python and R

library(reticulate)
reticulate::use_condaenv("dl")

## Sampling from the Posterior

Here we will implement “Prediction using Noisy Observations” because the Noise-free version can be understood as a special case of the noisy one with $$\sigma_n = 0$$.

### Implementations in Python

def sample_posterior(
x, y, x_star, mean_func, cov_func, cov_args = {},
sigma_n = 0.1, random_seed = -1, n_samples = 5
):
k_xx = cov_func(get_r(x, x), **cov_args)
k_xxs = cov_func(get_r(x, x_star), **cov_args)
k_xsx = cov_func(get_r(x_star, x), **cov_args)
k_xsxs = cov_func(get_r(x_star, x_star), **cov_args)

I = np.identity(k_xx.shape[1])
k_xx_noise = np.linalg.inv(k_xx + sigma ** 2 * I)
kxsx_kxxNoise = np.matmul(k_xsx, k_xx_noise)
# Eq.2.23, 24
fsb = np.matmul(kxsx_kxxNoise, y_train_N)
cov_fs = k_xsxs - np.matmul(kxsx_kxxNoise, k_xxs)
random_seed = int(random_seed)
if random_seed < 0:
prng = np.random
else:
prng = np.random.RandomState(random_seed)
out = prng.multivariate_normal(fsb, cov_fs, n_samples)
return out, fsb, cov_fs

### Implementations in R

sample_posterior <- function(
x, y, x_star, mean_func, cov_func, ...,
sigma_n = 0.1, random_seed = -1, n_samples = 5
) {
k_xx = cov_func(get_r(x, x), ...)
k_xxs = cov_func(get_r(x, x_star), ...)
k_xsx = cov_func(get_r(x_star, x), ...)
k_xsxs = cov_func(get_r(x_star, x_star), ...)

I = diag(1, dim(k_xx)[1])
k_xx_noise = solve(k_xx + sigma_n ^ 2 * I)
kxsx_kxxNoise = k_xsx %*% k_xx_noise
# Eq.2.23, 24
fsb = kxsx_kxxNoise %*% y
cov_fs = k_xsxs - kxsx_kxxNoise %*% k_xxs
random_seed <- as.integer(random_seed)
if (random_seed > 0) set.seed(random_seed)
return(MASS::mvrnorm(n_samples, fsb, cov_fs))
}

This time let’s try to fit some points in R.

x_train <- c(-10, -8, -5, 3, 7, 15)
y_train <- c(-5, -2, 3, 3, 2, 5)
x_star <- seq(-20, 20, length = 200)

par(mfrow = c(3, 3))
l <- c(0.25, 1, 4)
v <- c(0.5, 2, 8)
for (i in seq(length(l))) {
for (j in seq(length(v))) {
dt <- sample_posterior(
x_train, y_train, x_star, mean_func, matern_kernel,
l = l[i], v = v[j], n_samples = 10, random_seed = 100
)
matplot(x_star, t(dt), "l", col = rgb(0.3, 0.3, 0.3, 0.3), lty = 1,
main = paste("l=", l[i], ",v=", v[j]), ylim = c(-8, 8),
xlab = "x", ylab = "y")
points(x = x_train, y = y_train, pch = 16, col = "red")
}
}

Note, that when $$\ell$$ is small, it is easier for the predicted posterior to return to normal (prior), which is the mean function, 0 (see the points around x = 0). As $$\ell$$ increases, it becomes more and more likely the predicted $$y_{x=0}$$ to stay at the “local” value, which is provided by the nearest neighbor in y.