if(!require("devtools")) install.packages("devtools")
devtools::install_github("emanuelhuber/GauProMod")
devtools::install_local("/media/huber/Elements/UNIBAS/software/codeR/package_GauProMod/GauProMod")
library(GauProMod)
library(plot3D)
library(RColorBrewer)
#------------------------------------------------------------------------------#
#observations
obs <- list(x=c(-4, -3, -1, 0, 4),
y=c(-2, 0, 1, 2, 0))
# targets
targ <- list("x"=seq(-4.5, 4.5, len = 200))
covModel <- list(kernel="gaussian",
l = 0.7, # correlation length
h = 1.) # std. deviation
# linear mean function
op <- 2
sigma <- 0.2
GP <- gpCond(obs = obs, targ = targ, covModels=list(pos=covModel),
sigma = sigma, op = op)
names(GP)
# GP$mean = mean value at location xstar
# GP$cov = covariance matrix of the conditioned GP
# GP$logLik = log-likelihood of the conditioned GP
# GP$xstar = x-coordinates at which the GP is simulated
#--- plot mean +/- sd
xp <-(GP$mean + sqrt(diag(GP$cov))) # mean + sd
xm <-(GP$mean - sqrt(diag(GP$cov))) # mean - sd
# initialise the plot
plot(cbind(obs$x, obs$y), type="p", xlab="x", ylab="y",
xlim = range(c(obs$x, targ$x)), ylim = range(c(xp, xm, obs$y)),
pch = 20, col = "black", main = GP$logLik)
lines(GP$xstar, GP$mean,col="red") # mean
lines(GP$xstar, xm,lty=3) # + sd
lines(GP$xstar, xp,lty=3) # - sd
legend("topleft", legend = c("obs", "mean", "sd"), lty = c(NA, 1, 3),
pch = c(20, NA, NA), col=c("black", "red", "black"), bty="n")
#------------------------------------------------------------------------------#
M <- matrix(c(1L, 0L, 0L, 1L), nrow = 2)
cholfac(M)
M1 <- M
storage.mode(M1) <- "numeric"
cholfac(M1)
obs <- list(x = cbind(c(2.17, 7.92, 8.98, 7.77, 2.79, 5.36, 4.27, 3.07, 6.31),
c(1.33, 7.24, 4.26, 2.67, 6.17, 8.04, 3.18, 5.63, 8.33)),
y = c(2.60, 1.48, 1.36, 8.61, 1.00, 1.58, 8.42, 8.39, 1.50))
targ <- list(x = cbind(c(2.17, 7.92, 8.98, 7.77, 2.79, 5.36, 4.27, 3.07, 6.31,
3.74, 5.93, 7.19, 6.61, 5.54, 2.27, 1.61, 4.02, 1.06),
c(1.33, 7.24, 4.26, 2.67, 6.17, 8.04, 3.18, 5.63, 8.33,
6.34, 3.68, 6.82, 1.79, 8.60, 7.73, 5.35, 2.45, 4.92))
)
# Matern kernel
covModel <- list(kernel="matern",
l = 5, # correlation length
v = 1, # smoothness
h = 2.45 # std. deviation
)
# 2D linear mean function
op <- 2
# Gaussian likelihood
sigma <- 0.2
GP <- gpCond(obs = obs, targ = targ, covModels=list(pos=covModel),
sigma = sigma, op = op)
names(GP)
# GP$mean = mean value at location xstar
# GP$cov = covariance matrix of the conditioned GP
# GP$logLik = log-likelihood of the conditioned GP
# GP$xstar = x-coordinates at which the GP is simulated
# mean
Ymean <- GP$mean
# standard deviation
YSD <- sqrt(diag(GP$cov))
ysdplus <- Ymean - 1.95* YSD
ysdminus <- Ymean + 1.95* YSD
par(mfrow = c(1,1))
ylim <- range(Ymean, obs$y)
plot3D::scatter3D(x = targ$x[,1], y = targ$x[,2], z = Ymean, clim = ylim,
pch = 20)
plot3D::arrows3D(x0 = targ$x[,1], y0 = targ$x[,2], z0 = ysdminus,
x1 = targ$x[,1], y1 = targ$x[,2], z1 = ysdplus,
length=0.05, angle=90, code=3, add = TRUE, col = "black")
# large dots = observations
plot3D::scatter3D(x = obs$x[,1], y = obs$x[,2], z = obs$y, add = TRUE,
pch = 20, cex = 3, clim = ylim)
par(mfrow = c(1, 2))
ylim <- range(ysdplus, ysdminus, obs$y)
plot(targ$x[,1], Ymean, type = "p", ylim = ylim, pch = 3, cex = 0.5)
arrows(targ$x[,1], ysdminus, targ$x[,1], ysdplus, length=0.05, angle=90, code=3)
points(obs$x[,1], obs$y, col = "dodgerblue", pch = 20)
plot(targ$x[,2], Ymean, type = "p", ylim = ylim, pch = 3, cex = 0.5)
arrows(targ$x[,2], ysdminus, targ$x[,2], ysdplus, length=0.05, angle=90, code=3)
points(obs$x[,2], obs$y, col = "dodgerblue", pch = 20)
L <- cholfac(GP$cov)
ystar <- gpSim(GP , L = L)
par(mfrow = c(1,1))
ylim <- range(ystar[,3], obs$y)
plot3D::scatter3D(x = targ$x[,1], y = targ$x[,2], z = ystar[,3], clim = ylim,
pch = 18)
# dots = observations
plot3D::scatter3D(x = obs$x[,1], y = obs$x[,2], z = obs$y, add = TRUE,
pch = 20, cex = 3, clim = ylim)
#observations
obs <- list(x = cbind(c(2.17, 7.92, 8.98, 7.77, 2.79, 5.36, 4.27, 3.07, 6.31,
3.74, 5.93, 7.19, 6.61, 5.54, 2.27, 1.61, 4.02, 1.06),
c(1.33, 7.24, 4.26, 2.67, 6.17, 8.04, 3.18, 5.63, 8.33,
6.34, 3.68, 6.82, 1.79, 8.60, 7.73, 5.35, 2.45, 4.92)),
y = c(2.60, 1.48, 1.36, 8.61, 1.00, 1.58, 8.42, 8.39, 1.50,
9.05, 1.14, 1.49, 9.19, 1.32, 1.03, 6.41, 6.16, 5.42))
# targets (=2D mesh)
vx <- seq(0, 10, by = 0.5)
vy <- seq(0, 10, by = 0.5)
targ <- list(x = vecGrid(vx, vy))
# linear kernel
covModel <- list(kernel="linear",
b = 1, # slope
h = 0.5, # std. deviation
c = 1) # constant
# Matern kernel
covModel <- list(kernel="matern",
l = 1, # correlation length
v = 2.5, # smoothness
h = 2.45) # std. deviation
# squared exponential kernel (Gaussian)
covModel <- list(kernel="gaussian",
l = 0.5, # correlation length
h = 0.25) # std. deviation
# 2D quadratic mean function
op <- 5
# zero-mean function (no trend)
op <- 0
# 2D linear mean function
op <- 2
sigma <- 0.2
GP <- gpCond(obs = obs, targ = targ, covModels=list(pos=covModel),
sigma = sigma, op = op)
names(GP)
# GP$mean = mean value at location xstar
# GP$cov = covariance matrix of the conditioned GP
# GP$logLik = log-likelihood of the conditioned GP
#
# mean
Ymean <- matrix(GP$mean, nrow = length(vx), ncol = length(vy), byrow = TRUE)
# standard deviation
YSD <- matrix(sqrt(diag(GP$cov)), nrow = length(vx), ncol = length(vy),
byrow = TRUE)
par(mfrow = c(2,2))
plot3D::image2D(x = vx, y = vy, z = Ymean, asp=1)
points(obs$x, col="white",pch=3)
title(main = "mean")
plot3D::contour2D(x = vx, y = vy, Ymean, asp=1)
points(obs$x, col="black",pch=3)
rect(vx[1], vy[1], vx[length(vx)], vy[length(vy)])
title(main = "mean")
plot3D::image2D(x = vx, y = vy, z = YSD, asp=1)
points(obs$x, col="white",pch=3)
title(main = "standard deviation")
plot3D::contour2D(x = vx, y = vy, YSD, asp=1)
points(obs$x, col="black",pch=3)
rect(vx[1], vy[1], vx[length(vx)], vy[length(vy)])
title(main = "standard deviation")
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.