Генерация функциональных данных из гауссовского процесса в R

Question

Модель:

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

t € [0; 1]

e(t) - это гауссовский процесс с нулевым средним и ковариационной функцией f(s, t) = exp( -|t - s| )

Окончательный результат за 100 прогонов (= 100 серых линий ) с 50 отобранными точками каждая должна быть похожа на серую область на рисунке.

Зеленая линия - это то, что я получаю из кода ниже.

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)

Answer 1

Поскольку у вас set.seed в draw_samples, вы получаете одни и те же случайные числа при каждом розыгрыше. Если вы удалите его, то вы можете сделать:

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

, чтобы получить