# R/predict.binaryGP.R In binaryGP: Fit and Predict a Gaussian Process Model with (Time-Series) Binary Response

#### Documented in predict.binaryGP

#' Predictions of Binary Gaussian Process
#'
#' @description The function computes the predicted response and its variance as well as its confidence interval.
#'
#' @param object a class binaryGP object estimated by \code{binaryGP_fit}.
#' @param xnew a testing matrix with dimension \code{n_new} by \code{d} in which each row corresponds to a predictive location.
#' @param conf.level a value from 0 to 1 specifying the level of confidence interval. The default is 0.95.
#' @param sim.number a positive integer specifying the simulation number for Monte-Carlo method. The default is 101.
#' @param ... for compatibility with generic method \code{predict}.
#'
#' @return
#' \item{mean}{a matrix with dimension \code{n_new} by \code{T} displaying predicted responses at locations \code{xnew}.}
#' \item{var}{a matrix with dimension \code{n_new} by \code{T} displaying predictive variances at locations \code{xnew}.}
#' \item{upper.bound}{a matrix with dimension \code{n_new} by \code{T} displaying upper bounds with \code{conf.level} confidence level.}
#' \item{lower.bound}{a matrix with dimension \code{n_new} by \code{T} displaying lower bounds with \code{conf.level} confidence level.}
#' \item{y_pred}{a matrix with dimension \code{n_new} by \code{T} displaying predicted binary responses at locations \code{xnew}.}
#'
#' @seealso \code{\link{binaryGP_fit}} for estimation of the binary Gaussian process.
#' @author Chih-Li Sung <[email protected]>
#'
#' @import logitnorm
#' @examples
#'
#' library(binaryGP)
#'
#' #####      Testing function: cos(x1 + x2) * exp(x1*x2) with TT sequences      #####
#' #####   Thanks to Sonja Surjanovic and Derek Bingham, Simon Fraser University #####
#' test_function <- function(X, TT)
#' {
#'   x1 <- X[,1]
#'   x2 <- X[,2]
#'
#'   eta_1 <- cos(x1 + x2) * exp(x1*x2)
#'
#'   p_1 <- exp(eta_1)/(1+exp(eta_1))
#'   y_1 <- rep(NA, length(p_1))
#'   for(i in 1:length(p_1)) y_1[i] <- rbinom(1,1,p_1[i])
#'   Y <- y_1
#'   P <- p_1
#'   if(TT > 1){
#'     for(tt in 2:TT){
#'       eta_2 <- 0.3 * y_1 + eta_1
#'       p_2 <- exp(eta_2)/(1+exp(eta_2))
#'       y_2 <- rep(NA, length(p_2))
#'       for(i in 1:length(p_2)) y_2[i] <- rbinom(1,1,p_2[i])
#'       Y <- cbind(Y, y_2)
#'       P <- cbind(P, p_2)
#'       y_1 <- y_2
#'     }
#'   }
#'
#'   return(list(Y = Y, P = P))
#' }
#'
#' set.seed(1)
#' n <- 30
#' n.test <- 10
#' d <- 2
#' X <- matrix(runif(d * n), ncol = d)
#' X.test <- matrix(runif(d * n.test), ncol = d)
#'
#' ##### without time-series #####
#' Y <- test_function(X, 1)$Y ## Y is a vector #' test.out <- test_function(X.test, 1) #' Y.test <- test.out$Y
#' P.true <- test.out$P #' \donttest{ #' # fitting #' binaryGP.model <- binaryGP_fit(X = X, Y = Y) #' #' # prediction #' binaryGP.prediction <- predict(binaryGP.model, xnew = X.test) #' print(binaryGP.prediction$mean)
#' print(binaryGP.prediction$var) #' print(binaryGP.prediction$upper.bound)
#' print(binaryGP.prediction$lower.bound) #' #' ##### with time-series ##### #' Y <- test_function(X, 10)$Y  ## Y is a matrix with 10 columns
#' test.out <- test_function(X.test, 10)
#' Y.test <- test.out$Y #' P.true <- test.out$P
#'
#' # fitting
#' binaryGP.model <- binaryGP_fit(X = X, Y = Y, R = 1)
#'
#' # prediction
#' binaryGP.prediction <- predict(binaryGP.model, xnew = X.test)
#' print(binaryGP.prediction$mean) #' print(binaryGP.prediction$var)
#' print(binaryGP.prediction$upper.bound) #' print(binaryGP.prediction$lower.bound)
#' }
#' @export predict.binaryGP
#' @export
#'

