Generate functional data from Gaussian Process in R

Question

Model:

X(t) = 4*t + e(t); 

t € [0; 1]

e(t) is a Gaussian process with zero mean and covariance function f(s, t) = exp( -|t - s| )

The final result over 100 runs (=100 gray lines) with 50 sampled points each should be like the gray area in the picture.

The green line is what I get from the code below.

library(MASS)

kernel_1 <- function(x, y){
    exp(- abs(x - y))
}
cov_matrix <- function(x, kernel_fn, ...) {
    outer(x, x, function(a, b) kernel_fn(a, b, ...))
}
draw_samples <- function(x, N=1, kernel_fn, ...) {
    set.seed(100)
    Y <- matrix(NA, nrow = length(x), ncol = N)
    for (n in 1:N) {
        K <- cov_matrix(x, kernel_fn, ...)
        Y[, n] <- mvrnorm(1, mu = rep(0, times = length(x)), Sigma = K)
    }
    Y
}

x <- seq(0, 1, length.out = 51)  # x-coordinates

model1 <- function(obs, x) {
    model1_data <- matrix(NA, nrow = obs, ncol = length(x))
    for(i in 1:obs){
        e <- draw_samples(x, 1, kernel_fn = kernel_1)
        X <- c()
        for (p in 1:length(x)){
            t <- x[p]
            val <- (4*t) + e[p,]
            X = c(X, val)
        }
        model1_data[i,] <- X
  }
    model1_data
}

# model1(100, x)

Surely you just need to remove set.seed from draw.samples to prevent you from drawing the same samples every time the function is called? — Allan Cameron
Indeed @AllanCameron, that's what I was going to say. I don't understand why the 100 trajectories are not exactly the same. — Stéphane Laurent
@StéphaneLaurent yes, confirmed - see below. The OP seemes to have plotted the data plus some sort of smoothed line, but the code for the plot isn't provided. — Allan Cameron

Allan Cameron Allan Cameron · Accepted Answer · 2020-03-06T22:46:06

Because you have set.seed in draw_samples, you are getting the same random numbers with each draw. If you remove it, then you can do:

a <- model1(100, x)
matplot(t(a), type = "l", col = 'gray')

to get

Generate functional data from Gaussian Process in R

1 Answers