从 R 中的高斯过程生成函数数据

Generate functional data from Gaussian Process in R

型号:

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

t € [0; 1]

e(t) 是具有零均值和协方差函数的高斯过程 f(s, t) = exp( -|t - s| )

每次50个采样点运行100次(=100条灰线)的最终结果应该像图中的灰色区域。

绿线是我从下面的代码中得到的。

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)

因为您在 draw_samples 中有 set.seed,所以每次抽奖都会得到相同的随机数。如果你删除它,那么你可以这样做:

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

得到