R/gibbs.R In deepgp: Deep Gaussian Processes using MCMC

Defines functions gibbs_three_layergibbs_two_layergibbs_one_layer_sepgibbs_one_layer

```# Function Contents -----------------------------------------------------------
# Internal:
#   gibbs_one_layer
#   gibbs_two_layer
#   gibbs_three_layer

# One layer Gibbs -------------------------------------------------------------

gibbs_one_layer <- function(x, y, nmcmc, verb, initial, true_g, settings, v) {

dx <- sq_dist(x)
g <- vector(length = nmcmc)
if (is.null(true_g)) g[1] <- initial\$g else g[1] <- true_g
theta <- vector(length = nmcmc)
theta[1] <- initial\$theta
tau2 <- vector(length = nmcmc)
tau2[1] <- initial\$tau2
ll <- NULL

for (j in 2:nmcmc) {

if(verb) if(j %% 500 == 0) cat(j, '\n')

# Sample nugget (g)
if (is.null(true_g)) {
samp <- sample_g(y, dx, g[j - 1], theta[j - 1], alpha = settings\$alpha\$g,
beta = settings\$beta\$g, l = settings\$l, u = settings\$u,
ll_prev = ll, v = v)
g[j] <- samp\$g
ll <- samp\$ll
} else g[j] <- true_g

# Sample length scale (theta)
samp <- sample_theta(y, dx, g[j], theta[j - 1],
alpha = settings\$alpha\$theta,
beta = settings\$beta\$theta, l = settings\$l,
u = settings\$u, outer = TRUE, ll_prev = ll, v = v,
tau2 = TRUE)
theta[j] <- samp\$theta
ll <- samp\$ll
if (is.null(samp\$tau2)) tau2[j] <- tau2[j - 1] else tau2[j] <- samp\$tau2
} # end of j for loop

return(list(g = g, theta = theta, tau2 = tau2))
}

# One layer Gibbs SEPARABLE ---------------------------------------------------

gibbs_one_layer_sep <- function(x, y, nmcmc, verb, initial, true_g, settings, v) {

d <- ncol(x)
g <- vector(length = nmcmc)
if (is.null(true_g)) g[1] <- initial\$g else g[1] <- true_g
theta <- matrix(nrow = nmcmc, ncol = d)
if (length(initial\$theta) == 1) initial\$theta <- rep(initial\$theta, d)
theta[1, ] <- initial\$theta
tau2 <- vector(length = nmcmc)
tau2[1] <- initial\$tau2
ll <- NULL

for (j in 2:nmcmc) {

if(verb) if(j %% 500 == 0) cat(j, '\n')

# Sample nugget (g)
if (is.null(true_g)) {
samp <- sample_g_sep(y, x, g[j - 1], theta[j - 1, ], alpha = settings\$alpha\$g,
beta = settings\$beta\$g, l = settings\$l, u = settings\$u,
ll_prev = ll, v = v)
g[j] <- samp\$g
ll <- samp\$ll
} else g[j] <- true_g

# Sample length scale (theta)
for (i in 1:d) {
samp <- sample_theta_sep(y, x, g[j], theta[j - 1, ], index = i,
alpha = settings\$alpha\$theta,
beta = settings\$beta\$theta, l = settings\$l,
u = settings\$u, ll_prev = ll, v = v,
tau2 = (i == d))
theta[j, i] <- samp\$theta
ll <- samp\$ll
if (i == 1) { # update tau2 (repeat original value if nothing was accepted)
if (is.null(samp\$tau2)) tau2[j] <- tau2[j - 1] else tau2[j] <- samp\$tau2
} else { # only update tau2 if there was an acceptance
if (!is.null(samp\$tau2)) tau2[j] <- samp\$tau2
}
}
} # end of j for loop

return(list(g = g, theta = theta, tau2 = tau2))
}

# Two layer Gibbs -------------------------------------------------------------

gibbs_two_layer <- function(x, y, nmcmc, D, verb, initial, true_g,
settings, v) {

dx <- sq_dist(x)
dw <- sq_dist(initial\$w)

g <- vector(length = nmcmc)
if (is.null(true_g)) g[1] <- initial\$g else g[1] <- true_g
theta_y <- vector(length = nmcmc)
theta_y[1] <- initial\$theta_y
theta_w <- matrix(nrow = nmcmc, ncol = D)
theta_w[1, ] <- initial\$theta_w
w <- list()
w[[1]] <- initial\$w
tau2 <- vector(length = nmcmc)
tau2[1] <- initial\$tau2
ll_outer <- NULL

for (j in 2:nmcmc) {

if(verb) if(j %% 500 == 0) cat(j, '\n')

# Sample nugget (g)
if (is.null(true_g)) {
samp <- sample_g(y, dw, g[j - 1], theta_y[j - 1],
alpha = settings\$alpha\$g, beta = settings\$beta\$g,
l = settings\$l, u = settings\$u, ll_prev = ll_outer,
v = v)
g[j] <- samp\$g
ll_outer <- samp\$ll
} else g[j] <- true_g

# Sample outer length scale (theta_y)
samp <- sample_theta(y, dw, g[j], theta_y[j - 1],
alpha = settings\$alpha\$theta_y,
beta = settings\$beta\$theta_y, l = settings\$l,
u = settings\$u, outer = TRUE, ll_prev = ll_outer,
v = v, tau2 = TRUE)
theta_y[j] <- samp\$theta
ll_outer <- samp\$ll
if (is.null(samp\$tau2)) tau2[j] <- tau2[j - 1] else tau2[j] <- samp\$tau2

# Sample inner length scale (theta_w) - separately for each dimension
for (i in 1:D) {
samp <- sample_theta(w[[j - 1]][, i], dx, g = eps, theta_w[j - 1, i],
alpha = settings\$alpha\$theta_w,
beta = settings\$beta\$theta_w, l = settings\$l,
u = settings\$u, outer = FALSE, v = v)
theta_w[j, i] <- samp\$theta
}

# Sample hidden Gaussian layer (w)
samp <- sample_w(y, w[[j - 1]], dw, dx, g[j], theta_y[j], theta_w[j, ],
ll_prev = ll_outer, v = v)
w[[j]] <- samp\$w
ll_outer <- samp\$ll
dw <- samp\$dw
} # end of j for loop

return(list(g = g, theta_y = theta_y, theta_w = theta_w, w = w, tau2 = tau2))
}

# Three layer Gibbs -----------------------------------------------------------

gibbs_three_layer <- function(x, y, nmcmc, D, verb, initial, true_g,
settings, v) {

dx <- sq_dist(x)
dz <- sq_dist(initial\$z)
dw <- sq_dist(initial\$w)

g <- vector(length = nmcmc)
if (is.null(true_g)) g[1] <- initial\$g else g[1] <- true_g
theta_y <- vector(length = nmcmc)
theta_y[1] <- initial\$theta_y
theta_w <- matrix(nrow = nmcmc, ncol = D)
theta_w[1, ] <- initial\$theta_w
theta_z <- matrix(nrow = nmcmc, ncol = D)
theta_z[1, ] <- initial\$theta_z
w <- list()
w[[1]] <- initial\$w
z <- list()
z[[1]] <- initial\$z
tau2 <- vector(length = nmcmc)
tau2[1] <- initial\$tau2
ll_outer <- NULL

for (j in 2:nmcmc) {

if(verb) if(j %% 500 == 0) cat(j, '\n')

# Sample nugget (g)
if (is.null(true_g)) {
samp <- sample_g(y, dw, g[j - 1], theta_y[j - 1],
alpha = settings\$alpha\$g, beta = settings\$beta\$g,
l = settings\$l, u = settings\$u, ll_prev = ll_outer,
v = v)
g[j] <- samp\$g
ll_outer <- samp\$ll
} else g[j] <- true_g

# Sample outer length scale (theta_y)
samp <- sample_theta(y, dw, g[j], theta_y[j - 1],
alpha = settings\$alpha\$theta_y,
beta = settings\$beta\$theta_y, l = settings\$l,
u = settings\$u, outer = TRUE, ll_prev = ll_outer,
v = v, tau2 = TRUE)
theta_y[j] <- samp\$theta
ll_outer <- samp\$ll
if (is.null(samp\$tau2)) tau2[j] <- tau2[j - 1] else tau2[j] <- samp\$tau2

# Sample middle length scale (theta_w)
ll_mid <- 0 # re-calculated each time since we have a new z
for (i in 1:D) {
samp <- sample_theta(w[[j - 1]][, i], dz, g = eps, theta_w[j - 1, i],
alpha = settings\$alpha\$theta_w,
beta = settings\$beta\$theta_w, l = settings\$l,
u = settings\$u, outer = FALSE, v = v)
theta_w[j, i] <- samp\$theta
ll_mid <- ll_mid + samp\$ll
}

# Sample inner length scale (theta_z)
for (i in 1:D) {
samp <- sample_theta(z[[j - 1]][, i], dx, g = eps, theta_z[j - 1, i],
alpha = settings\$alpha\$theta_z,
beta = settings\$beta\$theta_z, l = settings\$l,
u = settings\$u, outer = FALSE, v = v)
theta_z[j, i] <- samp\$theta
}

# Sample inner hidden Gaussian layer (z)
samp <- sample_z(w[[j - 1]], z[[j - 1]], dz, dx, g = eps, theta_w[j, ],
theta_z[j, ], ll_prev = ll_mid, v = v)
z[[j]] <- samp\$z
dz <- samp\$dz

# Sample middle hidden Gaussian layer (w)
samp <- sample_w(y, w[[j - 1]], dw, dz, g[j], theta_y[j], theta_w[j, ],
ll_prev = ll_outer, v = v)
w[[j]] <- samp\$w
ll_outer <- samp\$ll
dw <- samp\$dw
} # end of j for loop

return(list(g = g, theta_y = theta_y, theta_w = theta_w, theta_z = theta_z,
w = w, z = z, tau2 = tau2))
}
```

Try the deepgp package in your browser

Any scripts or data that you put into this service are public.

deepgp documentation built on Dec. 28, 2022, 1:32 a.m.