Title: | Bayesian Predictive Stacking for Scalable Geospatial Transfer Learning |
---|---|
Description: | Provides functions for Bayesian Predictive Stacking within the Bayesian transfer learning framework for geospatial artificial systems, as introduced in "Bayesian Transfer Learning for Artificially Intelligent Geospatial Systems: A Predictive Stacking Approach" (Presicce and Banerjee, 2024) <doi:10.48550/arXiv.2410.09504>. This methodology enables efficient Bayesian geostatistical modeling, utilizing predictive stacking to improve inference across spatial datasets. The core functions leverage 'C++' for high-performance computation, making the framework well-suited for large-scale spatial data analysis in parallel and distributed computing environments. Designed for scalability, it allows seamless application in computationally demanding scenarios. |
Authors: | Luca Presicce [aut, cre] , Sudipto Banerjee [aut] |
Maintainer: | Luca Presicce <[email protected]> |
License: | GPL (>= 3) |
Version: | 0.0-4 |
Built: | 2024-11-25 06:29:45 UTC |
Source: | CRAN |
Compute the Euclidean distance matrix
arma_dist(X)
arma_dist(X)
X |
matrix (tipically of |
matrix distance matrix of the elements of
## Compute the Distance matrix of dimension (n x n) n <- 100 p <- 2 X <- matrix(runif(n*p), nrow = n, ncol = p) distance.matrix <- arma_dist(X)
## Compute the Distance matrix of dimension (n x n) n <- 100 p <- 2 X <- matrix(runif(n*p), nrow = n, ncol = p) distance.matrix <- arma_dist(X)
Gibbs sampler for Conjugate Bayesian Multivariate Linear Models
bayesMvLMconjugate(Y, X, mu_B, V_B, nu, Psi, n_iter = 1000, burn_in = 500)
bayesMvLMconjugate(Y, X, mu_B, V_B, nu, Psi, n_iter = 1000, burn_in = 500)
Y |
matrix |
X |
matrix |
mu_B |
matrix |
V_B |
matrix |
nu |
double prior parameter for |
Psi |
matrix prior parameter for |
n_iter |
integer iteration number for Gibbs sampler |
burn_in |
integer number of burn-in iteration |
B_samples array of posterior sample for
Sigma_samples array of posterior samples for
## Generate data n <- 100 p <- 3 q <- 2 Y <- matrix(rnorm(n*q), nrow = n, ncol = q) X <- matrix(rnorm(n*p), nrow = n, ncol = p) ## Prior parameters mu_B <- matrix(0, p, q) V_B <- diag(10, p) nu <- 3 Psi <- diag(q) ## Samples from posteriors n_iter <- 1000 burn_in <- 500 set.seed(1234) samples <- spBPS::bayesMvLMconjugate(Y, X, mu_B, V_B, nu, Psi, n_iter, burn_in)
## Generate data n <- 100 p <- 3 q <- 2 Y <- matrix(rnorm(n*q), nrow = n, ncol = q) X <- matrix(rnorm(n*p), nrow = n, ncol = p) ## Prior parameters mu_B <- matrix(0, p, q) V_B <- diag(10, p) nu <- 3 Psi <- diag(q) ## Samples from posteriors n_iter <- 1000 burn_in <- 500 set.seed(1234) samples <- spBPS::bayesMvLMconjugate(Y, X, mu_B, V_B, nu, Psi, n_iter, burn_in)
Combine subset models wiht BPS
BPS_combine(fit_list, K, rp)
BPS_combine(fit_list, K, rp)
fit_list |
list K fitted model outputs composed by two elements each: first named |
K |
integer number of folds |
rp |
double percentage of observations to take into account for optimization ( |
matrix posterior predictive density evaluations (each columns represent a different model)
## Generate subsets of data n <- 100 p <- 3 X <- matrix(rnorm(n*p), nrow = n, ncol = p) Y <- matrix(rnorm(n), nrow = n, ncol = 1) crd <- matrix(runif(n*2), nrow = n, ncol = 2) data_part <- subset_data(data = list(Y = Y, X = X, crd = crd), K = 10) ## Select competitive set of values for hyperparameters delta_seq <- c(0.1, 0.2, 0.3) phi_seq <- c(3, 4, 5) ## Perform Bayesian Predictive Stacking within subsets fit_list <- vector(length = 10, mode = "list") for (i in 1:10) { Yi <- data_part$Y_list[[i]] Xi <- data_part$X_list[[i]] crd_i <- data_part$crd_list[[i]] p <- ncol(Xi) bps <- spBPS::BPS_weights(data = list(Y = Yi, X = Xi), priors = list(mu_b = matrix(rep(0, p)), V_b = diag(10, p), a = 2, b = 2), coords = crd_i, hyperpar = list(delta = delta_seq, phi = phi_seq), K = 5) w_hat <- bps$W epd <- bps$epd fit_list[[i]] <- list(epd, w_hat) } ## Combination weights between partitions using Bayesian Predictive Stacking comb_bps <- BPS_combine(fit_list = fit_list, K = 10, rp = 1)
## Generate subsets of data n <- 100 p <- 3 X <- matrix(rnorm(n*p), nrow = n, ncol = p) Y <- matrix(rnorm(n), nrow = n, ncol = 1) crd <- matrix(runif(n*2), nrow = n, ncol = 2) data_part <- subset_data(data = list(Y = Y, X = X, crd = crd), K = 10) ## Select competitive set of values for hyperparameters delta_seq <- c(0.1, 0.2, 0.3) phi_seq <- c(3, 4, 5) ## Perform Bayesian Predictive Stacking within subsets fit_list <- vector(length = 10, mode = "list") for (i in 1:10) { Yi <- data_part$Y_list[[i]] Xi <- data_part$X_list[[i]] crd_i <- data_part$crd_list[[i]] p <- ncol(Xi) bps <- spBPS::BPS_weights(data = list(Y = Yi, X = Xi), priors = list(mu_b = matrix(rep(0, p)), V_b = diag(10, p), a = 2, b = 2), coords = crd_i, hyperpar = list(delta = delta_seq, phi = phi_seq), K = 5) w_hat <- bps$W epd <- bps$epd fit_list[[i]] <- list(epd, w_hat) } ## Combination weights between partitions using Bayesian Predictive Stacking comb_bps <- BPS_combine(fit_list = fit_list, K = 10, rp = 1)
Perform the BPS sampling from posterior and posterior predictive given a set of stacking weights
BPS_post(data, X_u, priors, coords, crd_u, hyperpar, W, R)
BPS_post(data, X_u, priors, coords, crd_u, hyperpar, W, R)
data |
list two elements: first named |
X_u |
matrix unobserved instances covariate matrix |
priors |
list priors: named |
coords |
matrix sample coordinates for X and Y |
crd_u |
matrix unboserved instances coordinates |
hyperpar |
list two elemets: first named |
W |
matrix set of stacking weights |
R |
integer number of desired samples |
list BPS posterior predictive samples
## Generate subsets of data n <- 100 p <- 3 X <- matrix(rnorm(n*p), nrow = n, ncol = p) Y <- matrix(rnorm(n), nrow = n, ncol = 1) crd <- matrix(runif(n*2), nrow = n, ncol = 2) data_part <- subset_data(data = list(Y = Y, X = X, crd = crd), K = 10) ## Select competetive set of values for hyperparameters delta_seq <- c(0.1, 0.2, 0.3) phi_seq <- c(3, 4, 5) ## Fit local models fit_list <- vector(length = 10, mode = "list") for (i in 1:10) { Yi <- data_part$Y_list[[i]] Xi <- data_part$X_list[[i]] crd_i <- data_part$crd_list[[i]] p <- ncol(Xi) bps <- spBPS::BPS_weights(data = list(Y = Yi, X = Xi), priors = list(mu_b = matrix(rep(0, p)), V_b = diag(10, p), a = 2, b = 2), coords = crd_i, hyperpar = list(delta = delta_seq, phi = phi_seq), K = 5) w_hat <- bps$W epd <- bps$epd fit_list[[i]] <- list(epd, w_hat) } ## Model combination weights between partitions using Bayesian Predictive Stacking comb_bps <- BPS_combine(fit_list = fit_list, K = 10, rp = 1) Wbps <- comb_bps$W W_list <- comb_bps$W_list ## Generate prediction points m <- 100 X_new <- matrix(rnorm(m*p), nrow = m, ncol = p) crd_new <- matrix(runif(m*2), nrow = m, ncol = 2) ## Perform posterior and posterior predictive sampling R <- 250 subset_ind <- sample(1:10, R, TRUE, Wbps) postsmp_and_pred <- vector(length = R, mode = "list") for (r in 1:R) { ind_s <- subset_ind[r] Ys <- matrix(data_part$Y_list[[ind_s]]) Xs <- data_part$X_list[[ind_s]] crds <- data_part$crd_list[[ind_s]] Ws <- W_list[[ind_s]] result <- spBPS::BPS_post(data = list(Y = Ys, X = Xs), coords = crds, X_u = X_new, crd_u = crd_new, priors = list(mu_b = matrix(rep(0, p)), V_b = diag(10, p), a = 2, b = 2), hyperpar = list(delta = delta_seq, phi = phi_seq), W = Ws, R = 1) postsmp_and_pred[[r]] <- result}
## Generate subsets of data n <- 100 p <- 3 X <- matrix(rnorm(n*p), nrow = n, ncol = p) Y <- matrix(rnorm(n), nrow = n, ncol = 1) crd <- matrix(runif(n*2), nrow = n, ncol = 2) data_part <- subset_data(data = list(Y = Y, X = X, crd = crd), K = 10) ## Select competetive set of values for hyperparameters delta_seq <- c(0.1, 0.2, 0.3) phi_seq <- c(3, 4, 5) ## Fit local models fit_list <- vector(length = 10, mode = "list") for (i in 1:10) { Yi <- data_part$Y_list[[i]] Xi <- data_part$X_list[[i]] crd_i <- data_part$crd_list[[i]] p <- ncol(Xi) bps <- spBPS::BPS_weights(data = list(Y = Yi, X = Xi), priors = list(mu_b = matrix(rep(0, p)), V_b = diag(10, p), a = 2, b = 2), coords = crd_i, hyperpar = list(delta = delta_seq, phi = phi_seq), K = 5) w_hat <- bps$W epd <- bps$epd fit_list[[i]] <- list(epd, w_hat) } ## Model combination weights between partitions using Bayesian Predictive Stacking comb_bps <- BPS_combine(fit_list = fit_list, K = 10, rp = 1) Wbps <- comb_bps$W W_list <- comb_bps$W_list ## Generate prediction points m <- 100 X_new <- matrix(rnorm(m*p), nrow = m, ncol = p) crd_new <- matrix(runif(m*2), nrow = m, ncol = 2) ## Perform posterior and posterior predictive sampling R <- 250 subset_ind <- sample(1:10, R, TRUE, Wbps) postsmp_and_pred <- vector(length = R, mode = "list") for (r in 1:R) { ind_s <- subset_ind[r] Ys <- matrix(data_part$Y_list[[ind_s]]) Xs <- data_part$X_list[[ind_s]] crds <- data_part$crd_list[[ind_s]] Ws <- W_list[[ind_s]] result <- spBPS::BPS_post(data = list(Y = Ys, X = Xs), coords = crds, X_u = X_new, crd_u = crd_new, priors = list(mu_b = matrix(rep(0, p)), V_b = diag(10, p), a = 2, b = 2), hyperpar = list(delta = delta_seq, phi = phi_seq), W = Ws, R = 1) postsmp_and_pred[[r]] <- result}
Perform the BPS sampling from posterior and posterior predictive given a set of stacking weights
BPS_post_MvT(data, X_u, priors, coords, crd_u, hyperpar, W, R)
BPS_post_MvT(data, X_u, priors, coords, crd_u, hyperpar, W, R)
data |
list two elements: first named |
X_u |
matrix unobserved instances covariate matrix |
priors |
list priors: named |
coords |
matrix sample coordinates for X and Y |
crd_u |
matrix unboserved instances coordinates |
hyperpar |
list two elemets: first named |
W |
matrix set of stacking weights |
R |
integer number of desired samples |
list BPS posterior predictive samples
## Generate subsets of data n <- 100 p <- 3 q <- 2 X <- matrix(rnorm(n*p), nrow = n, ncol = p) Y <- matrix(rnorm(n*q), nrow = n, ncol = q) crd <- matrix(runif(n*2), nrow = n, ncol = 2) data_part <- subset_data(data = list(Y = Y, X = X, crd = crd), K = 10) ## Select competitive set of values for hyperparameters alfa_seq <- c(0.7, 0.8, 0.9) phi_seq <- c(3, 4, 5) ## Fit local models fit_list <- vector(length = 10, mode = "list") for (i in 1:10) { Yi <- data_part$Y_list[[i]] Xi <- data_part$X_list[[i]] crd_i <- data_part$crd_list[[i]] bps <- spBPS::BPS_weights_MvT(data = list(Y = Yi, X = Xi), priors = list(mu_B = matrix(0, nrow = p, ncol = q), V_r = diag(10, p), Psi = diag(1, q), nu = 3), coords = crd_i, hyperpar = list(alpha = alfa_seq, phi = phi_seq), K = 5) w_hat <- bps$W epd <- bps$epd fit_list[[i]] <- list(epd, w_hat) } ## Model combination weights between partitions using Bayesian Predictive Stacking comb_bps <- BPS_combine(fit_list = fit_list, K = 10, rp = 1) Wbps <- comb_bps$W W_list <- comb_bps$W_list ## Generate prediction points m <- 100 X_new <- matrix(rnorm(m*p), nrow = m, ncol = p) crd_new <- matrix(runif(m*2), nrow = m, ncol = 2) ## Perform posterior and posterior predictive sampling R <- 250 subset_ind <- sample(1:10, R, TRUE, Wbps) postsmp_and_pred <- vector(length = R, mode = "list") for (r in 1:R) { ind_s <- subset_ind[r] Ys <- data_part$Y_list[[ind_s]] Xs <- data_part$X_list[[ind_s]] crds <- data_part$crd_list[[ind_s]] Ws <- W_list[[ind_s]] result <- spBPS::BPS_post_MvT(data = list(Y = Ys, X = Xs), coords = crds, X_u = X_new, crd_u = crd_new, priors = list(mu_B = matrix(0, nrow = p, ncol = q), V_r = diag(10, p), Psi = diag(1, q), nu = 3), hyperpar = list(alpha = alfa_seq, phi = phi_seq), W = Ws, R = 1) postsmp_and_pred[[r]] <- result}
## Generate subsets of data n <- 100 p <- 3 q <- 2 X <- matrix(rnorm(n*p), nrow = n, ncol = p) Y <- matrix(rnorm(n*q), nrow = n, ncol = q) crd <- matrix(runif(n*2), nrow = n, ncol = 2) data_part <- subset_data(data = list(Y = Y, X = X, crd = crd), K = 10) ## Select competitive set of values for hyperparameters alfa_seq <- c(0.7, 0.8, 0.9) phi_seq <- c(3, 4, 5) ## Fit local models fit_list <- vector(length = 10, mode = "list") for (i in 1:10) { Yi <- data_part$Y_list[[i]] Xi <- data_part$X_list[[i]] crd_i <- data_part$crd_list[[i]] bps <- spBPS::BPS_weights_MvT(data = list(Y = Yi, X = Xi), priors = list(mu_B = matrix(0, nrow = p, ncol = q), V_r = diag(10, p), Psi = diag(1, q), nu = 3), coords = crd_i, hyperpar = list(alpha = alfa_seq, phi = phi_seq), K = 5) w_hat <- bps$W epd <- bps$epd fit_list[[i]] <- list(epd, w_hat) } ## Model combination weights between partitions using Bayesian Predictive Stacking comb_bps <- BPS_combine(fit_list = fit_list, K = 10, rp = 1) Wbps <- comb_bps$W W_list <- comb_bps$W_list ## Generate prediction points m <- 100 X_new <- matrix(rnorm(m*p), nrow = m, ncol = p) crd_new <- matrix(runif(m*2), nrow = m, ncol = 2) ## Perform posterior and posterior predictive sampling R <- 250 subset_ind <- sample(1:10, R, TRUE, Wbps) postsmp_and_pred <- vector(length = R, mode = "list") for (r in 1:R) { ind_s <- subset_ind[r] Ys <- data_part$Y_list[[ind_s]] Xs <- data_part$X_list[[ind_s]] crds <- data_part$crd_list[[ind_s]] Ws <- W_list[[ind_s]] result <- spBPS::BPS_post_MvT(data = list(Y = Ys, X = Xs), coords = crds, X_u = X_new, crd_u = crd_new, priors = list(mu_B = matrix(0, nrow = p, ncol = q), V_r = diag(10, p), Psi = diag(1, q), nu = 3), hyperpar = list(alpha = alfa_seq, phi = phi_seq), W = Ws, R = 1) postsmp_and_pred[[r]] <- result}
Compute the BPS posterior samples given a set of stacking weights
BPS_postdraws(data, priors, coords, hyperpar, W, R)
BPS_postdraws(data, priors, coords, hyperpar, W, R)
data |
list two elements: first named |
priors |
list priors: named |
coords |
matrix sample coordinates for X and Y |
hyperpar |
list two elemets: first named |
W |
matrix set of stacking weights |
R |
integer number of desired samples |
matrix BPS posterior samples
Compute the BPS posterior samples given a set of stacking weights
BPS_postdraws_MvT(data, priors, coords, hyperpar, W, R, par)
BPS_postdraws_MvT(data, priors, coords, hyperpar, W, R, par)
data |
list two elements: first named |
priors |
list priors: named |
coords |
matrix sample coordinates for X and Y |
hyperpar |
list two elemets: first named |
W |
matrix set of stacking weights |
R |
integer number of desired samples |
par |
if |
matrix BPS posterior samples
Compute the BPS spatial prediction given a set of stacking weights
BPS_pred(data, X_u, priors, coords, crd_u, hyperpar, W, R)
BPS_pred(data, X_u, priors, coords, crd_u, hyperpar, W, R)
data |
list two elements: first named |
X_u |
matrix unobserved instances covariate matrix |
priors |
list priors: named |
coords |
matrix sample coordinates for X and Y |
crd_u |
matrix unboserved instances coordinates |
hyperpar |
list two elemets: first named |
W |
matrix set of stacking weights |
R |
integer number of desired samples |
list BPS posterior predictive samples
## Generate subsets of data n <- 100 p <- 3 X <- matrix(rnorm(n*p), nrow = n, ncol = p) Y <- matrix(rnorm(n), nrow = n, ncol = 1) crd <- matrix(runif(n*2), nrow = n, ncol = 2) data_part <- subset_data(data = list(Y = Y, X = X, crd = crd), K = 10) ## Select competetive set of values for hyperparameters delta_seq <- c(0.1, 0.2, 0.3) phi_seq <- c(3, 4, 5) ## Fit local models fit_list <- vector(length = 10, mode = "list") for (i in 1:10) { Yi <- data_part$Y_list[[i]] Xi <- data_part$X_list[[i]] crd_i <- data_part$crd_list[[i]] p <- ncol(Xi) bps <- spBPS::BPS_weights(data = list(Y = Yi, X = Xi), priors = list(mu_b = matrix(rep(0, p)), V_b = diag(10, p), a = 2, b = 2), coords = crd_i, hyperpar = list(delta = delta_seq, phi = phi_seq), K = 5) w_hat <- bps$W epd <- bps$epd fit_list[[i]] <- list(epd, w_hat) } ## Model combination weights between partitions using Bayesian Predictive Stacking comb_bps <- BPS_combine(fit_list = fit_list, K = 10, rp = 1) Wbps <- comb_bps$W W_list <- comb_bps$W_list ## Generate prediction points m <- 50 X_new <- matrix(rnorm(m*p), nrow = m, ncol = p) crd_new <- matrix(runif(m*2), nrow = m, ncol = 2) ## Perform posterior predictive sampling R <- 250 subset_ind <- sample(1:10, R, TRUE, Wbps) predictions <- vector(length = R, mode = "list") for (r in 1:R) { ind_s <- subset_ind[r] Ys <- matrix(data_part$Y_list[[ind_s]]) Xs <- data_part$X_list[[ind_s]] crds <- data_part$crd_list[[ind_s]] Ws <- W_list[[ind_s]] result <- spBPS::BPS_pred(data = list(Y = Ys, X = Xs), coords = crds, X_u = X_new, crd_u = crd_new, priors = list(mu_b = matrix(rep(0, p)), V_b = diag(10, p), a = 2, b = 2), hyperpar = list(delta = delta_seq, phi = phi_seq), W = Ws, R = 1) predictions[[r]] <- result}
## Generate subsets of data n <- 100 p <- 3 X <- matrix(rnorm(n*p), nrow = n, ncol = p) Y <- matrix(rnorm(n), nrow = n, ncol = 1) crd <- matrix(runif(n*2), nrow = n, ncol = 2) data_part <- subset_data(data = list(Y = Y, X = X, crd = crd), K = 10) ## Select competetive set of values for hyperparameters delta_seq <- c(0.1, 0.2, 0.3) phi_seq <- c(3, 4, 5) ## Fit local models fit_list <- vector(length = 10, mode = "list") for (i in 1:10) { Yi <- data_part$Y_list[[i]] Xi <- data_part$X_list[[i]] crd_i <- data_part$crd_list[[i]] p <- ncol(Xi) bps <- spBPS::BPS_weights(data = list(Y = Yi, X = Xi), priors = list(mu_b = matrix(rep(0, p)), V_b = diag(10, p), a = 2, b = 2), coords = crd_i, hyperpar = list(delta = delta_seq, phi = phi_seq), K = 5) w_hat <- bps$W epd <- bps$epd fit_list[[i]] <- list(epd, w_hat) } ## Model combination weights between partitions using Bayesian Predictive Stacking comb_bps <- BPS_combine(fit_list = fit_list, K = 10, rp = 1) Wbps <- comb_bps$W W_list <- comb_bps$W_list ## Generate prediction points m <- 50 X_new <- matrix(rnorm(m*p), nrow = m, ncol = p) crd_new <- matrix(runif(m*2), nrow = m, ncol = 2) ## Perform posterior predictive sampling R <- 250 subset_ind <- sample(1:10, R, TRUE, Wbps) predictions <- vector(length = R, mode = "list") for (r in 1:R) { ind_s <- subset_ind[r] Ys <- matrix(data_part$Y_list[[ind_s]]) Xs <- data_part$X_list[[ind_s]] crds <- data_part$crd_list[[ind_s]] Ws <- W_list[[ind_s]] result <- spBPS::BPS_pred(data = list(Y = Ys, X = Xs), coords = crds, X_u = X_new, crd_u = crd_new, priors = list(mu_b = matrix(rep(0, p)), V_b = diag(10, p), a = 2, b = 2), hyperpar = list(delta = delta_seq, phi = phi_seq), W = Ws, R = 1) predictions[[r]] <- result}
Compute the BPS spatial prediction given a set of stacking weights
BPS_pred_MvT(data, X_u, priors, coords, crd_u, hyperpar, W, R)
BPS_pred_MvT(data, X_u, priors, coords, crd_u, hyperpar, W, R)
data |
list two elements: first named |
X_u |
matrix unobserved instances covariate matrix |
priors |
list priors: named |
coords |
matrix sample coordinates for X and Y |
crd_u |
matrix unboserved instances coordinates |
hyperpar |
list two elemets: first named |
W |
matrix set of stacking weights |
R |
integer number of desired samples |
list BPS posterior predictive samples
## Generate subsets of data n <- 100 p <- 3 q <- 2 X <- matrix(rnorm(n*p), nrow = n, ncol = p) Y <- matrix(rnorm(n*q), nrow = n, ncol = q) crd <- matrix(runif(n*2), nrow = n, ncol = 2) data_part <- subset_data(data = list(Y = Y, X = X, crd = crd), K = 10) ## Select competitive set of values for hyperparameters alfa_seq <- c(0.7, 0.8, 0.9) phi_seq <- c(3, 4, 5) ## Fit local models fit_list <- vector(length = 10, mode = "list") for (i in 1:10) { Yi <- data_part$Y_list[[i]] Xi <- data_part$X_list[[i]] crd_i <- data_part$crd_list[[i]] bps <- spBPS::BPS_weights_MvT(data = list(Y = Yi, X = Xi), priors = list(mu_B = matrix(0, nrow = p, ncol = q), V_r = diag(10, p), Psi = diag(1, q), nu = 3), coords = crd_i, hyperpar = list(alpha = alfa_seq, phi = phi_seq), K = 5) w_hat <- bps$W epd <- bps$epd fit_list[[i]] <- list(epd, w_hat) } ## Model combination weights between partitions using Bayesian Predictive Stacking comb_bps <- BPS_combine(fit_list = fit_list, K = 10, rp = 1) Wbps <- comb_bps$W W_list <- comb_bps$W_list ## Generate prediction points m <- 100 X_new <- matrix(rnorm(m*p), nrow = m, ncol = p) crd_new <- matrix(runif(m*2), nrow = m, ncol = 2) ## Perform posterior predictive sampling R <- 250 subset_ind <- sample(1:10, R, TRUE, Wbps) predictions <- vector(length = R, mode = "list") for (r in 1:R) { ind_s <- subset_ind[r] Ys <- data_part$Y_list[[ind_s]] Xs <- data_part$X_list[[ind_s]] crds <- data_part$crd_list[[ind_s]] Ws <- W_list[[ind_s]] result <- spBPS::BPS_pred_MvT(data = list(Y = Ys, X = Xs), coords = crds, X_u = X_new, crd_u = crd_new, priors = list(mu_B = matrix(0, nrow = p, ncol = q), V_r = diag(10, p), Psi = diag(1, q), nu = 3), hyperpar = list(alpha = alfa_seq, phi = phi_seq), W = Ws, R = 1) predictions[[r]] <- result}
## Generate subsets of data n <- 100 p <- 3 q <- 2 X <- matrix(rnorm(n*p), nrow = n, ncol = p) Y <- matrix(rnorm(n*q), nrow = n, ncol = q) crd <- matrix(runif(n*2), nrow = n, ncol = 2) data_part <- subset_data(data = list(Y = Y, X = X, crd = crd), K = 10) ## Select competitive set of values for hyperparameters alfa_seq <- c(0.7, 0.8, 0.9) phi_seq <- c(3, 4, 5) ## Fit local models fit_list <- vector(length = 10, mode = "list") for (i in 1:10) { Yi <- data_part$Y_list[[i]] Xi <- data_part$X_list[[i]] crd_i <- data_part$crd_list[[i]] bps <- spBPS::BPS_weights_MvT(data = list(Y = Yi, X = Xi), priors = list(mu_B = matrix(0, nrow = p, ncol = q), V_r = diag(10, p), Psi = diag(1, q), nu = 3), coords = crd_i, hyperpar = list(alpha = alfa_seq, phi = phi_seq), K = 5) w_hat <- bps$W epd <- bps$epd fit_list[[i]] <- list(epd, w_hat) } ## Model combination weights between partitions using Bayesian Predictive Stacking comb_bps <- BPS_combine(fit_list = fit_list, K = 10, rp = 1) Wbps <- comb_bps$W W_list <- comb_bps$W_list ## Generate prediction points m <- 100 X_new <- matrix(rnorm(m*p), nrow = m, ncol = p) crd_new <- matrix(runif(m*2), nrow = m, ncol = 2) ## Perform posterior predictive sampling R <- 250 subset_ind <- sample(1:10, R, TRUE, Wbps) predictions <- vector(length = R, mode = "list") for (r in 1:R) { ind_s <- subset_ind[r] Ys <- data_part$Y_list[[ind_s]] Xs <- data_part$X_list[[ind_s]] crds <- data_part$crd_list[[ind_s]] Ws <- W_list[[ind_s]] result <- spBPS::BPS_pred_MvT(data = list(Y = Ys, X = Xs), coords = crds, X_u = X_new, crd_u = crd_new, priors = list(mu_B = matrix(0, nrow = p, ncol = q), V_r = diag(10, p), Psi = diag(1, q), nu = 3), hyperpar = list(alpha = alfa_seq, phi = phi_seq), W = Ws, R = 1) predictions[[r]] <- result}
Combine subset models wiht Pseudo-BMA
BPS_PseudoBMA(fit_list)
BPS_PseudoBMA(fit_list)
fit_list |
list K fitted model outputs composed by two elements each: first named |
matrix posterior predictive density evaluations (each columns represent a different model)
## Generate subsets of data n <- 100 p <- 3 X <- matrix(rnorm(n*p), nrow = n, ncol = p) Y <- matrix(rnorm(n), nrow = n, ncol = 1) crd <- matrix(runif(n*2), nrow = n, ncol = 2) data_part <- subset_data(data = list(Y = Y, X = X, crd = crd), K = 10) ## Select competitive set of values for hyperparameters delta_seq <- c(0.1, 0.2, 0.3) phi_seq <- c(3, 4, 5) ## Perform Bayesian Predictive Stacking within subsets fit_list <- vector(length = 10, mode = "list") for (i in 1:10) { Yi <- data_part$Y_list[[i]] Xi <- data_part$X_list[[i]] crd_i <- data_part$crd_list[[i]] p <- ncol(Xi) bps <- spBPS::BPS_weights(data = list(Y = Yi, X = Xi), priors = list(mu_b = matrix(rep(0, p)), V_b = diag(10, p), a = 2, b = 2), coords = crd_i, hyperpar = list(delta = delta_seq, phi = phi_seq), K = 5) w_hat <- bps$W epd <- bps$epd fit_list[[i]] <- list(epd, w_hat) } ## Combination weights between partitions using Pseudo Bayesian Model Averaging comb_bps <- BPS_PseudoBMA(fit_list = fit_list)
## Generate subsets of data n <- 100 p <- 3 X <- matrix(rnorm(n*p), nrow = n, ncol = p) Y <- matrix(rnorm(n), nrow = n, ncol = 1) crd <- matrix(runif(n*2), nrow = n, ncol = 2) data_part <- subset_data(data = list(Y = Y, X = X, crd = crd), K = 10) ## Select competitive set of values for hyperparameters delta_seq <- c(0.1, 0.2, 0.3) phi_seq <- c(3, 4, 5) ## Perform Bayesian Predictive Stacking within subsets fit_list <- vector(length = 10, mode = "list") for (i in 1:10) { Yi <- data_part$Y_list[[i]] Xi <- data_part$X_list[[i]] crd_i <- data_part$crd_list[[i]] p <- ncol(Xi) bps <- spBPS::BPS_weights(data = list(Y = Yi, X = Xi), priors = list(mu_b = matrix(rep(0, p)), V_b = diag(10, p), a = 2, b = 2), coords = crd_i, hyperpar = list(delta = delta_seq, phi = phi_seq), K = 5) w_hat <- bps$W epd <- bps$epd fit_list[[i]] <- list(epd, w_hat) } ## Combination weights between partitions using Pseudo Bayesian Model Averaging comb_bps <- BPS_PseudoBMA(fit_list = fit_list)
Compute the BPS weights by convex optimization
BPS_weights(data, priors, coords, hyperpar, K)
BPS_weights(data, priors, coords, hyperpar, K)
data |
list two elements: first named |
priors |
list priors: named |
coords |
matrix sample coordinates for X and Y |
hyperpar |
list two elemets: first named |
K |
integer number of folds |
matrix posterior predictive density evaluations (each columns represent a different model)
## Generate subsets of data n <- 100 p <- 3 X <- matrix(rnorm(n*p), nrow = n, ncol = p) Y <- matrix(rnorm(n), nrow = n) crd <- matrix(runif(n*2), nrow = n, ncol = 2) ## Select competitive set of values for hyperparameters delta_seq <- c(0.1, 0.2, 0.3) phi_seq <- c(3, 4, 5) ## Perform Bayesian Predictive Stacking within subsets bps <- spBPS::BPS_weights(data = list(Y = Y, X = X), priors = list(mu_b = matrix(rep(0, p)), V_b = diag(10, p), a = 2, b = 2), coords = crd, hyperpar = list(delta = delta_seq, phi = phi_seq), K = 5)
## Generate subsets of data n <- 100 p <- 3 X <- matrix(rnorm(n*p), nrow = n, ncol = p) Y <- matrix(rnorm(n), nrow = n) crd <- matrix(runif(n*2), nrow = n, ncol = 2) ## Select competitive set of values for hyperparameters delta_seq <- c(0.1, 0.2, 0.3) phi_seq <- c(3, 4, 5) ## Perform Bayesian Predictive Stacking within subsets bps <- spBPS::BPS_weights(data = list(Y = Y, X = X), priors = list(mu_b = matrix(rep(0, p)), V_b = diag(10, p), a = 2, b = 2), coords = crd, hyperpar = list(delta = delta_seq, phi = phi_seq), K = 5)
Compute the BPS weights by convex optimization
BPS_weights_MvT(data, priors, coords, hyperpar, K)
BPS_weights_MvT(data, priors, coords, hyperpar, K)
data |
list two elements: first named |
priors |
list priors: named |
coords |
matrix sample coordinates for X and Y |
hyperpar |
list two elemets: first named |
K |
integer number of folds |
matrix posterior predictive density evaluations (each columns represent a different model)
## Generate subsets of data n <- 100 p <- 3 q <- 2 X <- matrix(rnorm(n*p), nrow = n, ncol = p) Y <- matrix(rnorm(n*q), nrow = n, ncol = q) crd <- matrix(runif(n*2), nrow = n, ncol = 2) ## Select competitive set of values for hyperparameters alfa_seq <- c(0.7, 0.8, 0.9) phi_seq <- c(3, 4, 5) ## Perform Bayesian Predictive Stacking within subsets bps <- spBPS::BPS_weights_MvT(data = list(Y = Y, X = X), priors = list(mu_B = matrix(0, nrow = p, ncol = q), V_r = diag(10, p), Psi = diag(1, q), nu = 3), coords = crd, hyperpar = list(alpha = alfa_seq, phi = phi_seq), K = 5)
## Generate subsets of data n <- 100 p <- 3 q <- 2 X <- matrix(rnorm(n*p), nrow = n, ncol = p) Y <- matrix(rnorm(n*q), nrow = n, ncol = q) crd <- matrix(runif(n*2), nrow = n, ncol = 2) ## Select competitive set of values for hyperparameters alfa_seq <- c(0.7, 0.8, 0.9) phi_seq <- c(3, 4, 5) ## Perform Bayesian Predictive Stacking within subsets bps <- spBPS::BPS_weights_MvT(data = list(Y = Y, X = X), priors = list(mu_B = matrix(0, nrow = p, ncol = q), V_r = diag(10, p), Psi = diag(1, q), nu = 3), coords = crd, hyperpar = list(alpha = alfa_seq, phi = phi_seq), K = 5)
Solver for Bayesian Predictive Stacking of Predictive densities convex optimization problem
conv_opt(scores)
conv_opt(scores)
scores |
matrix |
W matrix of Bayesian Predictive Stacking weights for the K models considered
## Generate (randomly) K predictive scores for n observations n <- 50 K <- 5 scores <- matrix(runif(n*K), nrow = n, ncol = K) ## Find Bayesian Predictive Stacking weights opt_weights <- conv_opt(scores)
## Generate (randomly) K predictive scores for n observations n <- 50 K <- 5 scores <- matrix(runif(n*K), nrow = n, ncol = K) ## Find Bayesian Predictive Stacking weights opt_weights <- conv_opt(scores)
Compute the BPS weights by convex optimization
CVXR_opt(scores)
CVXR_opt(scores)
scores |
matrix |
conv_opt function to perform convex optimiazion with CVXR R package
Evaluate the density of a set of unobserved response with respect to the conditional posterior predictive
d_pred_cpp(data, X_u, Y_u, d_u, d_us, hyperpar, poster)
d_pred_cpp(data, X_u, Y_u, d_u, d_us, hyperpar, poster)
data |
list two elements: first named |
X_u |
matrix unobserved instances covariate matrix |
Y_u |
matrix unobserved instances response matrix |
d_u |
matrix unobserved instances distance matrix |
d_us |
matrix cross-distance between unobserved and observed instances matrix |
hyperpar |
list two elemets: first named |
poster |
list output from |
vector posterior predictive density evaluations
Evaluate the density of a set of unobserved response with respect to the conditional posterior predictive
d_pred_cpp_MvT(data, X_u, Y_u, d_u, d_us, hyperpar, poster)
d_pred_cpp_MvT(data, X_u, Y_u, d_u, d_us, hyperpar, poster)
data |
list two elements: first named |
X_u |
matrix unobserved instances covariate matrix |
Y_u |
matrix unobserved instances response matrix |
d_u |
matrix unobserved instances distance matrix |
d_us |
matrix cross-distance between unobserved and observed instances matrix |
hyperpar |
list two elemets: first named |
poster |
list output from |
double posterior predictive density evaluation
Compute the KCV of the density evaluations for fixed values of the hyperparameters
dens_kcv(data, priors, coords, hyperpar, K)
dens_kcv(data, priors, coords, hyperpar, K)
data |
list two elements: first named |
priors |
list priors: named |
coords |
matrix sample coordinates for X and Y |
hyperpar |
list two elemets: first named |
K |
integer number of folds |
vector posterior predictive density evaluations
Compute the KCV of the density evaluations for fixed values of the hyperparameters
dens_kcv_MvT(data, priors, coords, hyperpar, K)
dens_kcv_MvT(data, priors, coords, hyperpar, K)
data |
list two elements: first named |
priors |
list priors: named |
coords |
matrix sample coordinates for X and Y |
hyperpar |
list two elemets: first named |
K |
integer number of folds |
vector posterior predictive density evaluations
Compute the LOOCV of the density evaluations for fixed values of the hyperparameters
dens_loocv(data, priors, coords, hyperpar)
dens_loocv(data, priors, coords, hyperpar)
data |
list two elements: first named |
priors |
list priors: named |
coords |
matrix sample coordinates for X and Y |
hyperpar |
list two elemets: first named |
vector posterior predictive density evaluations
Compute the LOOCV of the density evaluations for fixed values of the hyperparameters
dens_loocv_MvT(data, priors, coords, hyperpar)
dens_loocv_MvT(data, priors, coords, hyperpar)
data |
list two elements: first named |
priors |
list priors: named |
coords |
matrix sample coordinates for X and Y |
hyperpar |
list two elemets: first named |
vector posterior predictive density evaluations
expand.grid()
in R
)Build a grid from two vector (i.e. equivalent to expand.grid()
in R
)
expand_grid_cpp(x, y)
expand_grid_cpp(x, y)
x |
vector first vector of numeric elements |
y |
vector second vector of numeric elements |
matrix expanded grid of combinations
## Create a matrix from all combination of vectors x <- seq(0, 10, length.out = 100) y <- seq(-1, 1, length.out = 20) grid <- expand_grid_cpp(x = x, y = y)
## Create a matrix from all combination of vectors x <- seq(0, 10, length.out = 100) y <- seq(-1, 1, length.out = 20) grid <- expand_grid_cpp(x = x, y = y)
and
(i.e. updated parameters)Compute the parameters for the posteriors distribution of and
(i.e. updated parameters)
fit_cpp(data, priors, coords, hyperpar)
fit_cpp(data, priors, coords, hyperpar)
data |
list two elements: first named |
priors |
list priors: named |
coords |
matrix sample coordinates for X and Y |
hyperpar |
list two elemets: first named |
list posterior update parameters
and
(i.e. updated parameters)Compute the parameters for the posteriors distribution of and
(i.e. updated parameters)
fit_cpp_MvT(data, priors, coords, hyperpar)
fit_cpp_MvT(data, priors, coords, hyperpar)
data |
list two elements: first named |
priors |
list priors: named |
coords |
matrix sample coordinates for X and Y |
hyperpar |
list two elemets: first named |
list posterior update parameters
Function to subset data for meta-analysis
forceSymmetry_cpp(mat)
forceSymmetry_cpp(mat)
mat |
matrix not-symmetric matrix |
matrix symmetric matrix (lower triangular of mat
is used)
## Force matrix to be symmetric (avoiding numerical problems) n <- 4 X <- matrix(runif(n*n), nrow = n, ncol = n) X <- forceSymmetry_cpp(mat = X)
## Force matrix to be symmetric (avoiding numerical problems) n <- 4 X <- matrix(runif(n*n), nrow = n, ncol = n) X <- forceSymmetry_cpp(mat = X)
Return the CV predictive density evaluations for all the model combinations
models_dens(data, priors, coords, hyperpar, useKCV, K)
models_dens(data, priors, coords, hyperpar, useKCV, K)
data |
list two elements: first named |
priors |
list priors: named |
coords |
matrix sample coordinates for X and Y |
hyperpar |
list two elemets: first named |
useKCV |
if |
K |
integer number of folds |
matrix posterior predictive density evaluations (each columns represent a different model)
Return the CV predictive density evaluations for all the model combinations
models_dens_MvT(data, priors, coords, hyperpar, useKCV, K)
models_dens_MvT(data, priors, coords, hyperpar, useKCV, K)
data |
list two elements: first named |
priors |
list priors: named |
coords |
matrix sample coordinates for X and Y |
hyperpar |
list two elemets: first named |
useKCV |
if |
K |
integer number of folds |
matrix posterior predictive density evaluations (each columns represent a different model)
Sample R draws from the posterior distributions
post_draws(poster, R, par, p)
post_draws(poster, R, par, p)
poster |
list output from |
R |
integer number of posterior samples |
par |
if |
p |
integer if |
list posterior samples
Sample R draws from the posterior distributions
post_draws_MvT(poster, R, par, p)
post_draws_MvT(poster, R, par, p)
poster |
list output from |
R |
integer number of posterior samples |
par |
if |
p |
integer if |
list posterior samples
Predictive sampler for Conjugate Bayesian Multivariate Linear Models
pred_bayesMvLMconjugate(X_new, B_samples, Sigma_samples)
pred_bayesMvLMconjugate(X_new, B_samples, Sigma_samples)
X_new |
matrix |
B_samples |
array of posterior sample for |
Sigma_samples |
array of posterior samples for |
Y_pred matrix of posterior mean for response matrix Y predictions
Y_pred_samples array of posterior predictive sample for response matrix Y
## Generate data n <- 100 p <- 3 q <- 2 Y <- matrix(rnorm(n*q), nrow = n, ncol = q) X <- matrix(rnorm(n*p), nrow = n, ncol = p) ## Prior parameters mu_B <- matrix(0, p, q) V_B <- diag(10, p) nu <- 3 Psi <- diag(q) ## Samples from posteriors n_iter <- 1000 burn_in <- 500 set.seed(1234) samples <- spBPS::bayesMvLMconjugate(Y, X, mu_B, V_B, nu, Psi, n_iter, burn_in) ## Extract posterior samples B_samples <- samples$B_samples Sigma_samples <- samples$Sigma_samples ## Samples from predictive posterior (based posterior samples) m <- 50 X_new <- matrix(rnorm(m*p), nrow = m, ncol = p) pred <- spBPS::pred_bayesMvLMconjugate(X_new, B_samples, Sigma_samples)
## Generate data n <- 100 p <- 3 q <- 2 Y <- matrix(rnorm(n*q), nrow = n, ncol = q) X <- matrix(rnorm(n*p), nrow = n, ncol = p) ## Prior parameters mu_B <- matrix(0, p, q) V_B <- diag(10, p) nu <- 3 Psi <- diag(q) ## Samples from posteriors n_iter <- 1000 burn_in <- 500 set.seed(1234) samples <- spBPS::bayesMvLMconjugate(Y, X, mu_B, V_B, nu, Psi, n_iter, burn_in) ## Extract posterior samples B_samples <- samples$B_samples Sigma_samples <- samples$Sigma_samples ## Samples from predictive posterior (based posterior samples) m <- 50 X_new <- matrix(rnorm(m*p), nrow = m, ncol = p) pred <- spBPS::pred_bayesMvLMconjugate(X_new, B_samples, Sigma_samples)
Draw from the conditional posterior predictive for a set of unobserved covariates
r_pred_cond(data, X_u, d_u, d_us, hyperpar, poster, post)
r_pred_cond(data, X_u, d_u, d_us, hyperpar, poster, post)
data |
list two elements: first named |
X_u |
matrix unobserved instances covariate matrix |
d_u |
matrix unobserved instances distance matrix |
d_us |
matrix cross-distance between unobserved and observed instances matrix |
hyperpar |
list two elemets: first named |
poster |
list output from |
post |
list output from |
list posterior predictive samples
Draw from the conditional posterior predictive for a set of unobserved covariates
r_pred_cond_MvT(data, X_u, d_u, d_us, hyperpar, poster, post)
r_pred_cond_MvT(data, X_u, d_u, d_us, hyperpar, poster, post)
data |
list two elements: first named |
X_u |
matrix unobserved instances covariate matrix |
d_u |
matrix unobserved instances distance matrix |
d_us |
matrix cross-distance between unobserved and observed instances matrix |
hyperpar |
list two elemets: first named |
poster |
list output from |
post |
list output from |
list posterior predictive samples
Draw from the joint posterior predictive for a set of unobserved covariates
r_pred_joint(data, X_u, d_u, d_us, hyperpar, poster, R)
r_pred_joint(data, X_u, d_u, d_us, hyperpar, poster, R)
data |
list two elements: first named |
X_u |
matrix unobserved instances covariate matrix |
d_u |
matrix unobserved instances distance matrix |
d_us |
matrix cross-distance between unobserved and observed instances matrix |
hyperpar |
list two elemets: first named |
poster |
list output from |
R |
integer number of posterior predictive samples |
list posterior predictive samples
Draw from the joint posterior predictive for a set of unobserved covariates
r_pred_joint_MvT(data, X_u, d_u, d_us, hyperpar, poster, R)
r_pred_joint_MvT(data, X_u, d_u, d_us, hyperpar, poster, R)
data |
list two elements: first named |
X_u |
matrix unobserved instances covariate matrix |
d_u |
matrix unobserved instances distance matrix |
d_us |
matrix cross-distance between unobserved and observed instances matrix |
hyperpar |
list two elemets: first named |
poster |
list output from |
R |
integer number of posterior predictive samples |
list posterior predictive samples
Draw from the marginals posterior predictive for a set of unobserved covariates
r_pred_marg(data, X_u, d_u, d_us, hyperpar, poster, R)
r_pred_marg(data, X_u, d_u, d_us, hyperpar, poster, R)
data |
list two elements: first named |
X_u |
matrix unobserved instances covariate matrix |
d_u |
matrix unobserved instances distance matrix |
d_us |
matrix cross-distance between unobserved and observed instances matrix |
hyperpar |
list two elemets: first named |
poster |
list output from |
R |
integer number of posterior predictive samples |
list posterior predictive samples
Draw from the joint posterior predictive for a set of unobserved covariates
r_pred_marg_MvT(data, X_u, d_u, d_us, hyperpar, poster, R)
r_pred_marg_MvT(data, X_u, d_u, d_us, hyperpar, poster, R)
data |
list two elements: first named |
X_u |
matrix unobserved instances covariate matrix |
d_u |
matrix unobserved instances distance matrix |
d_us |
matrix cross-distance between unobserved and observed instances matrix |
hyperpar |
list two elemets: first named |
poster |
list output from |
R |
integer number of posterior predictive samples |
list posterior predictive samples
Function to sample integers (index)
sample_index(size, length, p)
sample_index(size, length, p)
size |
integer dimension of the set to sample |
length |
integer number of elements to sample |
p |
vector sampling probabilities |
vector sample of integers
Perform prediction for BPS accelerated models - loop over prediction set
spPredict_BPS(data, X_u, priors, coords, crd_u, hyperpar, W, R, J)
spPredict_BPS(data, X_u, priors, coords, crd_u, hyperpar, W, R, J)
data |
list two elements: first named |
X_u |
matrix unobserved instances covariate matrix |
priors |
list priors: named |
coords |
matrix sample coordinates for X and Y |
crd_u |
matrix unboserved instances coordinates |
hyperpar |
list two elemets: first named |
W |
matrix set of stacking weights |
R |
integer number of desired samples |
J |
integer number of desired partition of prediction set |
list BPS posterior predictive samples
Function to subset data for meta-analysis
subset_data(data, K)
subset_data(data, K)
data |
list three elements: first named |
K |
integer number of desired subsets |
list subsets of data, and the set of indexes
## Create a list of K random subsets given a list with Y, X, and crd n <- 100 p <- 3 q <- 2 X <- matrix(rnorm(n*p), nrow = n, ncol = p) Y <- matrix(rnorm(n*q), nrow = n, ncol = q) crd <- matrix(runif(n*2), nrow = n, ncol = 2) subsets <- subset_data(data = list(Y = Y, X = X, crd = crd), K = 10)
## Create a list of K random subsets given a list with Y, X, and crd n <- 100 p <- 3 q <- 2 X <- matrix(rnorm(n*p), nrow = n, ncol = p) Y <- matrix(rnorm(n*q), nrow = n, ncol = q) crd <- matrix(runif(n*2), nrow = n, ncol = 2) subsets <- subset_data(data = list(Y = Y, X = X, crd = crd), K = 10)