markjrieke / riekelib

A collection of functions I use regulary
https://markjrieke.github.io/riekelib/
Other
2 stars 0 forks source link

fix condition Gaussian process #25

Open markjrieke opened 4 weeks ago

markjrieke commented 4 weeks ago

IIRC, I ran into something the other day where the function throws an error when you pass a single value to x_new

Idk --- try to break it again then fix it

markjrieke commented 1 week ago
# this fails for non-conformable arguments
condition_gaussian_process(
  x = 0,
  y = 0,
  x_new = 1:100/100,
  amplitude = 0.2,
  length_scale = 0.4
)

# this works (basically rewrite the internals)
n <- 1000

x <- c(0, 0.5, 0.7) * 0.002
y <- c(0, 0.3, 0.1) * 0.002
x_new <- 1:100/100
amplitude <- 0.2 * 0.002
length_scale <- 0.4
delta <- 1e-9

x1 <- x_new
x2 <- x
x_combined <- c(x1, x2)
S <- cov_exp_quad(x_combined, amplitude, length_scale, delta)
S11 <- matrix(S[1:length(x1), 1:length(x1)], nrow = length(x1), ncol = length(x1))
S22 <- matrix(S[(length(x1) + 1):nrow(S), (length(x1) + 1):ncol(S)], nrow = length(x), ncol = length(x))
S21 <- matrix(S[(length(x1) + 1):nrow(S), 1:length(x1)], nrow = length(x), ncol = length(x1))
S12 <- t(S21)
mu_bar <- S12 %*% MASS::ginv(S22) %*% y
Sigma_bar <- S11 - S12 %*% MASS::ginv(S22) %*% S21
y_new <- matrix(nrow = length(x_new), ncol = n)
for (i in 1:n) {
  y_new[, i] <- MASS::mvrnorm(1, mu_bar, Sigma_bar)
}