Модель регрессии гауссовского процесса с несколькими переменными в R

Question

Я работал над GP с несколькими переменными, я не уверен, должен ли я использовать тот же код, что и одномерный GP, или есть другие методы. Я пытаюсь использовать массив как результат функции ядра, однако он не работает. Я показываю код ниже, пожалуйста, скажите мне, где я делаю ошибки.

## Univariate GP
##   
## a covariance matrix
calcSigma <- function(X1,X2,l=1) {
      Sigma <- matrix(rep(0, length(X1)*length(X2)), nrow=length(X1))
      for (i in 1:nrow(Sigma)) {
        for (j in 1:ncol(Sigma)) {
          Sigma[i,j] <- exp(-0.5*(abs(X1[i]-X2[j])/l)^2)
        }
      }
      return(Sigma)
    }

# Define the points at which we want to define the functions
x.star <- seq(-5,5,len=50)
# with unknown data
Sigma <- calcSigma(x.star,x.star)
numSamps<-5

####
library(gptk)
value.gp<- gaussSamp(mu=matrix(0,nrow=dim(Sigma)[1]), Sigma, numSamps)
value.gp<-t(value.gp)
value.gp<- cbind(x=x.star,as.data.frame(value.gp))
value.gp <- melt(value.gp,id="x")

# Plot the result
fig2a <- ggplot(value.gp,aes(x=x,y=value)) +
  geom_rect(xmin=-Inf, xmax=Inf, ymin=-2, ymax=2, fill="grey80") +
  geom_line(aes(group=variable)) +
  theme_bw() +
  scale_y_continuous(lim=c(-2.5,2.5), name="output, f(x)") +
  xlab("input, x")
fig2a
############################################################################
##
##  With known data
##
f <- data.frame(x=c(-4,-3,-1,0,2),
                y=c(-2,0,1,2,-1))

# Calculate the covariance matrices
# using the same x.star values as above
x <- f$x
k.xx <- calcSigma(x,x)
k.xxs <- calcSigma(x,x.star)
k.xsx <- calcSigma(x.star,x)
k.xsxs <- calcSigma(x.star,x.star)

# These matrix calculations correspond to equation (2.19)
# in the book.
f.star.bar <- k.xsx%*%solve(k.xx)%*%f$y
cov.f.star <- k.xsxs - k.xsx%*%solve(k.xx)%*%k.xxs

# This time we'll plot more samples.  We could of course
# simply plot a +/- 2 standard deviation confidence interval
# as in the book but I want to show the samples explicitly here.
numSamps.01 <- 5
value.gp1<- gaussSamp(mu=f.star.bar , cov.f.star, numSamps.01)
value.gp1<-t(value.gp1)
value.gp1<- cbind(x=x.star,as.data.frame(value.gp1))
value.gp1 <- melt(value.gp1,id="x")
###########
fig2b <- ggplot(value.gp1,aes(x=x,y=value)) +
  geom_line(aes(group=variable), colour="grey80") +
  geom_point(data=f,aes(x=x,y=y),pch = 3,cex = 3) +
  theme_bw() + geom_line(aes(colour=variable)) +
  scale_y_continuous(lim=c(-3,3), name="output, f(x)") +
  xlab("input, x")
fig2b
############################################################################
##
##  With noisy observations
##
sigma.n <- 0.1

# Recalculate the mean and covariance functions
f.bar.star <- k.xsx%*%solve(k.xx + sigma.n^2*diag(1, ncol(k.xx)))%*%f$y
cov.f.star <- k.xsxs - k.xsx%*%solve(k.xx + sigma.n^2*diag(1, ncol(k.xx)))%*%k.xxs
value.gp2<- gaussSamp(mu=f.star.bar , cov.f.star, numSamps.01)
value.gp2<-t(value.gp2)
value.gp2<- cbind(x=x.star,as.data.frame(value.gp2))
value.gp2 <- melt(value.gp2,id="x")

fig2c <- ggplot(value.gp2,aes(x=x,y=value)) +
  geom_line(aes(group=variable), colour="grey80") +
  geom_point(data=f,aes(x=x,y=y),pch = 3,cex = 3) +
  theme_bw() + geom_line(aes(colour=variable)) +
  scale_y_continuous(lim=c(-3,3), name="output, f(x)") +
  geom_errorbar(data=f,aes(x=x,y=NULL,ymin=y-2*sigma.n, ymax=y+2*sigma.n), width=0.2) +
  xlab("input, x")
##############################################################################
##
## multivariate GP
##
##############################################################################
calcSigma.multi <- function(X1,X2,l=1) {
  Sigma <- array(NA, dim = c(1,ncol(X1), nrow(X1),nrow(X2)))
  for (i in 1:nrow(X1)) {
    for (j in 1:nrow(X2)) {
      Sigma[,,i,j] <- exp(-0.5*(abs(X1[i,]-X2[j,])/l)^2)
    }
  }
  return(Sigma)
}

x.star <- matrix(seq(-5,5,len=60),nrow=6,ncol=10)
##nrow = the number of data set(observations) we have
##ncol = the number of variables(x1,x2,x3,...etc.)
Sigma.multi <- calcSigma.multi(x.star,x.star)
numSamps<-nrow(x.star)
####
value.gp<- gaussSamp(mu=matrix(0,nrow=dim(Sigma.multi)[3]), Sigma.multi, numSamps)  ##stuck here!!!!