From 06673a48f68b5be9e6ba480e8482db0dbd9d0c3a Mon Sep 17 00:00:00 2001 From: sciome-bot Date: Fri, 12 Jan 2024 19:23:11 -0500 Subject: [PATCH] Fix all indentation errors --- .lintr | 1 + R/Log_Likelihood.R | 9 +- R/PrestoGP_CreateU_Multivariate.R | 12 +- R/PrestoGP_Model.R | 18 +- R/PrestoGP_Multivariate_Vecchia.R | 137 ++-- R/PrestoGP_Util_Functions.R | 3 +- R/PrestoGP_Vecchia.R | 116 +-- R/RcppExports.R | 4 +- tests/testthat/sim_multivariate.R | 30 +- tests/testthat/sim_multivariate_big.R | 70 +- tests/testthat/sim_multivariate_big_pred.R | 90 ++- tests/testthat/sim_multivariate_big_st.R | 78 +-- tests/testthat/sim_multivariate_lik.R | 64 +- tests/testthat/sim_multivariate_small.R | 70 +- tests/testthat/sim_multivariate_small_pred.R | 86 ++- tests/testthat/sim_vecchia.R | 10 +- tests/testthat/sim_vecchia_pred.R | 10 +- tests/testthat/sim_vecchia_small.R | 10 +- tests/testthat/sim_vecchia_small_pred.R | 10 +- tests/testthat/sim_vecchia_st.R | 10 +- tests/testthat/test-Log_Likelihood.R | 663 +++++++++--------- .../test-PrestoGP_CreateU_Multivariate.R | 270 +++---- 22 files changed, 857 insertions(+), 914 deletions(-) diff --git a/.lintr b/.lintr index 9139bc8..6c1cc72 100644 --- a/.lintr +++ b/.lintr @@ -6,6 +6,7 @@ linters: linters_with_defaults( object_length_linter = NULL, indentation_linter( indent = 2L, + hanging_indent_style = "never", assignment_as_infix = FALSE ) ) diff --git a/R/Log_Likelihood.R b/R/Log_Likelihood.R index a7260c3..6b9fce0 100755 --- a/R/Log_Likelihood.R +++ b/R/Log_Likelihood.R @@ -250,15 +250,12 @@ create.cov.upper.flex <- function(P, marg.var, marg.range, marg.smooth, j <- combs[iter, 2] smoothness.mat[i, j] <- (marg.smooth[i] + marg.smooth[j]) / 2 - range.mat[i, j] <- 1 / sqrt(((1 / marg.range[i])^2 + - (1 / marg.range[j])^2) / 2) + range.mat[i, j] <- 1 / sqrt(((1 / marg.range[i])^2 + (1 / marg.range[j])^2) / 2) s1 <- sqrt(marg.var[i] * marg.var[j]) - s2 <- ((1 / marg.range[i])^marg.smooth[i] * - (1 / marg.range[j])^marg.smooth[j]) / + s2 <- ((1 / marg.range[i])^marg.smooth[i] * (1 / marg.range[j])^marg.smooth[j]) / ((1 / range.mat[i, j])^(2 * smoothness.mat[i, j])) - s3 <- gamma(smoothness.mat[i, j]) / (sqrt(gamma(marg.smooth[i])) * - sqrt(gamma(marg.smooth[j]))) + s3 <- gamma(smoothness.mat[i, j]) / (sqrt(gamma(marg.smooth[i])) * sqrt(gamma(marg.smooth[j]))) s4 <- R.corr[iter] sig2.mat[i, j] <- s1 * s2 * s3 * s4 } diff --git a/R/PrestoGP_CreateU_Multivariate.R b/R/PrestoGP_CreateU_Multivariate.R index 1ec3e7b..c486541 100755 --- a/R/PrestoGP_CreateU_Multivariate.R +++ b/R/PrestoGP_CreateU_Multivariate.R @@ -112,8 +112,7 @@ sparseNN <- function(ordered_locs, n_neighbors, for (row in 1:n_neighbors) { # for the locations from 1 to n_neighbors, use the entire locs list to find the neighbors nn <- knn_indices( - ordered_locs[1: - (n_neighbors + 1), , drop = FALSE][-row, , + ordered_locs[1:(n_neighbors + 1), , drop = FALSE][-row, , drop = FALSE ], ordered_locs[row, , drop = FALSE], n_neighbors, @@ -365,8 +364,7 @@ createUMultivariate <- function(vec.approx, params, cov_func = NULL) { ajj <- 1 / rangep[ondx[2]] aij <- sqrt((aii^2 + ajj^2) / 2) K1 <- rho.mat[ondx[1], ondx[2]] * sqrt(sig2[ondx[1]]) * sqrt(sig2[ondx[2]]) * - aii^vii * ajj^vjj * gamma(vij) / (aij^(2 * vij) * sqrt(gamma(vii) * - gamma(vjj))) * + aii^vii * ajj^vjj * gamma(vij) / (aij^(2 * vij) * sqrt(gamma(vii) * gamma(vjj))) * cov_func(dist_func(olocs[1, , drop = FALSE], olocs[2, , drop = FALSE], ), smoothness = vij, alpha = aij ) @@ -421,8 +419,7 @@ createUMultivariate <- function(vec.approx, params, cov_func = NULL) { ajj <- 1 / rangep[ondx[2]] aij <- sqrt((aii^2 + ajj^2) / 2) K1 <- rho.mat[ondx[1], ondx[2]] * sqrt(sig2[ondx[1]]) * sqrt(sig2[ondx[2]]) * - aii^vii * ajj^vjj * gamma(vij) / (aij^(2 * vij) * sqrt(gamma(vii) * - gamma(vjj))) * + aii^vii * ajj^vjj * gamma(vij) / (aij^(2 * vij) * sqrt(gamma(vii) * gamma(vjj))) * cov_func(dist_func(olocs[1, , drop = FALSE], olocs[2, , drop = FALSE], ), smoothness = vij, alpha = aij ) @@ -461,8 +458,7 @@ createUMultivariate <- function(vec.approx, params, cov_func = NULL) { # positive definite. See equation (9) in Apanasovich (2011). K1[j] <- rho.mat[ondx[i], ondx[cur.q[j]]] * sqrt(sig2[ondx[i]]) * sqrt(sig2[ondx[cur.q[j]]]) * - aii^vii * ajj^vjj * gamma(vij) / (aij^(2 * vij) * - sqrt(gamma(vii) * gamma(vjj))) * + aii^vii * ajj^vjj * gamma(vij) / (aij^(2 * vij) * sqrt(gamma(vii) * gamma(vjj))) * cov_func( dist_func( olocs[i, , drop = FALSE], diff --git a/R/PrestoGP_Model.R b/R/PrestoGP_Model.R index 8f2e865..dab605f 100755 --- a/R/PrestoGP_Model.R +++ b/R/PrestoGP_Model.R @@ -521,8 +521,7 @@ setMethod("scale_locs", "PrestoGPModel", function(model, locs) { for (j in 1:model@nscale) { locs.out[[i]][, model@scaling == j] <- locs[[i]][, model@scaling == j] / - model@covparams[model@param_sequence[2, 1] + - model@nscale * (i - 1) + j - 1] + model@covparams[model@param_sequence[2, 1] + model@nscale * (i - 1) + j - 1] } } return(locs.out) @@ -535,25 +534,20 @@ setMethod("transform_covariance_parameters", "PrestoGPModel", function(model) { model@covparams <- c( exp(model@logparams[1:model@param_sequence[2, 2]]), gtools::inv.logit( - model@logparams[model@param_sequence[3, 1]: - model@param_sequence[3, 2]], + model@logparams[model@param_sequence[3, 1]:model@param_sequence[3, 2]], 0, 2.5 ), - exp(model@logparams[model@param_sequence[4, 1]: - model@param_sequence[4, 2]]), - tanh(model@logparams[model@param_sequence[5, 1]: - model@param_sequence[5, 2]]) + exp(model@logparams[model@param_sequence[4, 1]:model@param_sequence[4, 2]]), + tanh(model@logparams[model@param_sequence[5, 1]:model@param_sequence[5, 2]]) ) } else { model@covparams <- c( exp(model@logparams[1:model@param_sequence[2, 2]]), gtools::inv.logit( - model@logparams[model@param_sequence[3, 1]: - model@param_sequence[3, 2]], + model@logparams[model@param_sequence[3, 1]:model@param_sequence[3, 2]], 0, 2.5 ), - exp(model@logparams[model@param_sequence[4, 1]: - model@param_sequence[4, 2]]), 1 + exp(model@logparams[model@param_sequence[4, 1]:model@param_sequence[4, 2]]), 1 ) } invisible(model) diff --git a/R/PrestoGP_Multivariate_Vecchia.R b/R/PrestoGP_Multivariate_Vecchia.R index 388ee07..37823f4 100755 --- a/R/PrestoGP_Multivariate_Vecchia.R +++ b/R/PrestoGP_Multivariate_Vecchia.R @@ -45,78 +45,76 @@ setMethod("initialize", "MultivariateVecchiaModel", function(.Object, n_neighbor #' prediction <- prestogp_predict(model, X.test, locs.test) #' Vec.mean <- prediction[[1]] #' Vec.sds <- prediction[[2]] -setMethod("prestogp_predict", - "MultivariateVecchiaModel", +setMethod("prestogp_predict", "MultivariateVecchiaModel", function(model, X, locs, m = NULL, ordering.pred = c("obspred", "general"), pred.cond = c("independent", "general"), return.values = c("mean", "meanvar")) { - # validate parameters - ordering.pred <- match.arg(ordering.pred) - pred.cond <- match.arg(pred.cond) - return.values <- match.arg(return.values) - X <- check_input_pred(model, X, locs) - if (is.null(m)) { # m defaults to the value used for training - m <- model@n_neighbors - } - if (m < model@min_m) { - stop(paste("m must be at least ", model@min_m, sep = "")) - } - if (m >= nrow(model@X_train)) { - warning("Conditioning set size m chosen to be >=n. Changing to m=n-1") - m <- nrow(model@X_train) - 1 - } - - # Vecchia prediction at new locations - Vecchia.Pred <- predict(model@linear_model, newx = X, s = model@linear_model$lambda[model@lambda_1se_idx]) - # Vecchia trend prediction at observed data - Vecchia.hat <- predict(model@linear_model, newx = model@X_train, s = model@linear_model$lambda[model@lambda_1se_idx]) + # validate parameters + ordering.pred <- match.arg(ordering.pred) + pred.cond <- match.arg(pred.cond) + return.values <- match.arg(return.values) + X <- check_input_pred(model, X, locs) + if (is.null(m)) { # m defaults to the value used for training + m <- model@n_neighbors + } + if (m < model@min_m) { + stop(paste("m must be at least ", model@min_m, sep = "")) + } + if (m >= nrow(model@X_train)) { + warning("Conditioning set size m chosen to be >=n. Changing to m=n-1") + m <- nrow(model@X_train) - 1 + } - # Test set prediction - res <- model@Y_train - Vecchia.hat + # Vecchia prediction at new locations + Vecchia.Pred <- predict(model@linear_model, newx = X, s = model@linear_model$lambda[model@lambda_1se_idx]) + # Vecchia trend prediction at observed data + Vecchia.hat <- predict(model@linear_model, newx = model@X_train, s = model@linear_model$lambda[model@lambda_1se_idx]) - locs.train.scaled <- scale_locs(model, model@locs_train) - locs.scaled <- scale_locs(model, locs) - vec.approx.test <- vecchia_Mspecify(locs.train.scaled, m, - locs.list.pred = locs.scaled, - ordering.pred = ordering.pred, - pred.cond = pred.cond - ) + # Test set prediction + res <- model@Y_train - Vecchia.hat - ## carry out prediction - if (!model@apanasovich) { - params <- model@covparams - param.seq <- model@param_sequence - pred <- vecchia_Mprediction(res, vec.approx.test, - c( - params[1:param.seq[1, 2]], - rep(1, param.seq[2, 2] - param.seq[2, 1] + 1), - params[param.seq[3, 1]: - param.seq[5, 2]] - ), - return.values = return.values - ) - } else { - pred <- vecchia_Mprediction(res, vec.approx.test, model@covparams, - return.values = return.values + locs.train.scaled <- scale_locs(model, model@locs_train) + locs.scaled <- scale_locs(model, locs) + vec.approx.test <- vecchia_Mspecify(locs.train.scaled, m, + locs.list.pred = locs.scaled, + ordering.pred = ordering.pred, + pred.cond = pred.cond ) - } - # prediction function can return both mean and sds - # returns a list with elements mu.pred,mu.obs,var.pred,var.obs,V.ord - Vec.mean <- pred$mu.pred + Vecchia.Pred # residual + mean trend - if (return.values == "mean") { - return.list <- list(means = Vec.mean) - } else { - warning("Variance estimates do not include model fitting variance and are anticonservative. Use with caution.") - vec.var <- pred$var.pred - ndx.out <- NULL - for (i in seq_along(length(locs))) { - vec.sds[ndx.out == i] <- sqrt(vec.var[ndx.out == i] + - model@covparams[model@param_sequence[4, i]]) + ## carry out prediction + if (!model@apanasovich) { + params <- model@covparams + param.seq <- model@param_sequence + pred <- vecchia_Mprediction(res, vec.approx.test, + c( + params[1:param.seq[1, 2]], + rep(1, param.seq[2, 2] - param.seq[2, 1] + 1), + params[param.seq[3, 1]:param.seq[5, 2]] + ), + return.values = return.values + ) + } else { + pred <- vecchia_Mprediction(res, vec.approx.test, model@covparams, + return.values = return.values + ) } - return.list <- list(means = Vec.mean, sds = vec.sds) - } - return(return.list) -}) + # prediction function can return both mean and sds + # returns a list with elements mu.pred,mu.obs,var.pred,var.obs,V.ord + Vec.mean <- pred$mu.pred + Vecchia.Pred # residual + mean trend + if (return.values == "mean") { + return.list <- list(means = Vec.mean) + } else { + warning("Variance estimates do not include model fitting variance and are anticonservative. Use with caution.") + vec.var <- pred$var.pred + ndx.out <- NULL + for (i in seq_along(length(locs))) { + vec.sds[ndx.out == i] <- sqrt(vec.var[ndx.out == i] + model@covparams[model@param_sequence[4, i]]) + } + return.list <- list(means = Vec.mean, sds = vec.sds) + } + + return(return.list) + } +) setMethod("check_input", "MultivariateVecchiaModel", function(model, Y, X, locs) { if (!is.list(locs)) { @@ -222,8 +220,7 @@ setMethod("specify", "MultivariateVecchiaModel", function(model) { for (j in 1:model@nscale) { olocs.scaled[model@vecchia_approx$ondx == i, model@scaling == j] <- olocs.scaled[model@vecchia_approx$ondx == i, model@scaling == j] * - model@covparams[model@param_sequence[2, 1] + - model@nscale * (i - 1) + j - 1] + model@covparams[model@param_sequence[2, 1] + model@nscale * (i - 1) + j - 1] } } model@vecchia_approx$locsord <- olocs.scaled @@ -294,12 +291,8 @@ setMethod("transform_data", "MultivariateVecchiaModel", function(model, Y, X) { vecchia.approx = vecchia.approx, c( params[1:param.seq[1, 2]], - rep( - 1, - param.seq[2, 2] - param.seq[2, 1] + 1 - ), - params[param.seq[3, 1]: - param.seq[5, 2]] + rep(1, param.seq[2, 2] - param.seq[2, 1] + 1), + params[param.seq[3, 1]:param.seq[5, 2]] ) ) } else { diff --git a/R/PrestoGP_Util_Functions.R b/R/PrestoGP_Util_Functions.R index 26590e2..3213bb0 100644 --- a/R/PrestoGP_Util_Functions.R +++ b/R/PrestoGP_Util_Functions.R @@ -232,8 +232,7 @@ revMat <- function(mat) { } #' @export -vecchia_Mprediction <- function(z, vecchia.approx, covparms, var.exact = NULL, - return.values = "mean") { +vecchia_Mprediction <- function(z, vecchia.approx, covparms, var.exact = NULL, return.values = "mean") { GPvecchia:::removeNAs() U.obj <- createUMultivariate(vecchia.approx, covparms) V.ord <- U2V(U.obj) diff --git a/R/PrestoGP_Vecchia.R b/R/PrestoGP_Vecchia.R index 5c76dca..eda0117 100755 --- a/R/PrestoGP_Vecchia.R +++ b/R/PrestoGP_Vecchia.R @@ -42,71 +42,71 @@ setMethod("initialize", "VecchiaModel", function(.Object, n_neighbors = 25, ...) #' prediction <- prestogp_predict(model, X.test, locs.test) #' Vec.mean <- prediction[[1]] #' Vec.sds <- prediction[[2]] -setMethod("prestogp_predict", - "VecchiaModel", +setMethod("prestogp_predict", "VecchiaModel", function(model, X, locs, m = NULL, ordering.pred = c("obspred", "general"), pred.cond = c("independent", "general"), return.values = c("mean", "meanvar")) { - # validate parameters - ordering.pred <- match.arg(ordering.pred) - pred.cond <- match.arg(pred.cond) - return.values <- match.arg(return.values) - model <- check_input_pred(model, X, locs) - if (is.null(m)) { # m defaults to the value used for training - m <- model@n_neighbors - } - if (m < model@min_m) { - stop(paste("m must be at least ", model@min_m, sep = "")) - } - if (m >= nrow(model@X_train)) { - warning("Conditioning set size m chosen to be >=n. Changing to m=n-1") - m <- nrow(model@X_train) - 1 - } + # validate parameters + ordering.pred <- match.arg(ordering.pred) + pred.cond <- match.arg(pred.cond) + return.values <- match.arg(return.values) + model <- check_input_pred(model, X, locs) + if (is.null(m)) { # m defaults to the value used for training + m <- model@n_neighbors + } + if (m < model@min_m) { + stop(paste("m must be at least ", model@min_m, sep = "")) + } + if (m >= nrow(model@X_train)) { + warning("Conditioning set size m chosen to be >=n. Changing to m=n-1") + m <- nrow(model@X_train) - 1 + } - # Vecchia prediction at new locations - # Vecchia.Pred <- predict(model@Vecchia_SCAD_fit[[1]], X = X, which = model@lambda_1se_idx[[1]]) - Vecchia.Pred <- predict(model@linear_model, newx = X, s = model@linear_model$lambda[model@lambda_1se_idx]) - # Vecchia trend prediction at observed data - # Vecchia.hat <- predict(model@Vecchia_SCAD_fit[[1]], X = model@X_train, which = model@lambda_1se_idx[[1]]) - Vecchia.hat <- predict(model@linear_model, newx = model@X_train, s = model@linear_model$lambda[model@lambda_1se_idx]) + # Vecchia prediction at new locations + # Vecchia.Pred <- predict(model@Vecchia_SCAD_fit[[1]], X = X, which = model@lambda_1se_idx[[1]]) + Vecchia.Pred <- predict(model@linear_model, newx = X, s = model@linear_model$lambda[model@lambda_1se_idx]) + # Vecchia trend prediction at observed data + # Vecchia.hat <- predict(model@Vecchia_SCAD_fit[[1]], X = model@X_train, which = model@lambda_1se_idx[[1]]) + Vecchia.hat <- predict(model@linear_model, newx = model@X_train, s = model@linear_model$lambda[model@lambda_1se_idx]) - # Test set prediction - res <- model@Y_train - Vecchia.hat + # Test set prediction + res <- model@Y_train - Vecchia.hat - locs.train.scaled <- scale_locs(model, model@locs_train)[[1]] - locs.scaled <- scale_locs(model, list(locs))[[1]] - vec.approx.test <- vecchia_specify(locs.train.scaled, m, locs.pred = locs.scaled, ordering.pred = ordering.pred, pred.cond = pred.cond) + locs.train.scaled <- scale_locs(model, model@locs_train)[[1]] + locs.scaled <- scale_locs(model, list(locs))[[1]] + vec.approx.test <- vecchia_specify(locs.train.scaled, m, locs.pred = locs.scaled, ordering.pred = ordering.pred, pred.cond = pred.cond) - ## carry out prediction - if (!model@apanasovich) { - pred <- vecchia_prediction( - res, - vec.approx.test, - c(model@covparams[1], 1, model@covparams[3]), - model@covparams[4], - return.values = return.values - ) - } else { - pred <- vecchia_prediction( - res, - vec.approx.test, - c(model@covparams[1], model@covparams[2], model@covparams[3]), - model@covparams[4], return.values = return.values - ) - } + ## carry out prediction + if (!model@apanasovich) { + pred <- vecchia_prediction( + res, + vec.approx.test, + c(model@covparams[1], 1, model@covparams[3]), + model@covparams[4], + return.values = return.values + ) + } else { + pred <- vecchia_prediction( + res, + vec.approx.test, + c(model@covparams[1], model@covparams[2], model@covparams[3]), + model@covparams[4], return.values = return.values + ) + } - # prediction function can return both mean and sds - # returns a list with elements mu.pred,mu.obs,var.pred,var.obs,V.ord - Vec.mean <- pred$mu.pred + Vecchia.Pred # residual + mean trend - if (return.values == "mean") { - return.list <- list(means = Vec.mean) - } else { - warning("Variance estimates do not include model fitting variance and are anticonservative. Use with caution.") - # TODO: @Eric.Bair is this a typo/bug? Capital 'V' in Vec.sds but 'vec.sds' in return.list - Vec.sds <- sqrt(pred$var.pred + model@covparams[4]) - return.list <- list(means = Vec.mean, sds = vec.sds) - } + # prediction function can return both mean and sds + # returns a list with elements mu.pred,mu.obs,var.pred,var.obs,V.ord + Vec.mean <- pred$mu.pred + Vecchia.Pred # residual + mean trend + if (return.values == "mean") { + return.list <- list(means = Vec.mean) + } else { + warning("Variance estimates do not include model fitting variance and are anticonservative. Use with caution.") + # TODO: @Eric.Bair is this a typo/bug? Capital 'V' in Vec.sds but 'vec.sds' in return.list + Vec.sds <- sqrt(pred$var.pred + model@covparams[4]) + return.list <- list(means = Vec.mean, sds = vec.sds) + } - return(return.list) -}) + return(return.list) + } +) setMethod("check_input", "VecchiaModel", function(model, Y, X, locs) { if (!is.matrix(locs)) { diff --git a/R/RcppExports.R b/R/RcppExports.R index 0e90b6b..293cbae 100644 --- a/R/RcppExports.R +++ b/R/RcppExports.R @@ -2,9 +2,9 @@ # Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393 na_omit_c <- function(x) { - .Call("_PrestoGP_na_omit_c", PACKAGE = "PrestoGP", x) + .Call("_PrestoGP_na_omit_c", PACKAGE = "PrestoGP", x) } createU_helper_mat <- function(olocs, ondx, curqys, curqzs, vijs, aijs, full_const, nugget, sig2, U_beginning) { - .Call("_PrestoGP_createU_helper_mat", PACKAGE = "PrestoGP", olocs, ondx, curqys, curqzs, vijs, aijs, full_const, nugget, sig2, U_beginning) + .Call("_PrestoGP_createU_helper_mat", PACKAGE = "PrestoGP", olocs, ondx, curqys, curqzs, vijs, aijs, full_const, nugget, sig2, U_beginning) } diff --git a/tests/testthat/sim_multivariate.R b/tests/testthat/sim_multivariate.R index f1bf46f..6763784 100644 --- a/tests/testthat/sim_multivariate.R +++ b/tests/testthat/sim_multivariate.R @@ -102,20 +102,20 @@ d22 <- fields::rdist(expand.grid(space.dim2, space.dim2, time.dim2)) d33 <- fields::rdist(expand.grid(space.dim3, space.dim3, time.dim3)) d12 <- fields::rdist( - expand.grid(space.dim1, space.dim1, time.dim1), - expand.grid(space.dim2, space.dim2, time.dim2) + expand.grid(space.dim1, space.dim1, time.dim1), + expand.grid(space.dim2, space.dim2, time.dim2) ) d21 <- t(d12) d13 <- fields::rdist( - expand.grid(space.dim1, space.dim1, time.dim1), - expand.grid(space.dim3, space.dim3, time.dim3) + expand.grid(space.dim1, space.dim1, time.dim1), + expand.grid(space.dim3, space.dim3, time.dim3) ) d31 <- t(d13) d23 <- fields::rdist( - expand.grid(space.dim2, space.dim2, time.dim2), - expand.grid(space.dim3, space.dim3, time.dim3) + expand.grid(space.dim2, space.dim2, time.dim2), + expand.grid(space.dim3, space.dim3, time.dim3) ) d32 <- t(d23) @@ -132,8 +132,8 @@ aii <- 1 / ranges[1] ajj <- 1 / ranges[2] aij <- sqrt((aii^2 + ajj^2) / 2) Sigma12 <- rho[1, 2] * sqrt(marg.var[1]) * sqrt(marg.var[2]) * aii^vii * - ajj^vjj * gamma(vij) / (aij^(2 * vij) * sqrt(gamma(vii) * gamma(vjj))) * - fields::Matern(d12, smoothness = vij, alpha = aij) + ajj^vjj * gamma(vij) / (aij^(2 * vij) * sqrt(gamma(vii) * gamma(vjj))) * + fields::Matern(d12, smoothness = vij, alpha = aij) Sigma21 <- t(Sigma12) vii <- marg.smoothness[1] @@ -143,8 +143,8 @@ aii <- 1 / ranges[1] ajj <- 1 / ranges[3] aij <- sqrt((aii^2 + ajj^2) / 2) Sigma13 <- rho[1, 3] * sqrt(marg.var[1]) * sqrt(marg.var[3]) * aii^vii * - ajj^vjj * gamma(vij) / (aij^(2 * vij) * sqrt(gamma(vii) * gamma(vjj))) * - fields::Matern(d13, smoothness = vij, alpha = aij) + ajj^vjj * gamma(vij) / (aij^(2 * vij) * sqrt(gamma(vii) * gamma(vjj))) * + fields::Matern(d13, smoothness = vij, alpha = aij) Sigma31 <- t(Sigma13) vii <- marg.smoothness[2] @@ -154,15 +154,15 @@ aii <- 1 / ranges[2] ajj <- 1 / ranges[3] aij <- sqrt((aii^2 + ajj^2) / 2) Sigma23 <- rho[2, 3] * sqrt(marg.var[2]) * sqrt(marg.var[3]) * aii^vii * - ajj^vjj * gamma(vij) / (aij^(2 * vij) * sqrt(gamma(vii) * gamma(vjj))) * - fields::Matern(d23, smoothness = vij, alpha = aij) + ajj^vjj * gamma(vij) / (aij^(2 * vij) * sqrt(gamma(vii) * gamma(vjj))) * + fields::Matern(d23, smoothness = vij, alpha = aij) Sigma32 <- t(Sigma23) # Combine into the super Multivariate covariance matrix Sigma.All <- rbind( - cbind(Sigma11, Sigma12, Sigma13), - cbind(Sigma21, Sigma22, Sigma23), - cbind(Sigma31, Sigma32, Sigma33) + cbind(Sigma11, Sigma12, Sigma13), + cbind(Sigma21, Sigma22, Sigma23), + cbind(Sigma31, Sigma32, Sigma33) ) # Cholesky decomposition diff --git a/tests/testthat/sim_multivariate_big.R b/tests/testthat/sim_multivariate_big.R index ce9b9da..aed9017 100755 --- a/tests/testthat/sim_multivariate_big.R +++ b/tests/testthat/sim_multivariate_big.R @@ -12,9 +12,9 @@ beta.all <- rep(beta1, ny) X.st <- list() for (i in 1:ny) { - set.seed(1212 + i) - Sigma.X <- exp(-rdist(sample(1:p)) / 3) - X.st[[i]] <- mvrnorm(n.spatial.xy^2, rep(0, p), Sigma.X) + set.seed(1212 + i) + Sigma.X <- exp(-rdist(sample(1:p)) / 3) + X.st[[i]] <- mvrnorm(n.spatial.xy^2, rep(0, p), Sigma.X) } X.all <- superMatrix(X.st) mean.trend.st <- X.all %*% beta.all @@ -35,39 +35,33 @@ params.all <- c(x.variance, ranges, marg.smoothness, nuggets, rho.vec) locs.list <- list() for (i in 1:ny) { - loc1 <- seq(0, 1, length.out = n.spatial.xy) + rnorm(n.spatial.xy, 0, 0.001) - loc2 <- seq(0, 1, length.out = n.spatial.xy) + rnorm(n.spatial.xy, 0, 0.001) - locs.list[[i]] <- as.matrix(expand.grid(loc1, loc2)) + loc1 <- seq(0, 1, length.out = n.spatial.xy) + rnorm(n.spatial.xy, 0, 0.001) + loc2 <- seq(0, 1, length.out = n.spatial.xy) + rnorm(n.spatial.xy, 0, 0.001) + locs.list[[i]] <- as.matrix(expand.grid(loc1, loc2)) } Sigma.All <- matrix(nrow = ny * n.spatial.xy^2, ncol = ny * n.spatial.xy^2) for (i in 1:ny) { - for (j in i:ny) { - if (i == j) { - ndx1 <- ((i - 1) * n.spatial.xy^2 + 1):(n.spatial.xy^2 * i) - dij <- fields::rdist(locs.list[[i]]) - Sigma.All[ndx1, ndx1] <- marg.var[i] * - fields::Matern(dij, - range = ranges[i], - smoothness = marg.smoothness[i] - ) - } else { - ndx1 <- ((i - 1) * n.spatial.xy^2 + 1):(n.spatial.xy^2 * i) - ndx2 <- ((j - 1) * n.spatial.xy^2 + 1):(n.spatial.xy^2 * j) - dij <- fields::rdist(locs.list[[i]], locs.list[[j]]) - vii <- marg.smoothness[i] - vjj <- marg.smoothness[j] - vij <- (vii + vjj) / 2 - aii <- 1 / ranges[i] - ajj <- 1 / ranges[j] - aij <- sqrt((aii^2 + ajj^2) / 2) - Sigma.All[ndx1, ndx2] <- rho[i, j] * sqrt(marg.var[i]) * - sqrt(marg.var[j]) * aii^vii * ajj^vjj * gamma(vij) / - (aij^(2 * vij) * sqrt(gamma(vii) * gamma(vjj))) * - fields::Matern(dij, smoothness = vij, alpha = aij) - Sigma.All[ndx2, ndx1] <- t(Sigma.All[ndx1, ndx2]) - } + for (j in i:ny) { + if (i == j) { + ndx1 <- ((i - 1) * n.spatial.xy^2 + 1):(n.spatial.xy^2 * i) + dij <- fields::rdist(locs.list[[i]]) + Sigma.All[ndx1, ndx1] <- marg.var[i] * fields::Matern(dij, range = ranges[i], smoothness = marg.smoothness[i]) + } else { + ndx1 <- ((i - 1) * n.spatial.xy^2 + 1):(n.spatial.xy^2 * i) + ndx2 <- ((j - 1) * n.spatial.xy^2 + 1):(n.spatial.xy^2 * j) + dij <- fields::rdist(locs.list[[i]], locs.list[[j]]) + vii <- marg.smoothness[i] + vjj <- marg.smoothness[j] + vij <- (vii + vjj) / 2 + aii <- 1 / ranges[i] + ajj <- 1 / ranges[j] + aij <- sqrt((aii^2 + ajj^2) / 2) + Sigma.All[ndx1, ndx2] <- rho[i, j] * sqrt(marg.var[i]) * sqrt(marg.var[j]) * aii^vii * ajj^vjj * gamma(vij) / + (aij^(2 * vij) * sqrt(gamma(vii) * gamma(vjj))) * fields::Matern(dij, smoothness = vij, alpha = aij) + Sigma.All[ndx2, ndx1] <- t(Sigma.All[ndx1, ndx2]) } + } } L.C <- chol(Sigma.All) @@ -76,18 +70,18 @@ st.error <- rnorm(n.spatial.xy^2 * ny) %*% L.C nug.error <- NULL for (i in 1:ny) { - nug.error <- c(nug.error, nuggets[i] * rnorm(n.spatial.xy^2)) + nug.error <- c(nug.error, nuggets[i] * rnorm(n.spatial.xy^2)) } y.list <- list() for (i in 1:ny) { - ndx1 <- ((i - 1) * n.spatial.xy^2 + 1):(n.spatial.xy^2 * i) - y.list[[i]] <- mean.trend.st[ndx1] + st.error[ndx1] + nug.error[ndx1] + ndx1 <- ((i - 1) * n.spatial.xy^2 + 1):(n.spatial.xy^2 * i) + y.list[[i]] <- mean.trend.st[ndx1] + st.error[ndx1] + nug.error[ndx1] } rm( - ny, p, p.nz, n.spatial.xy, beta1, i, j, Sigma.X, mean.trend.st, n.rho, - loc1, loc2, ndx1, ndx2, dij, vii, vjj, vij, aii, ajj, aij, L.C, - st.error, nug.error, X.all, rho, rho.vec, ranges, Sigma.All, nuggets, - marg.smoothness, marg.var + ny, p, p.nz, n.spatial.xy, beta1, i, j, Sigma.X, mean.trend.st, n.rho, + loc1, loc2, ndx1, ndx2, dij, vii, vjj, vij, aii, ajj, aij, L.C, + st.error, nug.error, X.all, rho, rho.vec, ranges, Sigma.All, nuggets, + marg.smoothness, marg.var ) diff --git a/tests/testthat/sim_multivariate_big_pred.R b/tests/testthat/sim_multivariate_big_pred.R index 6b6a972..2df4cb6 100755 --- a/tests/testthat/sim_multivariate_big_pred.R +++ b/tests/testthat/sim_multivariate_big_pred.R @@ -14,8 +14,8 @@ beta.all <- rep(beta1, ny) X.st <- list() for (i in 1:ny) { - Sigma.X <- exp(-rdist(sample(1:p)) / 3) - X.st[[i]] <- mvrnorm(n.spatial.xy^2, rep(0, p), Sigma.X) + Sigma.X <- exp(-rdist(sample(1:p)) / 3) + X.st[[i]] <- mvrnorm(n.spatial.xy^2, rep(0, p), Sigma.X) } X.all <- superMatrix(X.st) mean.trend.st <- X.all %*% beta.all @@ -36,39 +36,37 @@ params.all <- c(x.variance, ranges, marg.smoothness, nuggets, rho.vec) locs.list <- list() for (i in 1:ny) { - loc1 <- seq(0, 1, length.out = n.spatial.xy) + rnorm(n.spatial.xy, 0, 0.001) - loc2 <- seq(0, 1, length.out = n.spatial.xy) + rnorm(n.spatial.xy, 0, 0.001) - locs.list[[i]] <- as.matrix(expand.grid(loc1, loc2)) + loc1 <- seq(0, 1, length.out = n.spatial.xy) + rnorm(n.spatial.xy, 0, 0.001) + loc2 <- seq(0, 1, length.out = n.spatial.xy) + rnorm(n.spatial.xy, 0, 0.001) + locs.list[[i]] <- as.matrix(expand.grid(loc1, loc2)) } Sigma.All <- matrix(nrow = ny * n.spatial.xy^2, ncol = ny * n.spatial.xy^2) for (i in 1:ny) { - for (j in i:ny) { - if (i == j) { - ndx1 <- ((i - 1) * n.spatial.xy^2 + 1):(n.spatial.xy^2 * i) - dij <- fields::rdist(locs.list[[i]]) - Sigma.All[ndx1, ndx1] <- marg.var[i] * - fields::Matern(dij, - range = ranges[i], - smoothness = marg.smoothness[i] - ) - } else { - ndx1 <- ((i - 1) * n.spatial.xy^2 + 1):(n.spatial.xy^2 * i) - ndx2 <- ((j - 1) * n.spatial.xy^2 + 1):(n.spatial.xy^2 * j) - dij <- fields::rdist(locs.list[[i]], locs.list[[j]]) - vii <- marg.smoothness[i] - vjj <- marg.smoothness[j] - vij <- (vii + vjj) / 2 - aii <- 1 / ranges[i] - ajj <- 1 / ranges[j] - aij <- sqrt((aii^2 + ajj^2) / 2) - Sigma.All[ndx1, ndx2] <- rho[i, j] * sqrt(marg.var[i]) * - sqrt(marg.var[j]) * aii^vii * ajj^vjj * gamma(vij) / - (aij^(2 * vij) * sqrt(gamma(vii) * gamma(vjj))) * - fields::Matern(dij, smoothness = vij, alpha = aij) - Sigma.All[ndx2, ndx1] <- t(Sigma.All[ndx1, ndx2]) - } + for (j in i:ny) { + if (i == j) { + ndx1 <- ((i - 1) * n.spatial.xy^2 + 1):(n.spatial.xy^2 * i) + dij <- fields::rdist(locs.list[[i]]) + Sigma.All[ndx1, ndx1] <- marg.var[i] * fields::Matern(dij, + range = ranges[i], + smoothness = marg.smoothness[i] + ) + } else { + ndx1 <- ((i - 1) * n.spatial.xy^2 + 1):(n.spatial.xy^2 * i) + ndx2 <- ((j - 1) * n.spatial.xy^2 + 1):(n.spatial.xy^2 * j) + dij <- fields::rdist(locs.list[[i]], locs.list[[j]]) + vii <- marg.smoothness[i] + vjj <- marg.smoothness[j] + vij <- (vii + vjj) / 2 + aii <- 1 / ranges[i] + ajj <- 1 / ranges[j] + aij <- sqrt((aii^2 + ajj^2) / 2) + Sigma.All[ndx1, ndx2] <- rho[i, j] * sqrt(marg.var[i]) * sqrt(marg.var[j]) * aii^vii * ajj^vjj * gamma(vij) / + (aij^(2 * vij) * sqrt(gamma(vii) * gamma(vjj))) * + fields::Matern(dij, smoothness = vij, alpha = aij) + Sigma.All[ndx2, ndx1] <- t(Sigma.All[ndx1, ndx2]) } + } } L.C <- chol(Sigma.All) @@ -77,7 +75,7 @@ st.error <- rnorm(n.spatial.xy^2 * ny) %*% L.C nug.error <- NULL for (i in 1:ny) { - nug.error <- c(nug.error, nuggets[i] * rnorm(n.spatial.xy^2)) + nug.error <- c(nug.error, nuggets[i] * rnorm(n.spatial.xy^2)) } y.list <- list() @@ -88,22 +86,22 @@ y.list.otst <- list() X.st.otr <- list() X.st.otst <- list() for (i in 1:ny) { - ndx1 <- ((i - 1) * n.spatial.xy^2 + 1):(n.spatial.xy^2 * i) - y.list[[i]] <- mean.trend.st[ndx1] + st.error[ndx1] + nug.error[ndx1] - nn <- n.spatial.xy^2 - otr <- rep(FALSE, nn) - otr[sample(1:nn, size = floor(nn / 2))] <- TRUE - locs.list.otr[[i]] <- locs.list[[i]][otr, ] - locs.list.otst[[i]] <- locs.list[[i]][!otr, ] - y.list.otr[[i]] <- y.list[[i]][otr] - y.list.otst[[i]] <- y.list[[i]][!otr] - X.st.otr[[i]] <- X.st[[i]][otr, ] - X.st.otst[[i]] <- X.st[[i]][!otr, ] + ndx1 <- ((i - 1) * n.spatial.xy^2 + 1):(n.spatial.xy^2 * i) + y.list[[i]] <- mean.trend.st[ndx1] + st.error[ndx1] + nug.error[ndx1] + nn <- n.spatial.xy^2 + otr <- rep(FALSE, nn) + otr[sample(1:nn, size = floor(nn / 2))] <- TRUE + locs.list.otr[[i]] <- locs.list[[i]][otr, ] + locs.list.otst[[i]] <- locs.list[[i]][!otr, ] + y.list.otr[[i]] <- y.list[[i]][otr] + y.list.otst[[i]] <- y.list[[i]][!otr] + X.st.otr[[i]] <- X.st[[i]][otr, ] + X.st.otst[[i]] <- X.st[[i]][!otr, ] } rm( - ny, p, p.nz, n.spatial.xy, beta1, i, j, Sigma.X, mean.trend.st, n.rho, - loc1, loc2, ndx1, ndx2, dij, vii, vjj, vij, aii, ajj, aij, L.C, - st.error, nug.error, X.all, rho, rho.vec, ranges, Sigma.All, nuggets, - marg.smoothness, marg.var, x.variance, nn, otr + ny, p, p.nz, n.spatial.xy, beta1, i, j, Sigma.X, mean.trend.st, n.rho, + loc1, loc2, ndx1, ndx2, dij, vii, vjj, vij, aii, ajj, aij, L.C, + st.error, nug.error, X.all, rho, rho.vec, ranges, Sigma.All, nuggets, + marg.smoothness, marg.var, x.variance, nn, otr ) diff --git a/tests/testthat/sim_multivariate_big_st.R b/tests/testthat/sim_multivariate_big_st.R index 2408cb8..b0adc6c 100755 --- a/tests/testthat/sim_multivariate_big_st.R +++ b/tests/testthat/sim_multivariate_big_st.R @@ -14,8 +14,8 @@ beta.all <- rep(beta1, ny) X.st <- list() for (i in 1:ny) { - Sigma.X <- exp(-rdist(sample(1:p)) / 3) - X.st[[i]] <- mvrnorm(n.spatial.xy^3, rep(0, p), Sigma.X) + Sigma.X <- exp(-rdist(sample(1:p)) / 3) + X.st[[i]] <- mvrnorm(n.spatial.xy^3, rep(0, p), Sigma.X) } X.all <- superMatrix(X.st) mean.trend.st <- X.all %*% beta.all @@ -39,44 +39,38 @@ params.all <- c(x.variance, ranges, marg.smoothness, nuggets, rho.vec) locs.list <- list() locs.listb <- list() for (i in 1:ny) { - loc1 <- seq(0, 1, length.out = n.spatial.xy) + rnorm(n.spatial.xy, 0, 0.001) - loc2 <- seq(0, 1, length.out = n.spatial.xy) + rnorm(n.spatial.xy, 0, 0.001) - loc3 <- seq(0, 365, length.out = n.spatial.xy) + rnorm(n.spatial.xy, 0, 0.001) - loc1b <- loc1 / ranges[2 * i - 1] - loc2b <- loc2 / ranges[2 * i - 1] - loc3b <- loc3 / ranges[2 * i] - locs.list[[i]] <- as.matrix(expand.grid(loc1, loc2, loc3)) - locs.listb[[i]] <- as.matrix(expand.grid(loc1b, loc2b, loc3b)) + loc1 <- seq(0, 1, length.out = n.spatial.xy) + rnorm(n.spatial.xy, 0, 0.001) + loc2 <- seq(0, 1, length.out = n.spatial.xy) + rnorm(n.spatial.xy, 0, 0.001) + loc3 <- seq(0, 365, length.out = n.spatial.xy) + rnorm(n.spatial.xy, 0, 0.001) + loc1b <- loc1 / ranges[2 * i - 1] + loc2b <- loc2 / ranges[2 * i - 1] + loc3b <- loc3 / ranges[2 * i] + locs.list[[i]] <- as.matrix(expand.grid(loc1, loc2, loc3)) + locs.listb[[i]] <- as.matrix(expand.grid(loc1b, loc2b, loc3b)) } Sigma.All <- matrix(nrow = ny * n.spatial.xy^3, ncol = ny * n.spatial.xy^3) for (i in 1:ny) { - for (j in i:ny) { - if (i == j) { - ndx1 <- ((i - 1) * n.spatial.xy^3 + 1):(n.spatial.xy^3 * i) - dij <- fields::rdist(locs.listb[[i]]) - Sigma.All[ndx1, ndx1] <- marg.var[i] * - fields::Matern(dij, - range = 1, - smoothness = marg.smoothness[i] - ) - } else { - ndx1 <- ((i - 1) * n.spatial.xy^3 + 1):(n.spatial.xy^3 * i) - ndx2 <- ((j - 1) * n.spatial.xy^3 + 1):(n.spatial.xy^3 * j) - dij <- fields::rdist(locs.listb[[i]], locs.listb[[j]]) - vii <- marg.smoothness[i] - vjj <- marg.smoothness[j] - vij <- (vii + vjj) / 2 - aii <- 1 - ajj <- 1 - aij <- sqrt((aii^2 + ajj^2) / 2) - Sigma.All[ndx1, ndx2] <- rho[i, j] * sqrt(marg.var[i]) * - sqrt(marg.var[j]) * aii^vii * ajj^vjj * gamma(vij) / - (aij^(2 * vij) * sqrt(gamma(vii) * gamma(vjj))) * - fields::Matern(dij, smoothness = vij, alpha = aij) - Sigma.All[ndx2, ndx1] <- t(Sigma.All[ndx1, ndx2]) - } + for (j in i:ny) { + if (i == j) { + ndx1 <- ((i - 1) * n.spatial.xy^3 + 1):(n.spatial.xy^3 * i) + dij <- fields::rdist(locs.listb[[i]]) + Sigma.All[ndx1, ndx1] <- marg.var[i] * fields::Matern(dij, range = 1, smoothness = marg.smoothness[i]) + } else { + ndx1 <- ((i - 1) * n.spatial.xy^3 + 1):(n.spatial.xy^3 * i) + ndx2 <- ((j - 1) * n.spatial.xy^3 + 1):(n.spatial.xy^3 * j) + dij <- fields::rdist(locs.listb[[i]], locs.listb[[j]]) + vii <- marg.smoothness[i] + vjj <- marg.smoothness[j] + vij <- (vii + vjj) / 2 + aii <- 1 + ajj <- 1 + aij <- sqrt((aii^2 + ajj^2) / 2) + Sigma.All[ndx1, ndx2] <- rho[i, j] * sqrt(marg.var[i]) * sqrt(marg.var[j]) * aii^vii * ajj^vjj * gamma(vij) / + (aij^(2 * vij) * sqrt(gamma(vii) * gamma(vjj))) * fields::Matern(dij, smoothness = vij, alpha = aij) + Sigma.All[ndx2, ndx1] <- t(Sigma.All[ndx1, ndx2]) } + } } L.C <- chol(Sigma.All) @@ -85,18 +79,18 @@ st.error <- rnorm(n.spatial.xy^3 * ny) %*% L.C nug.error <- NULL for (i in 1:ny) { - nug.error <- c(nug.error, nuggets[i] * rnorm(n.spatial.xy^3)) + nug.error <- c(nug.error, nuggets[i] * rnorm(n.spatial.xy^3)) } y.list <- list() for (i in 1:ny) { - ndx1 <- ((i - 1) * n.spatial.xy^3 + 1):(n.spatial.xy^3 * i) - y.list[[i]] <- mean.trend.st[ndx1] + st.error[ndx1] + nug.error[ndx1] + ndx1 <- ((i - 1) * n.spatial.xy^3 + 1):(n.spatial.xy^3 * i) + y.list[[i]] <- mean.trend.st[ndx1] + st.error[ndx1] + nug.error[ndx1] } rm( - ny, p, p.nz, n.spatial.xy, beta1, i, j, Sigma.X, mean.trend.st, n.rho, - loc1, loc2, ndx1, ndx2, dij, vii, vjj, vij, aii, ajj, aij, L.C, - st.error, nug.error, X.all, rho, rho.vec, ranges, Sigma.All, nuggets, - marg.smoothness, marg.var, locs.listb, loc3, loc1b, loc2b, loc3b + ny, p, p.nz, n.spatial.xy, beta1, i, j, Sigma.X, mean.trend.st, n.rho, + loc1, loc2, ndx1, ndx2, dij, vii, vjj, vij, aii, ajj, aij, L.C, + st.error, nug.error, X.all, rho, rho.vec, ranges, Sigma.All, nuggets, + marg.smoothness, marg.var, locs.listb, loc3, loc1b, loc2b, loc3b ) diff --git a/tests/testthat/sim_multivariate_lik.R b/tests/testthat/sim_multivariate_lik.R index 7bac2f5..d14dd05 100755 --- a/tests/testthat/sim_multivariate_lik.R +++ b/tests/testthat/sim_multivariate_lik.R @@ -23,39 +23,33 @@ params.all <- c(x.variance, ranges, marg.smoothness, nuggets, rho.vec) locs.list <- list() for (i in 1:ny) { - loc1 <- seq(0, 1, length.out = n.spatial.xy) + rnorm(n.spatial.xy, 0, 0.001) - loc2 <- seq(0, 1, length.out = n.spatial.xy) + rnorm(n.spatial.xy, 0, 0.001) - locs.list[[i]] <- as.matrix(expand.grid(loc1, loc2)) + loc1 <- seq(0, 1, length.out = n.spatial.xy) + rnorm(n.spatial.xy, 0, 0.001) + loc2 <- seq(0, 1, length.out = n.spatial.xy) + rnorm(n.spatial.xy, 0, 0.001) + locs.list[[i]] <- as.matrix(expand.grid(loc1, loc2)) } Sigma.All <- matrix(nrow = ny * n.spatial.xy^2, ncol = ny * n.spatial.xy^2) for (i in 1:ny) { - for (j in i:ny) { - if (i == j) { - ndx1 <- ((i - 1) * n.spatial.xy^2 + 1):(n.spatial.xy^2 * i) - dij <- fields::rdist(locs.list[[i]]) - Sigma.All[ndx1, ndx1] <- marg.var[i] * - fields::Matern(dij, - range = ranges[i], - smoothness = marg.smoothness[i] - ) - } else { - ndx1 <- ((i - 1) * n.spatial.xy^2 + 1):(n.spatial.xy^2 * i) - ndx2 <- ((j - 1) * n.spatial.xy^2 + 1):(n.spatial.xy^2 * j) - dij <- fields::rdist(locs.list[[i]], locs.list[[j]]) - vii <- marg.smoothness[i] - vjj <- marg.smoothness[j] - vij <- (vii + vjj) / 2 - aii <- 1 / ranges[i] - ajj <- 1 / ranges[j] - aij <- sqrt((aii^2 + ajj^2) / 2) - Sigma.All[ndx1, ndx2] <- rho[i, j] * sqrt(marg.var[i]) * - sqrt(marg.var[j]) * aii^vii * ajj^vjj * gamma(vij) / - (aij^(2 * vij) * sqrt(gamma(vii) * gamma(vjj))) * - fields::Matern(dij, smoothness = vij, alpha = aij) - Sigma.All[ndx2, ndx1] <- t(Sigma.All[ndx1, ndx2]) - } + for (j in i:ny) { + if (i == j) { + ndx1 <- ((i - 1) * n.spatial.xy^2 + 1):(n.spatial.xy^2 * i) + dij <- fields::rdist(locs.list[[i]]) + Sigma.All[ndx1, ndx1] <- marg.var[i] * fields::Matern(dij, range = ranges[i], smoothness = marg.smoothness[i]) + } else { + ndx1 <- ((i - 1) * n.spatial.xy^2 + 1):(n.spatial.xy^2 * i) + ndx2 <- ((j - 1) * n.spatial.xy^2 + 1):(n.spatial.xy^2 * j) + dij <- fields::rdist(locs.list[[i]], locs.list[[j]]) + vii <- marg.smoothness[i] + vjj <- marg.smoothness[j] + vij <- (vii + vjj) / 2 + aii <- 1 / ranges[i] + ajj <- 1 / ranges[j] + aij <- sqrt((aii^2 + ajj^2) / 2) + Sigma.All[ndx1, ndx2] <- rho[i, j] * sqrt(marg.var[i]) * sqrt(marg.var[j]) * aii^vii * ajj^vjj * gamma(vij) / + (aij^(2 * vij) * sqrt(gamma(vii) * gamma(vjj))) * fields::Matern(dij, smoothness = vij, alpha = aij) + Sigma.All[ndx2, ndx1] <- t(Sigma.All[ndx1, ndx2]) } + } } L.C <- chol(Sigma.All) @@ -64,18 +58,18 @@ st.error <- rnorm(n.spatial.xy^2 * ny) %*% L.C nug.error <- NULL for (i in 1:ny) { - nug.error <- c(nug.error, nuggets[i] * rnorm(n.spatial.xy^2)) + nug.error <- c(nug.error, nuggets[i] * rnorm(n.spatial.xy^2)) } y.list <- list() for (i in 1:ny) { - ndx1 <- ((i - 1) * n.spatial.xy^2 + 1):(n.spatial.xy^2 * i) - y.list[[i]] <- st.error[ndx1] + nug.error[ndx1] + ndx1 <- ((i - 1) * n.spatial.xy^2 + 1):(n.spatial.xy^2 * i) + y.list[[i]] <- st.error[ndx1] + nug.error[ndx1] } rm( - ny, n.spatial.xy, i, j, n.rho, - loc1, loc2, ndx1, ndx2, dij, vii, vjj, vij, aii, ajj, aij, L.C, - st.error, nug.error, rho, rho.vec, ranges, Sigma.All, nuggets, - marg.smoothness, marg.var + ny, n.spatial.xy, i, j, n.rho, + loc1, loc2, ndx1, ndx2, dij, vii, vjj, vij, aii, ajj, aij, L.C, + st.error, nug.error, rho, rho.vec, ranges, Sigma.All, nuggets, + marg.smoothness, marg.var ) diff --git a/tests/testthat/sim_multivariate_small.R b/tests/testthat/sim_multivariate_small.R index 0fbcd93..aefab32 100755 --- a/tests/testthat/sim_multivariate_small.R +++ b/tests/testthat/sim_multivariate_small.R @@ -12,9 +12,9 @@ beta.all <- rep(beta1, ny) X.st <- list() for (i in 1:ny) { - set.seed(1212 + i) - Sigma.X <- exp(-rdist(sample(1:p)) / 3) - X.st[[i]] <- mvrnorm(n.spatial.xy^2, rep(0, p), Sigma.X) + set.seed(1212 + i) + Sigma.X <- exp(-rdist(sample(1:p)) / 3) + X.st[[i]] <- mvrnorm(n.spatial.xy^2, rep(0, p), Sigma.X) } X.all <- superMatrix(X.st) mean.trend.st <- X.all %*% beta.all @@ -35,39 +35,33 @@ params.all <- c(x.variance, ranges, marg.smoothness, nuggets, rho.vec) locs.list <- list() for (i in 1:ny) { - loc1 <- seq(0, 1, length.out = n.spatial.xy) + rnorm(n.spatial.xy, 0, 0.001) - loc2 <- seq(0, 1, length.out = n.spatial.xy) + rnorm(n.spatial.xy, 0, 0.001) - locs.list[[i]] <- as.matrix(expand.grid(loc1, loc2)) + loc1 <- seq(0, 1, length.out = n.spatial.xy) + rnorm(n.spatial.xy, 0, 0.001) + loc2 <- seq(0, 1, length.out = n.spatial.xy) + rnorm(n.spatial.xy, 0, 0.001) + locs.list[[i]] <- as.matrix(expand.grid(loc1, loc2)) } Sigma.All <- matrix(nrow = ny * n.spatial.xy^2, ncol = ny * n.spatial.xy^2) for (i in 1:ny) { - for (j in i:ny) { - if (i == j) { - ndx1 <- ((i - 1) * n.spatial.xy^2 + 1):(n.spatial.xy^2 * i) - dij <- fields::rdist(locs.list[[i]]) - Sigma.All[ndx1, ndx1] <- marg.var[i] * - fields::Matern(dij, - range = ranges[i], - smoothness = marg.smoothness[i] - ) - } else { - ndx1 <- ((i - 1) * n.spatial.xy^2 + 1):(n.spatial.xy^2 * i) - ndx2 <- ((j - 1) * n.spatial.xy^2 + 1):(n.spatial.xy^2 * j) - dij <- fields::rdist(locs.list[[i]], locs.list[[j]]) - vii <- marg.smoothness[i] - vjj <- marg.smoothness[j] - vij <- (vii + vjj) / 2 - aii <- 1 / ranges[i] - ajj <- 1 / ranges[j] - aij <- sqrt((aii^2 + ajj^2) / 2) - Sigma.All[ndx1, ndx2] <- rho[i, j] * sqrt(marg.var[i]) * - sqrt(marg.var[j]) * aii^vii * ajj^vjj * gamma(vij) / - (aij^(2 * vij) * sqrt(gamma(vii) * gamma(vjj))) * - fields::Matern(dij, smoothness = vij, alpha = aij) - Sigma.All[ndx2, ndx1] <- t(Sigma.All[ndx1, ndx2]) - } + for (j in i:ny) { + if (i == j) { + ndx1 <- ((i - 1) * n.spatial.xy^2 + 1):(n.spatial.xy^2 * i) + dij <- fields::rdist(locs.list[[i]]) + Sigma.All[ndx1, ndx1] <- marg.var[i] * fields::Matern(dij, range = ranges[i], smoothness = marg.smoothness[i]) + } else { + ndx1 <- ((i - 1) * n.spatial.xy^2 + 1):(n.spatial.xy^2 * i) + ndx2 <- ((j - 1) * n.spatial.xy^2 + 1):(n.spatial.xy^2 * j) + dij <- fields::rdist(locs.list[[i]], locs.list[[j]]) + vii <- marg.smoothness[i] + vjj <- marg.smoothness[j] + vij <- (vii + vjj) / 2 + aii <- 1 / ranges[i] + ajj <- 1 / ranges[j] + aij <- sqrt((aii^2 + ajj^2) / 2) + Sigma.All[ndx1, ndx2] <- rho[i, j] * sqrt(marg.var[i]) * sqrt(marg.var[j]) * aii^vii * ajj^vjj * gamma(vij) / + (aij^(2 * vij) * sqrt(gamma(vii) * gamma(vjj))) * fields::Matern(dij, smoothness = vij, alpha = aij) + Sigma.All[ndx2, ndx1] <- t(Sigma.All[ndx1, ndx2]) } + } } L.C <- chol(Sigma.All) @@ -76,18 +70,18 @@ st.error <- rnorm(n.spatial.xy^2 * ny) %*% L.C nug.error <- NULL for (i in 1:ny) { - nug.error <- c(nug.error, nuggets[i] * rnorm(n.spatial.xy^2)) + nug.error <- c(nug.error, nuggets[i] * rnorm(n.spatial.xy^2)) } y.list <- list() for (i in 1:ny) { - ndx1 <- ((i - 1) * n.spatial.xy^2 + 1):(n.spatial.xy^2 * i) - y.list[[i]] <- mean.trend.st[ndx1] + st.error[ndx1] + nug.error[ndx1] + ndx1 <- ((i - 1) * n.spatial.xy^2 + 1):(n.spatial.xy^2 * i) + y.list[[i]] <- mean.trend.st[ndx1] + st.error[ndx1] + nug.error[ndx1] } rm( - ny, p, p.nz, n.spatial.xy, beta1, i, j, Sigma.X, mean.trend.st, n.rho, - loc1, loc2, ndx1, ndx2, dij, vii, vjj, vij, aii, ajj, aij, L.C, - st.error, nug.error, X.all, rho, rho.vec, ranges, Sigma.All, nuggets, - marg.smoothness, marg.var + ny, p, p.nz, n.spatial.xy, beta1, i, j, Sigma.X, mean.trend.st, n.rho, + loc1, loc2, ndx1, ndx2, dij, vii, vjj, vij, aii, ajj, aij, L.C, + st.error, nug.error, X.all, rho, rho.vec, ranges, Sigma.All, nuggets, + marg.smoothness, marg.var ) diff --git a/tests/testthat/sim_multivariate_small_pred.R b/tests/testthat/sim_multivariate_small_pred.R index 0fda2db..961b63b 100755 --- a/tests/testthat/sim_multivariate_small_pred.R +++ b/tests/testthat/sim_multivariate_small_pred.R @@ -14,8 +14,8 @@ beta.all <- rep(beta1, ny) X.st <- list() for (i in 1:ny) { - Sigma.X <- exp(-rdist(sample(1:p)) / 3) - X.st[[i]] <- mvrnorm(n.spatial.xy^2, rep(0, p), Sigma.X) + Sigma.X <- exp(-rdist(sample(1:p)) / 3) + X.st[[i]] <- mvrnorm(n.spatial.xy^2, rep(0, p), Sigma.X) } X.all <- superMatrix(X.st) mean.trend.st <- X.all %*% beta.all @@ -36,39 +36,33 @@ params.all <- c(x.variance, ranges, marg.smoothness, nuggets, rho.vec) locs.list <- list() for (i in 1:ny) { - loc1 <- seq(0, 1, length.out = n.spatial.xy) + rnorm(n.spatial.xy, 0, 0.001) - loc2 <- seq(0, 1, length.out = n.spatial.xy) + rnorm(n.spatial.xy, 0, 0.001) - locs.list[[i]] <- as.matrix(expand.grid(loc1, loc2)) + loc1 <- seq(0, 1, length.out = n.spatial.xy) + rnorm(n.spatial.xy, 0, 0.001) + loc2 <- seq(0, 1, length.out = n.spatial.xy) + rnorm(n.spatial.xy, 0, 0.001) + locs.list[[i]] <- as.matrix(expand.grid(loc1, loc2)) } Sigma.All <- matrix(nrow = ny * n.spatial.xy^2, ncol = ny * n.spatial.xy^2) for (i in 1:ny) { - for (j in i:ny) { - if (i == j) { - ndx1 <- ((i - 1) * n.spatial.xy^2 + 1):(n.spatial.xy^2 * i) - dij <- fields::rdist(locs.list[[i]]) - Sigma.All[ndx1, ndx1] <- marg.var[i] * - fields::Matern(dij, - range = ranges[i], - smoothness = marg.smoothness[i] - ) - } else { - ndx1 <- ((i - 1) * n.spatial.xy^2 + 1):(n.spatial.xy^2 * i) - ndx2 <- ((j - 1) * n.spatial.xy^2 + 1):(n.spatial.xy^2 * j) - dij <- fields::rdist(locs.list[[i]], locs.list[[j]]) - vii <- marg.smoothness[i] - vjj <- marg.smoothness[j] - vij <- (vii + vjj) / 2 - aii <- 1 / ranges[i] - ajj <- 1 / ranges[j] - aij <- sqrt((aii^2 + ajj^2) / 2) - Sigma.All[ndx1, ndx2] <- rho[i, j] * sqrt(marg.var[i]) * - sqrt(marg.var[j]) * aii^vii * ajj^vjj * gamma(vij) / - (aij^(2 * vij) * sqrt(gamma(vii) * gamma(vjj))) * - fields::Matern(dij, smoothness = vij, alpha = aij) - Sigma.All[ndx2, ndx1] <- t(Sigma.All[ndx1, ndx2]) - } + for (j in i:ny) { + if (i == j) { + ndx1 <- ((i - 1) * n.spatial.xy^2 + 1):(n.spatial.xy^2 * i) + dij <- fields::rdist(locs.list[[i]]) + Sigma.All[ndx1, ndx1] <- marg.var[i] * fields::Matern(dij, range = ranges[i], smoothness = marg.smoothness[i]) + } else { + ndx1 <- ((i - 1) * n.spatial.xy^2 + 1):(n.spatial.xy^2 * i) + ndx2 <- ((j - 1) * n.spatial.xy^2 + 1):(n.spatial.xy^2 * j) + dij <- fields::rdist(locs.list[[i]], locs.list[[j]]) + vii <- marg.smoothness[i] + vjj <- marg.smoothness[j] + vij <- (vii + vjj) / 2 + aii <- 1 / ranges[i] + ajj <- 1 / ranges[j] + aij <- sqrt((aii^2 + ajj^2) / 2) + Sigma.All[ndx1, ndx2] <- rho[i, j] * sqrt(marg.var[i]) * sqrt(marg.var[j]) * aii^vii * ajj^vjj * gamma(vij) / + (aij^(2 * vij) * sqrt(gamma(vii) * gamma(vjj))) * fields::Matern(dij, smoothness = vij, alpha = aij) + Sigma.All[ndx2, ndx1] <- t(Sigma.All[ndx1, ndx2]) } + } } L.C <- chol(Sigma.All) @@ -77,7 +71,7 @@ st.error <- rnorm(n.spatial.xy^2 * ny) %*% L.C nug.error <- NULL for (i in 1:ny) { - nug.error <- c(nug.error, nuggets[i] * rnorm(n.spatial.xy^2)) + nug.error <- c(nug.error, nuggets[i] * rnorm(n.spatial.xy^2)) } y.list <- list() @@ -88,22 +82,22 @@ y.list.otst <- list() X.st.otr <- list() X.st.otst <- list() for (i in 1:ny) { - ndx1 <- ((i - 1) * n.spatial.xy^2 + 1):(n.spatial.xy^2 * i) - y.list[[i]] <- mean.trend.st[ndx1] + st.error[ndx1] + nug.error[ndx1] - nn <- n.spatial.xy^2 - otr <- rep(FALSE, nn) - otr[sample(1:nn, size = floor(nn / 2))] <- TRUE - locs.list.otr[[i]] <- locs.list[[i]][otr, ] - locs.list.otst[[i]] <- locs.list[[i]][!otr, ] - y.list.otr[[i]] <- y.list[[i]][otr] - y.list.otst[[i]] <- y.list[[i]][!otr] - X.st.otr[[i]] <- X.st[[i]][otr, ] - X.st.otst[[i]] <- X.st[[i]][!otr, ] + ndx1 <- ((i - 1) * n.spatial.xy^2 + 1):(n.spatial.xy^2 * i) + y.list[[i]] <- mean.trend.st[ndx1] + st.error[ndx1] + nug.error[ndx1] + nn <- n.spatial.xy^2 + otr <- rep(FALSE, nn) + otr[sample(1:nn, size = floor(nn / 2))] <- TRUE + locs.list.otr[[i]] <- locs.list[[i]][otr, ] + locs.list.otst[[i]] <- locs.list[[i]][!otr, ] + y.list.otr[[i]] <- y.list[[i]][otr] + y.list.otst[[i]] <- y.list[[i]][!otr] + X.st.otr[[i]] <- X.st[[i]][otr, ] + X.st.otst[[i]] <- X.st[[i]][!otr, ] } rm( - ny, p, p.nz, n.spatial.xy, beta1, i, j, Sigma.X, mean.trend.st, n.rho, - loc1, loc2, ndx1, ndx2, dij, vii, vjj, vij, aii, ajj, aij, L.C, - st.error, nug.error, X.all, rho, rho.vec, ranges, Sigma.All, nuggets, - marg.smoothness, marg.var, x.variance, nn, otr + ny, p, p.nz, n.spatial.xy, beta1, i, j, Sigma.X, mean.trend.st, n.rho, + loc1, loc2, ndx1, ndx2, dij, vii, vjj, vij, aii, ajj, aij, L.C, + st.error, nug.error, X.all, rho, rho.vec, ranges, Sigma.All, nuggets, + marg.smoothness, marg.var, x.variance, nn, otr ) diff --git a/tests/testthat/sim_vecchia.R b/tests/testthat/sim_vecchia.R index 63c77cf..ecfb02c 100755 --- a/tests/testthat/sim_vecchia.R +++ b/tests/testthat/sim_vecchia.R @@ -27,8 +27,8 @@ loc2 <- seq(0, 1, length.out = n.spatial.xy) + rnorm(n.spatial.xy, 0, 0.001) locs <- as.matrix(expand.grid(loc1, loc2)) Sigma.All <- marg.var * Matern(rdist(locs), - range = ranges, - smoothness = marg.smoothness + range = ranges, + smoothness = marg.smoothness ) L.C <- chol(Sigma.All) @@ -37,7 +37,7 @@ nug.error <- nuggets * rnorm(n.spatial.xy^2) y <- mean.trend.st + st.error + nug.error rm( - p, p.nz, n.spatial.xy, beta1, Sigma.X, mean.trend.st, loc1, loc2, L.C, - st.error, nug.error, ranges, Sigma.All, nuggets, marg.smoothness, marg.var, - x.variance + p, p.nz, n.spatial.xy, beta1, Sigma.X, mean.trend.st, loc1, loc2, L.C, + st.error, nug.error, ranges, Sigma.All, nuggets, marg.smoothness, marg.var, + x.variance ) diff --git a/tests/testthat/sim_vecchia_pred.R b/tests/testthat/sim_vecchia_pred.R index dbda626..566ff72 100755 --- a/tests/testthat/sim_vecchia_pred.R +++ b/tests/testthat/sim_vecchia_pred.R @@ -27,8 +27,8 @@ loc2 <- seq(0, 1, length.out = n.spatial.xy) + rnorm(n.spatial.xy, 0, 0.001) locs <- as.matrix(expand.grid(loc1, loc2)) Sigma.All <- marg.var * Matern(rdist(locs), - range = ranges, - smoothness = marg.smoothness + range = ranges, + smoothness = marg.smoothness ) L.C <- chol(Sigma.All) @@ -48,7 +48,7 @@ X.otr <- X[otr, ] X.otst <- X[!otr, ] rm( - p, p.nz, n.spatial.xy, beta1, Sigma.X, mean.trend.st, loc1, loc2, L.C, - st.error, nug.error, ranges, Sigma.All, nuggets, marg.smoothness, marg.var, - x.variance, nn, otr + p, p.nz, n.spatial.xy, beta1, Sigma.X, mean.trend.st, loc1, loc2, L.C, + st.error, nug.error, ranges, Sigma.All, nuggets, marg.smoothness, marg.var, + x.variance, nn, otr ) diff --git a/tests/testthat/sim_vecchia_small.R b/tests/testthat/sim_vecchia_small.R index 39905dc..1bb9325 100755 --- a/tests/testthat/sim_vecchia_small.R +++ b/tests/testthat/sim_vecchia_small.R @@ -27,8 +27,8 @@ loc2 <- seq(0, 1, length.out = n.spatial.xy) + rnorm(n.spatial.xy, 0, 0.001) locs <- as.matrix(expand.grid(loc1, loc2)) Sigma.All <- marg.var * Matern(rdist(locs), - range = ranges, - smoothness = marg.smoothness + range = ranges, + smoothness = marg.smoothness ) L.C <- chol(Sigma.All) @@ -37,7 +37,7 @@ nug.error <- nuggets * rnorm(n.spatial.xy^2) y <- mean.trend.st + st.error + nug.error rm( - p, p.nz, n.spatial.xy, beta1, Sigma.X, mean.trend.st, loc1, loc2, L.C, - st.error, nug.error, ranges, Sigma.All, nuggets, marg.smoothness, marg.var, - x.variance + p, p.nz, n.spatial.xy, beta1, Sigma.X, mean.trend.st, loc1, loc2, L.C, + st.error, nug.error, ranges, Sigma.All, nuggets, marg.smoothness, marg.var, + x.variance ) diff --git a/tests/testthat/sim_vecchia_small_pred.R b/tests/testthat/sim_vecchia_small_pred.R index 36ba071..c082605 100755 --- a/tests/testthat/sim_vecchia_small_pred.R +++ b/tests/testthat/sim_vecchia_small_pred.R @@ -27,8 +27,8 @@ loc2 <- seq(0, 1, length.out = n.spatial.xy) + rnorm(n.spatial.xy, 0, 0.001) locs <- as.matrix(expand.grid(loc1, loc2)) Sigma.All <- marg.var * Matern(rdist(locs), - range = ranges, - smoothness = marg.smoothness + range = ranges, + smoothness = marg.smoothness ) L.C <- chol(Sigma.All) @@ -48,7 +48,7 @@ X.otr <- X[otr, ] X.otst <- X[!otr, ] rm( - p, p.nz, n.spatial.xy, beta1, Sigma.X, mean.trend.st, loc1, loc2, L.C, - st.error, nug.error, ranges, Sigma.All, nuggets, marg.smoothness, marg.var, - x.variance, nn, otr + p, p.nz, n.spatial.xy, beta1, Sigma.X, mean.trend.st, loc1, loc2, L.C, + st.error, nug.error, ranges, Sigma.All, nuggets, marg.smoothness, marg.var, + x.variance, nn, otr ) diff --git a/tests/testthat/sim_vecchia_st.R b/tests/testthat/sim_vecchia_st.R index da65aa4..401bbee 100755 --- a/tests/testthat/sim_vecchia_st.R +++ b/tests/testthat/sim_vecchia_st.R @@ -32,8 +32,8 @@ locss[, 1:2] <- locss[, 1:2] / ranges[1] locss[, 3] <- locs[, 3] / ranges[2] Sigma.All <- marg.var * Matern(rdist(locss), - range = 1, - smoothness = marg.smoothness + range = 1, + smoothness = marg.smoothness ) L.C <- chol(Sigma.All) @@ -42,7 +42,7 @@ nug.error <- nuggets * rnorm(n.spatial.xy^3) y <- mean.trend.st + st.error + nug.error rm( - p, p.nz, n.spatial.xy, beta1, Sigma.X, mean.trend.st, loc1, loc2, loc3, L.C, - st.error, nug.error, ranges, Sigma.All, nuggets, marg.smoothness, marg.var, - x.variance, locss + p, p.nz, n.spatial.xy, beta1, Sigma.X, mean.trend.st, loc1, loc2, loc3, L.C, + st.error, nug.error, ranges, Sigma.All, nuggets, marg.smoothness, marg.var, + x.variance, locss ) diff --git a/tests/testthat/test-Log_Likelihood.R b/tests/testthat/test-Log_Likelihood.R index 943d9ac..9326acb 100755 --- a/tests/testthat/test-Log_Likelihood.R +++ b/tests/testthat/test-Log_Likelihood.R @@ -1,371 +1,366 @@ test_that("negloglik_vecchia", { - set.seed(1234) - - y.orig <- rnorm(100) - - locs.dim1 <- seq(1, 10, length = 10) + rnorm(10, 0, 0.1) - locs.dim2 <- seq(1, 10, length = 10) + rnorm(10, 0, 0.1) - locs <- as.matrix(expand.grid(locs.dim1, locs.dim2)) - - covmat.true <- Matern(rdist(locs), range = 1, smoothness = 0.5) + diag(100) - - y <- y.orig %*% chol(covmat.true) - y <- as.vector(y) - - logparams <- create.initial.values.flex(0.9, 0.5, 0.5, 0.1, 1, 1) - pseq <- create.param.sequence(1) - vec.approx <- vecchia_specify(locs, 5) - LL.vecchia <- negloglik_vecchia(logparams, y, vec.approx, pseq) - expect_equal(194.75, LL.vecchia, tolerance = 1e-2) - - vec.mapprox <- vecchia_Mspecify(list(locs), 5) - LL.mvecchia <- mvnegloglik(logparams, vec.mapprox, y, pseq, 1) - # Univariate likelihood should equal multivariate likelihood - expect_equal(LL.vecchia, LL.mvecchia, tolerance = 1e-1) - - logparams2 <- create.initial.values.flex(0.9, 1, 0.5, 0.1, 1, 1) - vec.approx2 <- vec.approx - vec.approx2$locsord <- vec.approx$locsord / 0.5 - LL.vecchia2 <- negloglik_vecchia(logparams2, y, vec.approx2, pseq) - LL.vecchia.st <- negloglik_vecchia_ST(logparams, y, vec.approx, pseq, 1, 1) - vec.mapprox2 <- vec.mapprox - vec.mapprox2$locsord <- vec.mapprox$locsord / 0.5 - LL.mvecchia2 <- mvnegloglik_ST(logparams2, vec.mapprox2, y, pseq, 1, 1, 1) - - # Likelihood should equal spatiotemporal likelihood after scaling locs - expect_equal(LL.vecchia, LL.vecchia2, tolerance = 1e-3) - expect_equal(LL.vecchia, LL.vecchia.st, tolerance = 1e-3) - expect_equal(LL.vecchia, LL.mvecchia2, tolerance = 1e-1) + set.seed(1234) + + y.orig <- rnorm(100) + locs.dim1 <- seq(1, 10, length = 10) + rnorm(10, 0, 0.1) + locs.dim2 <- seq(1, 10, length = 10) + rnorm(10, 0, 0.1) + locs <- as.matrix(expand.grid(locs.dim1, locs.dim2)) + covmat.true <- Matern(rdist(locs), range = 1, smoothness = 0.5) + diag(100) + + y <- y.orig %*% chol(covmat.true) + y <- as.vector(y) + + logparams <- create.initial.values.flex(0.9, 0.5, 0.5, 0.1, 1, 1) + pseq <- create.param.sequence(1) + vec.approx <- vecchia_specify(locs, 5) + LL.vecchia <- negloglik_vecchia(logparams, y, vec.approx, pseq) + expect_equal(194.75, LL.vecchia, tolerance = 1e-2) + + vec.mapprox <- vecchia_Mspecify(list(locs), 5) + LL.mvecchia <- mvnegloglik(logparams, vec.mapprox, y, pseq, 1) + # Univariate likelihood should equal multivariate likelihood + expect_equal(LL.vecchia, LL.mvecchia, tolerance = 1e-1) + + logparams2 <- create.initial.values.flex(0.9, 1, 0.5, 0.1, 1, 1) + vec.approx2 <- vec.approx + vec.approx2$locsord <- vec.approx$locsord / 0.5 + LL.vecchia2 <- negloglik_vecchia(logparams2, y, vec.approx2, pseq) + LL.vecchia.st <- negloglik_vecchia_ST(logparams, y, vec.approx, pseq, 1, 1) + vec.mapprox2 <- vec.mapprox + vec.mapprox2$locsord <- vec.mapprox$locsord / 0.5 + LL.mvecchia2 <- mvnegloglik_ST(logparams2, vec.mapprox2, y, pseq, 1, 1, 1) + + # Likelihood should equal spatiotemporal likelihood after scaling locs + expect_equal(LL.vecchia, LL.vecchia2, tolerance = 1e-3) + expect_equal(LL.vecchia, LL.vecchia.st, tolerance = 1e-3) + expect_equal(LL.vecchia, LL.mvecchia2, tolerance = 1e-1) }) test_that("negloglik_vecchia_ST", { - set.seed(1234) - - y.orig <- rnorm(125) - - locs.dim1 <- seq(1, 5, length = 5) + rnorm(5, 0, 0.1) - locs.dim2 <- seq(1, 5, length = 5) + rnorm(5, 0, 0.1) - locs.dim3 <- seq(1, 5, length = 5) + rnorm(5, 0, 0.1) - locs <- as.matrix(expand.grid(locs.dim1, locs.dim2, locs.dim3)) - - covmat.true <- Matern(rdist(locs), range = 1, smoothness = 0.5) + diag(125) - - y <- y.orig %*% chol(covmat.true) - y <- as.vector(y) - - logparams <- create.initial.values.flex(0.9, c(2, 3), 0.5, 0.1, 1, 1) - pseq <- create.param.sequence(1, 2) - vec.approx <- vecchia_specify(locs, 5) - LL.vecchia <- negloglik_vecchia_ST( - logparams, y, vec.approx, pseq, - c(1, 1, 2), 2 - ) - expect_equal(284.73, LL.vecchia, tolerance = 1e-2) - - logparams2 <- create.initial.values.flex(0.9, 1, 0.5, 0.1, 1, 1) - pseq2 <- create.param.sequence(1) - vec.approx2 <- vec.approx - vec.approx2$locsord[, 1:2] <- vec.approx$locsord[, 1:2] / 2 - vec.approx2$locsord[, 3] <- vec.approx$locsord[, 3] / 3 - LL.vecchia2 <- negloglik_vecchia(logparams2, y, vec.approx2, pseq2) - expect_equal(LL.vecchia, LL.vecchia2, tolerance = 1e-3) + set.seed(1234) + + y.orig <- rnorm(125) + + locs.dim1 <- seq(1, 5, length = 5) + rnorm(5, 0, 0.1) + locs.dim2 <- seq(1, 5, length = 5) + rnorm(5, 0, 0.1) + locs.dim3 <- seq(1, 5, length = 5) + rnorm(5, 0, 0.1) + locs <- as.matrix(expand.grid(locs.dim1, locs.dim2, locs.dim3)) + + covmat.true <- Matern(rdist(locs), range = 1, smoothness = 0.5) + diag(125) + + y <- y.orig %*% chol(covmat.true) + y <- as.vector(y) + + logparams <- create.initial.values.flex(0.9, c(2, 3), 0.5, 0.1, 1, 1) + pseq <- create.param.sequence(1, 2) + vec.approx <- vecchia_specify(locs, 5) + LL.vecchia <- negloglik_vecchia_ST( + logparams, y, vec.approx, pseq, + c(1, 1, 2), 2 + ) + expect_equal(284.73, LL.vecchia, tolerance = 1e-2) + + logparams2 <- create.initial.values.flex(0.9, 1, 0.5, 0.1, 1, 1) + pseq2 <- create.param.sequence(1) + vec.approx2 <- vec.approx + vec.approx2$locsord[, 1:2] <- vec.approx$locsord[, 1:2] / 2 + vec.approx2$locsord[, 3] <- vec.approx$locsord[, 3] / 3 + LL.vecchia2 <- negloglik_vecchia(logparams2, y, vec.approx2, pseq2) + expect_equal(LL.vecchia, LL.vecchia2, tolerance = 1e-3) }) test_that("negloglik.full", { - set.seed(1234) + set.seed(1234) - y.orig <- rnorm(100) + y.orig <- rnorm(100) - locs.dim <- seq(1, 10, length = sqrt(length(y.orig))) - locs <- as.matrix(expand.grid(locs.dim, locs.dim)) + locs.dim <- seq(1, 10, length = sqrt(length(y.orig))) + locs <- as.matrix(expand.grid(locs.dim, locs.dim)) - covmat.true <- Matern(rdist(locs), range = 1, smoothness = 0.5) + diag(100) + covmat.true <- Matern(rdist(locs), range = 1, smoothness = 0.5) + diag(100) - y <- y.orig %*% chol(covmat.true) - y <- as.vector(y) + y <- y.orig %*% chol(covmat.true) + y <- as.vector(y) - params.init <- rep(NA, 4) - params.init[1] <- 0.9 * var(y) - params.init[2] <- 1 - params.init[3] <- 0.5 - params.init[4] <- 0.1 * var(y) + params.init <- rep(NA, 4) + params.init[1] <- 0.9 * var(y) + params.init[2] <- 1 + params.init[3] <- 0.5 + params.init[4] <- 0.1 * var(y) - params.init <- create.initial.values.flex( - params.init[1], - params.init[2], - params.init[3], - params.init[4], - 1, 1 - ) + params.init <- create.initial.values.flex( + params.init[1], + params.init[2], + params.init[3], + params.init[4], + 1, 1 + ) - d <- rdist(locs) - pseq <- create.param.sequence(1) + d <- rdist(locs) + pseq <- create.param.sequence(1) - res.optim.NM <- optim( - par = params.init, fn = negloglik.full, d = d, y = y, - param.seq = pseq, control = list(maxit = 5000) - ) + res.optim.NM <- optim( + par = params.init, fn = negloglik.full, d = d, y = y, + param.seq = pseq, control = list(maxit = 5000) + ) - LL.full <- negloglik.full(res.optim.NM$par, d, y, pseq) + LL.full <- negloglik.full(res.optim.NM$par, d, y, pseq) - params.final <- c( - exp(res.optim.NM$par[1:2]), - gtools::inv.logit(res.optim.NM$par[3], 0, 2.5), - exp(res.optim.NM$par[4]) - ) + params.final <- c( + exp(res.optim.NM$par[1:2]), + gtools::inv.logit(res.optim.NM$par[3], 0, 2.5), + exp(res.optim.NM$par[4]) + ) - pgp.params <- create.initial.values.flex( - params.final[1], - params.final[2], - params.final[3], - params.final[4], - 1, 1 - ) + pgp.params <- create.initial.values.flex( + params.final[1], + params.final[2], + params.final[3], + params.final[4], + 1, 1 + ) - LL.full.pgp <- mvnegloglik.full(pgp.params, list(locs), y, pseq) + LL.full.pgp <- mvnegloglik.full(pgp.params, list(locs), y, pseq) - vec.approx <- vecchia_specify(locs, 99) + vec.approx <- vecchia_specify(locs, 99) - LL.vecchia <- -1 * vecchia_likelihood( - y, vec.approx, params.final[1:3], - params.final[4] - ) + LL.vecchia <- -1 * vecchia_likelihood( + y, vec.approx, params.final[1:3], + params.final[4] + ) - vec.approx.pgp <- vecchia_Mspecify(list(locs), 99) - vec.U.pgp <- createUMultivariate(vec.approx.pgp, c(params.final, 1)) + vec.approx.pgp <- vecchia_Mspecify(list(locs), 99) + vec.U.pgp <- createUMultivariate(vec.approx.pgp, c(params.final, 1)) - LL.pgp <- -1 * GPvecchia:::vecchia_likelihood_U(y, vec.U.pgp) + LL.pgp <- -1 * GPvecchia:::vecchia_likelihood_U(y, vec.U.pgp) - expect_equal(173.315, LL.full, tolerance = 1e-3) - # Univariate likelihood should equal the multivariate likelihood - expect_equal(LL.full, LL.full.pgp, tolerance = 1e-3) - # Full likelihood should equal both the univariate and multivariate - # Vecchia approximations - expect_equal(LL.full, LL.vecchia, tolerance = 1e-3) - expect_equal(LL.full, LL.pgp, tolerance = 1e-3) + expect_equal(173.315, LL.full, tolerance = 1e-3) + # Univariate likelihood should equal the multivariate likelihood + expect_equal(LL.full, LL.full.pgp, tolerance = 1e-3) + # Full likelihood should equal both the univariate and multivariate + # Vecchia approximations + expect_equal(LL.full, LL.vecchia, tolerance = 1e-3) + expect_equal(LL.full, LL.pgp, tolerance = 1e-3) }) test_that("negloglik_full_ST", { - set.seed(1234) - - y.orig <- rnorm(125) - - locs.dim1 <- seq(1, 5, length = 5) + rnorm(5, 0, 0.1) - locs.dim2 <- seq(1, 5, length = 5) + rnorm(5, 0, 0.1) - locs.dim3 <- seq(1, 5, length = 5) + rnorm(5, 0, 0.1) - locs <- as.matrix(expand.grid(locs.dim1, locs.dim2, locs.dim3)) - - covmat.true <- Matern(rdist(locs), range = 1, smoothness = 0.5) + diag(125) - - y <- y.orig %*% chol(covmat.true) - y <- as.vector(y) - - logparams <- create.initial.values.flex(0.9, c(2, 3), 0.5, 0.1, 1, 1) - pseq <- create.param.sequence(1, 2) - LL.full <- negloglik_full_ST( - logparams, locs, y, pseq, - c(1, 1, 2), 2 - ) - expect_equal(286.84, LL.full, tolerance = 1e-2) - - logparams2 <- create.initial.values.flex(0.9, 1, 0.5, 0.1, 1, 1) - pseq2 <- create.param.sequence(1) - locs2 <- locs - locs2[, 1:2] <- locs[, 1:2] / 2 - locs2[, 3] <- locs[, 3] / 3 - LL.full2 <- negloglik.full(logparams2, rdist(locs2), y, pseq2) - expect_equal(LL.full, LL.full2, tolerance = 1e-3) + set.seed(1234) + + y.orig <- rnorm(125) + + locs.dim1 <- seq(1, 5, length = 5) + rnorm(5, 0, 0.1) + locs.dim2 <- seq(1, 5, length = 5) + rnorm(5, 0, 0.1) + locs.dim3 <- seq(1, 5, length = 5) + rnorm(5, 0, 0.1) + locs <- as.matrix(expand.grid(locs.dim1, locs.dim2, locs.dim3)) + + covmat.true <- Matern(rdist(locs), range = 1, smoothness = 0.5) + diag(125) + + y <- y.orig %*% chol(covmat.true) + y <- as.vector(y) + + logparams <- create.initial.values.flex(0.9, c(2, 3), 0.5, 0.1, 1, 1) + pseq <- create.param.sequence(1, 2) + LL.full <- negloglik_full_ST( + logparams, locs, y, pseq, + c(1, 1, 2), 2 + ) + expect_equal(286.84, LL.full, tolerance = 1e-2) + + logparams2 <- create.initial.values.flex(0.9, 1, 0.5, 0.1, 1, 1) + pseq2 <- create.param.sequence(1) + locs2 <- locs + locs2[, 1:2] <- locs[, 1:2] / 2 + locs2[, 3] <- locs[, 3] / 3 + LL.full2 <- negloglik.full(logparams2, rdist(locs2), y, pseq2) + expect_equal(LL.full, LL.full2, tolerance = 1e-3) }) test_that("mvnegloglik", { - load("sim_multivariate_big.RData") - set.seed(1234) - P <- 3 - Y <- cbind(runif(10), runif(10), runif(10)) - cor.matrix <- cor(Y) - cov_mat <- c(cor.matrix[upper.tri(cor.matrix)]) - logparams <- create.initial.values.flex( - rep(0.9, P), # marginal variance - rep(0.5, P), # range - rep(0.5, P), # smoothness - rep(0.1, P), # nuggets - cov_mat, - P - ) - pseq <- create.param.sequence(P) - vec.approx <- vecchia_Mspecify(locs.list, 25) - neg_likelihood <- mvnegloglik( - logparams, vec.approx, - unlist(y.list), pseq, P - ) - expect_equal(34384.23, neg_likelihood, tolerance = 1e-2) + load("sim_multivariate_big.RData") + set.seed(1234) + P <- 3 + Y <- cbind(runif(10), runif(10), runif(10)) + cor.matrix <- cor(Y) + cov_mat <- c(cor.matrix[upper.tri(cor.matrix)]) + logparams <- create.initial.values.flex( + rep(0.9, P), # marginal variance + rep(0.5, P), # range + rep(0.5, P), # smoothness + rep(0.1, P), # nuggets + cov_mat, + P + ) + pseq <- create.param.sequence(P) + vec.approx <- vecchia_Mspecify(locs.list, 25) + neg_likelihood <- mvnegloglik( + logparams, vec.approx, + unlist(y.list), pseq, P + ) + expect_equal(34384.23, neg_likelihood, tolerance = 1e-2) }) test_that("mvnegloglik_ST", { - source("sim_multivariate_big_st.R") - P <- 3 - Y <- cbind(runif(10), runif(10), runif(10)) - cor.matrix <- cor(Y) - cov_mat <- c(cor.matrix[upper.tri(cor.matrix)]) - logparams <- create.initial.values.flex( - rep(0.9, P), # marginal variance - rep(c(2, 3), P), # range - rep(0.5, P), # smoothness - rep(0.1, P), # nuggets - cov_mat, - P - ) - pseq <- create.param.sequence(P, 2) - vec.approx <- vecchia_Mspecify(locs.list, 25) - neg_likelihood <- mvnegloglik_ST( - logparams, vec.approx, - unlist(y.list), pseq, P, c(1, 1, 2), 2 - ) - expect_equal(35106.73, neg_likelihood, tolerance = 1e-2) - - vec.approx2 <- vec.approx - for (i in 1:P) { - vec.approx2$locsord[, 1:2] <- vec.approx$locsord[, 1:2] / 2 - vec.approx2$locsord[, 3] <- vec.approx$locsord[, 3] / 3 - } - logparams2 <- logparams - logparams2[pseq[2, 1]:pseq[2, 2]] <- 0 - neg_likelihood2 <- mvnegloglik_ST( - logparams2, vec.approx2, - unlist(y.list), pseq, P, c(1, 1, 2), 2 - ) - expect_equal(neg_likelihood, neg_likelihood2, tolerance = 1e-3) - - logparams <- create.initial.values.flex( - rep(0.9, P), # marginal variance - 2:7, # range - rep(0.5, P), # smoothness - rep(0.1, P), # nuggets - cov_mat, - P - ) - neg_likelihood <- mvnegloglik_ST( - logparams, vec.approx, - unlist(y.list), pseq, P, c(1, 1, 2), 2 - ) - expect_equal(36354.9, neg_likelihood, tolerance = 1e-2) - - vec.approx2 <- vec.approx - vec.approx2$locsord[vec.approx$ondx == 1, 1:2] <- - vec.approx$locsord[vec.approx$ondx == 1, 1:2] / 2 - vec.approx2$locsord[vec.approx$ondx == 1, 3] <- - vec.approx$locsord[vec.approx$ondx == 1, 3] / 3 - vec.approx2$locsord[vec.approx$ondx == 2, 1:2] <- - vec.approx$locsord[vec.approx$ondx == 2, 1:2] / 4 - vec.approx2$locsord[vec.approx$ondx == 2, 3] <- - vec.approx$locsord[vec.approx$ondx == 2, 3] / 5 - vec.approx2$locsord[vec.approx$ondx == 3, 1:2] <- - vec.approx$locsord[vec.approx$ondx == 3, 1:2] / 6 - vec.approx2$locsord[vec.approx$ondx == 3, 3] <- - vec.approx$locsord[vec.approx$ondx == 3, 3] / 7 - - neg_likelihood2 <- mvnegloglik_ST( - logparams2, vec.approx2, - unlist(y.list), pseq, P, c(1, 1, 2), 2 - ) - expect_equal(neg_likelihood, neg_likelihood2, tolerance = 1e-3) + source("sim_multivariate_big_st.R") + P <- 3 + Y <- cbind(runif(10), runif(10), runif(10)) + cor.matrix <- cor(Y) + cov_mat <- c(cor.matrix[upper.tri(cor.matrix)]) + logparams <- create.initial.values.flex( + rep(0.9, P), # marginal variance + rep(c(2, 3), P), # range + rep(0.5, P), # smoothness + rep(0.1, P), # nuggets + cov_mat, + P + ) + pseq <- create.param.sequence(P, 2) + vec.approx <- vecchia_Mspecify(locs.list, 25) + neg_likelihood <- mvnegloglik_ST( + logparams, vec.approx, + unlist(y.list), pseq, P, c(1, 1, 2), 2 + ) + expect_equal(35106.73, neg_likelihood, tolerance = 1e-2) + + vec.approx2 <- vec.approx + for (i in 1:P) { + vec.approx2$locsord[, 1:2] <- vec.approx$locsord[, 1:2] / 2 + vec.approx2$locsord[, 3] <- vec.approx$locsord[, 3] / 3 + } + logparams2 <- logparams + logparams2[pseq[2, 1]:pseq[2, 2]] <- 0 + neg_likelihood2 <- mvnegloglik_ST( + logparams2, vec.approx2, + unlist(y.list), pseq, P, c(1, 1, 2), 2 + ) + expect_equal(neg_likelihood, neg_likelihood2, tolerance = 1e-3) + + logparams <- create.initial.values.flex( + rep(0.9, P), # marginal variance + 2:7, # range + rep(0.5, P), # smoothness + rep(0.1, P), # nuggets + cov_mat, + P + ) + neg_likelihood <- mvnegloglik_ST( + logparams, vec.approx, + unlist(y.list), pseq, P, c(1, 1, 2), 2 + ) + expect_equal(36354.9, neg_likelihood, tolerance = 1e-2) + + vec.approx2 <- vec.approx + vec.approx2$locsord[vec.approx$ondx == 1, 1:2] <- + vec.approx$locsord[vec.approx$ondx == 1, 1:2] / 2 + vec.approx2$locsord[vec.approx$ondx == 1, 3] <- + vec.approx$locsord[vec.approx$ondx == 1, 3] / 3 + vec.approx2$locsord[vec.approx$ondx == 2, 1:2] <- + vec.approx$locsord[vec.approx$ondx == 2, 1:2] / 4 + vec.approx2$locsord[vec.approx$ondx == 2, 3] <- + vec.approx$locsord[vec.approx$ondx == 2, 3] / 5 + vec.approx2$locsord[vec.approx$ondx == 3, 1:2] <- + vec.approx$locsord[vec.approx$ondx == 3, 1:2] / 6 + vec.approx2$locsord[vec.approx$ondx == 3, 3] <- + vec.approx$locsord[vec.approx$ondx == 3, 3] / 7 + + neg_likelihood2 <- mvnegloglik_ST( + logparams2, vec.approx2, + unlist(y.list), pseq, P, c(1, 1, 2), 2 + ) + expect_equal(neg_likelihood, neg_likelihood2, tolerance = 1e-3) }) test_that("mvnegloglik.full", { - source("sim_multivariate_lik.R") - - pseq <- create.param.sequence(3) - - param.marg.var <- 0.9 * unlist(lapply(y.list, var)) - param.marg.scale <- rep(1, 3) - param.marg.smooth <- rep(0.5, 3) - param.marg.nugget <- 0.1 * unlist(lapply(y.list, var)) - param.rho <- rep(0.5, 3) - params.init <- create.initial.values.flex( - param.marg.var, - param.marg.scale, - param.marg.smooth, - param.marg.nugget, - param.rho, 3 - ) - - res.optim.NM <- optim( - par = params.init, fn = mvnegloglik.full, - locs = locs.list, y = unlist(y.list), - param.seq = pseq, - method = "Nelder-Mead", - control = list(trace = 0, maxit = 10000, reltol = 1e-4) - ) - - LL.full.mv <- mvnegloglik.full( - res.optim.NM$par, locs.list, - unlist(y.list), pseq - ) - - param.seq.begin <- pseq[, 1] - param.seq.end <- pseq[, 2] - params.init.final <- res.optim.NM$par - params.init.final.t <- c( - exp(params.init.final[param.seq.begin[1]: - param.seq.end[2]]), - gtools::inv.logit(params.init.final[ - param.seq.begin[3]: - param.seq.end[3] - ], 0, 2.5), - exp(params.init.final[param.seq.begin[4]: - param.seq.end[4]]), - tanh(params.init.final[ - param.seq.begin[5]:param.seq.end[5] - ]) - ) - - params.final.flex.test <- create.initial.values.flex( - params.init.final.t[param.seq.begin[1]:param.seq.end[1]], - params.init.final.t[param.seq.begin[2]:param.seq.end[2]], - params.init.final.t[param.seq.begin[3]:param.seq.end[3]], - params.init.final.t[param.seq.begin[4]:param.seq.end[4]], - params.init.final.t[param.seq.begin[5]:param.seq.end[5]], 3 - ) - - cov.list <- create.cov.upper.flex( - 3, - params.init.final.t[param.seq.begin[1]:param.seq.end[1]], - params.init.final.t[param.seq.begin[2]:param.seq.end[2]], - params.init.final.t[param.seq.begin[3]:param.seq.end[3]], - params.init.final.t[param.seq.begin[4]:param.seq.end[4]], - params.init.final.t[param.seq.begin[5]:param.seq.end[5]] - ) - - cov.mat <- cat.covariances( - locs.list, cov.list$variance, - cov.list$range, - cov.list$smoothness, cov.list$nugget - ) - - LL.full.calc <- -1 * mvtnorm::dmvnorm(unlist(y.list), - rep(0, length(unlist(y.list))), - cov.mat, - log = TRUE - ) - - vec.mapprox <- vecchia_Mspecify(locs.list, length(unlist(y.list)) - 1) - U.mobj <- createUMultivariate(vec.mapprox, params.init.final.t) - - LL.vecchia.mv <- -1 * GPvecchia:::vecchia_likelihood_U(unlist(y.list), U.mobj) - - expect_equal(541.31, LL.full.mv, tolerance = 1e-3) - expect_equal(LL.full.calc, LL.full.mv, tolerance = 1e-3) - # Full likelihood should equal the Vecchia likelihood - expect_equal(LL.full.mv, LL.vecchia.mv, tolerance = 1e-3) - # Check create.initial.values.flex: - expect_equal(params.init.final, params.final.flex.test, tolerance = 1e-3) - # Check create.cov.upper.flex: - cov.list.true <- readRDS("covlist.rds") - expect_equal(cov.list$variance, cov.list.true$variance, tolerance = 1e-3) - expect_equal(cov.list$range, cov.list.true$range, tolerance = 1e-3) - expect_equal(cov.list$smoothness, cov.list.true$smoothness, tolerance = 1e-3) - expect_equal(cov.list$nugget, cov.list.true$nugget, tolerance = 1e-3) - # Check cat.covariances: - cov.mat.true <- readRDS("covmat.rds") - expect_equal(cov.mat, cov.mat.true, tolerance = 1e-3) + source("sim_multivariate_lik.R") + + pseq <- create.param.sequence(3) + + param.marg.var <- 0.9 * unlist(lapply(y.list, var)) + param.marg.scale <- rep(1, 3) + param.marg.smooth <- rep(0.5, 3) + param.marg.nugget <- 0.1 * unlist(lapply(y.list, var)) + param.rho <- rep(0.5, 3) + params.init <- create.initial.values.flex( + param.marg.var, + param.marg.scale, + param.marg.smooth, + param.marg.nugget, + param.rho, 3 + ) + + res.optim.NM <- optim( + par = params.init, fn = mvnegloglik.full, + locs = locs.list, y = unlist(y.list), + param.seq = pseq, + method = "Nelder-Mead", + control = list(trace = 0, maxit = 10000, reltol = 1e-4) + ) + + LL.full.mv <- mvnegloglik.full( + res.optim.NM$par, locs.list, + unlist(y.list), pseq + ) + + param.seq.begin <- pseq[, 1] + param.seq.end <- pseq[, 2] + params.init.final <- res.optim.NM$par + params.init.final.t <- c( + exp(params.init.final[param.seq.begin[1]:param.seq.end[2]]), + gtools::inv.logit(params.init.final[ + param.seq.begin[3]:param.seq.end[3] + ], 0, 2.5), + exp(params.init.final[param.seq.begin[4]:param.seq.end[4]]), + tanh(params.init.final[ + param.seq.begin[5]:param.seq.end[5] + ]) + ) + + params.final.flex.test <- create.initial.values.flex( + params.init.final.t[param.seq.begin[1]:param.seq.end[1]], + params.init.final.t[param.seq.begin[2]:param.seq.end[2]], + params.init.final.t[param.seq.begin[3]:param.seq.end[3]], + params.init.final.t[param.seq.begin[4]:param.seq.end[4]], + params.init.final.t[param.seq.begin[5]:param.seq.end[5]], 3 + ) + + cov.list <- create.cov.upper.flex( + 3, + params.init.final.t[param.seq.begin[1]:param.seq.end[1]], + params.init.final.t[param.seq.begin[2]:param.seq.end[2]], + params.init.final.t[param.seq.begin[3]:param.seq.end[3]], + params.init.final.t[param.seq.begin[4]:param.seq.end[4]], + params.init.final.t[param.seq.begin[5]:param.seq.end[5]] + ) + + cov.mat <- cat.covariances( + locs.list, cov.list$variance, + cov.list$range, + cov.list$smoothness, cov.list$nugget + ) + + LL.full.calc <- -1 * mvtnorm::dmvnorm(unlist(y.list), + rep(0, length(unlist(y.list))), + cov.mat, + log = TRUE + ) + + vec.mapprox <- vecchia_Mspecify(locs.list, length(unlist(y.list)) - 1) + U.mobj <- createUMultivariate(vec.mapprox, params.init.final.t) + + LL.vecchia.mv <- -1 * GPvecchia:::vecchia_likelihood_U(unlist(y.list), U.mobj) + + expect_equal(541.31, LL.full.mv, tolerance = 1e-3) + expect_equal(LL.full.calc, LL.full.mv, tolerance = 1e-3) + # Full likelihood should equal the Vecchia likelihood + expect_equal(LL.full.mv, LL.vecchia.mv, tolerance = 1e-3) + # Check create.initial.values.flex: + expect_equal(params.init.final, params.final.flex.test, tolerance = 1e-3) + # Check create.cov.upper.flex: + cov.list.true <- readRDS("covlist.rds") + expect_equal(cov.list$variance, cov.list.true$variance, tolerance = 1e-3) + expect_equal(cov.list$range, cov.list.true$range, tolerance = 1e-3) + expect_equal(cov.list$smoothness, cov.list.true$smoothness, tolerance = 1e-3) + expect_equal(cov.list$nugget, cov.list.true$nugget, tolerance = 1e-3) + # Check cat.covariances: + cov.mat.true <- readRDS("covmat.rds") + expect_equal(cov.mat, cov.mat.true, tolerance = 1e-3) }) diff --git a/tests/testthat/test-PrestoGP_CreateU_Multivariate.R b/tests/testthat/test-PrestoGP_CreateU_Multivariate.R index c50fc24..9c98a09 100755 --- a/tests/testthat/test-PrestoGP_CreateU_Multivariate.R +++ b/tests/testthat/test-PrestoGP_CreateU_Multivariate.R @@ -1,148 +1,148 @@ test_that("create.param.sequence", { - seq <- create.param.sequence(1) - colnames(seq) <- NULL - expect_equal(2, ncol(seq)) - expect_equal(5, nrow(seq)) - expect_equal(c(1, 1), seq[1, ]) - expect_equal(c(2, 2), seq[2, ]) - expect_equal(c(3, 3), seq[3, ]) - expect_equal(c(4, 4), seq[4, ]) - expect_equal(c(5, 5), seq[5, ]) - - seq <- create.param.sequence(3) - colnames(seq) <- NULL - expect_equal(2, ncol(seq)) - expect_equal(5, nrow(seq)) - expect_equal(c(1, 3), seq[1, ]) - expect_equal(c(4, 6), seq[2, ]) - expect_equal(c(7, 9), seq[3, ]) - expect_equal(c(10, 12), seq[4, ]) - expect_equal(c(13, 15), seq[5, ]) - - seq <- create.param.sequence(3, 2) - colnames(seq) <- NULL - expect_equal(2, ncol(seq)) - expect_equal(5, nrow(seq)) - expect_equal(c(1, 3), seq[1, ]) - expect_equal(c(4, 9), seq[2, ]) - expect_equal(c(10, 12), seq[3, ]) - expect_equal(c(13, 15), seq[4, ]) - expect_equal(c(16, 18), seq[5, ]) + seq <- create.param.sequence(1) + colnames(seq) <- NULL + expect_equal(2, ncol(seq)) + expect_equal(5, nrow(seq)) + expect_equal(c(1, 1), seq[1, ]) + expect_equal(c(2, 2), seq[2, ]) + expect_equal(c(3, 3), seq[3, ]) + expect_equal(c(4, 4), seq[4, ]) + expect_equal(c(5, 5), seq[5, ]) + + seq <- create.param.sequence(3) + colnames(seq) <- NULL + expect_equal(2, ncol(seq)) + expect_equal(5, nrow(seq)) + expect_equal(c(1, 3), seq[1, ]) + expect_equal(c(4, 6), seq[2, ]) + expect_equal(c(7, 9), seq[3, ]) + expect_equal(c(10, 12), seq[4, ]) + expect_equal(c(13, 15), seq[5, ]) + + seq <- create.param.sequence(3, 2) + colnames(seq) <- NULL + expect_equal(2, ncol(seq)) + expect_equal(5, nrow(seq)) + expect_equal(c(1, 3), seq[1, ]) + expect_equal(c(4, 9), seq[2, ]) + expect_equal(c(10, 12), seq[3, ]) + expect_equal(c(13, 15), seq[4, ]) + expect_equal(c(16, 18), seq[5, ]) }) test_that("max_min_ordering", { - set.seed(7919) - load("multivariate_sim_spatial3.Rdata") - order <- max_min_ordering(locs_train, fields::rdist) - exp_order <- c( - 63, 21, 50, 30, 76, 78, 40, 36, 23, 69, 9, 67, 32, 2, 20, 62, - 22, 31, 74, 39, 35, 58, 68, 54, 41, 3, 80, 46, 6, 88, 12, 47, - 72, 42, 13, 83, 25, 52, 11, 60, 24, 1, 28, 84, 29, 64, 66, 81, - 82, 55, 61, 87, 17, 33, 43, 45, 10, 79, 53, 75, 89, 51, 73, 27, - 26, 77, 44, 38, 65, 16, 19, 37, 57, 70, 15, 4, 5, 86, 14, 49, - 85, 34, 48, 59, 18, 8, 71, 7, 56, 90 - ) - order.gpv <- order_maxmin_exact(locs_train) - expect_equal(exp_order, order) - expect_equal(order, order.gpv) + set.seed(7919) + load("multivariate_sim_spatial3.Rdata") + order <- max_min_ordering(locs_train, fields::rdist) + exp_order <- c( + 63, 21, 50, 30, 76, 78, 40, 36, 23, 69, 9, 67, 32, 2, 20, 62, + 22, 31, 74, 39, 35, 58, 68, 54, 41, 3, 80, 46, 6, 88, 12, 47, + 72, 42, 13, 83, 25, 52, 11, 60, 24, 1, 28, 84, 29, 64, 66, 81, + 82, 55, 61, 87, 17, 33, 43, 45, 10, 79, 53, 75, 89, 51, 73, 27, + 26, 77, 44, 38, 65, 16, 19, 37, 57, 70, 15, 4, 5, 86, 14, 49, + 85, 34, 48, 59, 18, 8, 71, 7, 56, 90 + ) + order.gpv <- order_maxmin_exact(locs_train) + expect_equal(exp_order, order) + expect_equal(order, order.gpv) }) test_that("sparseNN", { - set.seed(1212) - - locs <- matrix(nrow = 100, ncol = 2) - locs[1, ] <- rep(0, 2) - for (i in 2:nrow(locs)) { - cur.r <- rnorm(1, 5) - cur.t <- runif(1, 0, 2 * pi) - locs[i, ] <- locs[i - 1, ] + c(cur.r * cos(cur.t), cur.r * sin(cur.t)) + set.seed(1212) + + locs <- matrix(nrow = 100, ncol = 2) + locs[1, ] <- rep(0, 2) + for (i in 2:nrow(locs)) { + cur.r <- rnorm(1, 5) + cur.t <- runif(1, 0, 2 * pi) + locs[i, ] <- locs[i - 1, ] + c(cur.r * cos(cur.t), cur.r * sin(cur.t)) + } + + mm.order <- order_maxmin_exact(locs) + olocs <- locs[mm.order, ] + pgp.nn <- sparseNN(olocs, 5, fields::rdist, "rdist") + gpv.nn <- GpGp:::find_ordered_nn(olocs, 5) + + indices <- matrix(nrow = nrow(olocs), ncol = 5) + distances <- indices + for (i in seq_len(nrow(olocs))) { + if (i <= 5) { + cur.dist <- fields::rdist( + olocs[(1:(5 + 1)), ][-i, ], + olocs[i, , drop = FALSE] + ) + indices[i, ] <- order(cur.dist) + } else { + cur.dist <- fields::rdist(olocs[(1:(i - 1)), ], olocs[i, , drop = FALSE]) + indices[i, ] <- order(cur.dist)[1:5] } - - mm.order <- order_maxmin_exact(locs) - olocs <- locs[mm.order, ] - pgp.nn <- sparseNN(olocs, 5, fields::rdist, "rdist") - gpv.nn <- GpGp:::find_ordered_nn(olocs, 5) - - indices <- matrix(nrow = nrow(olocs), ncol = 5) - distances <- indices - for (i in 1:nrow(olocs)) { - if (i <= 5) { - cur.dist <- fields::rdist( - olocs[(1:(5 + 1)), ][-i, ], - olocs[i, , drop = FALSE] - ) - indices[i, ] <- order(cur.dist) - } else { - cur.dist <- fields::rdist(olocs[(1:(i - 1)), ], olocs[i, , drop = FALSE]) - indices[i, ] <- order(cur.dist)[1:5] - } - distances[i, ] <- cur.dist[indices[i, ]] - } - - # Should produce the same nearest neighbors as GPvecchia - expect_equal(pgp.nn$indices[-(1:5), ], gpv.nn[-(1:5), -1]) - # Should obtain the same nearest neighbors and distances when we calculate - # the neighbors by brute force. - expect_equal(pgp.nn$indices[-(1:5), ], indices[-(1:5), ]) - expect_equal(pgp.nn$distances[-(1:5), ], distances[-(1:5), ], tolerance = 1e-2) + distances[i, ] <- cur.dist[indices[i, ]] + } + + # Should produce the same nearest neighbors as GPvecchia + expect_equal(pgp.nn$indices[-(1:5), ], gpv.nn[-(1:5), -1]) + # Should obtain the same nearest neighbors and distances when we calculate + # the neighbors by brute force. + expect_equal(pgp.nn$indices[-(1:5), ], indices[-(1:5), ]) + expect_equal(pgp.nn$distances[-(1:5), ], distances[-(1:5), ], tolerance = 1e-2) }) test_that("createUMultivariate", { - set.seed(1212) - - dim1 <- (0:9)^2 + rnorm(10, 0, 1e-2) - dim2 <- (1:10)^2 + rnorm(10, 0, 1e-2) - - locs <- as.matrix(expand.grid(dim1, dim2)) - - params <- c(3, 1.5, 0.6, 2) - - vec.approx <- vecchia_specify(locs, m = 5) - U.obj <- createU(vec.approx, covparms = params[1:3], nuggets = params[4]) - - vec.mapprox <- vecchia_Mspecify(list(locs), 5) - U.mobj <- createUMultivariate(vec.mapprox, c(params, 1)) - # Should produce the same ordered locs as GPvecchia in the univariate case - expect_equal(vec.approx$locsord, vec.mapprox$locsord) - expect_equal(vec.approx$ord, vec.mapprox$ord) - # Should produce the same U matrix as GPvecchia in the univariate case - expect_equal(sum(abs(U.obj$U - U.mobj$U)), 0, tolerance = 1e-3) - - # Should identify the same nearest neighbors and latents as GPvecchia in the - # univariate case - NNarray <- vec.approx$U.prep$revNNarray[, -6] - cond <- vec.approx$U.prep$revCond[, -6] - - for (i in 6:100) { - cur.y <- NNarray[i, cond[i, ]] - expect_equal(sort(cur.y), sort(vec.mapprox$q.list$q.y[[i]])) - if (sum(!cond[i, ]) > 0) { - cur.z <- NNarray[i, !cond[i, ]] - expect_equal(sort(cur.z), sort(vec.mapprox$q.list$q.z[[i]])) - } - } + set.seed(1212) + + dim1 <- (0:9)^2 + rnorm(10, 0, 1e-2) + dim2 <- (1:10)^2 + rnorm(10, 0, 1e-2) + + locs <- as.matrix(expand.grid(dim1, dim2)) - source("sim_multivariate.R") - - vec.mapprox <- vecchia_Mspecify(locs.list, 25) - U.mobj <- createUMultivariate(vec.mapprox, c( - marg.var, ranges, - marg.smoothness, - nuggets, rho.vec - )) - Umat.tst <- readRDS("Umat.rds") - expect_equal(sum(abs(U.mobj$U - Umat.tst)), 0, tolerance = 1e-4) - - vec.mapprox.full <- vecchia_Mspecify(locs.list, 374) - U.mobj.full <- createUMultivariate(vec.mapprox.full, c( - marg.var, ranges, - marg.smoothness, - nuggets, rho.vec - )) - - Sigma.hat <- solve(U.mobj.full$U %*% t(U.mobj.full$U)) - Sigma.hat <- Sigma.hat[U.mobj.full$latent, U.mobj.full$latent] - Sigma.hat <- Sigma.hat[order(U.mobj.full$ord), order(U.mobj.full$ord)] - # Vecchia approximation of covariance should equal true covariance when m=n-1 - expect_equal(sum(abs(Sigma.All - Sigma.hat)), 0, tolerance = 1e-4) + params <- c(3, 1.5, 0.6, 2) + + vec.approx <- vecchia_specify(locs, m = 5) + U.obj <- createU(vec.approx, covparms = params[1:3], nuggets = params[4]) + + vec.mapprox <- vecchia_Mspecify(list(locs), 5) + U.mobj <- createUMultivariate(vec.mapprox, c(params, 1)) + # Should produce the same ordered locs as GPvecchia in the univariate case + expect_equal(vec.approx$locsord, vec.mapprox$locsord) + expect_equal(vec.approx$ord, vec.mapprox$ord) + # Should produce the same U matrix as GPvecchia in the univariate case + expect_equal(sum(abs(U.obj$U - U.mobj$U)), 0, tolerance = 1e-3) + + # Should identify the same nearest neighbors and latents as GPvecchia in the + # univariate case + NNarray <- vec.approx$U.prep$revNNarray[, -6] + cond <- vec.approx$U.prep$revCond[, -6] + + for (i in 6:100) { + cur.y <- NNarray[i, cond[i, ]] + expect_equal(sort(cur.y), sort(vec.mapprox$q.list$q.y[[i]])) + if (sum(!cond[i, ]) > 0) { + cur.z <- NNarray[i, !cond[i, ]] + expect_equal(sort(cur.z), sort(vec.mapprox$q.list$q.z[[i]])) + } + } + + source("sim_multivariate.R") + + vec.mapprox <- vecchia_Mspecify(locs.list, 25) + U.mobj <- createUMultivariate(vec.mapprox, c( + marg.var, ranges, + marg.smoothness, + nuggets, rho.vec + )) + Umat.tst <- readRDS("Umat.rds") + expect_equal(sum(abs(U.mobj$U - Umat.tst)), 0, tolerance = 1e-4) + + vec.mapprox.full <- vecchia_Mspecify(locs.list, 374) + U.mobj.full <- createUMultivariate(vec.mapprox.full, c( + marg.var, ranges, + marg.smoothness, + nuggets, rho.vec + )) + + Sigma.hat <- solve(U.mobj.full$U %*% t(U.mobj.full$U)) + Sigma.hat <- Sigma.hat[U.mobj.full$latent, U.mobj.full$latent] + Sigma.hat <- Sigma.hat[order(U.mobj.full$ord), order(U.mobj.full$ord)] + # Vecchia approximation of covariance should equal true covariance when m=n-1 + expect_equal(sum(abs(Sigma.All - Sigma.hat)), 0, tolerance = 1e-4) })