predict.binaryGP <- function(object, xnew, conf.level = 0.95, sim.number = 101, ...){

if (is.binaryGP(object) == FALSE) {
stop("The object in question is not of class binaryGP \n")
}
if (is.matrix(xnew) == FALSE) {
xnew <- as.matrix(xnew)
}

X            <- object$X Y <- object$Y
R            <- object$R L <- object$L
corr         <- object$corr constantMean <- object$constantMean
alpha_hat    <- object$alpha_hat varphi_hat <- object$varphi_hat
gamma_hat    <- object$gamma_hat theta_hat <- object$theta_hat
rho_hat      <- -log10(theta_hat)
sigma_hat    <- object$sigma_hat nugget <- object$nugget_hat
eta_hat      <- object$eta_hat M.mx <- object$Model.mat
standardize  <- object$standardize orthogonalGP <- object$orthogonalGP
mean.x       <- object$mean.x range.x <- object$range.x

###################       Setting       ###################
n <- nrow(X)
n.test <- nrow(xnew)
d <- ncol(X)
T_val <- ncol(Y)
if(is.null(T_val)) {
T_val <- 1
Y <- matrix(Y, ncol = 1)
}
N <- n * T_val

if(standardize){
xnew <- t((t(xnew) - mean.x)/range.x)
}

R.mx <- corr_matrix(X, rho_hat, corr)
R.mx <- R.mx + diag(nugget/(sigma_hat)^2, n)  # add a nugget term to make the variance-covariance matrix more stable
if(orthogonalGP) R.mx <- R.mx - orthogonalize(X, rho_hat, corr)
R_inv <- solve(R.mx)

p.vt[p.vt > 1-1e-8] <- 1-1e-8
p.vt[p.vt < 1e-8] <- 1e-8
p.mx <- matrix(p.vt, ncol = T_val)  ## initial p

beta_hat <- c(alpha_hat, varphi_hat, gamma_hat)
mu <- c(M.mx %*% beta_hat)
mu.mx <- matrix(mu, ncol = T_val)

P.new <- P_mean.new <- P_var.new <-
P_quantile_u.new <- P_quantile_l.new <- Y.new <- array(NA, c(n.test, sim.number, T_val))

for(j in 1:sim.number){
###################      MH algorithm     ###################
for(tt in 1:T_val){
for(k in 1:n){
Q_kk <- R_inv[k,k]
Q_kj <- R_inv[k,-k]
m_D <- mu.mx[k,tt] - (sum(Q_kj * (log(p.mx[-k,tt]) - log(1-p.mx[-k,tt]) - mu.mx[-k,tt])))/Q_kk
v_D <- sigma_hat^2/Q_kk

p_star <- rlogitnorm(n = 1, mu = m_D, sigma = sqrt(v_D))
p_star <- min(c(p_star, 1-1e-8))
p_star <- max(c(p_star, 1e-8))
U <- runif(1, 0, 1)
f_1 <- dbinom(Y[k,tt], size = 1, prob = p_star)
f_2 <- dbinom(Y[k,tt], size = 1, prob = p.mx[k,tt])
if(U < min(c(1, f_1/f_2))){
p.mx[k,tt] <- p_star
}
}
}
###################      w/o time-series: monte carlo method                                           ###################
###################      w/ time-series : draw conditional distribution and then monte carlo method    ###################
for(i in 1:n.test){
x.new <- xnew[i,,drop=FALSE]
p.new <- p_mean.new <- p_var.new <-
p_quantile_u.new  <- p_quantile_l.new <- y.new <- rep(NA, T_val)
r.vt <- corr_vec(x.new, X, rho_hat, corr)
if(orthogonalGP) r.vt <- r.vt - orthogonalize_vec(x.new, X, rho_hat, corr)
k.vt <- r.vt %*% R_inv
rRinvr <- c(k.vt %*% t(r.vt))
for(tt in 1:T_val){
if(constantMean)  M.mx.new <- matrix(1, ncol = 1, nrow = 1) else M.mx.new <- cbind(1, x.new)
if(R > 0){
if(tt <= R){
ytmp <- rep(0,R)
ytmp[0:length(y.new[(tt-1):max(tt-R, 0)])] <- y.new[(tt-1):max(tt-R, 0)]
}else{
ytmp <- y.new[(tt-1):(tt-R)]
}
M.mx.new <- cbind(M.mx.new, matrix(ytmp, nrow=1))
}
if(L > 0){
if(tt <= L){
ytmp <- rep(0,L)
ytmp[0:length(y.new[(tt-1):max(tt-L, 0)])] <- y.new[(tt-1):max(tt-L, 0)]
}else{
ytmp <- y.new[(tt-1):(tt-L)]
}
xy.tmp <- matrix(0, nrow = 1, ncol = d * L)
for(l in 1:L) xy.tmp[1,(d *(l-1) + 1):(d * l)] <- x.new * ytmp[l]
#xy.tmp <- xy.tmp[,c(ncol(xy.tmp)-1), drop=FALSE] #temp
M.mx.new <- cbind(M.mx.new, xy.tmp)
}

m_D <- c(M.mx.new %*% beta_hat) + k.vt %*% (log(p.mx[,tt])-log(1-p.mx[,tt])-mu.mx[,tt])
v_D <- sigma_hat^2 * max((1 - rRinvr), 0)

p_mean.new[tt] <- momentsLogitnorm(mu = m_D, sigma = sqrt(v_D), stop.on.error = FALSE)[1]
p_var.new[tt] <- momentsLogitnorm(mu = m_D, sigma = sqrt(v_D), stop.on.error = FALSE)[2]

if(sqrt(v_D) < 1e-3 | p_mean.new[tt] < 1e-4){ ### bug from momentsLogitnorm package due to sigma is too small
p_mean.new[tt] <- rlogitnorm(n = 1, mu = m_D, sigma = sqrt(v_D))
}

p_quantile_u.new[tt] <- qlogitnorm(p = 1 - (1 - conf.level)/2, mu = m_D, sigma = sqrt(v_D))
p_quantile_l.new[tt] <- qlogitnorm(p = (1 - conf.level)/2, mu = m_D, sigma = sqrt(v_D))

if(T_val > 1){
### update if time-series is included ###
p.new[tt] <- rlogitnorm(n = 1, mu = m_D, sigma = sqrt(v_D))
y.new[tt] <- rbinom(n = 1, size = 1, prob = p.new[tt])
}else{
p.new[tt] <- rlogitnorm(n = 1, mu = m_D, sigma = sqrt(v_D))
y.new[tt] <- rbinom(n = 1, size = 1, prob = p_mean.new[tt])
}
}
P_mean.new[i, j, ] <- p_mean.new
P_var.new[i, j, ] <- p_var.new
P_quantile_u.new[i, j, ] <- p_quantile_u.new
P_quantile_l.new[i, j, ] <- p_quantile_l.new
P.new[i, j, ] <- p.new
Y.new[i, j, ] <- y.new
}
}

if(T_val > 1){
p_new <- apply(P.new, c(1,3), median)
p_var_new <- apply(P.new, c(1,3), var)
p_quantile_u_new <- apply(P.new, c(1,3), quantile, probs = 1 - (1 - conf.level)/2)
p_quantile_l_new <- apply(P.new, c(1,3), quantile, probs = (1 - conf.level)/2)
}else{
p_new <- apply(P_mean.new, c(1,3), median)
p_var_new <- apply(P_var.new, c(1,3), median) + apply(P_mean.new, c(1,3), var)
p_quantile_u_new <- apply(P.new, c(1,3), quantile, probs = 1 - (1 - conf.level)/2)
p_quantile_l_new <- apply(P.new, c(1,3), quantile, probs = (1 - conf.level)/2)
}

y_new <- apply(Y.new, c(1,3), median)

return(list(mean = p_new, var = p_var_new, upper.bound = p_quantile_u_new, lower.bound = p_quantile_l_new, y_pred = y_new))
}


## Try the binaryGP package in your browser

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

binaryGP documentation built on Sept. 19, 2017, 9:02 a.m.