diff --git a/.Rbuildignore b/.Rbuildignore index c4defe2..b70af4e 100644 --- a/.Rbuildignore +++ b/.Rbuildignore @@ -2,4 +2,8 @@ ^\.Rproj\.user$ ^figures ^data -^simulations \ No newline at end of file +^simulations +^application +^other_code +^testing +^\.vscode \ No newline at end of file diff --git a/DESCRIPTION b/DESCRIPTION index dfe8c69..a5c20b9 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,11 +1,11 @@ Package: spiox Type: Package -Title: Inside-Out Cross-covariance for spatial multivariate data -Version: 0.3 -Date: 2024-10-11 +Title: Inside-out cross-covariance (IOX) for spatial multivariate data +Version: 0.5 +Date: 2024-12-02 Author: Michele Peruzzi Maintainer: Michele Peruzzi -Description: Fits Bayesian models for multivariate spatial data using Inside-Out Cross-covariance and Vecchia approximations of GPs. +Description: Fits Bayesian models for multivariate spatial data using IOX and Vecchia approximations of GPs. License: GPL (>= 2) Encoding: UTF-8 RoxygenNote: 7.3.1 diff --git a/R/RcppExports.R b/R/RcppExports.R index 74acc9a..94f8adf 100644 --- a/R/RcppExports.R +++ b/R/RcppExports.R @@ -57,19 +57,23 @@ run_spf_model <- function(Y, n_factors, delta_gamma_shape, delta_gamma_rate, dl_ .Call(`_spiox_run_spf_model`, Y, n_factors, delta_gamma_shape, delta_gamma_rate, dl_dirichlet_a, Lambda_start, Delta_start, mcmc, print_every, seq_lambda) } -spiox_wishart <- function(Y, X, coords, custom_dag, theta_opts, Sigma_start, mvreg_B_start, mcmc = 1000L, print_every = 100L, matern = TRUE, sample_iwish = TRUE, sample_mvr = TRUE, sample_theta_gibbs = TRUE, upd_theta_opts = TRUE, num_threads = 1L) { +spiox_wishart <- function(Y, X, coords, custom_dag, theta_opts, Sigma_start, mvreg_B_start, mcmc = 1000L, print_every = 100L, matern = 1L, sample_iwish = TRUE, sample_mvr = TRUE, sample_theta_gibbs = TRUE, upd_theta_opts = TRUE, num_threads = 1L) { .Call(`_spiox_spiox_wishart`, Y, X, coords, custom_dag, theta_opts, Sigma_start, mvreg_B_start, mcmc, print_every, matern, sample_iwish, sample_mvr, sample_theta_gibbs, upd_theta_opts, num_threads) } -spiox_latent <- function(Y, X, coords, custom_dag, theta_opts, Sigma_start, mvreg_B_start, mcmc = 1000L, print_every = 100L, matern = TRUE, sample_iwish = TRUE, sample_mvr = TRUE, sample_theta_gibbs = TRUE, upd_theta_opts = TRUE, num_threads = 1L, sampling = 2L) { +spiox_latent <- function(Y, X, coords, custom_dag, theta_opts, Sigma_start, mvreg_B_start, mcmc = 1000L, print_every = 100L, matern = 1L, sample_iwish = TRUE, sample_mvr = TRUE, sample_theta_gibbs = TRUE, upd_theta_opts = TRUE, num_threads = 1L, sampling = 2L) { .Call(`_spiox_spiox_latent`, Y, X, coords, custom_dag, theta_opts, Sigma_start, mvreg_B_start, mcmc, print_every, matern, sample_iwish, sample_mvr, sample_theta_gibbs, upd_theta_opts, num_threads, sampling) } -spiox_predict <- function(X_new, coords_new, Y, X, coords, dag, B, Sigma, theta, matern = TRUE, num_threads = 1L) { +spiox_predict <- function(X_new, coords_new, Y, X, coords, dag, B, Sigma, theta, matern = 1L, num_threads = 1L) { .Call(`_spiox_spiox_predict`, X_new, coords_new, Y, X, coords, dag, B, Sigma, theta, matern, num_threads) } -spiox_latent_predict <- function(X_new, coords_new, coords, dag, W, B, Sigma, Dvec, theta, matern = TRUE, num_threads = 1L) { +spiox_predict_part <- function(Y_new, X_new, coords_new, Y, X, coords, dag, B, Sigma, theta, matern = 1L, num_threads = 1L) { + .Call(`_spiox_spiox_predict_part`, Y_new, X_new, coords_new, Y, X, coords, dag, B, Sigma, theta, matern, num_threads) +} + +spiox_latent_predict <- function(X_new, coords_new, coords, dag, W, B, Sigma, Dvec, theta, matern = 1L, num_threads = 1L) { .Call(`_spiox_spiox_latent_predict`, X_new, coords_new, coords, dag, W, B, Sigma, Dvec, theta, matern, num_threads) } diff --git a/README.md b/README.md index b4b76d6..905713e 100644 --- a/README.md +++ b/README.md @@ -1,11 +1,11 @@ -# Inside-Out Cross-covariance for multivariate Gaussian processes +# Inside-out cross-covariance for multivariate Gaussian processes -![Simulated data from IOX](figures/prior_sample_3.png) +![Simulated data from IOX](https://github.com/mkln/spiox-paper/figures/prior_sample_3.png) Noise-free trivariate GP using IOX. -![Cross-covariance](figures/cij_plot.png) +![Cross-covariance](https://github.com/mkln/spiox-paper/figures/cij_plot.png) Cross-covariance between two variables, for varying spatial range (left) or smoothness (right) of one of them relative to the other. diff --git a/figures/cij_plot.pdf b/figures/cij_plot.pdf deleted file mode 100644 index 38cf4bb..0000000 Binary files a/figures/cij_plot.pdf and /dev/null differ diff --git a/figures/cij_plot.png b/figures/cij_plot.png deleted file mode 100644 index 8431b15..0000000 Binary files a/figures/cij_plot.png and /dev/null differ diff --git a/figures/cij_plot.r b/figures/cij_plot.r deleted file mode 100644 index a554e5c..0000000 --- a/figures/cij_plot.r +++ /dev/null @@ -1,89 +0,0 @@ -rm(list=ls()) -library(tidyverse) -library(magrittr) -library(Matrix) -library(latex2exp) -library(scico) - -set.seed(4) - -image.matrix <- \(x) { - image(as(x, "dgCMatrix")) -} - -iox_spcor <- function(outcome_i, outcome_j, h_values, - x_test_grid, - cx, theta_mat, - n_angles=5, diag_only=T, at_limit=T){ - - all_combos_setup <- expand.grid(x_test_grid, x_test_grid, - h_values, seq(0,2*pi,length.out=n_angles)) %>% as.matrix() - colnames(all_combos_setup) <- c("s1x", "s1y", "h", "angle") - - all_combos_setup %<>% data.frame() %>% - mutate(s2x = s1x + h*cos(angle), s2y = s1y + h*sin(angle)) - - testx <- all_combos_setup %>% dplyr::select(s1x, s1y) %>% as.matrix() - testy <- all_combos_setup %>% dplyr::select(s2x, s2y) %>% as.matrix() - - cov_computed <- iox(testx, testy, outcome_i, outcome_j, - cx, theta_mat, - matern = TRUE, - diag_only=diag_only, at_limit=at_limit) - - all_combos_setup %<>% mutate(spcov = cov_computed) - - spcov_by_h <- all_combos_setup %>% group_by(h) %>% summarise(spcov = mean(spcov)) - - return( spcov_by_h ) -} - - -q <- 2 - -cor_pair <- function(phi, nu){ - - theta_mat <- matrix(1, ncol=2, nrow=4) - theta_mat[1,] <- c(10, phi) - theta_mat[2,] <- c(1, 1) - theta_mat[3,] <- c(1, nu) - theta_mat[4,] <- rep(1e-10, 2) - - set.seed(1) - xx <- seq(0, 1, length.out=10) - coords <- expand.grid(xx, xx) - cx <- as.matrix(coords) - n_g <- 5 - - xg <- xx[ seq(1, length(xx), length(xx)/(n_g+1)) ] %>% tail(-1) - hvec <- c(10^seq(-5, -2, length.out=25), 10^seq(-2, -0.5, length.out=25)) - - spcor <- iox_spcor(1, 2, hvec, xg, cx, theta_mat, n_angles=10, diag_only=T, at_limit=F) - return(spcor %>% mutate(phi = phi, nu = nu)) -} - -phinu <- expand.grid(phi = seq(1, 40, 2), - nu = 1) %>% bind_rows( - expand.grid(phi = 10, nu = seq(0.5, 1.9, length.out=25))) - -results <- phinu %>% as.matrix() %>% plyr::alply(1, \(x) cor_pair(x[1], x[2])) - -phi_corr <- results %>% bind_rows() %>% - filter(nu==1) %>% - ggplot(aes(h, spcov, group=phi, color=phi)) + - scale_color_scico(palette="managua") + - geom_line() + - theme_minimal() + - labs(color=TeX("$\\phi$"), x="Distance", y="Correlation") - -nu_corr <- results %>% bind_rows() %>% - filter(phi==10) %>% - ggplot(aes(h, spcov, group=nu, color=nu)) + - scale_color_scico(palette="managua", direction = -1) + - geom_line() + - theme_minimal() + - labs(color=TeX("$\\nu$"), x="Distance", y="Correlation") - -cij_plot <- grid.arrange(phi_corr, nu_corr, nrow=1) -ggsave("figures/cij_plot.pdf", plot=cij_plot, width=8.75, height=3) -ggsave("figures/cij_plot.png", plot=cij_plot, width=8.75, height=3) diff --git a/figures/prior_sample_3.png b/figures/prior_sample_3.png deleted file mode 100644 index cb1d7d6..0000000 Binary files a/figures/prior_sample_3.png and /dev/null differ diff --git a/figures/prior_sample_3.r b/figures/prior_sample_3.r deleted file mode 100644 index af8f277..0000000 --- a/figures/prior_sample_3.r +++ /dev/null @@ -1,93 +0,0 @@ -rm(list=ls()) -library(tidyverse) -library(magrittr) -library(Matrix) -library(latex2exp) - -set.seed(4) - -image.matrix <- \(x) { - image(as(x, "dgCMatrix")) -} - - - - - - -q <- 3 - -theta_mat <- matrix(1, ncol=q, nrow=4) -theta_mat[1,] <- c(5, 15, 30) -theta_mat[2,] <- c(1, 1, 1) -theta_mat[3,] <- c(0.5, 1.2, 1.9) -theta_mat[4,] <- rep(1e-3, q) - -sds <- c(1, 1, 1) -Omega <- matrix(c(1,-.9,0.7,-.9,1,-.5,0.7,-.5,1), ncol=3) -Sigma <- diag(sds) %*% Omega %*% diag(sds) - -St <- chol(Sigma) -S <- t(St) - -# spatial -xx <- seq(0, 1, length.out=20) -coords <- expand.grid(xx, xx) -nr <- nrow(coords) -cx <- as.matrix(coords) -A <- matrix(runif(nr), ncol=1) - -m <- 20 - -custom_dag <- dag_vecchia(cx, m) - - -m <- 50 - -custom_dag <- dag_vecchia(cx, m) - -Hlist <- 1:q %>% sapply(\(j) - spiox::daggp_build(cx, custom_dag, - phi=theta_mat[1,j], sigmasq=theta_mat[2,j], - nu=theta_mat[3,j], tausq=theta_mat[4,j], - matern=T, 16)$H -) - - -V <- matrix(rnorm(nr * q), ncol=q) %*% St - -Y <- V - -save_which <- rep(0, q) - -# make q spatial outcomes -which_process <- rep(0, q) -for(i in 1:q){ - cat(i, "\n") - Y[,i] <- Matrix::solve(Hlist[[i]], V[,i], sparse=T) -} - -df <- data.frame(coords, y=Y) %>% - pivot_longer(cols=-c(Var1, Var2)) %>% - mutate(name = ifelse(name == "y.1", "Outcome 1", ifelse(name == "y.2", "Outcome 2", "Outcome 3"))) - -( p2 <- ggplot(df, - aes(Var1, Var2, fill=value)) + - geom_raster() + - scico::scale_fill_scico(palette="vik") + # bam, broc, cork, managua, vik - facet_wrap(~name, ncol=q) + - theme_minimal() + - labs(x=NULL, y=NULL, fill="Value") + - scale_x_continuous(breaks=c(0.5, 1), expand=c(0,0)) + - scale_y_continuous(breaks=c(0, 0.5, 1), expand=c(0,0)) + - theme( - panel.grid = element_blank(), - panel.spacing.x = unit(10, "pt"), - panel.border = element_rect(fill=NA, color="black"), - axis.text.x = element_text(margin = margin(t = 0), hjust=1), - - axis.text.y = element_text(margin = margin(r = 0), vjust=1) ) ) - -#ggsave("figures/prior_sample_3.pdf", plot=p2, width=8.7, height=3) -#ggsave("figures/prior_sample_3.png", plot=p2, width=8.7, height=3) - diff --git a/figures/prior_sample_3_diff.r b/figures/prior_sample_3_diff.r deleted file mode 100644 index b41f12f..0000000 --- a/figures/prior_sample_3_diff.r +++ /dev/null @@ -1,103 +0,0 @@ -rm(list=ls()) -library(tidyverse) -library(magrittr) -library(Matrix) -library(latex2exp) - -set.seed(2024) - -image.matrix <- \(x) { - image(as(x, "dgCMatrix")) -} - -q <- 4 - -theta_mat <- matrix(1, ncol=q, nrow=4) -theta_mat[1,] <- c(15, 15, 15, 15) # phi -theta_mat[2,] <- c(1, 1, 1, 1) # sigmasq -theta_mat[3,] <- c(1, 14.9, 1.25, 1.5) # nu -theta_mat[4,] <- rep(1e-3, q) - -sds <- c(1, 1, 1, 1) -Omega <- cbind(c(1,-.9,0.7, 0.8), c(-.9,1,-.5, -0.7), c(0.7,-.5,1, 0.8), c(0.8, -0.7, 0.8, 1)) -Sigma <- diag(sds) %*% Omega %*% diag(sds) - -St <- chol(Sigma) -S <- t(St) - -# spatial -xx <- seq(0, 1, length.out=200) -coords <- expand.grid(xx, xx) -nr <- nrow(coords) -cx <- as.matrix(coords) -A <- matrix(runif(nr), ncol=1) - -m <- 40 - -custom_dag <- dag_vecchia(cx, m) - - -# warping for the 3rd outcome -warpx <- spiox::daggp_build(cx, custom_dag, phi=1, sigmasq=.03, - nu=1, tausq=0, matern=1, 16)$H %>% solve(rnorm(nr)) -warpy <- spiox::daggp_build(cx, custom_dag, phi=1, sigmasq=.03, - nu=1, tausq=0, matern=1, 16)$H %>% solve(rnorm(nr)) -cx_warp <- cx + cbind(warpx, warpy) - -cxlist <- list(cx, cx, cx_warp, cx) - -cov_choice <- c(1, 2, 1, 1) -Hlist <- 1:q %>% sapply(\(j) {cat(j, "\n") - spiox::daggp_build(cxlist[[j]], custom_dag, - phi=theta_mat[1,j], sigmasq=theta_mat[2,j], - nu=theta_mat[3,j], tausq=theta_mat[4,j], - matern=cov_choice[j], 16)$H -}) - -# warping for the 3rd outcome - -# sill for the 4th outcome -wsdelta <- spiox::daggp_build(cx, custom_dag, phi=30, sigmasq=.5, - nu=1.5, tausq=0, matern=1, 16)$H %>% solve(rnorm(nr)) - -V <- matrix(rnorm(nr * q), ncol=q) %*% St -V[,4] <- (wsdelta+1) * V[,4] - -Y <- V - -save_which <- rep(0, q) - -# make q spatial outcomes -which_process <- rep(0, q) -for(i in 1:q){ - cat(i, "\n") - Y[,i] <- Matrix::solve(Hlist[[i]], V[,i], sparse=T) -} - -df <- data.frame(coords, y=Y) %>% - pivot_longer(cols=-c(Var1, Var2)) %>% - mutate(name = - ifelse(name == "y.1", "Outcome 1", - ifelse(name == "y.2", "Outcome 2", - ifelse(name == "y.3", "Outcome 3", "Outcome 4")))) - -( p2 <- ggplot(df, - aes(Var1, Var2, fill=value)) + - geom_raster() + - scico::scale_fill_scico(palette="vik") + # bam, broc, cork, managua, vik - facet_wrap(~name, ncol=2) + - theme_minimal() + - labs(x=NULL, y=NULL, fill="Value") + - scale_x_continuous(breaks=c(0.5, 1), expand=c(0,0)) + - scale_y_continuous(breaks=c(0, 0.5, 1), expand=c(0,0)) + - theme( - panel.grid = element_blank(), - panel.spacing.x = unit(10, "pt"), - panel.border = element_rect(fill=NA, color="black"), - axis.text.x = element_text(margin = margin(t = 0), hjust=1), - - axis.text.y = element_text(margin = margin(r = 0), vjust=1) ) ) - -ggsave("figures/prior_sample_3_diff2.pdf", plot=p2, width=7, height=6.7) -ggsave("figures/prior_sample_3_diff2.png", plot=p2, width=7, height=6.7) - diff --git a/figures/prior_sample_Sc.r b/figures/prior_sample_Sc.r deleted file mode 100644 index 45c96d3..0000000 --- a/figures/prior_sample_Sc.r +++ /dev/null @@ -1,107 +0,0 @@ -rm(list=ls()) -library(tidyverse) -library(magrittr) -library(Matrix) -library(latex2exp) -library(spiox) - -set.seed(4) - -image.matrix <- \(x) { - image(as(x, "dgCMatrix")) -} - -q <- 3 - -theta_mat <- matrix(1, ncol=q, nrow=4) -theta_mat[1,] <- c(5, 15, 30) -theta_mat[2,] <- c(1, 1, 1) -theta_mat[3,] <- c(0.5, 1.2, 1.9) -theta_mat[4,] <- rep(1e-3, q) - -sds <- c(1, 1, 1) -Omega <- matrix(c(1,-.9,0.7,-.9,1,-.5,0.7,-.5,1), ncol=3) -Sigma <- diag(sds) %*% Omega %*% diag(sds) - -St <- chol(Sigma) -S <- t(St) - -# referejce location -xS <- seq(0, 1, length.out=50) -cooS <- expand.grid(xS, xS) -nrS <- nrow(cooS) -cxS <- as.matrix(cooS) + matrix(runif(nrS*2, -1e-5, 1e-5), ncol=2) - -m <- 50 - -custom_dag_S <- dag_vecchia(cxS, m) - -# sample at ref locations using vecchia dag GP -Hlist <- 1:q %>% sapply(\(j) - spiox::daggp_build(cxS, custom_dag_S, - phi=theta_mat[1,j], sigmasq=theta_mat[2,j], - nu=theta_mat[3,j], tausq=theta_mat[4,j], - matern=T, 16)$H -) -nr <- nrow(cooS) -V <- matrix(rnorm(nr * q), ncol=q) %*% St -Ys <- V -save_which <- rep(0, q) -# make q spatial outcomes -for(i in 1:q){ - cat(i, "\n") - Ys[,i] <- Matrix::solve(Hlist[[i]], V[,i], sparse=T) -} - -# use spiox predictions to make the sample -# an alternative is to augment the dag with the U section and sample jointly at S and U -# since that does the same thing as making predictions based on Ys - - -# spiox_predict wants output from a response model -# here we pretend we have a sample of size 1 (the arrays below) -# and use a zero beta coefficient matrix - -# non reference locations -xU <- seq(0, 1, length.out=200) -cooU <- expand.grid(xU, xU) -cxU <- as.matrix(cooU) - -X <- matrix(0, nrow=nrow(cxS), ncol=1) -X_new <- matrix(0, nrow=nrow(cxU), ncol=1) -B_arr <- array(0, dim=c(1, q, 1)) -Sigma_arr <- array(0, dim=c(q,q,1)) -Sigma_arr[,,1] <- Sigma -theta_arr <- array(0, dim=c(4,q,1)) -theta_arr[,,1] <- theta_mat - - -custom_dag_U <- dag_vecchia_predict(cxS, cxU, m = m) - -Yu_out <- spiox::spiox_predict(X_new, cxU, Ys, X, cxS, custom_dag_U, B_arr, Sigma_arr, theta_arr, 8) - -# plot - -df <- cooU %>% cbind(data.frame(y=Yu_out$Y, type="U")) %>% - pivot_longer(cols=-c(Var1, Var2, type)) %>% - mutate(name = ifelse(name == "y.1", "Outcome 1", ifelse(name == "y.2", "Outcome 2", "Outcome 3"))) - - -( p2 <- ggplot(df, - aes(Var1, Var2, fill=value)) + - geom_raster() + - scico::scale_fill_scico(palette="vik") + # bam, broc, cork, managua, vik - facet_wrap(~name, ncol=q) + - theme_minimal() + - labs(x=NULL, y=NULL, fill="Value") + - scale_x_continuous(breaks=c(0.5, 1), expand=c(0,0)) + - scale_y_continuous(breaks=c(0, 0.5, 1), expand=c(0,0)) + - theme( - panel.grid = element_blank(), - panel.spacing.x = unit(10, "pt"), - panel.border = element_rect(fill=NA, color="black"), - axis.text.x = element_text(margin = margin(t = 0), hjust=1), - - axis.text.y = element_text(margin = margin(r = 0), vjust=1) ) ) - -ggsave("figures/prior_sample_Sc.pdf", plot=p2, width=8.7, height=3) diff --git a/figures/visualize_iox.r b/figures/visualize_iox.r deleted file mode 100644 index e9c5ff5..0000000 --- a/figures/visualize_iox.r +++ /dev/null @@ -1,82 +0,0 @@ -library(tidyverse) -library(Matrix) -library(scico) -set.seed(1) -nulist <- c(.5, .7, 1.5, 1) -q <- length(nulist) - -Sigma <- solve(rWishart(1, q+1, diag(q))[,,1]) - -x <- seq(0,1, length.out=4) -coords <- expand.grid(x,x) -cx <- as.matrix(coords) -nr <- nrow(coords) -In <- diag(nr) - -C <- 1:q %>% lapply(\(j) - spiox::Correlationc(cx,cx,theta = c(1,1,j,0), covar = 1, T)) - - -# iox -L <- C %>% lapply(\(Cm) t(chol(Cm))) -Lbig <- Matrix::bdiag(L) -C_iox <- as(Lbig %*% (Sigma %x% In) %*% t(Lbig), "CsparseMatrix") -image(C_iox) - -# lmc -A <- t(chol(Sigma)) -C_lmc <- as((A %x% In) %*% Matrix::bdiag(C) %*% (t(A) %x% In), "CsparseMatrix") - - -ggplot_matrix_image <- function(M, label=NULL, full_symm=F){ - M <- M %>% as("CsparseMatrix") - M_df <- as.data.frame(summary(M)) - - colnames(M_df) <- c("row", "col", "value") # Renaming the columns for clarity - - if(full_symm){ - M_df <- bind_rows(M_df, M_df %>% - rename(rowx=row, colx=col) %>% - rename(row=colx, col=rowx)) %>% unique() - } - # Create the ggplot - plotmade <- ggplot(M_df, aes(x = col, y = row, fill = value)) + - geom_raster() + - scale_fill_scico(palette="batlowK") + - scale_y_reverse() + # Reverse the y-axis to match matrix convention - theme_void() + - labs(title = label, - fill = "Value") + - theme(legend.position="none", - panel.border = element_rect(fill=NA, color="black")) - return(plotmade) -} - -iox_overall <- ggplot_matrix_image(C_iox) -iox_left <- ggplot_matrix_image(Lbig) -iox_center <- ggplot_matrix_image(Sigma %x% In, NULL, T) -iox_right <- ggplot_matrix_image(t(Lbig)) - -iox_visual <- gridExtra::grid.arrange(iox_overall, iox_left, iox_center, iox_right, nrow=1) - - -lmc_overall <- ggplot_matrix_image(C_lmc) -lmc_left <- ggplot_matrix_image(A %x% In) -lmc_center <- ggplot_matrix_image(Matrix::bdiag(C), NULL, T) -lmc_right <- ggplot_matrix_image(t(A) %x% In) - -lmc_visual <- gridExtra::grid.arrange(lmc_overall, lmc_left, lmc_center, lmc_right, nrow=1) - - -iox_vs_lmc <- gridExtra::grid.arrange(iox_visual, lmc_visual, ncol=1) - -ggsave(filename="figures/iox_overall.pdf", plot=iox_overall, width=3, height=3) -ggsave(filename="figures/iox_left.pdf", plot=iox_left, width=3, height=3) -ggsave(filename="figures/iox_center.pdf", plot=iox_center, width=3, height=3) -ggsave(filename="figures/iox_right.pdf", plot=iox_right, width=3, height=3) - -ggsave(filename="figures/lmc_overall.pdf", plot=lmc_overall, width=3, height=3) -ggsave(filename="figures/lmc_left.pdf", plot=lmc_left, width=3, height=3) -ggsave(filename="figures/lmc_center.pdf", plot=lmc_center, width=3, height=3) -ggsave(filename="figures/lmc_right.pdf", plot=lmc_right, width=3, height=3) - diff --git a/iox_testing.r b/iox_testing.r deleted file mode 100644 index 41518df..0000000 --- a/iox_testing.r +++ /dev/null @@ -1,149 +0,0 @@ -rm(list=ls()) -library(tidyverse) -library(magrittr) -library(Matrix) - -set.seed(2) - -sumabs <- \(x) sum(abs(x)) -trroot <- \(x) sum( sqrt(svd(x)$d )) -svdroot <- \(x) with(svd(x), u %*% diag(sqrt(d)) %*% t(u)) -image.matrix <- \(x) { - image(as(x, "dgCMatrix")) -} - -# spatial -xx <- seq(0, 1, length.out=14) -coords <- #cbind(runif(100), runif(100))# - expand.grid(xx, xx) -cx <- as.matrix(coords) -nr <- nrow(cx) - -# multivariate -q <- 3 -Q <- rWishart(1, q+2, diag(q))[,,1] # -Sigma <- solve(Q) - - -theta_mat <- matrix(1, ncol=q, nrow=4) -theta_mat[1,] <- c(5, 20, 30)[1:q] -theta_mat[2,] <- c(1, 1, 1)[1:q] -theta_mat[3,] <- c(0.5, 1.2, 1.9)[1:q] -theta_mat[4,] <- rep(0, q) - - -Clist <- 1:q %>% lapply(\(j) spiox::Correlationc(cx, cx, theta_mat[,j], 1, 1)) -Llist <- Clist %>% lapply(\(C) t(chol(C))) -Lilist <- Llist %>% lapply(\(L) solve(L)) - - - - - -set.seed(1) -new1 <- matrix(runif(2), ncol=2) -new2 <- matrix(runif(2), ncol=2) - - - -Clist1 <- 1:q %>% lapply(\(j) spiox::Correlationc(rbind(cx, new1), rbind(cx, new1), theta_mat[,j], 1, 1)) -Llist1 <- Clist1 %>% lapply(\(C) t(chol(C))) -Lilist1 <- Llist1 %>% lapply(\(L) solve(L)) - -Lsparse <- Matrix::bdiag(Llist %>% lapply(\(L) { L[sample(x=1:prod(dim(L)), size=round(0.9*prod(dim(L))), replace=F)] <- 0; return(L)})) -Cbig <- Matrix::bdiag(Llist1) %*% (Sigma %x% diag(nr+1)) %*% t(Matrix::bdiag(Llist1)) -outx <- c(197, 394, 591) - - - -H <- Cbig[outx, -outx] %*% solve(Cbig[-outx, -outx]) -R <- Cbig[outx, outx] - H %*% Cbig[-outx, outx] - -Hmake <- 1:q %>% lapply(\(i) Clist1[[i]][outx[1], -outx[1]] %*% solve(Clist1[[i]][-outx[1], -outx[1]]) ) - -1:q - -Clist2 <- 1:q %>% lapply(\(j) spiox::Correlationc(rbind(cx, new2), rbind(cx, new2), theta_mat[,j], 1, 1)) -Llist2 <- Clist2 %>% lapply(\(C) t(chol(C))) -Lilist2 <- Llist2 %>% lapply(\(L) solve(L)) - - -Libig <- 1:q %>% lapply(\(j) { - Li <- matrix(0, nr + 2, nr + 2) - Li[1:nr, 1:nr] <- Lilist[[j]] - Li[nr+1, 1:(nr+1)] <- Lilist1[[j]][nr+1,] - Li[nr+2, c(1:nr, nr+2)] <- Lilist2[[j]][nr+1,] - return (Li) }) - -Lbig <- Libig %>% lapply(\(Li) solve(Li)) -ixnew <- c(nr+1, nr+2, 2*nr + 2 + 1, 2*nr +2 + 2) -(Matrix::bdiag(Lbig) %*% (Sigma %x% diag(nr+2)) %*% t(Matrix::bdiag(Lbig)))[ixnew, ixnew] - - -Lbig <- 1:q %>% lapply(\(j) { - L <- matrix(0, nr + 2, nr + 2) - L[1:nr, 1:nr] <- Llist[[j]] - L[nr+1, 1:(nr+1)] <- Llist1[[j]][nr+1,] - L[nr+2, c(1:nr, nr+2)] <- Llist2[[j]][nr+1,] - return (L) }) - -ixold <- -ixnew - -Cbig <- (Matrix::bdiag(Lbig) %*% ((Sigma) %x% diag(nr+2)) %*% - t(Matrix::bdiag(Lbig))) - -o1 <- 1:196 -o23 <- 197:594 - -Cond_1_given_23 <- Cbig[o1, o1] - Cbig[o1, o23] %*% solve(Cbig[o23, o23]) %*% Cbig[o23, o1] - - -Cbig2 <- (Sigma %x% matrix(1, nr+2, nr+2)) * - (do.call(rbind, Lbig) %*% t(do.call(rbind, Lbig))) - -Csave12_r <- (Cbig[ixnew, ixnew] - Cbig[ixnew, ixold] %*% solve(Cbig[ixold, ixold]) %*% Cbig[ixold, ixnew]) -Csave12_m <- Cbig[ixnew, ixnew] - -Lbig_sub <- Lbig %>% lapply(\(L) L %>% tail(2)) -(Sigma %x% matrix(1, 2, 2)) * (do.call(rbind, Lbig_sub) %*% t(do.call(rbind, Lbig_sub))) - -C <- function(x,y,theta){ - spiox::Correlationc(x,y,theta,T,F) -} -R <- function(x,theta){ - C(x,x,theta) - C(x,cx,theta) %*% solve(C(cx,cx,theta)) %*% C(cx,x,theta) -} - - -newx <- matrix(runif(10*2), ncol=2) -Cmat <- iox_mat(cx[1:2,], cx[1:2,], cx, theta_mat, matern = T, D_only = F) - - -for(i in 1:q){ - for(j in 1:q){ - cat( iox(rbind(new1, new2), rbind(new1, new2), i, j, cx, theta_mat) , "\n" ) - } -} - - - -i <- 1; j <- 1 -Sigma[i, j] * - (Lbig[[i]][nr+1,,drop=F] %*% t(Lbig[[j]][nr+1,,drop=F])) - -Sigma[i, j] * - iox(rbind(new1, new2), rbind(new1, new2), i, j, cx, theta_mat, matern = T, diag_only = F, limit = F) - -Sigma[i,j] * (C(new2, cx, theta_mat[,i]) %*% t(Lilist[[i]]) %*% Lilist[[j]] %*% C(cx, new2, theta_mat[,j]) + - sqrt( R(new2, theta_mat[,i]) * R(new2, theta_mat[,j]) ) ) - - -Ciox <- function(x, y, i, j){ - iox(x, y, i, j, cx, theta_mat, matern = T, diag_only = F, limit = F) -} - -Riox <- function(x, i, j){ - Ciox(x,x,i,j) - Ciox(x,cx,i,j) %*% solve(Ciox(cx,cx,i,j)) %*% Ciox(cx,x,i,j) -} -Ciox(rbind(new1, new2), rbind(new1, new2), 1, 2) -Sigma[2, 2] * Riox(rbind(new1, new2), 2, 2) diff --git a/simulations/multi_lmc.r b/simulations/multi_lmc.r deleted file mode 100644 index 09f3cc0..0000000 --- a/simulations/multi_lmc.r +++ /dev/null @@ -1,373 +0,0 @@ -args <- commandArgs(TRUE) - -starts <- args[1] -ends <- args[2] - -cat("Simulations ", starts:ends, "\n") - -library(tidyverse) -library(magrittr) -library(Matrix) -library(spiox) - -# many outcomes, parsimonious matern model - -image.matrix <- \(x) { - image(as(x, "dgCMatrix")) -} - -rtail <- \(x, nt){ - tail(x, nt) -} -tailh <- \(x, nt){ - tail(x, round(length(x)/2)) -} -rtailh <- \(x){ rtail(x, round(nrow(x)/2)) } -ctail <- \(x, nt){ - tail(x, c(NA,nt)) -} -ctailh <- \(x){ ctail(x, round(ncol(x)/2)) } -stail <- \(x, nt){ - tail(x, c(NA,NA,nt)) -} -stailh <- \(x){ stail(x, round(dim(x)[3]/2)) } - -perturb <- function(x, sd=1, symm=F){ - pert <- matrix(rnorm(prod(dim(x)), sd), ncol=ncol(x)) - if(symm){ - pert <- crossprod(pert) - } - - return(x + pert) -} -q <- 24 - -nthreads <- 7 - -for(oo in starts:ends){ - set.seed(oo) - cat(oo, "\n") - - optlist <- seq(15, 50, length.out=q) %>% sample(q, replace=T) - - # spatial - cx_in <- matrix(runif(2500*2), ncol=2) #3000 - colnames(cx_in) <- c("Var1","Var2") - n_in <- nrow(cx_in) - which_in <- 1:n_in - - xout <- seq(0, 1, length.out=20) #20 - coords_out <- expand.grid(xout, xout) - cx_out <- as.matrix(coords_out) - n_out <- nrow(cx_out) - which_out <- (n_in+1):(n_in+n_out) - - cx_all <- rbind(cx_in, cx_out) - nr_all <- nrow(cx_all) - - Clist <- optlist %>% lapply(\(phi) spiox::Correlationc(cx_all, cx_all, c(phi,1,1,1e-15), 1, TRUE) ) - Llist <- Clist %>% lapply(\(C) t(chol(C))) - - Q <- rWishart(1, q+1, 1/2 * diag(q))[,,1] # - Sigma <- solve(Q) - St <- chol(Sigma) - S <- t(St) - - V <- matrix(rnorm(nr_all * q), ncol=q) - V <- V - for(i in 1:q){ - V[,i] <- Llist[[i]] %*% V[,i] - } - Y_sp <- V %*% St - - # regression - p <- 2 - X <- matrix(1, ncol=1, nrow=nr_all) %>% cbind(matrix(rnorm(nr_all*(p-1)), ncol=p-1)) - - Beta <- matrix(rnorm(q * p), ncol=q) - - Y_regression <- X %*% Beta - - Y <- as.matrix(Y_sp + Y_regression) + matrix(rnorm(nr_all * q), ncol=q) %*% diag(rep(0.01, q)) - - Y_in <- Y[which_in,] - X_in <- X[which_in,] - - if(F){ - df <- data.frame(cx_out, y=as.matrix(Y_sp[which_out,])) %>% - pivot_longer(cols=-c(Var1, Var2)) - ggplot(df, - aes(Var1, Var2, fill=value)) + - geom_raster() + - scale_fill_viridis_c() + - facet_wrap(~name, ncol=5) - } - - simdata <- data.frame(coords=cx_all, Y_spatial=Y_sp, Y=Y, X=X) - #ggplot(simdata, aes(coords.Var1, coords.Var2, color=Y_spatial.1)) + geom_point() + scale_color_viridis_c() - - if(F){ - save(file=glue::glue("simulations/lmc_m/data_{oo}.RData"), - list=c("simdata", "oo", "simdata", "Sigma", "optlist")) - } - - ############################## - - set.seed(1) - - phitsq <- expand.grid(phi <- seq(10, 60, length.out=20), - tsq <- c(2*1e-6, 1e-5, 1e-4)) - - theta_opts <- rbind(phitsq$Var1, 1, 1, phitsq$Var2) - - theta_opts_metrop <- theta_opts[,1:q] - theta_opts_metrop[2,1] <- 2 - ############################################## - - m_nn <- 20 - mcmc <- 10000 - - if(F){ - custom_dag <- dag_vecchia(cx_in, m_nn) - - ############################################## - set.seed(1) - RhpcBLASctl::omp_set_num_threads(1) - RhpcBLASctl::blas_set_num_threads(1) - estim_time <- system.time({ - spiox_gibbs_out <- spiox::spiox_wishart(Y_in, X_in, cx_in, - custom_dag = custom_dag, - theta=theta_opts, - - Sigma_start = Sigma %>% perturb(.1, T), - mvreg_B_start = Beta %>% perturb(.1), - - mcmc = mcmc, - print_every = 100, - - sample_iwish=T, - sample_mvr=T, - sample_theta_gibbs=T, - upd_theta_opts=F, - num_threads = nthreads) - }) - - - - predict_dag <- dag_vecchia_predict(cx_in, cx_all[which_out,], m_nn) - - predict_time <- system.time({ - spiox_gibbs_predicts <- spiox::spiox_predict(X_new = X[which_out,], - coords_new = cx_all[which_out,], - - # data - Y_in, X_in, cx_in, - predict_dag, - spiox_gibbs_out$B %>% tail(c(NA, NA, round(mcmc/2))), - spiox_gibbs_out$S %>% tail(c(NA, NA, round(mcmc/2))), - spiox_gibbs_out$theta %>% tail(c(NA, NA, round(mcmc/2))), - num_threads = nthreads) - }) - - Ytest <- spiox_gibbs_predicts$Y %>% stailh() %>% apply(1:2, mean) - Ytrue <- Y[which_out,] - 1:q %>% sapply(\(j) cor(Ytest[,j], Ytrue[,j])) - - Y_spiox_sum_post_mean <- with(spiox_gibbs_predicts, apply(Y[,1,]+Y[,2,], 1, mean)) - sqrt(mean( (Y_spiox_sum_post_mean - Ytrue[,1]-Ytrue[,2])^2 )) - - total_time <- estim_time + predict_time - - save(file=glue::glue("simulations/lmc_m/spiox_gibbs_{oo}.RData"), - list=c("spiox_gibbs_out", "spiox_gibbs_predicts", "estim_time", "predict_time", "total_time")) - rm(list=c("spiox_gibbs_out", "spiox_gibbs_predicts")) - } - - if(F){ - custom_dag <- dag_vecchia(cx_in, m_nn) - - ############################################## - set.seed(1) - RhpcBLASctl::omp_set_num_threads(1) - RhpcBLASctl::blas_set_num_threads(1) - - estim_time <- system.time({ - spiox_metrop_out <- spiox::spiox_wishart(Y_in, X_in, cx_in, - custom_dag = custom_dag, - theta=theta_opts_metrop, - - Sigma_start = Sigma %>% perturb(.1, T), - mvreg_B_start = Beta %>% perturb(.1), - - mcmc = mcmc, - print_every = 100, - - sample_iwish=T, - sample_mvr=T, - sample_theta_gibbs=F, - upd_theta_opts=T, - num_threads = nthreads) - }) - - - - predict_dag <- dag_vecchia_predict(cx_in, cx_all[which_out,], m_nn) - - predict_time <- system.time({ - spiox_metrop_predicts <- spiox::spiox_predict(X_new = X[which_out,], - coords_new = cx_all[which_out,], - - # data - Y_in, X_in, cx_in, - predict_dag, - spiox_metrop_out$B %>% tail(c(NA, NA, round(mcmc/2))), - spiox_metrop_out$S %>% tail(c(NA, NA, round(mcmc/2))), - spiox_metrop_out$theta %>% tail(c(NA, NA, round(mcmc/2))), - num_threads = nthreads) - }) - - Ytest <- spiox_metrop_predicts$Y %>% stailh() %>% apply(1:2, mean) - Ytrue <- Y[which_out,] - 1:q %>% sapply(\(j) cor(Ytest[,j], Ytrue[,j])) - - Y_spiox_sum_post_mean <- with(spiox_metrop_predicts, apply(Y[,1,]+Y[,2,], 1, mean)) - sqrt(mean( (Y_spiox_sum_post_mean - Ytrue[,1]-Ytrue[,2])^2 )) - - total_time <- estim_time + predict_time - - save(file=glue::glue("simulations/lmc_m/spiox_metrop_{oo}.RData"), - list=c("spiox_metrop_out", "spiox_metrop_predicts", "estim_time", "predict_time", "total_time")) - - rm(list=c("spiox_metrop_out", "spiox_metrop_predicts")) - } - - if(F){ - custom_dag <- dag_vecchia(cx_in, m_nn) - - ############################################## - set.seed(1) - RhpcBLASctl::omp_set_num_threads(1) - RhpcBLASctl::blas_set_num_threads(1) - estim_time <- system.time({ - spiox_clust_out <- spiox::spiox_wishart(Y_in, X_in, cx_in, - custom_dag = custom_dag, - theta=theta_opts_metrop[,1:6], - - Sigma_start = Sigma %>% perturb(.1, T), - mvreg_B_start = Beta %>% perturb(.1), - - mcmc = mcmc, - print_every = 100, - - sample_iwish=T, - sample_mvr=T, - sample_theta_gibbs=T, - upd_theta_opts=T, - num_threads = nthreads) - }) - - - - predict_dag <- dag_vecchia_predict(cx_in, cx_all[which_out,], m_nn) - - predict_time <- system.time({ - spiox_clust_predicts <- spiox::spiox_predict(X_new = X[which_out,], - coords_new = cx_all[which_out,], - - # data - Y_in, X_in, cx_in, - predict_dag, - spiox_clust_out$B %>% tail(c(NA, NA, round(mcmc/2))), - spiox_clust_out$S %>% tail(c(NA, NA, round(mcmc/2))), - spiox_clust_out$theta %>% tail(c(NA, NA, round(mcmc/2))), - num_threads = nthreads) - }) - - Ytest <- spiox_clust_predicts$Y %>% stailh() %>% apply(1:2, mean) - Ytrue <- Y[which_out,] - 1:q %>% sapply(\(j) cor(Ytest[,j], Ytrue[,j])) - - Y_spiox_sum_post_mean <- with(spiox_clust_predicts, apply(Y[,1,]+Y[,2,], 1, mean)) - sqrt(mean( (Y_spiox_sum_post_mean - Ytrue[,1]-Ytrue[,2])^2 )) - - total_time <- estim_time + predict_time - - save(file=glue::glue("simulations/lmc_m/spiox_clust_{oo}.RData"), - list=c("spiox_clust_out", "spiox_clust_predicts", "estim_time", "predict_time", "total_time")) - rm(list=c("spiox_clust_out", "spiox_clust_predicts")) - } - - if(F){ - # meshed - library(meshed) - - Y_meshed <- Y - Y_meshed[which_out,] <- NA - - meshed_time <- system.time({ - spmeshed_out <- meshed::spmeshed(y=Y_meshed, x=X, coords=cx_all, k=6, family = "gaussian", - block_size=40, - n_samples = round(mcmc/2), n_thin = 1, n_burn = round(mcmc/2), - n_threads = nthreads, verbose = 10, - predict_everywhere = T, - prior = list(phi=c(10, 60), tausq=c(1e-4,1e-4), nu=c(0.999, 1.001))) - }) - - m_order <- order(spmeshed_out$savedata$coords_blocking$ix) - Ymesh_out <- spmeshed_out$yhat_mcmc %>% tailh() %>% abind::abind(along=3) %>% `[`(m_order[which_out],,) - - save(file=glue::glue("simulations/lmc_m/meshed_{oo}.RData"), - list=c("spmeshed_out", "meshed_time")) - - } - - if(T){ - # nngp - library(spNNGP) - - nngp_time <- system.time({ - - starting <- list("phi"=20, "sigma.sq"=1, "tau.sq"=1e-19, "nu"=1) - tuning <- list("phi"=0.1, "sigma.sq"=0.1, "tau.sq"=0.1, "nu"=0) - priors.1 <- list("beta.Norm"=list(rep(0,ncol(X_in)), diag(1e3,ncol(X_in))), - "phi.Unif"=c(10, 60), "sigma.sq.IG"=c(2, 1), - "nu.Unif"=c(0.999, 1.001), - "tau.sq.IG"=c(1e-3, 1e-3)) - - verbose <- TRUE - n.neighbors <- 10 - mcmc_nngp <- 10000 - burnin <- 1:round(mcmc_nngp/2) - - nngp_results <- list() - for(j in 1:q){ - cat("NNGP ", j, "\n") - m.s.1 <- spNNGP::spNNGP(Y_in[,j] ~ X_in - 1, - coords=cx_in, - starting=starting, method="response", - n.neighbors=n.neighbors, - tuning=tuning, priors=priors.1, cov.model="matern", - n.samples=mcmc_nngp, n.omp.threads=nthreads) - - m.s.1$p.beta.samples %<>% `[`(-burnin,,drop=F) - m.s.1$p.theta.samples %<>% `[`(-burnin,,drop=F) - - nngp_pred.1 <- predict(m.s.1, X[which_out,,drop=F], - coords=cx_all[which_out,], n.omp.threads=nthreads) - - nngp_results[[j]] <- list(i=j, - fitmodel = m.s.1, - predicts = nngp_pred.1) - - } - - }) - - save(file=glue::glue("simulations/lmc_m/nngp_{oo}.RData"), - list=c("burnin", "nngp_results", "nngp_time")) - - - } -} - - diff --git a/simulations/multi_lmc_postprocess.r b/simulations/multi_lmc_postprocess.r deleted file mode 100644 index 761471a..0000000 --- a/simulations/multi_lmc_postprocess.r +++ /dev/null @@ -1,194 +0,0 @@ -#args <- commandArgs(TRUE) - -#i <- args[1] - -library(tidyverse) -library(magrittr) -library(Matrix) -library(spiox) - -RhpcBLASctl::blas_set_num_threads(2) -RhpcBLASctl::omp_set_num_threads(2) - -# many outcomes, parsimonious matern model - -image.matrix <- \(x) { - image(as(x, "dgCMatrix")) -} - -rtail <- \(x, nt){ - tail(x, nt) -} -tailh <- \(x, nt){ - tail(x, round(length(x)/2)) -} -rtailh <- \(x){ rtail(x, round(nrow(x)/2)) } -ctail <- \(x, nt){ - tail(x, c(NA,nt)) -} -ctailh <- \(x){ ctail(x, round(ncol(x)/2)) } -stail <- \(x, nt){ - tail(x, c(NA,NA,nt)) -} -stailh <- \(x){ stail(x, round(dim(x)[3]/2)) } - -perturb <- function(x, sd=1){ - return(x + matrix(rnorm(prod(dim(x)), sd), ncol=ncol(x))) -} - -q <- 24 -g <- glue::glue - - -mcmc <- 10000 -mcmc_keep <- 5000 -nr <- 2500 -nout <- 400 - -Sigbuild <- function(spiox_out){ - result <- 1:10000 %>% sapply(\(i) with(spiox_out, sqrt(diag(theta[2,,i])) %*% Sigma[,,i] %*% sqrt(diag(theta[2,,i])) )) %>% - array(dim=c(q,q,10000)) - return( - result - ) -} - -for(i in 1:20){ - cat(i, "\n") - - load(g("simulations/lmc_m/data_{i}.RData")) - Y <- simdata %>% dplyr::select(contains("Y.")) - Y_out <- Y %>% tail(nout) - Y_out_arr <- 1:mcmc %>% lapply(\(i) Y_out) %>% abind::abind(along=3) - coords_all <- simdata %>% dplyr::select(contains("coords")) - cx <- coords_all %>% head(nr) %>% as.matrix() - - cat("load_results", "\n") - - load(g("simulations/lmc_m/spiox_metrop_{i}.RData")) - time_metrop <- total_time - load(g("simulations/lmc_m/spiox_gibbs_{i}.RData")) - time_gibbs <- total_time - load(g("simulations/lmc_m/spiox_clust_{i}.RData")) - time_clust <- total_time - load(g("simulations/lmc_m/meshed_{i}.RData")) - - meshed_outdata <- simdata %>% mutate(sample = c(rep("insample", nr), rep("outsample", nout))) %>% - left_join(spmeshed_out$coordsdata %>% mutate(ix=1:(nr+nout)), by=c("coords.Var1"="Var1", "coords.Var2"="Var2")) - - meshed_outsample <- meshed_outdata %>% - dplyr::filter(sample=="outsample") %$% ix - - Y_target_meshed <- meshed_outdata %>% dplyr::select(contains("Y.")) %>% `[`(meshed_outsample,) - meshed_predicts <- list(Y = spmeshed_out$yhat_mcmc %>% abind::abind(along=3) %>% `[`(meshed_outsample,,)) - - # meshed_time - load(g("simulations/lmc_m/nngp_{i}.RData")) - nngp_predicts <- list() - nngp_predicts$Y <- array(0, dim=c(400, 24, 5000)) - for(j in 1:24){ - cat(j, " ") - nngp_predicts$Y[,j,] <- nngp_results[[j]]$predicts$p.y.0 - } - nngp_Sigma <- nngp_predicts$Y %>% apply(3, \(x) cor(x)) %>% array(dim=c(q,q,5000)) - nngp_Cor0 <- nngp_Sigma %>% apply(1:2, mean) - cat("\n") - # nngp_time - - - cat("predictions", "\n") - - Y_out_gdiff <- rowSums(Y_out[,1:12]) - rowSums(Y_out[,13:24]) - target <- spiox_clust_predicts - - calc_performance <- function(target){ - perf1 <- 1:q %>% sapply(\(i) sqrt(mean((apply(target$Y[,i,], 1, mean) - Y_out[,i])^2)) ) - - target_gdiff <- apply(target$Y[,1:12,], c(1,3), sum) - apply(target$Y[,13:24,], c(1,3), sum) - perf2 <- sqrt( mean( ( apply(target_gdiff, 1, mean) - Y_out_gdiff)^2 ) ) - - return(list(marginal_rmspe = perf1, - gdiff_rmspe = perf2)) - } - - table_marginals <- \(listx){ - listx %>% lapply(\(x) x$marginal_rmspe) %>% do.call(cbind, .) - } - - table_gdiff <- \(listx){ - listx %>% lapply(\(x) x$gdiff_rmspe) %>% do.call(cbind, .) - } - - perf_spiox_gibbs <- calc_performance(spiox_gibbs_predicts) - perf_spiox_metrop <- calc_performance(spiox_metrop_predicts) - perf_spiox_clust <- calc_performance(spiox_clust_predicts) - perf_spmeshed <- calc_performance(meshed_predicts) - perf_nngp <- calc_performance(nngp_predicts) - - cat("cor at zero", "\n") - - spiox_cor_at_zero <- function(spiox_out){ - spiox_theta <- spiox_out$theta - - Sig <- Sigbuild(spiox_out) - - spiox_theta[2,,] <- 1 - spiox_scaling_factors <- scaling_factor_at_zero(cx, spiox_theta %>% tail(c(NA, NA, 5000)) %>% apply(1:2, mean)) - spiox_scaling_factors[upper.tri(spiox_scaling_factors)] <- spiox_scaling_factors[lower.tri(spiox_scaling_factors)] - - return( - Sig %>% apply(3, \(s) cov2cor(s * spiox_scaling_factors)) %>% array(dim=c(q,q,10000)) ) - } - - - theta_mat <- rbind(optlist, 1, 1, 1e-15) - - spiox_scaling_factors <- scaling_factor_at_zero(cx, theta_mat) - spiox_scaling_factors[upper.tri(spiox_scaling_factors)] <- spiox_scaling_factors[lower.tri(spiox_scaling_factors)] - Cor_at_zero <- - cov2cor(Sigma * spiox_scaling_factors) - - spiox_gibbs_Cor0 <- spiox_cor_at_zero(spiox_gibbs_out) %>% tail(c(NA, NA, 5000)) %>% apply(1:2, mean) - spiox_metrop_Cor0 <- spiox_cor_at_zero(spiox_metrop_out) %>% tail(c(NA, NA, 5000)) %>% apply(1:2, mean) - spiox_clust_Cor0 <- spiox_cor_at_zero(spiox_clust_out) %>% tail(c(NA, NA, 5000)) %>% apply(1:2, mean) - - omega_spmeshed <- 1:nr %>% lapply(\(m) with(spmeshed_out, - cov2cor(tcrossprod(lambda_mcmc[,,m]) + diag(tausq_mcmc[,m])))) %>% - abind::abind(along=3) %>% tail(c(NA, NA, 5000)) %>% apply(1:2, mean) - - - - results_gdiff <- table_gdiff(list(perf_spiox_gibbs, perf_spiox_metrop, perf_spiox_clust, perf_spmeshed, perf_nngp)) %>% - as.data.frame() %>% mutate(sim = i) - results_margs <- table_marginals(list(perf_spiox_gibbs, perf_spiox_metrop, perf_spiox_clust, perf_spmeshed, perf_nngp)) %>% - as.data.frame() %>% mutate(sim = i, variable=1:n()) - colnames(results_gdiff)[1:5] <- colnames(results_margs)[1:5] <- c("spiox_gibbs", "spiox_metrop", "spiox_clust", "spmeshed", "nngp") - - mat_lt <- function(mat){ - return( mat %>% `[`(lower.tri(.)) ) - } - - results_estcor <- data.frame(spiox_gibbs = mat_lt(spiox_gibbs_Cor0), - spiox_metrop = mat_lt(spiox_metrop_Cor0), - spiox_clust = mat_lt(spiox_clust_Cor0), - spmeshed = mat_lt(omega_spmeshed), - nngp = mat_lt(nngp_Cor0), - true = mat_lt(Cor_at_zero), - sim = i) - - - timings <- list(sim=i, spiox_gibbs = time_gibbs["elapsed"], - spiox_metrop = time_metrop["elapsed"], - spiox_clust = time_clust["elapsed"], - meshed = meshed_time["elapsed"], - nngp = nngp_time["elapsed"]) - - results <- list(sim=i, timings=timings, - gdiff=results_gdiff, - margs=results_margs, - cor0 = results_estcor) - - - save(file=g("simulations/lmc_m/results_{i}.RData"), list=c("results")) - -} \ No newline at end of file diff --git a/simulations/multi_spiox.r b/simulations/multi_spiox.r deleted file mode 100644 index b4c08c7..0000000 --- a/simulations/multi_spiox.r +++ /dev/null @@ -1,374 +0,0 @@ -args <- commandArgs(TRUE) - -starts <- args[1] -ends <- args[2] - - -cat("Simulations ", starts:ends, "\n") - -library(tidyverse) -library(magrittr) -library(Matrix) -library(spiox) - -# many outcomes, parsimonious matern model - -image.matrix <- \(x) { - image(as(x, "dgCMatrix")) -} - -rtail <- \(x, nt){ - tail(x, nt) -} -tailh <- \(x, nt){ - tail(x, round(length(x)/2)) -} -rtailh <- \(x){ rtail(x, round(nrow(x)/2)) } -ctail <- \(x, nt){ - tail(x, c(NA,nt)) -} -ctailh <- \(x){ ctail(x, round(ncol(x)/2)) } -stail <- \(x, nt){ - tail(x, c(NA,NA,nt)) -} -stailh <- \(x){ stail(x, round(dim(x)[3]/2)) } - -perturb <- function(x, sd=1, symm=F){ - pert <- matrix(rnorm(prod(dim(x)), sd), ncol=ncol(x)) - if(symm){ - pert <- crossprod(pert) - } - - return(x + pert) -} - -q <- 24 - -nthreads <- 14 -oo <- 1 - -for(oo in starts:ends){ - set.seed(oo) - cat(oo, "\n") - - optlist <- seq(0.5, 2, length.out=6) %>% sample(q, replace=T) - - # spatial - cx_in <- matrix(runif(2500*2), ncol=2) #2500 - colnames(cx_in) <- c("Var1","Var2") - n_in <- nrow(cx_in) - which_in <- 1:n_in - - xout <- seq(0, 1, length.out=20) #20 - coords_out <- expand.grid(xout, xout) - cx_out <- as.matrix(coords_out) - n_out <- nrow(cx_out) - which_out <- (n_in+1):(n_in+n_out) - - cx_all <- rbind(cx_in, cx_out) - nr_all <- nrow(cx_all) - - Clist <- optlist %>% lapply(\(nu) spiox::Correlationc(cx_all, cx_all, c(30,1,nu,1e-3), 1, TRUE) ) - Llist <- Clist %>% lapply(\(C) t(chol(C))) - - Q <- rWishart(1, q+1, 1/2 * diag(q))[,,1] # - Sigma <- solve(Q) - St <- chol(Sigma) - S <- t(St) - - V <- matrix(rnorm(nr_all * q), ncol=q) %*% St - Y_sp <- V - for(i in 1:q){ - Y_sp[,i] <- Llist[[i]] %*% V[,i] - } - - # regression - p <- 2 - X <- matrix(1, ncol=1, nrow=nr_all) %>% cbind(matrix(rnorm(nr_all*(p-1)), ncol=p-1)) - - Beta <- 0*matrix(rnorm(q * p), ncol=q) - - Y_regression <- X %*% Beta - - Y <- as.matrix(Y_sp + Y_regression) - - Y_in <- Y[which_in,] - X_in <- X[which_in,] - - if(F){ - df <- data.frame(cx_out, y=as.matrix(Y_sp[which_out,])) %>% - pivot_longer(cols=-c(Var1, Var2)) - ggplot(df, - aes(Var1, Var2, fill=value)) + - geom_raster() + - scale_fill_viridis_c() + - facet_wrap(~name, ncol=5) - } - - simdata <- data.frame(coords=cx_all, Y_spatial=Y_sp, Y=Y, X=X) - #ggplot(simdata, aes(coords.Var1, coords.Var2, color=Y_spatial.1)) + geom_point() + scale_color_viridis_c() - - if(T){ - save(file=glue::glue("simulations/spiox_m/data_{oo}.RData"), - list=c("simdata", "oo", "Sigma", "optlist")) - } - - ############################## - set.seed(1) - - theta_opts <- rbind(30, 1, seq(0.5, 2, length.out=30), 1e-3) - - theta_opts_metrop <- rbind(30, 1, seq(0.5, 2, length.out=q), 1e-3) - theta_opts_metrop[2,1] <- 2 - - ############################################## - - m_nn <- 15 - mcmc <- 10000 - - if(T){ - custom_dag <- dag_vecchia(cx_in, m_nn) - - ############################################## - set.seed(1) - RhpcBLASctl::omp_set_num_threads(1) - RhpcBLASctl::blas_set_num_threads(1) - estim_time <- system.time({ - spiox_gibbs_out <- spiox::spiox_wishart(Y_in, X_in, cx_in, - custom_dag = custom_dag, - theta=theta_opts, - - Sigma_start = Sigma, - mvreg_B_start = Beta, - - mcmc = mcmc, - print_every = 100, - matern = TRUE, - sample_iwish=T, - sample_mvr=T, - sample_theta_gibbs=T, - upd_theta_opts=F, - num_threads = nthreads) - }) - - predict_dag <- dag_vecchia_predict(cx_in, cx_all[which_out,], m_nn) - - predict_time <- system.time({ - spiox_gibbs_predicts <- spiox::spiox_predict(X_new = X[which_out,], - coords_new = cx_all[which_out,], - - # data - Y_in, X_in, cx_in, - predict_dag, - spiox_gibbs_out$B %>% tail(c(NA, NA, round(mcmc/2))), - spiox_gibbs_out$Sigma %>% tail(c(NA, NA, round(mcmc/2))), - spiox_gibbs_out$theta %>% tail(c(NA, NA, round(mcmc/2))), - matern = TRUE, - num_threads = nthreads) - }) - - Ytest <- spiox_gibbs_predicts$Y %>% stailh() %>% apply(1:2, mean) - Ytrue <- Y[which_out,] - 1:q %>% sapply(\(j) cor(Ytest[,j], Ytrue[,j])) - - Y_spiox_sum_post_mean <- with(spiox_gibbs_predicts, apply(Y[,1,]+Y[,2,], 1, mean)) - sqrt(mean( (Y_spiox_sum_post_mean - Ytrue[,1]-Ytrue[,2])^2 )) - - total_time <- estim_time + predict_time - - save(file=glue::glue("simulations/spiox_m/spiox_gibbs_{oo}.RData"), - list=c("spiox_gibbs_out", "spiox_gibbs_predicts", "estim_time", "predict_time", "total_time")) - rm(list=c("spiox_gibbs_out", "spiox_gibbs_predicts")) - } - - if(F){ - custom_dag <- dag_vecchia(cx_in, m_nn) - - ############################################## - set.seed(1) - RhpcBLASctl::omp_set_num_threads(1) - RhpcBLASctl::blas_set_num_threads(1) - - estim_time <- system.time({ - spiox_metrop_out <- spiox::spiox_wishart(Y_in, X_in, cx_in, - custom_dag = custom_dag, - theta=theta_opts_metrop, - - Sigma_start = Sigma %>% perturb(.1, T), - mvreg_B_start = Beta %>% perturb(.1), - - mcmc = mcmc, - print_every = 100, - matern = TRUE, - sample_iwish=T, - sample_mvr=T, - sample_theta_gibbs=F, - upd_theta_opts=T, - num_threads = nthreads) - }) - - - - predict_dag <- dag_vecchia_predict(cx_in, cx_all[which_out,], m_nn) - - predict_time <- system.time({ - spiox_metrop_predicts <- spiox::spiox_predict(X_new = X[which_out,], - coords_new = cx_all[which_out,], - - # data - Y_in, X_in, cx_in, - predict_dag, - spiox_metrop_out$B %>% tail(c(NA, NA, round(mcmc/2))), - spiox_metrop_out$Sigma %>% tail(c(NA, NA, round(mcmc/2))), - spiox_metrop_out$theta %>% tail(c(NA, NA, round(mcmc/2))), - num_threads = nthreads) - }) - - Ytest <- spiox_metrop_predicts$Y %>% stailh() %>% apply(1:2, mean) - Ytrue <- Y[which_out,] - 1:q %>% sapply(\(j) cor(Ytest[,j], Ytrue[,j])) - - Y_spiox_sum_post_mean <- with(spiox_metrop_predicts, apply(Y[,1,]+Y[,2,], 1, mean)) - sqrt(mean( (Y_spiox_sum_post_mean - Ytrue[,1]-Ytrue[,2])^2 )) - - total_time <- estim_time + predict_time - - save(file=glue::glue("simulations/spiox_m/spiox_metrop_{oo}.RData"), - list=c("spiox_metrop_out", "spiox_metrop_predicts", "estim_time", "predict_time", "total_time")) - - rm(list=c("spiox_metrop_out", "spiox_metrop_predicts")) - } - - if(T){ - custom_dag <- dag_vecchia(cx_in, m_nn) - - ############################################## - set.seed(1) - RhpcBLASctl::omp_set_num_threads(1) - RhpcBLASctl::blas_set_num_threads(1) - estim_time <- system.time({ - spiox_clust_out <- spiox::spiox_wishart(Y_in, X_in, cx_in, - custom_dag = custom_dag, - theta=theta_opts_metrop[,1:6], - - Sigma_start = Sigma, - mvreg_B_start = Beta, - - mcmc = mcmc, - print_every = 100, - matern = TRUE, - sample_iwish = T, - sample_mvr = T, - sample_theta_gibbs = T, - upd_theta_opts = T, - num_threads = nthreads) - }) - - predict_dag <- dag_vecchia_predict(cx_in, cx_all[which_out,], m_nn) - - predict_time <- system.time({ - spiox_clust_predicts <- spiox::spiox_predict(X_new = X[which_out,], - coords_new = cx_all[which_out,], - - # data - Y_in, X_in, cx_in, - predict_dag, - spiox_clust_out$B %>% tail(c(NA, NA, round(mcmc/2))), - spiox_clust_out$S %>% tail(c(NA, NA, round(mcmc/2))), - spiox_clust_out$theta %>% tail(c(NA, NA, round(mcmc/2))), - num_threads = nthreads) - }) - - Ytest <- spiox_clust_predicts$Y %>% stailh() %>% apply(1:2, mean) - Ytrue <- Y[which_out,] - 1:q %>% sapply(\(j) cor(Ytest[,j], Ytrue[,j])) - - Y_spiox_sum_post_mean <- with(spiox_clust_predicts, apply(Y[,1,]+Y[,2,], 1, mean)) - sqrt(mean( (Y_spiox_sum_post_mean - Ytrue[,1]-Ytrue[,2])^2 )) - - total_time <- estim_time + predict_time - - save(file=glue::glue("simulations/spiox_m/spiox_clust_{oo}.RData"), - list=c("spiox_clust_out", "spiox_clust_predicts", "estim_time", "predict_time", "total_time")) - rm(list=c("spiox_clust_out", "spiox_clust_predicts")) - } - - if(T){ - # meshed - library(meshed) - - Y_meshed <- Y - Y_meshed[which_out,] <- NA - - meshed_time <- system.time({ - spmeshed_out <- meshed::spmeshed(y=Y_meshed, x=X, coords=cx_all, k=6, family = "gaussian", - block_size=40, - n_samples = round(mcmc/2), n_thin = 1, n_burn = round(mcmc/2), - n_threads = nthreads, verbose = 10, - predict_everywhere = T, - prior = list(phi=c(10, 50), tausq=c(1e-4,1e-4), nu=c(.5, 2))) - }) - - m_order <- order(spmeshed_out$savedata$coords_blocking$ix) - Ymesh_out <- spmeshed_out$yhat_mcmc %>% tailh() %>% abind::abind(along=3) %>% `[`(m_order[which_out],,) - - Ymesh_mean <- Ymesh_out %>% apply(1:2, mean) - 1:q %>% sapply(\(j) cor(Ymesh_mean[,j], Ytrue[,j])) - Y_meshed_sum_post_mean <- apply(Ymesh_out[,1,]+Ymesh_out[,2,], 1, mean) - sqrt(mean( (Y_meshed_sum_post_mean - Ytrue[,1]-Ytrue[,2])^2 )) - - save(file=glue::glue("simulations/spiox_m/meshed_{oo}.RData"), - list=c("spmeshed_out", "meshed_time")) - - } - - if(T){ - # nngp - library(spNNGP) - - nngp_time <- system.time({ - - starting <- list("phi"=30, "sigma.sq"=1, "tau.sq"=1e-19, "nu"=1) - tuning <- list("phi"=0, "sigma.sq"=0.1, "tau.sq"=0.1, "nu"=0.1) - priors.1 <- list("beta.Norm"=list(rep(0,ncol(X_in)), diag(1e3,ncol(X_in))), - "phi.Unif"=c(10, 60), "sigma.sq.IG"=c(2, 1), - "nu.Unif"=c(0.4, 2), - "tau.sq.IG"=c(2, 1)) - - verbose <- TRUE - n.neighbors <- 15 - mcmc_nngp <- 10000 - burnin <- 1:round(mcmc_nngp/2) - - nngp_results <- list() - for(j in 1:q){ - cat("NNGP ", j, "\n") - m.s.1 <- spNNGP::spNNGP(Y_in[,j] ~ X_in - 1, - coords=cx_in, - starting=starting, method="response", - n.neighbors=n.neighbors, - tuning=tuning, priors=priors.1, cov.model="matern", - n.samples=mcmc_nngp, n.omp.threads=nthreads) - - m.s.1$p.beta.samples %<>% `[`(-burnin,,drop=F) - m.s.1$p.theta.samples %<>% `[`(-burnin,,drop=F) - - nngp_pred.1 <- predict(m.s.1, X[which_out,,drop=F], - coords=cx_all[which_out,], n.omp.threads=nthreads) - - nngp_results[[j]] <- list(i=j, - fitmodel = m.s.1, - predicts = nngp_pred.1) - - } - - }) - - save(file=glue::glue("simulations/spiox_m/nngp_{oo}.RData"), - list=c("burnin", "nngp_results", "nngp_time")) - - - } -} - - diff --git a/simulations/multi_spiox_postprocess.r b/simulations/multi_spiox_postprocess.r deleted file mode 100644 index ca24203..0000000 --- a/simulations/multi_spiox_postprocess.r +++ /dev/null @@ -1,216 +0,0 @@ - -#args <- commandArgs(TRUE) - - -library(tidyverse) -library(magrittr) -library(Matrix) -library(spiox) - -# many outcomes, parsimonious matern model - -image.matrix <- \(x) { - image(as(x, "dgCMatrix")) -} - -rtail <- \(x, nt){ - tail(x, nt) -} -tailh <- \(x, nt){ - tail(x, round(length(x)/2)) -} -rtailh <- \(x){ rtail(x, round(nrow(x)/2)) } -ctail <- \(x, nt){ - tail(x, c(NA,nt)) -} -ctailh <- \(x){ ctail(x, round(ncol(x)/2)) } -stail <- \(x, nt){ - tail(x, c(NA,NA,nt)) -} -stailh <- \(x){ stail(x, round(dim(x)[3]/2)) } - -perturb <- function(x, sd=1){ - return(x + matrix(rnorm(prod(dim(x)), sd), ncol=ncol(x))) -} - -q <- 24 -g <- glue::glue - -mcmc <- 10000 -mcmc_keep <- 5000 -nr <- 2500 -nout <- 400 - -results <- list() - - -Sigbuild <- function(spiox_out){ - result <- 1:10000 %>% sapply(\(i) with(spiox_out, sqrt(diag(theta[2,,i])) %*% Sigma[,,i] %*% sqrt(diag(theta[2,,i])) )) %>% - array(dim=c(q,q,10000)) - return( - result - ) -} - -for(i in 1:20){ - - cat(i, "\n") - - load(g("simulations/spiox_m/data_{i}.RData")) - Y <- simdata %>% dplyr::select(contains("Y.")) - Y_out <- Y %>% tail(nout) - Y_out_arr <- 1:mcmc %>% lapply(\(i) Y_out) %>% abind::abind(along=3) - coords_all <- simdata %>% dplyr::select(contains("coords")) - cx <- coords_all %>% head(nr) %>% as.matrix() - - cat("load results", "\n") - #load(g("simulations/spiox_m/spiox_metrop_{i}.RData")) - #time_metrop <- total_time - load(g("simulations/spiox_m/spiox_gibbs_{i}.RData")) - time_gibbs <- total_time - load(g("simulations/spiox_m/spiox_clust_{i}.RData")) - time_clust <- total_time - load(g("simulations/spiox_m/meshed_{i}.RData")) - meshed_outdata <- simdata %>% mutate(sample = c(rep("insample", nr), rep("outsample", nout))) %>% - left_join(spmeshed_out$coordsdata %>% mutate(ix=1:(nr+nout)), by=c("coords.Var1"="Var1", "coords.Var2"="Var2")) - - meshed_outsample <- meshed_outdata %>% - dplyr::filter(sample=="outsample") %$% ix - - Y_target_meshed <- meshed_outdata %>% dplyr::select(contains("Y.")) %>% `[`(meshed_outsample,) - meshed_predicts <- list(Y = spmeshed_out$yhat_mcmc %>% abind::abind(along=3) %>% `[`(meshed_outsample,,)) - - load(g("simulations/spiox_m/nngp_{i}.RData")) - nngp_predicts <- list() - nngp_predicts$Y <- array(0, dim=c(400, 24, 5000)) - nngp_nu <- rep(0, q) - for(j in 1:q){ - cat(j, " ") - nngp_predicts$Y[,j,] <- nngp_results[[j]]$predicts$p.y.0 - nngp_nu[j] <- nngp_results[[j]]$fitmodel$p.theta.samples[,"nu"] %>% mean() - } - nngp_Sigma <- nngp_predicts$Y %>% apply(3, \(x) cor(x)) %>% array(dim=c(q,q,5000)) - nngp_Cor0 <- nngp_Sigma %>% apply(1:2, mean) - cat("\n") - - cat("predictions", "\n") - Y_out_gdiff <- rowSums(Y_out[,1:12]) - rowSums(Y_out[,13:24]) - target <- spiox_clust_predicts - - calc_performance <- function(target){ - perf1 <- 1:q %>% sapply(\(i) sqrt(mean((apply(target$Y[,i,], 1, mean) - Y_out[,i])^2)) ) - - target_gdiff <- apply(target$Y[,1:12,], c(1,3), sum) - apply(target$Y[,13:24,], c(1,3), sum) - perf2 <- sqrt( mean( ( apply(target_gdiff, 1, mean) - Y_out_gdiff)^2 ) ) - - return(list(marginal_rmspe = perf1, - gdiff_rmspe = perf2)) - } - - table_marginals <- \(listx){ - listx %>% lapply(\(x) x$marginal_rmspe) %>% do.call(cbind, .) - } - - table_gdiff <- \(listx){ - listx %>% lapply(\(x) x$gdiff_rmspe) %>% do.call(cbind, .) - } - - perf_spiox_gibbs <- calc_performance(spiox_gibbs_predicts) - #perf_spiox_metrop <- calc_performance(spiox_metrop_predicts) - perf_spiox_clust <- calc_performance(spiox_clust_predicts) - perf_nngp <- calc_performance(nngp_predicts) - perf_spmeshed <- calc_performance(meshed_predicts) - - cat("cor at zero", "\n") - spiox_cor_at_zero <- function(spiox_out){ - spiox_theta <- spiox_out$theta - - Sig <- Sigbuild(spiox_out) - - spiox_theta[2,,] <- 1 - spiox_scaling_factors <- scaling_factor_at_zero(cx, spiox_theta %>% tail(c(NA, NA, 5000)) %>% apply(1:2, mean)) - spiox_scaling_factors[upper.tri(spiox_scaling_factors)] <- spiox_scaling_factors[lower.tri(spiox_scaling_factors)] - - return( - Sig %>% apply(3, \(s) cov2cor(s * spiox_scaling_factors)) %>% array(dim=c(q,q,10000)) ) - } - - - theta_mat <- rbind(30, 1, optlist, 1e-3) - - spiox_scaling_factors <- scaling_factor_at_zero(cx, theta_mat) - spiox_scaling_factors[upper.tri(spiox_scaling_factors)] <- spiox_scaling_factors[lower.tri(spiox_scaling_factors)] - Cor_at_zero <- - cov2cor(Sigma * spiox_scaling_factors) - - spiox_gibbs_Cor0 <- spiox_cor_at_zero(spiox_gibbs_out) %>% tail(c(NA, NA, 5000)) %>% apply(1:2, mean) - #spiox_metrop_Cor0 <- spiox_cor_at_zero(spiox_metrop_out) %>% tail(c(NA, NA, 5000)) %>% apply(1:2, mean) - spiox_clust_Cor0 <- spiox_cor_at_zero(spiox_clust_out) %>% tail(c(NA, NA, 5000)) %>% apply(1:2, mean) - - omega_spmeshed <- 1:nr %>% lapply(\(m) with(spmeshed_out, - cov2cor(tcrossprod(lambda_mcmc[,,m]) + diag(tausq_mcmc[,m])))) %>% - abind::abind(along=3) %>% tail(c(NA, NA, 5000)) %>% apply(1:2, mean) - - results_gdiff <- table_gdiff(list(perf_spiox_gibbs, #perf_spiox_metrop, - perf_spiox_clust, perf_spmeshed, perf_nngp)) %>% - as.data.frame() %>% mutate(sim = i) - results_margs <- table_marginals(list(perf_spiox_gibbs, #perf_spiox_metrop, - perf_spiox_clust, perf_spmeshed, perf_nngp)) %>% - as.data.frame() %>% mutate(sim = i, variable=1:n()) - colnames(results_gdiff)[1:4] <- colnames(results_margs)[1:4] <- - c("spiox_gibbs", #"spiox_metrop", - "spiox_clust", "spmeshed", "nngp") - - mat_lt <- function(mat){ - return( mat %>% `[`(lower.tri(.)) ) - } - - extract_nu <- function(theta_mcmc){ - return( - theta_mcmc %>% tail(c(NA, NA, mcmc_keep)) %>% apply(1:2, mean) %>% `[`(3,) - ) - } - - results_estcor <- data.frame(spiox_gibbs = mat_lt(spiox_gibbs_Cor0), - #spiox_metrop = mat_lt(spiox_metrop_Cor0), - spiox_clust = mat_lt(spiox_clust_Cor0), - spmeshed = mat_lt(omega_spmeshed), - nngp = mat_lt(nngp_Cor0), - true = mat_lt(Cor_at_zero), - sim = i) - - results_estnu <- data.frame(spiox_gibbs = spiox_gibbs_out$theta %>% extract_nu(), - #spiox_metrop = spiox_metrop_out$theta %>% extract_nu(), - spiox_clust = spiox_clust_out$theta %>% extract_nu(), - nngp = nngp_nu, - spmeshed = NA, - true = optlist, - sim = i) - - timings <- list(sim=i, spiox_gibbs = time_gibbs["elapsed"], - #spiox_metrop = time_metrop["elapsed"], - spiox_clust = time_clust["elapsed"], - meshed = meshed_time["elapsed"], - nngp = nngp_time["elapsed"]) - - results <- list(sim=i, timings = timings, - gdiff=results_gdiff, - margs=results_margs, - cor0 = results_estcor, - nu = results_estnu) - - save(file=g("simulations/spiox_m/results_{i}.RData"), list=c("results")) - - -} - -for(i in 1:20){ - load(g("simulations/spiox_m/results_{i}.RData")) - colnames(results$gdiff)[1:4] <- colnames(results$margs)[1:4] <- - c("spiox_gibbs", #"spiox_metrop", - "spiox_clust", "spmeshed", "nngp") - save(file=g("simulations/spiox_m/results_{i}.RData"), list=c("results")) -} - - - diff --git a/simulations/trivariate_matern_toy.r b/simulations/trivariate_matern_toy.r deleted file mode 100644 index 8334935..0000000 --- a/simulations/trivariate_matern_toy.r +++ /dev/null @@ -1,491 +0,0 @@ -rm(list=ls()) -library(tidyverse) -library(magrittr) -library(Matrix) -library(spiox) - -image.matrix <- \(x) { - image(as(x, "dgCMatrix")) -} - -rtail <- \(x, nt){ - tail(x, nt) -} -tailh <- \(x, nt){ - tail(x, round(length(x)/2)) -} -rtailh <- \(x){ rtail(x, round(nrow(x)/2)) } -ctail <- \(x, nt){ - tail(x, c(NA,nt)) -} -ctailh <- \(x){ ctail(x, round(ncol(x)/2)) } -stail <- \(x, nt){ - tail(x, c(NA,NA,nt)) -} -stailh <- \(x){ stail(x, round(dim(x)[3]/2)) } - -perturb <- function(x, sd=1){ - return(x + matrix(rnorm(prod(dim(x)), sd), ncol=ncol(x))) -} - - -for(ii in 1:3){ - - set.seed(ii) - - q <- 3 - nu1 <- 0.5 - nu2 <- 0.8 - nu3 <- 1.2 - - sds <- c(1, 1, 1) - Omega <- matrix(c(1,-.9,0.7,-.9,1,-.5,0.7,-.5,1), ncol=3) - Sigma <- diag(sds) %*% Omega %*% diag(sds) - - St <- chol(Sigma) - S <- t(St) - - # spatial - cx_in <- matrix(runif(2500*2), ncol=2) #2500 - colnames(cx_in) <- c("Var1","Var2") - n_in <- nrow(cx_in) - which_in <- 1:n_in - - xout <- seq(0, 1, length.out=50) #50 - coords_out <- expand.grid(xout, xout) - cx_out <- as.matrix(coords_out) - n_out <- nrow(cx_out) - which_out <- (n_in+1):(n_in+n_out) - - cx_all <- rbind(cx_in, cx_out) - nr_all <- nrow(cx_all) - - - cx_12 <- rbind(cbind(cx_all, 1), cbind(cx_all, 2), cbind(cx_all, 3)) - - V <- Omega #cbind(c(1, rhocorr), c(rhocorr, 1)) - St <- chol(V) - - phi <- 30 - sigma11 <- sigma22 <- sigma33 <- 1 - sigma12 <- sqrt(sigma11 * sigma22) * V[2,1] * - sqrt(gamma(nu1+1))/sqrt(gamma(nu1)) * - sqrt(gamma(nu2+1))/sqrt(gamma(nu2)) * gamma(nu1/2 + nu2/2)/gamma(nu1/2 + nu2/2 + 1) - sigma13 <- sqrt(sigma11 * sigma33) * V[3,1] * - sqrt(gamma(nu1+1))/sqrt(gamma(nu1)) * - sqrt(gamma(nu3+1))/sqrt(gamma(nu3)) * gamma(nu1/2 + nu3/2)/gamma(nu1/2 + nu3/2 + 1) - sigma23 <- sqrt(sigma22 * sigma33) * V[3,2] * - sqrt(gamma(nu2+1))/sqrt(gamma(nu2)) * - sqrt(gamma(nu3+1))/sqrt(gamma(nu3)) * gamma(nu2/2 + nu3/2)/gamma(nu2/2 + nu3/2 + 1) - - if(T){ - - # parsimonious multi matern - Cparmat <- GpGpm::matern_multi(c(sigma11, sigma12, sigma22, sigma13, sigma23, sigma33, - 1/phi, 1/phi, 1/phi, 1/phi, 1/phi, 1/phi, - nu1, nu1/2+nu2/2, nu2, nu1/2+nu3/2, nu2/2+nu3/2, nu3, 1e-3, 1e-3, 1e-3, 1e-3, 1e-3, 1e-3), cx_12) - - Lparmat <- t(chol(Cparmat)) - - wvec <- Lparmat %*% rnorm(nr_all*q) - W <- matrix(wvec, ncol=q) - - # matern - Y_sp <- W - - # regression - p <- 1 - X <- matrix(1, ncol=1, nrow=nr_all) #%>% cbind(matrix(rnorm(nr_all*(p-1)), ncol=p-1)) - - Beta <- matrix(rnorm(q * p), ncol=q) - - Y_regression <- X %*% Beta - #Error <- matrix(rnorm(nrow(Y_regression) * q),ncol=q) %*% diag(D <- runif(q, 0, 0.1)) - Y <- as.matrix(Y_sp + Y_regression) # + Error - - Y_in <- Y[which_in,] - X_in <- X[which_in,,drop=F] - - if(F){ - df <- data.frame(cx_out, y=as.matrix(Y_sp[which_out,])) %>% - pivot_longer(cols=-c(Var1, Var2)) - ggplot(df, - aes(Var1, Var2, fill=value)) + - geom_raster() + - scale_fill_viridis_c() + - facet_wrap(~name, ncol=5) - } - - simdata <- data.frame(coords=cx_all, Y_spatial=Y_sp, Y=Y, X=X) - #ggplot(simdata, aes(coords.Var1, coords.Var2, color=Y_spatial.1)) + geom_point() + scale_color_viridis_c() - - save(file=glue::glue("simulations/trivariate_matern_toy_reps/data_{ii}.RData"), - list=c("simdata", - "Beta", "D", "Y_in", "X_in", "which_in", "which_out", - "Y_regression", #"Error", - "Y", "X", "W", "Lparmat", - "nu1", "nu2", "nu3", "phi", "Sigma")) - } else { - load(glue::glue("simulations/trivariate_matern_toy_reps/data_{ii}.RData")) - } - ############################## - - set.seed(1) - - #theta_opts <- cbind(c(10, 1, .5, 1e-19), c(20, 1, 1, 1e-19)) - theta_opts <- cbind(c(phi+1e-8, 1.001, nu1, 1e-3), c(phi-1e-8, 0.999, nu2, 0.0011), c(phi, 1, nu3, 1e-3)) - - ############################################## - - m_nn <- 20 - mcmc <- 10000 - RhpcBLASctl::blas_set_num_threads(1) - RhpcBLASctl::omp_set_num_threads(1) - - if(T){ - - custom_dag <- dag_vecchia(cx_in, m_nn) - - ############################################## - set.seed(1) - estim_time <- system.time({ - spiox_out <- spiox::spiox_wishart(Y_in, X_in, cx_in, - custom_dag = custom_dag, - theta=theta_opts, - - Sigma_start = Sigma, - mvreg_B_start = Beta,# %>% perturb(), - - mcmc = mcmc, - print_every = 100, - matern = TRUE, - sample_iwish=T, - sample_mvr=T, - sample_theta_gibbs=F, - upd_theta_opts=T, - num_threads = 8) - }) - - predict_dag <- dag_vecchia_predict(cx_in, cx_all[which_out,], m_nn) - - predict_time <- system.time({ - spiox_predicts <- spiox::spiox_predict(X_new = X[which_out,,drop=F], - coords_new = cx_all[which_out,], - - # data - Y_in, X_in, cx_in, - predict_dag, - spiox_out$B %>% tail(c(NA, NA, round(mcmc/2))), - spiox_out$Sigma %>% tail(c(NA, NA, round(mcmc/2))), - spiox_out$theta %>% tail(c(NA, NA, round(mcmc/2))), - matern = TRUE, - num_threads = 16) - }) - - - Y_out <- Y[which_out,] - Yhat_spiox <- spiox_predicts$Y %>% apply(1:2, mean) - - - - - total_time <- estim_time + predict_time - - save(file=glue::glue("simulations/trivariate_matern_toy_reps/spiox_{ii}.RData"), - list=c("spiox_out", "spiox_predicts", "estim_time", "predict_time", "total_time")) - } - - if(T){ - custom_dag <- dag_vecchia(cx_in, m_nn) - theta_opts_latent <- theta_opts - theta_opts_latent[2,1] <- 2 - theta_opts_latent[4,] <- 1e-10 - ############################################## - set.seed(1) - latentq_estim_time <- system.time({ - spiox_latentq_out <- spiox::spiox_latent(Y_in, X_in, cx_in, - custom_dag = custom_dag, - theta=theta_opts_latent, - - Sigma_start = Sigma, - mvreg_B_start = Beta,# %>% perturb(), - - mcmc = mcmc, - print_every = 100, - matern = TRUE, - sample_iwish=T, - sample_mvr=T, - sample_theta_gibbs=F, - upd_theta_opts=T, - num_threads = 8, - sampling = 3) - }) - - predict_dag <- dag_vecchia_predict(cx_in, cx_all[which_out,], m_nn) - - latentq_predict_time <- system.time({ - spiox_latentq_predicts <- - spiox_latentq_out %>% with( - spiox::spiox_latent_predict(X_new = X[which_out,,drop=F], - coords_new = cx_all[which_out,], - - # data - cx_in, - predict_dag, - W %>% tail(c(NA, NA, round(mcmc/2))), - B %>% tail(c(NA, NA, round(mcmc/2))), - Sigma %>% tail(c(NA, NA, round(mcmc/2))), - Ddiag %>% tail(c(NA, round(mcmc/2))), - theta %>% tail(c(NA, NA, round(mcmc/2))), - matern = TRUE, - num_threads = 16) ) - }) - - save(file=glue::glue("simulations/trivariate_matern_toy_reps/spiox_latentq_{ii}.RData"), - list=c("spiox_latentq_out", "spiox_latentq_predicts", "latentq_estim_time", "latentq_predict_time")) - - } - - if(T){ - - custom_dag <- dag_vecchia(cx_in, m_nn) - theta_opts_latent <- theta_opts - theta_opts_latent[2,1] <- 2 - theta_opts_latent[4,] <- 1e-10 - ############################################## - set.seed(1) - latentn_estim_time <- system.time({ - spiox_latentn_out <- spiox::spiox_latent(Y_in, X_in, cx_in, - custom_dag = custom_dag, - theta=theta_opts_latent, - - Sigma_start = Sigma, - mvreg_B_start = Beta,# %>% perturb(), - - mcmc = mcmc, - print_every = 100, - matern = TRUE, - sample_iwish=T, - sample_mvr=T, - sample_theta_gibbs=F, - upd_theta_opts=T, - num_threads = 8, - sampling = 2) - }) - - predict_dag <- dag_vecchia_predict(cx_in, cx_all[which_out,], m_nn) - - latentn_predict_time <- system.time({ - spiox_latentn_predicts <- - spiox_latentn_out %>% with( - spiox::spiox_latent_predict(X_new = X[which_out,,drop=F], - coords_new = cx_all[which_out,], - - # data - cx_in, - predict_dag, - W %>% tail(c(NA, NA, round(mcmc/2))), - B %>% tail(c(NA, NA, round(mcmc/2))), - Sigma %>% tail(c(NA, NA, round(mcmc/2))), - Ddiag %>% tail(c(NA, round(mcmc/2))), - theta %>% tail(c(NA, NA, round(mcmc/2))), - matern = TRUE, - num_threads = 16) ) - }) - - save(file=glue::glue("simulations/trivariate_matern_toy_reps/spiox_latentn_{ii}.RData"), - list=c("spiox_latentn_out", "spiox_latentn_predicts", "latentn_estim_time", "latentn_predict_time")) - - } - - ################################################################################ - # lmc via meshed - if(T){ - - # meshed - library(meshed) - - Y_meshed <- Y - Y_meshed[which_out,] <- NA - - keep <- round(mcmc/4) - burn <- round(mcmc/4*3) - - meshed_time <- system.time({ - spmeshed_out <- meshed::spmeshed(y=Y_meshed, x=X, coords=cx_all, k=q, family = "gaussian", - block_size=30, n_samples = keep, n_thin = 1, n_burn = burn, n_threads = 8, verbose = 10, - predict_everywhere = T, - prior = list(phi=c(3, 100), tausq=c(3, .1), nu=c(.5, 1.9))) - }) - - m_order <- order(spmeshed_out$savedata$coords_blocking$ix) - Ymesh_out <- spmeshed_out$yhat_mcmc %>% tailh() %>% abind::abind(along=3) %>% `[`(m_order[which_out],,) - Yhat_meshed <- Ymesh_out %>% apply(1:2, mean) - - save(file=glue::glue("simulations/trivariate_matern_toy_reps/meshed_{ii}.RData"), - list=c("spmeshed_out", "Ymesh_out", "meshed_time")) - - - - - } - - if(T){ - library(GpGpm) - source("~/GpGp_multi_paper/R/helper.R") - source("~/GpGp_multi_paper/R/fisher_multi.R") - source("~/GpGp_multi_paper/R/fit_multi.R") - source("~/GpGp_multi_paper/R/link_multi.R") - source("~/GpGp_multi_paper/R/check.R") - - RhpcBLASctl::blas_set_num_threads(8) - RhpcBLASctl::omp_set_num_threads(8) - - # fit GpGpm - locs <- rbind(cbind(cx_in, 1), cbind(cx_in, 2), cbind(cx_in, 3)) - y <- as.vector(Y_in[,1:3]) - X <- model.matrix(lm( y ~ -1 + as.factor(locs[,ncol(locs)]))) - - # some info - ncomp <- length(unique(locs[,ncol(locs)])) - neach <- ncomp*(ncomp+1)/2 - d <- ncol(locs) - 1 - M <- matrix(0.5, ncomp, ncomp) - diag(M) <- 1 - - # start marginal parms and logparms - start_parms <- get_start_parms(y, X , locs, "matern_multi")$start_parms - start_logparms <- log(start_parms) - start_logparms <- append(start_logparms, 0, 2*neach) - start_logparms <- append(start_logparms, 0, 3*neach+1) - - # some covparms indices - cross_nug_inds <- c() - for(j1 in 2:ncomp){ for(j2 in 1:(j1-1)){ - cross_nug_inds <- c( cross_nug_inds, multi_matern_parm_index(ncomp, j1, j2)$nugget ) - }} - inds <- matrix(FALSE, ncomp, ncomp) - diag(inds) <- TRUE - marginal_ran_inds <- neach + which(t(inds)[upper.tri(inds, diag = TRUE)] == TRUE) - - # some flex-a logparms indices - inds <- matrix(FALSE, ncomp, ncomp) - inds[upper.tri(inds, diag = FALSE)] <- TRUE - inds <- t(inds) - log_cross_var_inds <- which(t(inds)[upper.tri(inds, diag = TRUE)] == TRUE) - log_cross_ran_inds <- neach + log_cross_var_inds - log_delta_B_ind <- 2*neach + 1 - log_cross_smo_inds <- 2*neach + 1 + log_cross_var_inds - log_delta_A_ind <- 3*neach + 2 - log_cross_nug_inds <- 3*neach + 2 + log_cross_var_inds - log_nug_inds <- (3*neach + 3): length(start_logparms) - inds <- matrix(FALSE, ncomp, ncomp) - diag(inds) <- TRUE - log_marginal_ran_inds <- neach + which(t(inds)[upper.tri(inds, diag = TRUE)] == TRUE) - log_marginal_smo_inds <- 2*neach + 1 + which(t(inds)[upper.tri(inds, diag = TRUE)] == TRUE) - - # some flex-b logparms indices - log_beta_ind <- 4*neach + 3 - - # some pars logparms indices - pars_log_cross_nug_inds <- (neach + ncomp + 1) + log_cross_var_inds - - # fit sequence of models - - # Independent - - # using flex_a link, but does not matter which link since fitting independent - linkfuns <- list() - linkfuns$link <- link_flex_a - linkfuns$dlink <- dlink_flex_a - - start_logparms[log_cross_var_inds] <- 0 - start_logparms[log_cross_nug_inds] <- 0 - - active_logparms <- rep(TRUE, length(start_logparms)) - active_logparms[c(log_cross_var_inds, - log_cross_ran_inds, - log_delta_B_ind, - log_cross_smo_inds, - log_delta_A_ind, - log_cross_nug_inds)] <- FALSE - - print("starting independent fit") - t1 <- proc.time() - fit1 <- fit_multi( - y, locs, X, - neighbor_fun = nearest_multi_any, - order_fun = order_completely_random, - m = 40, - start_logparms = start_logparms, - linkfuns = linkfuns, - active_logparms = active_logparms - ) - t2 <- proc.time() - fit1$time_elapsed <- (t2-t1)[3] - fit1$valid <- TRUE - - # flex-e - - # flex-e link - linkfuns$link <- link_flex_e - linkfuns$dlink <- dlink_flex_e - - start_logparms <- fit1$logparms - start_logparms <- append(start_logparms, 0, length(start_logparms)) - start_logparms[log_cross_ran_inds] <- log(finv(M)) - start_logparms[log_delta_B_ind] <- log(0.01) - start_logparms[log_cross_smo_inds] <- log(finv(M)) - start_logparms[log_delta_A_ind] <- log(0.01) - start_logparms[log_beta_ind] <- log(1) - - - active_logparms <- rep(TRUE, length(start_logparms)) - if( FALSE ){ active_logparms[ log_cross_nug_inds ] <- FALSE } - if(ncomp==2){ - active_logparms[log_cross_ran_inds] <- FALSE - active_logparms[log_cross_smo_inds] <- FALSE - } - - print("starting flexible-e fit") - t1 <- proc.time() - fit4 <- fit_multi( - y, locs, X, - neighbor_fun = nearest_multi_any, - order_fun = order_completely_random, - m = 40, - start_logparms = start_logparms, - linkfuns = linkfuns, - active_logparms = active_logparms, - max_iter = 100 - ) - t2 <- proc.time() - fit4$time_elapsed <- (t2-t1)[3] - fit4$valid <- T - - - # 1:6 -- sigma11, sigma12, sigma22, sigma13, sigma23, sigma33, - # 7:12 -- 1/phi, 1/phi, 1/phi, 1/phi, 1/phi, 1/phi, - # 13:18 -- nu1, nu1/2+nu2/2, nu2, nu1/2+nu3/2, nu2/2+nu3/2, nu3, - # 19:24 --1e-3, 1e-3, 1e-3, 1e-3, 1e-3, 1e-3) - - # sigma, 1/phi, nu, tausq - maternify <- function(covparms){ - x <- matrix(0, q,q) - x[upper.tri(x, diag=T)] <- covparms - x[lower.tri(x)] <- x[upper.tri(x)] - return(x) - } - - Sigma_est <- maternify(fit4$covparms[1:6]) - phi_est <- maternify(1/fit4$covparms[7:12]) - nu_est <- maternify(fit4$covparms[13:18]) - tausq_est <- maternify(fit4$covparms[19:24]) - - gpgpm <- list(Sigma=Sigma_est, phi=phi_est, nu=nu_est, tausq=tausq_est) - - save(list=c("fit4", "gpgpm"),file = glue::glue("simulations/trivariate_matern_toy_reps/gpgpm_{ii}.RData")) - - } - -} \ No newline at end of file diff --git a/simulations/trivariate_matern_toy_postprocess.r b/simulations/trivariate_matern_toy_postprocess.r deleted file mode 100644 index 78b32e2..0000000 --- a/simulations/trivariate_matern_toy_postprocess.r +++ /dev/null @@ -1,236 +0,0 @@ -rm(list=ls()) -library(tidyverse) -library(magrittr) -library(Matrix) -library(spiox) - -image.matrix <- \(x) { - image(as(x, "dgCMatrix")) -} - -rtail <- \(x, nt){ - tail(x, nt) -} -tailh <- \(x, nt){ - tail(x, round(length(x)/2)) -} -rtailh <- \(x){ rtail(x, round(nrow(x)/2)) } -ctail <- \(x, nt){ - tail(x, c(NA,nt)) -} -ctailh <- \(x){ ctail(x, round(ncol(x)/2)) } -stail <- \(x, nt){ - tail(x, c(NA,NA,nt)) -} -stailh <- \(x){ stail(x, round(dim(x)[3]/2)) } - -perturb <- function(x, sd=1){ - return(x + matrix(rnorm(prod(dim(x)), sd), ncol=ncol(x))) -} - -set.seed(1) - - -q <- 3 - -nu1 <- 0.5 -nu2 <- .8 -nu3 <- 1.2 - -sds <- c(1, 1, 1) -Omega <- matrix(c(1,-.9,0.7,-.9,1,-.5,0.7,-.5,1), ncol=3) -Sigma <- diag(sds) %*% Omega %*% diag(sds) - -St <- chol(Sigma) -S <- t(St) - -# spatial -cx_in <- matrix(runif(2500*2), ncol=2) #2500 -colnames(cx_in) <- c("Var1","Var2") -n_in <- nrow(cx_in) -which_in <- 1:n_in - -xout <- seq(0, 1, length.out=50) #50 -coords_out <- expand.grid(xout, xout) -cx_out <- as.matrix(coords_out) -n_out <- nrow(cx_out) -which_out <- (n_in+1):(n_in+n_out) - -cx_all <- rbind(cx_in, cx_out) -nr_all <- nrow(cx_all) - - -cx_12 <- rbind(cbind(cx_all, 1), cbind(cx_all, 2), cbind(cx_all, 3)) - -V <- Omega #cbind(c(1, rhocorr), c(rhocorr, 1)) -St <- chol(V) - -phi <- 30 -sigma11 <- sigma22 <- sigma33 <- 1 -matern_scaling_factors <- diag(q) -matern_scaling_factors[1,2] <- matern_scaling_factors[2,1] <- - sqrt(sigma11 * sigma22) * - sqrt(gamma(nu1+1))/sqrt(gamma(nu1)) * - sqrt(gamma(nu2+1))/sqrt(gamma(nu2)) * gamma(nu1/2 + nu2/2)/gamma(nu1/2 + nu2/2 + 1) -matern_scaling_factors[1,3] <- matern_scaling_factors[3,1] <- - sqrt(sigma11 * sigma33) * - sqrt(gamma(nu1+1))/sqrt(gamma(nu1)) * - sqrt(gamma(nu3+1))/sqrt(gamma(nu3)) * gamma(nu1/2 + nu3/2)/gamma(nu1/2 + nu3/2 + 1) -matern_scaling_factors[2,3] <- matern_scaling_factors[3,2] <- - sqrt(sigma22 * sigma33) * - sqrt(gamma(nu2+1))/sqrt(gamma(nu2)) * - sqrt(gamma(nu3+1))/sqrt(gamma(nu3)) * gamma(nu2/2 + nu3/2)/gamma(nu2/2 + nu3/2 + 1) - -sigma12 <- V[2,1] * matern_scaling_factors[1,2] -sigma13 <- V[3,1] * matern_scaling_factors[1,3] -sigma23 <- V[3,2] * matern_scaling_factors[2,3] - - - -set.seed(1) - -#theta_opts <- cbind(c(10, 1, .5, 1e-19), c(20, 1, 1, 1e-19)) -theta_opts <- cbind(c(10, 2, nu1, 1e-1), c(20, 1, nu2, 1e-2), c(20, 1, 0.5, 1e-2)) -############################################## - -m_nn <- 20 -mcmc <- 10000 -RhpcBLASctl::blas_set_num_threads(1) -RhpcBLASctl::omp_set_num_threads(1) - -Sigbuild <- function(spiox_out){ - return( - 1:mcmc %>% sapply(\(i) with(spiox_out, sqrt(diag(theta[2,,i])) %*% Sigma[,,i] %*% sqrt(diag(theta[2,,i])) )) %>% - array(dim=c(q,q,mcmc)) - ) -} - -spiox_dens_plotter <- function(spiox_out, which_var, bw=.6){ - - targets <- spiox_out$theta %>% stailh() %>% `[`(which_var,,) %>% t() %>% as.data.frame() %>% - mutate(m = 1:n()) - colnames(targets) <- paste0("targets_", colnames(targets)) - targets %<>% pivot_longer(cols=-targets_m) - - ggplot(targets, aes(value)) + - geom_boxplot()+ #bw=bw) + - facet_grid( ~ name) + theme_minimal() - - -} - -spiox_cov_at_zero <- function(spiox_out){ - spiox_theta <- spiox_out$theta - - Sigbuild <- function(spiox_out){ - return( - 1:mcmc %>% sapply(\(i) with(spiox_out, sqrt(diag(theta[2,,i])) %*% Sigma[,,i] %*% sqrt(diag(theta[2,,i])) )) %>% - array(dim=c(q,q,mcmc)) - ) - } - Sig <- Sigbuild(spiox_out) - - spiox_theta[2,,] <- 1 - spiox_scaling_factors <- scaling_factor_at_zero(cx_in, spiox_theta %>% tail(c(NA, NA, 5000)) %>% apply(1:2, mean)) - spiox_scaling_factors[upper.tri(spiox_scaling_factors)] <- spiox_scaling_factors[lower.tri(spiox_scaling_factors)] - - return( - Sig %>% apply(3, \(s) cov2cor(s * spiox_scaling_factors)) %>% array(dim=c(q,q,mcmc)) ) -} - -# cov at zero -SS <- V * matern_scaling_factors # true - - -est_results <- list() - -for( rr in 1:60 ){ - cat(rr, "\n") - - # "spiox_out", "spiox_predicts", "estim_time", "predict_time", "total_time" - load(glue::glue("simulations/trivariate_matern_toy_reps/spiox_{rr}.RData")) - - # "spiox_latentq_out", "spiox_latentq_predicts", "latentq_estim_time", "latentq_predict_time" - load(glue::glue("simulations/trivariate_matern_toy_reps/spiox_latentq_{rr}.RData")) - - # "spiox_latentn_out", "spiox_latentn_predicts", "latentn_estim_time", "latentn_predict_time" - load(glue::glue("simulations/trivariate_matern_toy_reps/spiox_latentn_{rr}.RData")) - - # "spmeshed_out", "Ymesh_out", "meshed_time" - load(glue::glue("simulations/trivariate_matern_toy_reps/meshed_{rr}.RData")) - - # "fit2", "gpgpm" - load(glue::glue("simulations/trivariate_matern_toy_reps/gpgpm_{rr}.RData")) - - est_phi <- list( spiox_resp = spiox_out$theta %>% stail(5000) %>% `[`(1,,) %>% apply(1, mean), - spiox_latn = spiox_latentn_out$theta %>% stail(5000) %>% `[`(1,,) %>% apply(1, mean), - spiox_latq = spiox_latentq_out$theta %>% stail(5000) %>% `[`(1,,) %>% apply(1, mean), - meshed = spmeshed_out$theta_mcmc %>% `[`(1,,) %>% apply(1, mean), - gpgpm = diag(gpgpm$phi) ) - - est_nu <- list( spiox_resp = spiox_out$theta %>% stail(5000) %>% `[`(3,,) %>% apply(1, mean), - spiox_latn = spiox_latentn_out$theta %>% stail(5000) %>% `[`(3,,) %>% apply(1, mean), - spiox_latq = spiox_latentq_out$theta %>% stail(5000) %>% `[`(3,,) %>% apply(1, mean), - meshed = spmeshed_out$theta_mcmc %>% `[`(2,,) %>% apply(1, mean), - gpgpm = diag(gpgpm$nu) ) - - - est_tsq <- list( spiox_resp = spiox_out$theta %>% stail(5000) %>% `[`(4,,) %>% apply(1, mean), - spiox_latn = spiox_latentn_out$Ddiag %>% ctail(5000) %>% apply(1, mean), - spiox_latq = spiox_latentq_out$Ddiag %>% ctail(5000) %>% apply(1, mean), - meshed = spmeshed_out$tausq_mcmc %>% apply(1, mean), - gpgpm = diag(gpgpm$tausq) ) - - spiox_cor_at_zero <- function(spiox_out){ - spiox_theta <- spiox_out$theta - - Sigbuild <- function(spiox_out, mcmc=10000){ - return( - 1:mcmc %>% sapply(\(i) with(spiox_out, sqrt(diag(theta[2,,i])) %*% Sigma[,,i] %*% sqrt(diag(theta[2,,i])) )) %>% - array(dim=c(q,q,mcmc)) - ) - } - Sig <- Sigbuild(spiox_out) - - spiox_theta[2,,] <- 1 - spiox_scaling_factors <- scaling_factor_at_zero(cx_in, spiox_theta %>% tail(c(NA, NA, 5000)) %>% apply(1:2, mean)) - spiox_scaling_factors[upper.tri(spiox_scaling_factors)] <- spiox_scaling_factors[lower.tri(spiox_scaling_factors)] - - return( - Sig %>% apply(3, \(s) cov2cor(s * spiox_scaling_factors)) %>% array(dim=c(q,q,mcmc)) ) - } - - # covariance at zero: estimated - spiox_SS_resp <- spiox_out %>% spiox_cor_at_zero() %>% stail(5000) %>% apply(1:2, mean) - spiox_SS_latq <- spiox_latentq_out %>% spiox_cor_at_zero() %>% stail(5000) %>% apply(1:2, mean) - spiox_SS_latn <- spiox_latentn_out %>% spiox_cor_at_zero() %>% stail(5000) %>% apply(1:2, mean) - - meshed_SS <- spmeshed_out$lambda_mcmc %>% meshed:::cube_correl_from_lambda() %>% stail(5000) %>% apply(1:2, mean) - gpgpm_SS <- cov2cor(gpgpm$Sigma) - - - est_Corr0 <- list( spiox_resp = spiox_SS_resp[lower.tri(spiox_SS_resp)], - spiox_latn = spiox_SS_latn[lower.tri(spiox_SS_latn)], - spiox_latq = spiox_SS_latq[lower.tri(spiox_SS_latq)], - meshed = meshed_SS[lower.tri(meshed_SS)], - gpgpm = gpgpm_SS[lower.tri(gpgpm_SS)] ) %>% as.data.frame() %>% - mutate(parameter = c("corr21", "corr31", "corr32")) - - - est_results[[rr]] <- bind_rows( - as.data.frame(est_phi) %>% mutate(parameter = paste0("phi", 1:n())), - as.data.frame(est_nu) %>% mutate(parameter = paste0("nu", 1:n())), - as.data.frame(est_tsq) %>% mutate(parameter = paste0("tausq", 1:n())), - est_Corr0 - ) %>% mutate(rep = rr) %>% pivot_longer(cols=-c(parameter, rep), names_to = "model") - - -} - - -save(file="simulations/trivariate_matern_toy_reps/estimation_results_all.RData", list="est_results") - -est_results %>% bind_rows() %>% group_by(model, parameter) %>% summarise_all(mean) %>% pivot_wider(id_cols=model, names_from = parameter) - - - diff --git a/simulations/trivariate_spiox_toy.r b/simulations/trivariate_spiox_toy.r deleted file mode 100644 index 6aa6492..0000000 --- a/simulations/trivariate_spiox_toy.r +++ /dev/null @@ -1,482 +0,0 @@ - -library(tidyverse) -library(magrittr) -library(Matrix) -library(spiox) - -image.matrix <- \(x) { - image(as(x, "dgCMatrix")) -} - -rtail <- \(x, nt){ - tail(x, nt) -} -tailh <- \(x, nt){ - tail(x, round(length(x)/2)) -} -rtailh <- \(x){ rtail(x, round(nrow(x)/2)) } -ctail <- \(x, nt){ - tail(x, c(NA,nt)) -} -ctailh <- \(x){ ctail(x, round(ncol(x)/2)) } -stail <- \(x, nt){ - tail(x, c(NA,NA,nt)) -} -stailh <- \(x){ stail(x, round(dim(x)[3]/2)) } - -perturb <- function(x, sd=1){ - return(x + matrix(rnorm(prod(dim(x)), sd), ncol=ncol(x))) -} - -g <- glue::glue - -for(ii in 1:3){ - #set.seed(2024) - set.seed(ii) - - q <- 3 - - theta_mat <- matrix(1, ncol=q, nrow=4) - theta_mat[1,] <- c(phi <- 20, phi, phi) - theta_mat[2,] <- c(1, 1, 1) - theta_mat[3,] <- c(nu1 <- 0.5, nu2 <- .8, nu3 <- 1.2) - theta_mat[4,] <- rep(1e-13, q) - - sds <- c(1, 1, 1) - Omega <- matrix(c(1,-.9,0.7,-.9,1,-.5,0.7,-.5,1), ncol=3) - Sigma <- diag(sds) %*% Omega %*% diag(sds) - - St <- chol(Sigma) - S <- t(St) - - # spatial - cx_in <- matrix(runif(2500*2), ncol=2) #2500 - colnames(cx_in) <- c("Var1","Var2") - n_in <- nrow(cx_in) - which_in <- 1:n_in - - xout <- seq(0, 1, length.out=50) #50 - coords_out <- expand.grid(xout, xout) - cx_out <- as.matrix(coords_out) - n_out <- nrow(cx_out) - which_out <- (n_in+1):(n_in+n_out) - - cx_all <- rbind(cx_in, cx_out) - nr_all <- nrow(cx_all) - - - cx_12 <- rbind(cbind(cx_all, 1), cbind(cx_all, 2), cbind(cx_all, 3)) - - V <- Omega #cbind(c(1, rhocorr), c(rhocorr, 1)) - St <- chol(V) - - phi <- 30 - - if(T){ - D <- rep(1e-3, q) - nulist <- c(nu1, nu2, nu3) - - Clist <- 1:q %>% lapply(\(j) spiox::Correlationc(cx_all, cx_all, c(phi, 1, nulist[j], D[j]), TRUE, TRUE) ) - Llist <- Clist %>% lapply(\(C) t(chol(C))) - Lilist <- Llist %>% lapply(\(L) solve(L)) - - V <- matrix(rnorm(nr_all * q), ncol=q) %*% St - - Y_sp <- V - # make q spatial outcomes - for(i in 1:q){ - Y_sp[,i] <- Llist[[i]] %*% V[,i] - } - - # regression - p <- 1 - X <- matrix(1, ncol=1, nrow=nr_all) #%>% cbind(matrix(rnorm(nr_all*(p-1)), ncol=p-1)) - - Beta <- matrix(rnorm(q * p), ncol=q) - - Y_regression <- X %*% Beta - #Error <- matrix(rnorm(nrow(Y_regression) * q),ncol=q) %*% diag(D <- runif(q, 0, 0.1)) - Y <- as.matrix(Y_sp + Y_regression) # + Error - - Y_in <- Y[which_in,] - X_in <- X[which_in,,drop=F] - - if(F){ - df <- data.frame(cx_out, y=as.matrix(Y_sp[which_out,])) %>% - pivot_longer(cols=-c(Var1, Var2)) - ggplot(df, - aes(Var1, Var2, fill=value)) + - geom_raster() + - scale_fill_viridis_c() + - facet_wrap(~name, ncol=5) - } - - simdata <- data.frame(coords=cx_all, Y_spatial=Y_sp, Y=Y, X=X) - #ggplot(simdata, aes(coords.Var1, coords.Var2, color=Y_spatial.1)) + geom_point() + scale_color_viridis_c() - - save(file=glue::glue("simulations/trivariate_spiox_toy_reps/data_{ii}.RData"), - list=c("simdata", - "Beta", "D", "Y_in", "X_in", "which_in", "which_out", - "Y_regression", #"Error", - "Y", "X", - "nulist", "phi", "Sigma")) - } else { - load(g("simulations/trivariate_spiox_toy_reps/data_{ii}.RData")) - } - ############################## - - set.seed(1) - - #theta_opts <- cbind(c(10, 1, .5, 1e-19), c(20, 1, 1, 1e-19)) - theta_opts <- cbind(c(phi+1e-8, 1.001, nu1, 1e-3), c(phi-1e-8, 0.999, nu2, 0.0011), c(phi, 1, nu3, 1e-3)) - ############################################## - - m_nn <- 20 - mcmc <- 10000 - RhpcBLASctl::blas_set_num_threads(1) - RhpcBLASctl::omp_set_num_threads(1) - - if(T){ - custom_dag <- dag_vecchia(cx_in, m_nn) - - ############################################## - set.seed(1) - estim_time <- system.time({ - spiox_out <- spiox::spiox_wishart(Y_in, X_in, cx_in, - custom_dag = custom_dag, - theta=theta_opts, - - Sigma_start = Sigma, - mvreg_B_start = Beta,# %>% perturb(), - - mcmc = mcmc, - print_every = 100, - matern = TRUE, - sample_iwish=T, - sample_mvr=T, - sample_theta_gibbs=F, - upd_theta_opts=T, - num_threads = 8) - }) - - predict_dag <- dag_vecchia_predict(cx_in, cx_all[which_out,], m_nn) - - predict_time <- system.time({ - spiox_predicts <- spiox::spiox_predict(X_new = X[which_out,,drop=F], - coords_new = cx_all[which_out,], - - # data - Y_in, X_in, cx_in, - predict_dag, - spiox_out$B %>% tail(c(NA, NA, round(mcmc/2))), - spiox_out$Sigma %>% tail(c(NA, NA, round(mcmc/2))), - spiox_out$theta %>% tail(c(NA, NA, round(mcmc/2))), - matern = TRUE, - num_threads = 8) - }) - - - Y_out <- Y[which_out,] - Yhat_spiox <- spiox_predicts$Y %>% apply(1:2, mean) - - - - - total_time <- estim_time + predict_time - - save(file=glue::glue("simulations/trivariate_spiox_toy_reps/spiox_{ii}.RData"), - list=c("spiox_out", "spiox_predicts", "estim_time", "predict_time", "total_time")) - } - - if(T){ - custom_dag <- dag_vecchia(cx_in, m_nn) - theta_opts_latent <- theta_opts - theta_opts_latent[4,] <- 0 - ############################################## - set.seed(1) - latentq_estim_time <- system.time({ - spiox_latentq_out <- spiox::spiox_latent(Y_in, X_in, cx_in, - custom_dag = custom_dag, - theta=theta_opts_latent, - - Sigma_start = Sigma, - mvreg_B_start = Beta,# %>% perturb(), - - mcmc = mcmc, - print_every = 100, - matern = TRUE, - sample_iwish=T, - sample_mvr=T, - sample_theta_gibbs=F, - upd_theta_opts=T, - num_threads = 8, - sampling = 3) - }) - - predict_dag <- dag_vecchia_predict(cx_in, cx_all[which_out,], m_nn) - - latentq_predict_time <- system.time({ - spiox_latentq_predicts <- - spiox_latentq_out %>% with( - spiox::spiox_latent_predict(X_new = X[which_out,,drop=F], - coords_new = cx_all[which_out,], - - # data - cx_in, - predict_dag, - W %>% tail(c(NA, NA, round(mcmc/2))), - B %>% tail(c(NA, NA, round(mcmc/2))), - Sigma %>% tail(c(NA, NA, round(mcmc/2))), - Ddiag %>% tail(c(NA, round(mcmc/2))), - theta %>% tail(c(NA, NA, round(mcmc/2))), - matern = TRUE, - num_threads = 8) ) - }) - - save(file=glue::glue("simulations/trivariate_spiox_toy_reps/spiox_latentq_{ii}.RData"), - list=c("spiox_latentq_out", "spiox_latentq_predicts", "latentq_estim_time", "latentq_predict_time")) - - } - - if(T){ - - custom_dag <- dag_vecchia(cx_in, m_nn) - theta_opts_latent <- theta_opts - theta_opts_latent[2,1] <- 2 - theta_opts_latent[4,] <- 0 - ############################################## - set.seed(1) - latentn_estim_time <- system.time({ - spiox_latentn_out <- spiox::spiox_latent(Y_in, X_in, cx_in, - custom_dag = custom_dag, - theta=theta_opts_latent, - - Sigma_start = diag(q), - mvreg_B_start = 0*Beta,# %>% perturb(), - - mcmc = mcmc, - print_every = 100, - matern = TRUE, - sample_iwish=T, - sample_mvr=T, - sample_theta_gibbs=F, - upd_theta_opts=T, - num_threads = 8, - sampling = 2) - }) - - predict_dag <- dag_vecchia_predict(cx_in, cx_all[which_out,], m_nn) - - latentn_predict_time <- system.time({ - spiox_latentn_predicts <- - spiox_latentn_out %>% with( - spiox::spiox_latent_predict(X_new = X[which_out,,drop=F], - coords_new = cx_all[which_out,], - - # data - cx_in, - predict_dag, - W %>% tail(c(NA, NA, round(mcmc/2))), - B %>% tail(c(NA, NA, round(mcmc/2))), - Sigma %>% tail(c(NA, NA, round(mcmc/2))), - Ddiag %>% tail(c(NA, round(mcmc/2))), - theta %>% tail(c(NA, NA, round(mcmc/2))), - matern = TRUE, - num_threads = 8) ) - }) - - save(file=glue::glue("simulations/trivariate_spiox_toy_reps/spiox_latentn_{ii}.RData"), - list=c("spiox_latentn_out", "spiox_latentn_predicts", "latentn_estim_time", "latentn_predict_time")) - - } - - ################################################################################ - # lmc via meshed - if(T){ - - # meshed - library(meshed) - - Y_meshed <- Y - Y_meshed[which_out,] <- NA - - keep <- round(mcmc/4) - burn <- round(mcmc/4*3) - - meshed_time <- system.time({ - spmeshed_out <- meshed::spmeshed(y=Y_meshed, x=X, coords=cx_all, k=q, family = "gaussian", - block_size=30, n_samples = keep, n_thin = 1, n_burn = burn, - n_threads = 8, verbose = 10, - predict_everywhere = T, - prior = list(phi=c(3, 100), tausq=c(3, .1), nu=c(.5, 1.9))) - }) - - m_order <- order(spmeshed_out$savedata$coords_blocking$ix) - Ymesh_out <- spmeshed_out$yhat_mcmc %>% tailh() %>% abind::abind(along=3) %>% `[`(m_order[which_out],,) - Yhat_meshed <- Ymesh_out %>% apply(1:2, mean) - - save(file=glue::glue("simulations/trivariate_spiox_toy_reps/meshed_{ii}.RData"), - list=c("spmeshed_out", "Ymesh_out", "meshed_time")) - - - - - } - - if(T){ - library(GpGpm) - source("~/GpGp_multi_paper/R/helper.R") - source("~/GpGp_multi_paper/R/fisher_multi.R") - source("~/GpGp_multi_paper/R/fit_multi.R") - source("~/GpGp_multi_paper/R/link_multi.R") - source("~/GpGp_multi_paper/R/check.R") - - RhpcBLASctl::blas_set_num_threads(8) - RhpcBLASctl::omp_set_num_threads(8) - - # fit GpGpm - locs <- rbind(cbind(cx_in, 1), cbind(cx_in, 2), cbind(cx_in, 3)) - y <- as.vector(Y_in[,1:3]) - X <- model.matrix(lm( y ~ -1 + as.factor(locs[,ncol(locs)]))) - - # some info - ncomp <- length(unique(locs[,ncol(locs)])) - neach <- ncomp*(ncomp+1)/2 - d <- ncol(locs) - 1 - M <- matrix(0.5, ncomp, ncomp) - diag(M) <- 1 - - # start marginal parms and logparms - start_parms <- get_start_parms(y, X , locs, "matern_multi")$start_parms - start_logparms <- log(start_parms) - start_logparms <- append(start_logparms, 0, 2*neach) - start_logparms <- append(start_logparms, 0, 3*neach+1) - - # some covparms indices - cross_nug_inds <- c() - for(j1 in 2:ncomp){ for(j2 in 1:(j1-1)){ - cross_nug_inds <- c( cross_nug_inds, multi_matern_parm_index(ncomp, j1, j2)$nugget ) - }} - inds <- matrix(FALSE, ncomp, ncomp) - diag(inds) <- TRUE - marginal_ran_inds <- neach + which(t(inds)[upper.tri(inds, diag = TRUE)] == TRUE) - - # some flex-a logparms indices - inds <- matrix(FALSE, ncomp, ncomp) - inds[upper.tri(inds, diag = FALSE)] <- TRUE - inds <- t(inds) - log_cross_var_inds <- which(t(inds)[upper.tri(inds, diag = TRUE)] == TRUE) - log_cross_ran_inds <- neach + log_cross_var_inds - log_delta_B_ind <- 2*neach + 1 - log_cross_smo_inds <- 2*neach + 1 + log_cross_var_inds - log_delta_A_ind <- 3*neach + 2 - log_cross_nug_inds <- 3*neach + 2 + log_cross_var_inds - log_nug_inds <- (3*neach + 3): length(start_logparms) - inds <- matrix(FALSE, ncomp, ncomp) - diag(inds) <- TRUE - log_marginal_ran_inds <- neach + which(t(inds)[upper.tri(inds, diag = TRUE)] == TRUE) - log_marginal_smo_inds <- 2*neach + 1 + which(t(inds)[upper.tri(inds, diag = TRUE)] == TRUE) - - # some flex-b logparms indices - log_beta_ind <- 4*neach + 3 - - # some pars logparms indices - pars_log_cross_nug_inds <- (neach + ncomp + 1) + log_cross_var_inds - - # fit sequence of models - - # Independent - - # using flex_a link, but does not matter which link since fitting independent - linkfuns <- list() - linkfuns$link <- link_flex_a - linkfuns$dlink <- dlink_flex_a - - start_logparms[log_cross_var_inds] <- 0 - start_logparms[log_cross_nug_inds] <- 0 - - active_logparms <- rep(TRUE, length(start_logparms)) - active_logparms[c(log_cross_var_inds, - log_cross_ran_inds, - log_delta_B_ind, - log_cross_smo_inds, - log_delta_A_ind, - log_cross_nug_inds)] <- FALSE - - print("starting independent fit") - t1 <- proc.time() - fit1 <- fit_multi( - y, locs, X, - neighbor_fun = nearest_multi_any, - order_fun = order_completely_random, - m = 40, - start_logparms = start_logparms, - linkfuns = linkfuns, - active_logparms = active_logparms - ) - t2 <- proc.time() - fit1$time_elapsed <- (t2-t1)[3] - fit1$valid <- TRUE - - # flex-e - - # flex-e link - linkfuns$link <- link_flex_e - linkfuns$dlink <- dlink_flex_e - - start_logparms <- fit1$logparms - start_logparms <- append(start_logparms, 0, length(start_logparms)) - start_logparms[log_cross_ran_inds] <- log(finv(M)) - start_logparms[log_delta_B_ind] <- log(0.01) - start_logparms[log_cross_smo_inds] <- log(finv(M)) - start_logparms[log_delta_A_ind] <- log(0.01) - start_logparms[log_beta_ind] <- log(1) - - - active_logparms <- rep(TRUE, length(start_logparms)) - if( FALSE ){ active_logparms[ log_cross_nug_inds ] <- FALSE } - if(ncomp==2){ - active_logparms[log_cross_ran_inds] <- FALSE - active_logparms[log_cross_smo_inds] <- FALSE - } - - print("starting flexible-e fit") - t1 <- proc.time() - fit4 <- fit_multi( - y, locs, X, - neighbor_fun = nearest_multi_any, - order_fun = order_completely_random, - m = 40, - start_logparms = start_logparms, - linkfuns = linkfuns, - active_logparms = active_logparms - ) - t2 <- proc.time() - fit4$time_elapsed <- (t2-t1)[3] - fit4$valid <- T - - # 1:6 -- sigma11, sigma12, sigma22, sigma13, sigma23, sigma33, - # 7:12 -- 1/phi, 1/phi, 1/phi, 1/phi, 1/phi, 1/phi, - # 13:18 -- nu1, nu1/2+nu2/2, nu2, nu1/2+nu3/2, nu2/2+nu3/2, nu3, - # 19:24 --1e-3, 1e-3, 1e-3, 1e-3, 1e-3, 1e-3) - - # sigma, 1/phi, nu, tausq - maternify <- function(covparms){ - x <- matrix(0, q,q) - x[upper.tri(x, diag=T)] <- covparms - x[lower.tri(x)] <- x[upper.tri(x)] - return(x) - } - - Sigma_est <- maternify(fit4$covparms[1:6]) - phi_est <- maternify(1/fit4$covparms[7:12]) - nu_est <- maternify(fit4$covparms[13:18]) - tausq_est <- maternify(fit4$covparms[19:24]) - - gpgpm <- list(Sigma=Sigma_est, phi=phi_est, nu=nu_est, tausq=tausq_est) - - save(list=c("fit4", "gpgpm"),file = glue::glue("simulations/trivariate_spiox_toy_reps/gpgpm_{ii}.RData")) - - } - -} \ No newline at end of file diff --git a/simulations/trivariate_spiox_toy_postprocess.r b/simulations/trivariate_spiox_toy_postprocess.r deleted file mode 100644 index 181bd6c..0000000 --- a/simulations/trivariate_spiox_toy_postprocess.r +++ /dev/null @@ -1,206 +0,0 @@ -rm(list=ls()) -library(tidyverse) -library(magrittr) -library(Matrix) -library(spiox) - -image.matrix <- \(x) { - image(as(x, "dgCMatrix")) -} - -rtail <- \(x, nt){ - tail(x, nt) -} -tailh <- \(x, nt){ - tail(x, round(length(x)/2)) -} -rtailh <- \(x){ rtail(x, round(nrow(x)/2)) } -ctail <- \(x, nt){ - tail(x, c(NA,nt)) -} -ctailh <- \(x){ ctail(x, round(ncol(x)/2)) } -stail <- \(x, nt){ - tail(x, c(NA,NA,nt)) -} -stailh <- \(x){ stail(x, round(dim(x)[3]/2)) } - -perturb <- function(x, sd=1){ - return(x + matrix(rnorm(prod(dim(x)), sd), ncol=ncol(x))) -} - -set.seed(1) - - -q <- 3 - -theta_mat <- matrix(1, ncol=q, nrow=4) -theta_mat[1,] <- c(phi <- 30, phi, phi) -theta_mat[2,] <- c(1, 1, 1) -theta_mat[3,] <- c(nu1 <- 0.5, nu2 <- .8, nu3 <- 1.2) -theta_mat[4,] <- rep(0.001, q) -sds <- c(1, 1, 1) -Omega <- matrix(c(1,-.9,0.7,-.9,1,-.5,0.7,-.5,1), ncol=3) -Sigma <- diag(sds) %*% Omega %*% diag(sds) - -St <- chol(Sigma) -S <- t(St) - -# spatial -cx_in <- matrix(runif(2500*2), ncol=2) #2500 -colnames(cx_in) <- c("Var1","Var2") -n_in <- nrow(cx_in) -which_in <- 1:n_in - -xout <- seq(0, 1, length.out=50) #50 -coords_out <- expand.grid(xout, xout) -cx_out <- as.matrix(coords_out) -n_out <- nrow(cx_out) -which_out <- (n_in+1):(n_in+n_out) - -cx_all <- rbind(cx_in, cx_out) -nr_all <- nrow(cx_all) - - -cx_12 <- rbind(cbind(cx_all, 1), cbind(cx_all, 2), cbind(cx_all, 3)) - -V <- Omega #cbind(c(1, rhocorr), c(rhocorr, 1)) -St <- chol(V) - -phi <- 30 - -load("simulations/trivariate_spiox_toy/data.RData") - -############################## - -simdata %>% tail(2500) %>% - dplyr::select(contains("coords"), contains("Y_")) %>% - pivot_longer(cols=-c(coords.Var1, coords.Var2)) %>% - ggplot(aes(coords.Var1, coords.Var2, fill=value)) + - geom_raster() + - scale_fill_viridis_c() + - facet_grid(~name) + - theme_minimal() - - - -set.seed(1) - -#theta_opts <- cbind(c(10, 1, .5, 1e-19), c(20, 1, 1, 1e-19)) -theta_opts <- cbind(c(10, 2, nu1, 1e-1), c(20, 1, nu2, 1e-2), c(20, 1, 0.5, 1e-2)) -############################################## - -m_nn <- 20 -mcmc <- 10000 -RhpcBLASctl::blas_set_num_threads(1) -RhpcBLASctl::omp_set_num_threads(1) - - -spiox_dens_plotter <- function(spiox_out, which_var, bw=.6){ - - targets <- spiox_out$theta %>% stailh() %>% `[`(which_var,,) %>% t() %>% as.data.frame() %>% - mutate(m = 1:n()) - colnames(targets) <- paste0("targets_", colnames(targets)) - targets %<>% pivot_longer(cols=-targets_m) - - ggplot(targets, aes(value)) + - geom_boxplot()+ #bw=bw) + - facet_grid( ~ name) + theme_minimal() - - -} - - -# covariance at zero: true -spiox_scaling_factors <- scaling_factor_at_zero(cx_in, theta_mat) -spiox_scaling_factors[upper.tri(spiox_scaling_factors)] <- spiox_scaling_factors[lower.tri(spiox_scaling_factors)] -SS <- Sigma * spiox_scaling_factors - -est_results <- list() - -for( rr in 1:60 ){ - cat(rr, "\n") - - # "spiox_out", "spiox_predicts", "estim_time", "predict_time", "total_time" - load(glue::glue("simulations/trivariate_spiox_toy_reps/spiox_{rr}.RData")) - - # "spiox_latentq_out", "spiox_latentq_predicts", "latentq_estim_time", "latentq_predict_time" - load(glue::glue("simulations/trivariate_spiox_toy_reps/spiox_latentq_{rr}.RData")) - - # "spiox_latentn_out", "spiox_latentn_predicts", "latentn_estim_time", "latentn_predict_time" - load(glue::glue("simulations/trivariate_spiox_toy_reps/spiox_latentn_{rr}.RData")) - - # "spmeshed_out", "Ymesh_out", "meshed_time" - load(glue::glue("simulations/trivariate_spiox_toy_reps/meshed_{rr}.RData")) - - # "fit2", "gpgpm" - load(glue::glue("simulations/trivariate_spiox_toy_reps/gpgpm_{rr}.RData")) - - est_phi <- list( spiox_resp = spiox_out$theta %>% stail(5000) %>% `[`(1,,) %>% apply(1, mean), - spiox_latn = spiox_latentn_out$theta %>% stail(5000) %>% `[`(1,,) %>% apply(1, mean), - spiox_latq = spiox_latentq_out$theta %>% stail(5000) %>% `[`(1,,) %>% apply(1, mean), - meshed = spmeshed_out$theta_mcmc %>% `[`(1,,) %>% apply(1, mean), - gpgpm = diag(gpgpm$phi) ) - - est_nu <- list( spiox_resp = spiox_out$theta %>% stail(5000) %>% `[`(3,,) %>% apply(1, mean), - spiox_latn = spiox_latentn_out$theta %>% stail(5000) %>% `[`(3,,) %>% apply(1, mean), - spiox_latq = spiox_latentq_out$theta %>% stail(5000) %>% `[`(3,,) %>% apply(1, mean), - meshed = spmeshed_out$theta_mcmc %>% `[`(2,,) %>% apply(1, mean), - gpgpm = diag(gpgpm$nu) ) - - - est_tsq <- list( spiox_resp = spiox_out$theta %>% stail(5000) %>% `[`(4,,) %>% apply(1, mean), - spiox_latn = spiox_latentn_out$Ddiag %>% ctail(5000) %>% apply(1, mean), - spiox_latq = spiox_latentq_out$Ddiag %>% ctail(5000) %>% apply(1, mean), - meshed = spmeshed_out$tausq_mcmc %>% apply(1, mean), - gpgpm = diag(gpgpm$tausq) ) - - spiox_cor_at_zero <- function(spiox_out){ - spiox_theta <- spiox_out$theta - - Sigbuild <- function(spiox_out, mcmc=10000){ - return( - 1:mcmc %>% sapply(\(i) with(spiox_out, sqrt(diag(theta[2,,i])) %*% Sigma[,,i] %*% sqrt(diag(theta[2,,i])) )) %>% - array(dim=c(q,q,mcmc)) - ) - } - Sig <- Sigbuild(spiox_out) - - spiox_theta[2,,] <- 1 - spiox_scaling_factors <- scaling_factor_at_zero(cx_in, spiox_theta %>% tail(c(NA, NA, 5000)) %>% apply(1:2, mean)) - spiox_scaling_factors[upper.tri(spiox_scaling_factors)] <- spiox_scaling_factors[lower.tri(spiox_scaling_factors)] - - return( - Sig %>% apply(3, \(s) cov2cor(s * spiox_scaling_factors)) %>% array(dim=c(q,q,mcmc)) ) - } - - # covariance at zero: estimated - spiox_SS_resp <- spiox_out %>% spiox_cor_at_zero() %>% stail(5000) %>% apply(1:2, mean) - spiox_SS_latq <- spiox_latentq_out %>% spiox_cor_at_zero() %>% stail(5000) %>% apply(1:2, mean) - spiox_SS_latn <- spiox_latentn_out %>% spiox_cor_at_zero() %>% stail(5000) %>% apply(1:2, mean) - - meshed_SS <- spmeshed_out$lambda_mcmc %>% meshed:::cube_correl_from_lambda() %>% stail(5000) %>% apply(1:2, mean) - gpgpm_SS <- cov2cor(gpgpm$Sigma) - - - est_Corr0 <- list( spiox_resp = spiox_SS_resp[lower.tri(spiox_SS_resp)], - spiox_latn = spiox_SS_latn[lower.tri(spiox_SS_latn)], - spiox_latq = spiox_SS_latq[lower.tri(spiox_SS_latq)], - meshed = meshed_SS[lower.tri(meshed_SS)], - gpgpm = gpgpm_SS[lower.tri(gpgpm_SS)] ) %>% as.data.frame() %>% - mutate(parameter = c("corr21", "corr31", "corr32")) - - - est_results[[rr]] <- bind_rows( - as.data.frame(est_phi) %>% mutate(parameter = paste0("phi", 1:n())), - as.data.frame(est_nu) %>% mutate(parameter = paste0("nu", 1:n())), - as.data.frame(est_tsq) %>% mutate(parameter = paste0("tausq", 1:n())), - est_Corr0 - ) %>% mutate(rep = rr) %>% pivot_longer(cols=-c(parameter, rep), names_to = "model") - - -} - -save(file="simulations/trivariate_spiox_toy_reps/estimation_results_all.RData", list="est_results") - - - diff --git a/spiox_for_cij.r b/spiox_for_cij.r deleted file mode 100644 index d694900..0000000 --- a/spiox_for_cij.r +++ /dev/null @@ -1,182 +0,0 @@ -rm(list=ls()) -library(tidyverse) -library(magrittr) -library(Matrix) - -set.seed(2) - -image.matrix <- \(x) { - image(as(x, "dgCMatrix")) -} - -# spatial -xx <- seq(0, 1, length.out=30) -coords <- expand.grid(xx, xx) -cx <- as.matrix(coords) -nr <- nrow(cx) - -q <- 3 - -k <- 3 -optlist <- seq(3, 20, length.out=q) %>% sample(k, replace=T) - -Clist <- optlist %>% lapply(\(phi) (exp(-phi * as.matrix(dist(cx))^1) + 1e-16*diag(nr)) ) -#Clist <- optlist %>% lapply(\(nu) spiox::Correlationc(cx, cx, c(20,1,nu,1e-8), 1, TRUE) ) -Llist <- Clist %>% lapply(\(C) t(chol(C))) - -Q <- rWishart(1, q+2, diag(q))[,,1] # -Sigma <- solve(Q) - -St <- chol(Sigma) -S <- t(St) - -W <- 1:k %>% lapply(\(j) Llist[[j]] %*% rnorm(nr)) %>% abind::abind(along=2) - -Y_sp <- W %*% St - -# regression -p <- 2 -X <- matrix(1, ncol=1, nrow=nr) %>% cbind(matrix(rnorm(nr*(p-1)), ncol=p-1)) - -Beta <- matrix(rnorm(q * p), ncol=q) - -Y_regression <- X %*% Beta -Error <- matrix(rnorm(nr * q),ncol=q) # %*% diag(D <- runif(q, 0, 0.1)) -Y <- as.matrix(Y_sp + Y_regression) #+ Error - -cov_Y_cols <- cov(Y) %>% as("dgCMatrix") -cov_Y_rows <- cov(t(Y)) %>% as("dgCMatrix") - -df <- data.frame(coords, y=Y_sp) %>% - pivot_longer(cols=-c(Var1, Var2)) - -if(F){ - ggplot(df, - aes(Var1, Var2, fill=value)) + - geom_raster() + - scale_fill_viridis_c() + - facet_wrap(~name, ncol=5) -} - -############################## - -set.seed(1) - -radgp_rho <- .15 -test_radgp <- spiox::radgp_build(cx, radgp_rho, phi=10, sigmasq=1, nu=1.5, tausq=0, matern=F) -test_radgp$dag %>% sapply(\(x) nrow(x)) %>% summary() - -par_opts <- seq(5, 30, length.out=k) -theta_opts <- par_opts %>% sapply(\(phi) matrix( c(phi, 1, 1, 1e-16), ncol=1)) - -testset <- sample(1:nrow(Y), 10, replace=F) - -perturb <- function(x, sd=1){ - return(x + matrix(rnorm(prod(dim(x)), sd), ncol=ncol(x))) -} - -set.seed(1) -(total_time <- system.time({ - spiox_out <- spiox::spiox_wishart(Y[-testset,,drop=F], - X[-testset,,drop=F], - cx[-testset,], - radgp_rho = radgp_rho, theta=theta_opts, - - Sigma_start = diag(q), - mvreg_B_start = 0*Beta,# %>% perturb(), - - mcmc = mcmc <- 5000, - print_every = 50, - - sample_iwish=T, - sample_mvr=T, - sample_gp=T, - upd_opts=T, - num_threads = 16) -})) - -Sig_est <- 1:mcmc %>% sapply(\(m) with(spiox_out, diag(sqrt(theta[2,,m])) %*% - Sigma[,,m] %*% diag(sqrt(theta[2,,m])) ) ) %>% - array(c(q,q,mcmc)) %>% apply(1:2, mean) - -thetamean <- spiox_out$theta %>% apply(1:2, mean) -philist <- thetamean[1,] -siglist <- thetamean[2,] - -hvec <- c(seq(0, 0.008, length.out=10), seq(0.01, .5, length.out=10)) -xg <- seq(0.1, 0.9, length.out=5) -test_coords <- expand.grid(xg,xg) %>% as.matrix() - - -i <- 2; j <- 2 -testh <- Sig_est[i,j] * spiox::iox_cross_avg(hvec, i, j, - test_coords, cx, - philist, c(1,1,1), 12, 1, num_threads=15) - -truth <- hvec %>% sapply(\(h) (S %*% diag(exp(-optlist*h)) %*% t(S))[i,j]) - -mesh_total_time <- system.time({ - meshout <- meshed::spmeshed(y=Y, - x=X, - coords=cx, k = q, - block_size = 50, - n_samples = 10000, - n_burn = 1, - n_thin = 1, - n_threads = 16, - prior = list(phi=c(1, 20)), - verbose=50 - )}) -h <- 0 -mesh_mcmc <- dim(meshout$theta_mcmc)[3] - -mesh_h_f <- function(h){ - sig_mcmc <- 1:mesh_mcmc %>% sapply(\(m) with(meshout, - lambda_mcmc[,,m] %*% diag(exp(-theta_mcmc[1,,m]*h)) %*% t(lambda_mcmc[,,m]) + diag(tausq_mcmc[,m]) - )) %>% array(dim=c(3,3,mesh_mcmc)) - return(sig_mcmc %>% apply(1:2, mean)) -} - - - -phi_mesh_out <- meshout$theta_mcmc[1,,] %>% apply(1, mean) - -Sig_mesh_out <- meshout$lambda_mcmc %>% apply(3,\(l) tcrossprod(l)) %>% array(dim=c(q,q,mcmc)) -Sig_mesh_est <- Sig_mesh_out %>% apply(1:2, mean) -S_mesh <- t(chol(Sig_mesh_est)) - - - - -meshh <- hvec %>% sapply(\(h) mesh_h_f(h)[i,j]) - -plot(hvec, truth, type='l', ylim=c(-2,1)) -lines(hvec, meshh, col="red") -lines(hvec, testh, col="blue") - - - -set.seed(1) -(total_time2 <- system.time({ - spassso_out <- spassso::spassso(Y[-testset,], - X[-testset,,drop=F], - cx[-testset,], - radgp_rho = radgp_rho, theta=theta_opts, - - spf_k = q, spf_a_delta = .1, spf_b_delta = .1, spf_a_dl = 0.9, - - spf_Lambda_start = matrix(rnorm(q*q), ncol=q) %>% perturb(), - spf_Delta_start = runif(q),#diag(Delta), - mvreg_B_start = Beta,# %>% perturb(), - - mcmc = mcmc <- 1000, - print_every=100, - - sample_precision=1, - sample_mvr=T, - sample_gp=T) -})) - -Sig_lmc <- spassso_out$S %>% apply(3, \(s) crossprod(s)) %>% array(dim=c(q,q,mcmc)) -Sig_lmc_est <- Sig_lmc %>% apply(1:2, mean) - diff --git a/spiox_prior_sample.r b/spiox_prior_sample.r deleted file mode 100644 index e3ac7df..0000000 --- a/spiox_prior_sample.r +++ /dev/null @@ -1,260 +0,0 @@ -rm(list=ls()) -library(tidyverse) -library(magrittr) -library(Matrix) -library(latex2exp) - -set.seed(4) - -image.matrix <- \(x) { - image(as(x, "dgCMatrix")) -} - -q <- 3 - -theta_mat <- matrix(1, ncol=q, nrow=4) -theta_mat[1,] <- c(5, 15, 30) -theta_mat[2,] <- c(1, 1, 1) -theta_mat[3,] <- c(0.5, 1.2, 1.9) -theta_mat[4,] <- rep(1e-3, q) - -sds <- c(1, 1, 1) -Omega <- matrix(c(1,-.9,0.7,-.9,1,-.5,0.7,-.5,1), ncol=3) -Sigma <- diag(sds) %*% Omega %*% diag(sds) - -St <- chol(Sigma) -S <- t(St) - -# spatial -xx <- seq(0, 1, length.out=200) -coords <- expand.grid(xx, xx) -nr <- nrow(coords) -cx <- as.matrix(coords) - -j <- 1 -radius <- 0.025 -test <- spiox::radgp_build(cx, radius, - phi=5, sigmasq=1, - nu=1.5, tausq=1e-5, - matern=F, 16) -test$dag %>% sapply(\(x) prod(dim(x))) %>% summary() - -Hlist <- 1:q %>% sapply(\(j) - spiox::radgp_build(cx, radius, - phi=theta_mat[1,j], sigmasq=theta_mat[2,j], - nu=theta_mat[3,j], tausq=theta_mat[4,j], - matern=T, 16)$H - ) - - -V <- matrix(rnorm(nr * q), ncol=q) %*% St - -Y <- V - -save_which <- rep(0, q) - -# make q spatial outcomes -which_process <- rep(0, q) -for(i in 1:q){ - cat(i, "\n") - Y[,i] <- Matrix::solve(Hlist[[i]], V[,i], sparse=T) -} - -df <- data.frame(coords, y=Y) %>% - pivot_longer(cols=-c(Var1, Var2)) %>% - mutate(name = ifelse(name == "y.1", "Outcome 1", ifelse(name == "y.2", "Outcome 2", "Outcome 3"))) - -( p2 <- ggplot(df, - aes(Var1, Var2, fill=value)) + - geom_raster() + - scico::scale_fill_scico(palette="vik") + # bam, broc, cork, managua, vik - facet_wrap(~name, ncol=q) + - theme_minimal() + - labs(x=NULL, y=NULL) + - scale_x_continuous(breaks=c(0.5, 1), expand=c(0,0)) + - scale_y_continuous(breaks=c(0, 0.5, 1), expand=c(0,0)) + - theme( - panel.grid = element_blank(), - panel.spacing.x = unit(10, "pt"), - panel.border = element_rect(fill=NA, color="black"), - axis.text.x = element_text(margin = margin(t = 0), hjust=1), - - axis.text.y = element_text(margin = margin(r = 0), vjust=1) ) ) - -ggsave("figures/prior_sample_3.pdf", plot=p2, width=8.7, height=3) - - - - - -## compute correlation functions -# we only want to go linearly down testx and testy making only those comparisons, hence diag_only=T - -library(spiox) -set.seed(1) -xx <- seq(0, 1, length.out=10) -coords <- expand.grid(xx, xx) -cx <- as.matrix(coords) -n_g <- 4 -iox_spcor <- function(outcome_i, outcome_j, h_values, - x_test_grid, - cx, theta_mat, - n_angles=5, diag_only=T, at_limit=T){ - - all_combos_setup <- expand.grid(x_test_grid, x_test_grid, - h_values, seq(0,2*pi,length.out=n_angles)) %>% as.matrix() - colnames(all_combos_setup) <- c("s1x", "s1y", "h", "angle") - - all_combos_setup %<>% data.frame() %>% - mutate(s2x = s1x + h*cos(angle), s2y = s1y + h*sin(angle)) - - testx <- all_combos_setup %>% dplyr::select(s1x, s1y) %>% as.matrix() - testy <- all_combos_setup %>% dplyr::select(s2x, s2y) %>% as.matrix() - - cov_computed <- iox(testx, testy, outcome_i, outcome_j, - cx, theta_mat, - matern = TRUE, - diag_only=diag_only, at_limit=at_limit) - - all_combos_setup %<>% mutate(spcov = cov_computed) - - spcov_by_h <- all_combos_setup %>% group_by(h) %>% summarise(spcov = mean(spcov)) - - return( spcov_by_h ) -} - - -xg_base <- xx[ seq(1, length(xx), length(xx)/(n_g+1)) ] %>% tail(-1) - - -## compute IOX covariance between two locations both in S^c -xg <- xg_base + runif(n_g, 0, 1e-2) -# limit variances at 0 -marg_cor_lims <- 1:q %>% - lapply(\(j) iox_spcor(j, j, c(0), xg, cx, theta_mat, n_angles=8, diag_only=T, at_limit=T) %>% - mutate(outcome=j) ) %>% - bind_rows() %>% - rename(cov0 = spcov) - -hvec <- c(10^seq(-5, -1, length.out=25), 10^seq(-1, -0.5, length.out=25)) - #c(seq(0, 0.05, length.out=10), seq(0.06, .5, length.out=20)) - -combis <- cbind(c(1,1), c(2,2), c(3,3), c(1,2), c(1,3), c(2,3)) - -cor_all_combs_Sc <- plyr::alply(combis, 2, \(x) { - i <- x[1] - j <- x[2] - - spcov_by_h <- iox_spcor(i, j, hvec, xg, cx, theta_mat, 10) %>% - mutate(outcome_1 = i, outcome_2 = j) %>% - left_join(marg_cor_lims %>% dplyr::select(-h), by=c("outcome_1"="outcome")) %>% - left_join(marg_cor_lims %>% dplyr::select(-h), by=c("outcome_2"="outcome")) %>% - mutate(spcor = #Sigma[i,j]/sqrt(Sigma[i,i]*Sigma[j,j]) * spcov/sqrt(cov0.x * cov0.y)) - spcov) - return(spcov_by_h) -}) %>% - bind_rows() %>% - mutate(model="IOX_Sc") - -## compute IOX covariance between two locations, one in S -xg <- xg_base -# limit variances at 0 -marg_cor_lims <- 1:q %>% - lapply(\(j) iox_spcor(j, j, c(0), xg, cx, theta_mat, n_angles=8, diag_only=T, at_limit=T) %>% - mutate(outcome=j) ) %>% - bind_rows() %>% - rename(cov0 = spcov) - -hvec <- c(10^seq(-5, -1, length.out=25), 10^seq(-1, -0.5, length.out=25)) -#c(seq(0, 0.05, length.out=10), seq(0.06, .5, length.out=20)) - -combis <- cbind(c(1,1), c(2,2), c(3,3), c(1,2), c(1,3), c(2,3)) - -cor_all_combs_S <- plyr::alply(combis, 2, \(x) { - i <- x[1] - j <- x[2] - - spcov_by_h <- iox_spcor(i, j, hvec, xg, cx, theta_mat, 10) %>% - mutate(outcome_1 = i, outcome_2 = j) %>% - left_join(marg_cor_lims %>% dplyr::select(-h), by=c("outcome_1"="outcome")) %>% - left_join(marg_cor_lims %>% dplyr::select(-h), by=c("outcome_2"="outcome")) %>% - mutate(spcor = Sigma[i,j] * #/sqrt(Sigma[i,i]*Sigma[j,j]) * spcov/sqrt(cov0.x * cov0.y)) - spcov) - return(spcov_by_h) -}) %>% - bind_rows() %>% - mutate(model="IOX_S") - -cor_all_combs <- cor_all_combs_S #bind_rows(cor_all_combs_Sc, cor_all_combs_S) - -ggplot(cor_all_combs %>% dplyr::filter(outcome_1 == 1, outcome_2 == 2), aes(h, spcor)) + - geom_line(aes(color=model)) + - theme_minimal() + - scale_y_continuous(limits=c(-1,1)) - - -# correlation functions used to generate IOX -rho_corrs <- 1:q %>% lapply(\(i) { - thetai <- theta_mat[,i] - thetai[4] <- 0 - rho_corrs <- - hvec %>% sapply(\(h) spiox::Correlationc(matrix(c(0,0), nrow=1), matrix(c(h, 0), nrow=1), - theta = thetai, covar = T, same=F) ) - return(data.frame(h=hvec, outcome_1 = i, outcome_2 = i, spcor = rho_corrs, model="Matern")) -}) %>% bind_rows() - - -corr_plot <- cor_all_combs %>% bind_rows(rho_corrs) - -corr_plot$model <- factor(corr_plot$model) -levels(corr_plot$model) <- c(TeX("IOX $l \\in S$"), TeX("IOX $l,l' \\in S^c$"), "Matern" ) - - - - - -################ WIP - - -cross_plots <- ggplot(corr_plot %>% dplyr::filter(grepl("IOX", model), outcome_1 != outcome_2), - aes(h, spcor)) + - geom_line(aes(color=model)) + - theme_minimal() + - scale_color_manual(values=c("black", "red"), labels = scales::parse_format()) + - facet_wrap(outcome_1 ~ outcome_2) + - scale_y_continuous(limits=c(-1,1)) - -marg_plots <- ggplot(corr_plot %>% dplyr::filter(outcome_1 == outcome_2), - aes(h, spcor)) + - geom_line(aes(color=model)) + - theme_minimal() + - scale_color_manual(values=c("black", "red", "blue"), labels = scales::parse_format()) + - facet_wrap(outcome_1 ~ outcome_2) + - scale_y_continuous(limits=c(0,1)) - - - - -library(gridExtra) -library(grid) - -grid.arrange( - p11, p22, p33, - p12, p13, p23, nrow=2) - -i <- 1; j <- 1 -spcov1_by_h <- iox_spcor(i, j, hvec, xx, cx, theta_mat, 8) %>% - mutate(outcome_1 = i, outcome_2 = j) %>% - left_join(marg_cor_lims %>% dplyr::select(-h), by=c("outcome_1"="outcome")) %>% - left_join(marg_cor_lims %>% dplyr::select(-h), by=c("outcome_2"="outcome")) %>% - mutate(spcor = Sigma[i,j]/sqrt(Sigma[i,i]*Sigma[j,j]) * spcov/sqrt(cov0.x * cov0.y)) - - - -spcov1_by_h %>% ggplot(aes(x=h, y=spcor)) + - geom_line() + - - - - - diff --git a/spiox_response.r b/spiox_response.r deleted file mode 100644 index 1cf1f1e..0000000 --- a/spiox_response.r +++ /dev/null @@ -1,332 +0,0 @@ -rm(list=ls()) -library(tidyverse) -library(magrittr) -library(Matrix) - -set.seed(2) - -image.matrix <- \(x) { - image(as(x, "dgCMatrix")) -} - -# spatial -xx <- seq(0, 1, length.out=50) -coords <- expand.grid(xx, xx) -cx <- as.matrix(coords) -nr <- nrow(cx) - -test_radgp <- spiox::radgp_build(cx, 0.12, phi=20, sigmasq=1, nu=1.5, tausq=0, matern=T, 16) -w <- solve(test_radgp$H, rnorm(nr)) - -Ctest <- spiox::Correlationc(cx, cx, c(1,1,1.9,1-8), 1, TRUE) - -q <- 5 - -optlist <- seq(0.5, 1.9, length.out=q) #%>% sample(q, replace=T) - -#Clist <- optlist %>% lapply(\(phi) (exp(-phi * as.matrix(dist(cx))^1) + 1e-6*diag(nr)) ) -Clist <- optlist %>% lapply(\(nu) spiox::Correlationc(cx, cx, c(20,1,nu,1e-12), 1, TRUE) ) -Llist <- Clist %>% lapply(\(C) t(chol(C))) -Lilist <- Llist %>% lapply(\(L) solve(L)) - - -Q <- rWishart(1, q+1, 1/2 * diag(q))[,,1] # -Sigma <- solve(Q) - -St <- chol(Sigma) -S <- t(St) - - - -V <- matrix(rnorm(nr * q), ncol=q) %*% St - -Y_sp <- V - -save_which <- rep(0, q) - -# make q spatial outcomes -which_process <- rep(0, q) -for(i in 1:q){ - which_process[i] <- i #1:length(optlist) %>% sample(1) - Y_sp[,i] <- Llist[[which_process[i]]] %*% V[,i] -} -theta_true <- optlist[which_process] - -# regression -p <- 2 -X <- matrix(1, ncol=1, nrow=nr) %>% cbind(matrix(rnorm(nr*(p-1)), ncol=p-1)) - -Beta <- matrix(rnorm(q * p), ncol=q) - -Y_regression <- X %*% Beta -Error <- matrix(rnorm(nr * q),ncol=q) %*% diag(D <- runif(q, 0, 0.1)) -Y <- as.matrix(Y_sp + Y_regression) + Error - -cov_Y_cols <- cov(Y) %>% as("dgCMatrix") -cov_Y_rows <- cov(t(Y)) %>% as("dgCMatrix") - -df <- data.frame(coords, y=as.matrix(Y_sp)) %>% - pivot_longer(cols=-c(Var1, Var2)) - -if(F){ - ggplot(df, - aes(Var1, Var2, fill=value)) + - geom_raster() + - scale_fill_viridis_c() + - facet_wrap(~name, ncol=5) -} - -############################## - -set.seed(1) - -phi_opts <- seq(5, 30, length.out=5) -nu_opts <- seq(0.5, 1.9, length.out=20) -optsall <- expand.grid(phi_opts, nu_opts) %>% as.matrix() - - -#theta_opts <- cbind(c(20, 1, .51, 1e-15), c(20, 1, 1.5, 1e-15)) -theta_opts <- rbind(optsall[,1], 1, optsall[,2], 1e-8) -#theta_opts <- par_opts %>% sapply(\(nu) matrix( c(20, 1, nu, 1e-6), ncol=1)) -#theta_opts <- par_opts %>% sapply(\(phi) matrix( c(phi, 1, 1, 1e-16), ncol=1)) - -testset <- sample(1:nrow(Y), 10, replace=F) - -perturb <- function(x, sd=1){ - return(x + matrix(rnorm(prod(dim(x)), sd), ncol=ncol(x))) -} -m_nn <- 20 - -cx_in <- cx[-testset,] -custom_dag <- dag_vecchia(cx_in, m_nn) - -set.seed(1) -(total_time <- system.time({ - spiox_out <- spiox::spiox_wishart(Y[-testset,,drop=F], - X[-testset,,drop=F], - cx_in, - custom_dag = custom_dag, - theta=theta_opts, - - Sigma_start = Sigma, - mvreg_B_start = Beta,# %>% perturb(), - - mcmc = mcmc <- 1500, - print_every = 100, - - sample_iwish=T, - sample_mvr=T, - sample_theta_gibbs=T, - upd_theta_opts=F, - num_threads = 16) -})) - - - -boom2 <- spiox::spiox_predict(X_new = X[testset,,drop=F], - coords_new = cx[testset,], - - # data - Y[-testset,], - X[-testset,,drop=F], - cx[-testset,], - - radgp_rho, - spiox_out$B %>% tail(c(NA, NA, round(mcmc/2))), - spiox_out$S %>% tail(c(NA, NA, round(mcmc/2))), - spiox_out$theta %>% tail(c(NA, NA, round(mcmc/2))), - num_threads = 16) - -Ytest <- boom2$Y %>% apply(1:2, mean) -Ytrue <- Y[testset,] - -1:q %>% sapply(\(j) cor(Ytest[,j], Ytrue[,j])) - - -spiox_out$theta %>% tail(c(NA,NA,round(mcmc/2))) %>% apply(1:2, mean) - -Sig_est <- round(mcmc/2):mcmc %>% sapply(\(m) with(spiox_out, diag(sqrt(theta[2,,m])) %*% - Sigma[,,m] %*% diag(sqrt(theta[2,,m])) ) ) %>% - array(c(q,q,mcmc)) %>% apply(1:2, mean) - -cbind(spiox_out$theta %>% apply(1:2, mean), theta_true) %>% data.frame() %>% arrange(theta_true) - -Cor_mcmc <- spiox_out$Sigma %>% tail(c(NA, NA, round(mcmc/2))) %>% apply(3, \(x) cov2cor(x)) %>% - array(dim=c(q,q,round(mcmc/2))) - -image(abs(cov2cor(Sigma) - apply(Cor_mcmc, 1:2, mean))) - -# microergodics -diag(Sigma) * theta_true -j <- 2 -plot(spiox_out$theta[j,] * - spiox_out$Sigma[j,j,], type='l') - - -df_test <- data.frame(cx[testset,], y=Wtest) %>% - pivot_longer(cols=-c(Var1, Var2)) -ggplot(df_test, aes(Var1, Var2, color=value)) + - geom_point() + - facet_grid(~name) + - scale_color_viridis_c() + - theme_minimal() - - - - - - -mesh_total_time <- system.time({ - meshout <- meshed::spmeshed(y=Y, - x=X, - coords=cx, k = q, - block_size = 50, - n_samples = 5000, - n_burn = 1, - n_thin = 1, - n_threads = 16, - prior = list(phi=c(1, 20)), - verbose=50 - )}) -h <- 0 -mesh_mcmc <- dim(meshout$theta_mcmc)[3] - -mesh_h_f <- function(h){ - sig_mcmc <- 1:mesh_mcmc %>% sapply(\(m) with(meshout, - lambda_mcmc[,,m] %*% diag(exp(-theta_mcmc[1,,m]*h)) %*% t(lambda_mcmc[,,m]) + diag(tausq_mcmc[,m]) - )) %>% array(dim=c(3,3,mesh_mcmc)) - return(sig_mcmc %>% apply(1:2, mean)) -} - - - - -set.seed(1) -(total_time2 <- system.time({ - spassso_out <- spassso::spassso(Y[-testset,], - X[-testset,,drop=F], - cx[-testset,], - radgp_rho = radgp_rho, theta=theta_opts, - - spf_k = kfit, spf_a_delta = .1, spf_b_delta = .1, spf_a_dl = 0.5, - - spf_Lambda_start = Lambda[,1:kfit] %>% perturb(), - spf_Delta_start = runif(q),#diag(Delta), - mvreg_B_start = Beta,# %>% perturb(), - - mcmc = mcmc <- 1000, - print_every=100, - - sample_precision=1, - sample_mvr=T, - sample_gp=T) -})) - - - -spf_test <- spiox::run_spf_model(Y, kfit, .1, .1, 0.5, - Lambda[,1:kfit] %>% perturb(), runif(q), - mcmc = mcmc, 100, seq_lambda = F) - - - -Qc_samples <- 1:mcmc %>% - sapply(\(i) with(spf_test, cov2cor(tcrossprod(Lambda[,,i]) + diag(Delta[,i]) ))) %>% - array(dim=c(q, q, mcmc)) - -Qc_post_spf <- Qc_samples %>% apply(1:2, mean) - -mean((Qc_post_spf-cov2cor(Q))^2) - - - -Q_est_pattern <- Q_samples %>% tail(c(NA,NA,round(3*mcmc/4))) %>% - apply(1:2, \(x) 1*(sign(quantile(x, 0.05))==sign(quantile(x, 0.95)))) - -Q_true_pattern <- 1*(Q!=0) - -# percentage of correct off-diagonal nonzeros -check_pattern <- Q_true_pattern == Q_est_pattern -mean(check_pattern) -#(sum(check_pattern) - q)/(prod(dim(Q)) - q) - - -S_samples <- 1:mcmc %>% sapply(\(i) - with(spiox_out, solve(tcrossprod(Lambda[,,i]) + diag(Delta[,i])) ) ) %>% - array(dim=c(qstar,qstar,mcmc)) - -B_reg <- 1:mcmc %>% sapply(\(i) { - S <- with(spiox_out, solve(tcrossprod(Lambda[,,i]) + diag(Delta[,i])) ) - return( S[-(1:pstar),1:pstar] %*% solve(S[1:pstar,1:pstar]) ) -} -) %>% - array(dim=c(q,pstar,mcmc)) - -image(B_reg %>% apply(1:2, mean)) -image(Bstar) - - -B_true_pattern <- 1*(Bstar!=0) -B_est_pattern <- B_reg %>% tail(c(NA,NA,round(3*mcmc/4))) %>% - apply(1:2, \(x) 1*(sign(quantile(x, 0.05))==sign(quantile(x, 0.95)))) - -# percentage of correct nonzeros -check_pattern_B <- B_est_pattern == B_true_pattern -mean(check_pattern_B) - - -## setting zeros - - -eps <- 0.05 - -chandra_fdr <- function(eps, Qc_samples, beta=0.9){ - H1 <- Qc_samples %>% apply(1:2, \(x) abs(x)>eps) - H0 <- Qc_samples %>% apply(1:2, \(x) abs(x)% apply(2:3, mean) - post_H0 <- H0 %>% apply(2:3, mean) - - dij <- 1*(post_H1>beta) - - fdr <- sum(dij*post_H0)/max(c(sum(dij), 1)) - - return(fdr) -} - -eps_grid <- seq(0.05, 0.7, length.out=100) - -fdr_grid <- eps_grid %>% sapply(chandra_fdr, Qc_samples %>% tail(c(NA, NA, 1000))) - -plot(eps_grid, fdr_grid, type='l') - - -Qc_est_pattern <- abs(Qc_samples %>% tail(c(NA,NA,round(1000))) %>% - apply(1:2, mean)) > 0.3 - -Q_true_pattern <- 1*(Q!=0) - -image(Qc_est_pattern) -image(Q_true_pattern) -mean(Q_true_pattern == Qc_est_pattern) - - - - - - - - - - - - - - - - -preds <- ROCR::prediction(Q_est_pattern[lower.tri(Q_est_pattern)], Q_true_pattern[lower.tri(Q_true_pattern)]) -(perf <- ROCR::performance(preds, "fpr")@y.values[[1]][2]) - - - diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp index 0b0ac4f..019ba2b 100644 --- a/src/RcppExports.cpp +++ b/src/RcppExports.cpp @@ -215,7 +215,7 @@ BEGIN_RCPP END_RCPP } // spiox_wishart -Rcpp::List spiox_wishart(const arma::mat& Y, const arma::mat& X, const arma::mat& coords, const arma::field& custom_dag, arma::mat theta_opts, const arma::mat& Sigma_start, const arma::mat& mvreg_B_start, int mcmc, int print_every, bool matern, bool sample_iwish, bool sample_mvr, bool sample_theta_gibbs, bool upd_theta_opts, int num_threads); +Rcpp::List spiox_wishart(const arma::mat& Y, const arma::mat& X, const arma::mat& coords, const arma::field& custom_dag, arma::mat theta_opts, const arma::mat& Sigma_start, const arma::mat& mvreg_B_start, int mcmc, int print_every, int matern, bool sample_iwish, bool sample_mvr, bool sample_theta_gibbs, bool upd_theta_opts, int num_threads); RcppExport SEXP _spiox_spiox_wishart(SEXP YSEXP, SEXP XSEXP, SEXP coordsSEXP, SEXP custom_dagSEXP, SEXP theta_optsSEXP, SEXP Sigma_startSEXP, SEXP mvreg_B_startSEXP, SEXP mcmcSEXP, SEXP print_everySEXP, SEXP maternSEXP, SEXP sample_iwishSEXP, SEXP sample_mvrSEXP, SEXP sample_theta_gibbsSEXP, SEXP upd_theta_optsSEXP, SEXP num_threadsSEXP) { BEGIN_RCPP Rcpp::RObject rcpp_result_gen; @@ -229,7 +229,7 @@ BEGIN_RCPP Rcpp::traits::input_parameter< const arma::mat& >::type mvreg_B_start(mvreg_B_startSEXP); Rcpp::traits::input_parameter< int >::type mcmc(mcmcSEXP); Rcpp::traits::input_parameter< int >::type print_every(print_everySEXP); - Rcpp::traits::input_parameter< bool >::type matern(maternSEXP); + Rcpp::traits::input_parameter< int >::type matern(maternSEXP); Rcpp::traits::input_parameter< bool >::type sample_iwish(sample_iwishSEXP); Rcpp::traits::input_parameter< bool >::type sample_mvr(sample_mvrSEXP); Rcpp::traits::input_parameter< bool >::type sample_theta_gibbs(sample_theta_gibbsSEXP); @@ -240,7 +240,7 @@ BEGIN_RCPP END_RCPP } // spiox_latent -Rcpp::List spiox_latent(const arma::mat& Y, const arma::mat& X, const arma::mat& coords, const arma::field& custom_dag, arma::mat theta_opts, const arma::mat& Sigma_start, const arma::mat& mvreg_B_start, int mcmc, int print_every, bool matern, bool sample_iwish, bool sample_mvr, bool sample_theta_gibbs, bool upd_theta_opts, int num_threads, int sampling); +Rcpp::List spiox_latent(const arma::mat& Y, const arma::mat& X, const arma::mat& coords, const arma::field& custom_dag, arma::mat theta_opts, const arma::mat& Sigma_start, const arma::mat& mvreg_B_start, int mcmc, int print_every, int matern, bool sample_iwish, bool sample_mvr, bool sample_theta_gibbs, bool upd_theta_opts, int num_threads, int sampling); RcppExport SEXP _spiox_spiox_latent(SEXP YSEXP, SEXP XSEXP, SEXP coordsSEXP, SEXP custom_dagSEXP, SEXP theta_optsSEXP, SEXP Sigma_startSEXP, SEXP mvreg_B_startSEXP, SEXP mcmcSEXP, SEXP print_everySEXP, SEXP maternSEXP, SEXP sample_iwishSEXP, SEXP sample_mvrSEXP, SEXP sample_theta_gibbsSEXP, SEXP upd_theta_optsSEXP, SEXP num_threadsSEXP, SEXP samplingSEXP) { BEGIN_RCPP Rcpp::RObject rcpp_result_gen; @@ -254,7 +254,7 @@ BEGIN_RCPP Rcpp::traits::input_parameter< const arma::mat& >::type mvreg_B_start(mvreg_B_startSEXP); Rcpp::traits::input_parameter< int >::type mcmc(mcmcSEXP); Rcpp::traits::input_parameter< int >::type print_every(print_everySEXP); - Rcpp::traits::input_parameter< bool >::type matern(maternSEXP); + Rcpp::traits::input_parameter< int >::type matern(maternSEXP); Rcpp::traits::input_parameter< bool >::type sample_iwish(sample_iwishSEXP); Rcpp::traits::input_parameter< bool >::type sample_mvr(sample_mvrSEXP); Rcpp::traits::input_parameter< bool >::type sample_theta_gibbs(sample_theta_gibbsSEXP); @@ -266,7 +266,7 @@ BEGIN_RCPP END_RCPP } // spiox_predict -Rcpp::List spiox_predict(const arma::mat& X_new, const arma::mat& coords_new, const arma::mat& Y, const arma::mat& X, const arma::mat& coords, const arma::field& dag, const arma::cube& B, const arma::cube& Sigma, const arma::cube& theta, bool matern, int num_threads); +Rcpp::List spiox_predict(const arma::mat& X_new, const arma::mat& coords_new, const arma::mat& Y, const arma::mat& X, const arma::mat& coords, const arma::field& dag, const arma::cube& B, const arma::cube& Sigma, const arma::cube& theta, int matern, int num_threads); RcppExport SEXP _spiox_spiox_predict(SEXP X_newSEXP, SEXP coords_newSEXP, SEXP YSEXP, SEXP XSEXP, SEXP coordsSEXP, SEXP dagSEXP, SEXP BSEXP, SEXP SigmaSEXP, SEXP thetaSEXP, SEXP maternSEXP, SEXP num_threadsSEXP) { BEGIN_RCPP Rcpp::RObject rcpp_result_gen; @@ -280,14 +280,36 @@ BEGIN_RCPP Rcpp::traits::input_parameter< const arma::cube& >::type B(BSEXP); Rcpp::traits::input_parameter< const arma::cube& >::type Sigma(SigmaSEXP); Rcpp::traits::input_parameter< const arma::cube& >::type theta(thetaSEXP); - Rcpp::traits::input_parameter< bool >::type matern(maternSEXP); + Rcpp::traits::input_parameter< int >::type matern(maternSEXP); Rcpp::traits::input_parameter< int >::type num_threads(num_threadsSEXP); rcpp_result_gen = Rcpp::wrap(spiox_predict(X_new, coords_new, Y, X, coords, dag, B, Sigma, theta, matern, num_threads)); return rcpp_result_gen; END_RCPP } +// spiox_predict_part +Rcpp::List spiox_predict_part(const arma::mat& Y_new, const arma::mat& X_new, const arma::mat& coords_new, const arma::mat& Y, const arma::mat& X, const arma::mat& coords, const arma::field& dag, const arma::cube& B, const arma::cube& Sigma, const arma::cube& theta, int matern, int num_threads); +RcppExport SEXP _spiox_spiox_predict_part(SEXP Y_newSEXP, SEXP X_newSEXP, SEXP coords_newSEXP, SEXP YSEXP, SEXP XSEXP, SEXP coordsSEXP, SEXP dagSEXP, SEXP BSEXP, SEXP SigmaSEXP, SEXP thetaSEXP, SEXP maternSEXP, SEXP num_threadsSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< const arma::mat& >::type Y_new(Y_newSEXP); + Rcpp::traits::input_parameter< const arma::mat& >::type X_new(X_newSEXP); + Rcpp::traits::input_parameter< const arma::mat& >::type coords_new(coords_newSEXP); + Rcpp::traits::input_parameter< const arma::mat& >::type Y(YSEXP); + Rcpp::traits::input_parameter< const arma::mat& >::type X(XSEXP); + Rcpp::traits::input_parameter< const arma::mat& >::type coords(coordsSEXP); + Rcpp::traits::input_parameter< const arma::field& >::type dag(dagSEXP); + Rcpp::traits::input_parameter< const arma::cube& >::type B(BSEXP); + Rcpp::traits::input_parameter< const arma::cube& >::type Sigma(SigmaSEXP); + Rcpp::traits::input_parameter< const arma::cube& >::type theta(thetaSEXP); + Rcpp::traits::input_parameter< int >::type matern(maternSEXP); + Rcpp::traits::input_parameter< int >::type num_threads(num_threadsSEXP); + rcpp_result_gen = Rcpp::wrap(spiox_predict_part(Y_new, X_new, coords_new, Y, X, coords, dag, B, Sigma, theta, matern, num_threads)); + return rcpp_result_gen; +END_RCPP +} // spiox_latent_predict -Rcpp::List spiox_latent_predict(const arma::mat& X_new, const arma::mat& coords_new, const arma::mat& coords, const arma::field& dag, const arma::cube& W, const arma::cube& B, const arma::cube& Sigma, const arma::mat& Dvec, const arma::cube& theta, bool matern, int num_threads); +Rcpp::List spiox_latent_predict(const arma::mat& X_new, const arma::mat& coords_new, const arma::mat& coords, const arma::field& dag, const arma::cube& W, const arma::cube& B, const arma::cube& Sigma, const arma::mat& Dvec, const arma::cube& theta, int matern, int num_threads); RcppExport SEXP _spiox_spiox_latent_predict(SEXP X_newSEXP, SEXP coords_newSEXP, SEXP coordsSEXP, SEXP dagSEXP, SEXP WSEXP, SEXP BSEXP, SEXP SigmaSEXP, SEXP DvecSEXP, SEXP thetaSEXP, SEXP maternSEXP, SEXP num_threadsSEXP) { BEGIN_RCPP Rcpp::RObject rcpp_result_gen; @@ -301,7 +323,7 @@ BEGIN_RCPP Rcpp::traits::input_parameter< const arma::cube& >::type Sigma(SigmaSEXP); Rcpp::traits::input_parameter< const arma::mat& >::type Dvec(DvecSEXP); Rcpp::traits::input_parameter< const arma::cube& >::type theta(thetaSEXP); - Rcpp::traits::input_parameter< bool >::type matern(maternSEXP); + Rcpp::traits::input_parameter< int >::type matern(maternSEXP); Rcpp::traits::input_parameter< int >::type num_threads(num_threadsSEXP); rcpp_result_gen = Rcpp::wrap(spiox_latent_predict(X_new, coords_new, coords, dag, W, B, Sigma, Dvec, theta, matern, num_threads)); return rcpp_result_gen; @@ -326,6 +348,7 @@ static const R_CallMethodDef CallEntries[] = { {"_spiox_spiox_wishart", (DL_FUNC) &_spiox_spiox_wishart, 15}, {"_spiox_spiox_latent", (DL_FUNC) &_spiox_spiox_latent, 16}, {"_spiox_spiox_predict", (DL_FUNC) &_spiox_spiox_predict, 11}, + {"_spiox_spiox_predict_part", (DL_FUNC) &_spiox_spiox_predict_part, 12}, {"_spiox_spiox_latent_predict", (DL_FUNC) &_spiox_spiox_latent_predict, 11}, {NULL, NULL, 0} }; diff --git a/src/cholesky_lrupd.h b/src/cholesky_lrupd.h deleted file mode 100644 index 52db34e..0000000 --- a/src/cholesky_lrupd.h +++ /dev/null @@ -1,31 +0,0 @@ -#include -using namespace std; - -inline void c_uchol_upd_r1(int n, double *U, double *z) { - int i, j; - double r, c, s, temp; - double alpha=1; - - for (i = 0; i < n; ++i) { - r = sqrt(U[i * n + i] * U[i * n + i] + alpha * z[i] * z[i]); - c = r / U[i * n + i]; - s = z[i] / U[i * n + i]; - U[i * n + i] = r; - - for (j = i + 1; j < n; ++j) { - U[j * n + i] = (U[j * n + i] + alpha * s * z[j]) / c; - z[j] = c * z[j] - s * U[j * n + i]; - } - } -} - -inline void uchol_update(arma::mat& U, const arma::mat& V){ - // Goal: return upper cholesky of U^T U + VV^T - // where U is the upper chol of an (n x n) spd matrix - // where V is a (n x p) matrix - unsigned int n = U.n_rows; - for(unsigned int i=0; i(z.memptr())); - } -} diff --git a/src/daggp.cpp b/src/daggp.cpp index e0e5295..47f6ed9 100755 --- a/src/daggp.cpp +++ b/src/daggp.cpp @@ -20,7 +20,7 @@ DagGP::DagGP( oneuv = arma::ones(1); matern = covariance_matern; // 0 pexp or 1 matern or 2 wave - + //Rcpp::Rcout << "DagGP covariance choice " << matern << endl; //thread safe stuff n_threads = num_threads_in; diff --git a/src/dirichletlaplace.cpp b/src/dirichletlaplace.cpp deleted file mode 100644 index f4c8ad0..0000000 --- a/src/dirichletlaplace.cpp +++ /dev/null @@ -1,45 +0,0 @@ -#include "dirichletlaplace.h" - -using namespace std; - -arma::vec r_psi_cond(const arma::vec& phi, const arma::vec& theta, double tau){ - unsigned int p = phi.n_elem; - arma::vec psi = arma::zeros(p); - for(unsigned int i=0; i -//#include -//#include -#include - - -#define ZTOL (8.*DBL_EPSILON) - -/*****************************************************************************/ -/* Privat Functions */ -/*****************************************************************************/ - -double _gig_mode(double lambda, double omega) -/*---------------------------------------------------------------------------*/ -/* Compute mode of GIG distribution. */ -/* */ -/* Parameters: */ -/* lambda .. parameter for distribution */ -/* omega ... parameter for distribution */ -/* */ -/* Return: */ -/* mode */ -/*---------------------------------------------------------------------------*/ -{ - if (lambda >= 1.) - /* mode of fgig(x) */ - return (sqrt((lambda-1.)*(lambda-1.) + omega*omega)+(lambda-1.))/omega; - else - /* 0 <= lambda < 1: use reciprocal of mode of f(1/x) */ - return omega / (sqrt((1.-lambda)*(1.-lambda) + omega*omega)+(1.-lambda)); -} /* end of _gig_mode() */ - -/*---------------------------------------------------------------------------*/ - -void _rgig_ROU_noshift (arma::vec & res, int n, double lambda, double lambda_old, double omega, double alpha) -/*---------------------------------------------------------------------------*/ -/* Tpye 1: */ -/* Ratio-of-uniforms without shift. */ -/* Dagpunar (1988), Sect.~4.6.2 */ -/* Lehner (1989) */ -/*---------------------------------------------------------------------------*/ -{ - double xm, nc; /* location of mode; c=log(f(xm)) normalization constant */ - double ym, um; /* location of maximum of x*sqrt(f(x)); umax of MBR */ - double s, t; /* auxiliary variables */ - double U, V, X; /* random variables */ - - int i; /* loop variable (number of generated random variables) */ - int count = 0; /* counter for total number of iterations */ - - /* -- Setup -------------------------------------------------------------- */ - - /* shortcuts */ - t = 0.5 * (lambda-1.); - s = 0.25 * omega; - - /* mode = location of maximum of sqrt(f(x)) */ - xm = _gig_mode(lambda, omega); - - /* normalization constant: c = log(sqrt(f(xm))) */ - nc = t*log(xm) - s*(xm + 1./xm); - - /* location of maximum of x*sqrt(f(x)): */ - /* we need the positive root of */ - /* omega/2*y^2 - (lambda+1)*y - omega/2 = 0 */ - ym = ((lambda+1.) + sqrt((lambda+1.)*(lambda+1.) + omega*omega))/omega; - - /* boundaries of minmal bounding rectangle: */ - /* we us the "normalized" density f(x) / f(xm). hence */ - /* upper boundary: vmax = 1. */ - /* left hand boundary: umin = 0. */ - /* right hand boundary: umax = ym * sqrt(f(ym)) / sqrt(f(xm)) */ - um = exp(0.5*(lambda+1.)*log(ym) - s*(ym + 1./ym) - nc); - - /* -- Generate sample ---------------------------------------------------- */ - - for (i=0; i (t*log(X) - s*(X + 1./X) - nc))); - - /* store random point */ - res(i) = (lambda_old < 0.) ? (alpha / X) : (alpha * X); - } - - /* -- End ---------------------------------------------------------------- */ - - return; -} /* end of _rgig_ROU_noshift() */ - - -/*---------------------------------------------------------------------------*/ - -void _rgig_newapproach1 (arma::vec & res, int n, double lambda, double lambda_old, double omega, double alpha) -/*---------------------------------------------------------------------------*/ -/* Type 4: */ -/* New approach, constant hat in log-concave part. */ -/* Draw sample from GIG distribution. */ -/* */ -/* Case: 0 < lambda < 1, 0 < omega < 1 */ -/* */ -/* Parameters: */ -/* n ....... sample size (positive integer) */ -/* lambda .. parameter for distribution */ -/* omega ... parameter for distribution */ -/* */ -/* Return: */ -/* random sample of size 'n' */ -/*---------------------------------------------------------------------------*/ -{ - /* parameters for hat function */ - double A[3], Atot; /* area below hat */ - double k0; /* maximum of PDF */ - double k1, k2; /* multiplicative constant */ - - double xm; /* location of mode */ - double x0; /* splitting point T-concave / T-convex */ - double a; /* auxiliary variable */ - - double U, V, X; /* random numbers */ - double hx; /* hat at X */ - - int i; /* loop variable (number of generated random variables) */ - int count = 0; /* counter for total number of iterations */ - - /* -- Check arguments ---------------------------------------------------- */ - - if (lambda >= 1. || omega >1.) - error ("invalid parameters"); - - /* -- Setup -------------------------------------------------------------- */ - - /* mode = location of maximum of sqrt(f(x)) */ - xm = _gig_mode(lambda, omega); - - /* splitting point */ - x0 = omega/(1.-lambda); - - /* domain [0, x_0] */ - k0 = exp((lambda-1.)*log(xm) - 0.5*omega*(xm + 1./xm)); /* = f(xm) */ - A[0] = k0 * x0; - - /* domain [x_0, Infinity] */ - if (x0 >= 2./omega) { - k1 = 0.; - A[1] = 0.; - k2 = pow(x0, lambda-1.); - A[2] = k2 * 2. * exp(-omega*x0/2.)/omega; - } - - else { - /* domain [x_0, 2/omega] */ - k1 = exp(-omega); - A[1] = (lambda == 0.) - ? k1 * log(2./(omega*omega)) - : k1 / lambda * ( pow(2./omega, lambda) - pow(x0, lambda) ); - - /* domain [2/omega, Infinity] */ - k2 = pow(2/omega, lambda-1.); - A[2] = k2 * 2 * exp(-1.)/omega; - } - - /* total area */ - Atot = A[0] + A[1] + A[2]; - - /* -- Generate sample ---------------------------------------------------- */ - - for (i=0; i 2./omega) ? x0 : 2./omega; - X = -2./omega * log(exp(-omega/2. * a) - omega/(2.*k2) * V); - hx = k2 * exp(-omega/2. * X); - break; - - } while(0); - - /* accept or reject */ - U = unif_rand() * hx; - - if (log(U) <= (lambda-1.) * log(X) - omega/2. * (X+1./X)) { - /* store random point */ - res(i) = (lambda_old < 0.) ? (alpha / X) : (alpha * X); - break; - } - } while(1); - - } - - /* -- End ---------------------------------------------------------------- */ - - return; -} /* end of _rgig_newapproach1() */ - -/*---------------------------------------------------------------------------*/ - -void -_rgig_ROU_shift_alt (arma::vec & res, int n, double lambda, double lambda_old, double omega, double alpha) -/*---------------------------------------------------------------------------*/ -/* Type 8: */ -/* Ratio-of-uniforms with shift by 'mode', alternative implementation. */ -/* Dagpunar (1989) */ -/* Lehner (1989) */ -/*---------------------------------------------------------------------------*/ -{ - double xm, nc; /* location of mode; c=log(f(xm)) normalization constant */ - double s, t; /* auxiliary variables */ - double U, V, X; /* random variables */ - - int i; /* loop variable (number of generated random variables) */ - int count = 0; /* counter for total number of iterations */ - - double a, b, c; /* coefficent of cubic */ - double p, q; /* coefficents of depressed cubic */ - double fi, fak; /* auxiliary results for Cardano's rule */ - - double y1, y2; /* roots of (1/x)*sqrt(f((1/x)+m)) */ - - double uplus, uminus; /* maximum and minimum of x*sqrt(f(x+m)) */ - - /* -- Setup -------------------------------------------------------------- */ - - /* shortcuts */ - t = 0.5 * (lambda-1.); - s = 0.25 * omega; - - /* mode = location of maximum of sqrt(f(x)) */ - xm = _gig_mode(lambda, omega); - - /* normalization constant: c = log(sqrt(f(xm))) */ - nc = t*log(xm) - s*(xm + 1./xm); - - /* location of minimum and maximum of (1/x)*sqrt(f(1/x+m)): */ - - /* compute coeffients of cubic equation y^3+a*y^2+b*y+c=0 */ - a = -(2.*(lambda+1.)/omega + xm); /* < 0 */ - b = (2.*(lambda-1.)*xm/omega - 1.); - c = xm; - - /* we need the roots in (0,xm) and (xm,inf) */ - - /* substitute y=z-a/3 for depressed cubic equation z^3+p*z+q=0 */ - p = b - a*a/3.; - q = (2.*a*a*a)/27. - (a*b)/3. + c; - - /* use Cardano's rule */ - fi = acos(-q/(2.*sqrt(-(p*p*p)/27.))); - fak = 2.*sqrt(-p/3.); - y1 = fak * cos(fi/3.) - a/3.; - y2 = fak * cos(fi/3. + 4./3.*M_PI) - a/3.; - - /* boundaries of minmal bounding rectangle: */ - /* we us the "normalized" density f(x) / f(xm). hence */ - /* upper boundary: vmax = 1. */ - /* left hand boundary: uminus = (y2-xm) * sqrt(f(y2)) / sqrt(f(xm)) */ - /* right hand boundary: uplus = (y1-xm) * sqrt(f(y1)) / sqrt(f(xm)) */ - uplus = (y1-xm) * exp(t*log(y1) - s*(y1 + 1./y1) - nc); - uminus = (y2-xm) * exp(t*log(y2) - s*(y2 + 1./y2) - nc); - - /* -- Generate sample ---------------------------------------------------- */ - - for (i=0; i (t*log(X) - s*(X + 1./X) - nc))); - - /* store random point */ - res(i) = (lambda_old < 0.) ? (alpha / X) : (alpha * X); - } - - /* -- End ---------------------------------------------------------------- */ - - return; -} /* end of _rgig_ROU_shift_alt() */ - -/*---------------------------------------------------------------------------*/ - -double -_unur_bessel_k_nuasympt (double x, double nu, int islog, int expon_scaled) -/*---------------------------------------------------------------------------*/ -/* Asymptotic expansion of Bessel K_nu(x) function */ -/* when BOTH nu and x are large. */ -/* */ -/* parameters: */ -/* x ... argument for K_nu() */ -/* nu ... order or Bessel function */ -/* islog ... return logarithm of result TRUE and result when FALSE */ -/* expon_scaled ... return exp(-x)*K_nu(x) when TRUE and K_nu(x) when FALSE*/ -/* */ -/*---------------------------------------------------------------------------*/ -/* */ -/* references: */ -/* ## Abramowitz & Stegun , p.378, __ 9.7.8. __ */ -/* */ -/* ## K_nu(nu * z) ~ sqrt(pi/(2*nu)) * exp(-nu*eta)/(1+z^2)^(1/4) */ -/* ## * {1 - u_1(t)/nu + u_2(t)/nu^2 - ... } */ -/* */ -/* ## where t = 1 / sqrt(1 + z^2), */ -/* ## eta = sqrt(1 + z^2) + log(z / (1 + sqrt(1+z^2))) */ -/* ## */ -/* ## and u_k(t) from p.366 __ 9.3.9 __ */ -/* */ -/* ## u0(t) = 1 */ -/* ## u1(t) = (3*t - 5*t^3)/24 */ -/* ## u2(t) = (81*t^2 - 462*t^4 + 385*t^6)/1152 */ -/* ## ... */ -/* */ -/* ## with recursion 9.3.10 for k = 0, 1, .... : */ -/* ## */ -/* ## u_{k+1}(t) = t^2/2 * (1 - t^2) * u'_k(t) + */ -/* ## 1/8 \int_0^t (1 - 5*s^2)* u_k(s) ds */ -/*---------------------------------------------------------------------------*/ -/* */ -/* Original implementation in R code (R package "Bessel" v. 0.5-3) by */ -/* Martin Maechler, Date: 23 Nov 2009, 13:39 */ -/* */ -/* Translated into C code by Kemal Dingic, Oct. 2011. */ -/* */ -/* Modified by Josef Leydold on Tue Nov 1 13:22:09 CET 2011 */ -/* */ -/*---------------------------------------------------------------------------*/ -{ -#define M_LNPI 1.14472988584940017414342735135 /* ln(pi) */ - - double z; /* rescaled argument for K_nu() */ - double sz, t, t2, eta; /* auxiliary variables */ - double d, u1t,u2t,u3t,u4t; /* (auxiliary) results for Debye polynomials */ - double res; /* value of log(K_nu(x)) [= result] */ - - /* rescale: we comute K_nu(z * nu) */ - z = x / nu; - - /* auxiliary variables */ - sz = hypot(1,z); /* = sqrt(1+z^2) */ - t = 1. / sz; - t2 = t*t; - - eta = (expon_scaled) ? (1./(z + sz)) : sz; - eta += log(z) - log1p(sz); /* = log(z/(1+sz)) */ - - /* evaluate Debye polynomials u_j(t) */ - u1t = (t * (3. - 5.*t2))/24.; - u2t = t2 * (81. + t2*(-462. + t2 * 385.))/1152.; - u3t = t*t2 * (30375. + t2 * (-369603. + t2 * (765765. - t2 * 425425.)))/414720.; - u4t = t2*t2 * (4465125. - + t2 * (-94121676. - + t2 * (349922430. - + t2 * (-446185740. - + t2 * 185910725.)))) / 39813120.; - d = (-u1t + (u2t + (-u3t + u4t/nu)/nu)/nu)/nu; - - /* log(K_nu(x)) */ - res = log(1.+d) - nu*eta - 0.5*(log(2.*nu*sz) - M_LNPI); - - return (islog ? res : exp(res)); -} /* end of _unur_bessel_k_nuasympt() */ - -/*---------------------------------------------------------------------------*/ -/*---------------------------------------------------------------------------*/ -/* Draw sample from GIG distribution. */ -/* without calling GetRNGstate() ... PutRNGstate() */ -/* */ -/* Parameters: */ -/* n ....... sample size (positive integer) */ -/* lambda .. parameter for distribution */ -/* chi ... parameter for distribution */ -/* psi ... parameter for distribution */ -/* */ -/* Return: */ -/* random sample of size 'n' */ -/*---------------------------------------------------------------------------*/ - - -arma::vec c_rgigauss_vec(int n, double lambda, double chi, double psi) -{ - double omega, alpha; /* parameters of standard distribution */ - double abs_lambda; /* absolute value of lambda */ - arma::vec res = arma::zeros(n); /* results */ - int i; - - /* check sample size */ - if (n<0) { - error("sample size 'n' must be non-negative integer."); - } - - /* check GIG parameters: */ - if ( !(R_FINITE(lambda) && R_FINITE(chi) && R_FINITE(psi)) || - (chi < 0. || psi < 0) || - (chi == 0. && lambda <= 0.) || - (psi == 0. && lambda >= 0.) ) { - error("invalid parameters for GIG distribution: lambda=%g, chi=%g, psi=%g", - lambda, chi, psi); - } - - /* alternative parametrization */ - omega = sqrt(psi*chi); - alpha = sqrt(chi/psi); - abs_lambda = (lambda >= 0.) ? lambda : -lambda; - - /* run generator */ - do { - - if (omega < ZTOL) { - /* for very small values of omega we have serious round-off errors */ - if (lambda > 0.0) { - /* special cases which are basically Gamma distribution */ - for (i=0; i 2. || omega > 3.) { - /* Ratio-of-uniforms with shift by 'mode', alternative implementation */ - _rgig_ROU_shift_alt(res, n, abs_lambda, lambda, omega, alpha); - break; - } - - if (abs_lambda >= 1.-2.25*omega*omega || omega > 0.2) { - /* Ratio-of-uniforms without shift */ - _rgig_ROU_noshift(res, n, abs_lambda, lambda, omega, alpha); - break; - } - - if (omega > 0.) { /* remaining case */ - /* New approach, constant hat in log-concave part. */ - _rgig_newapproach1(res, n, abs_lambda, lambda, omega, alpha); - break; - } - - else - error("parameters must satisfy lambda>=0 and omega>0."); - - /* Remark: case lambda == 0.0 and very small omega may cause numerical problems */ - - } while (0); - - /* return result */ - //UNPROTECT(1); - return res; - -} /* end of c_rgigauss_vec() */ - - \ No newline at end of file diff --git a/src/rinvgauss.cpp b/src/rinvgauss.cpp deleted file mode 100644 index 4e3345f..0000000 --- a/src/rinvgauss.cpp +++ /dev/null @@ -1,36 +0,0 @@ -#include - -using namespace Rcpp; -using namespace std; - - -arma::vec c_rigauss_vec(int n, double mu, double lambda){ - arma::vec random_vector(n); - double z,y,x,u; - for(int i=0; i - -double c_rigauss(double mu, double lambda); -arma::vec c_rgigauss_vec(int n, double lambda, double chi, double psi); \ No newline at end of file diff --git a/src/spf.cpp b/src/spf.cpp deleted file mode 100644 index 220e2b3..0000000 --- a/src/spf.cpp +++ /dev/null @@ -1,136 +0,0 @@ -#include "spf.h" - -void SparsePrecisionFactor::compute_P(){ - P = Iq + Lambda.t() * arma::diagmat(1.0/Delta) * Lambda; - try { - Pchol = arma::chol(P, "lower"); - } catch (...) { - Rcpp::Rcout << "error at P" << endl - << P << endl; - Rcpp::Rcout << 1.0/Delta << endl; - Rcpp::stop("stop"); - } - arma::mat Pcholi = arma::inv(arma::trimatl(Pchol)); - Pi = Pcholi.t() * Pcholi; -} - -void SparsePrecisionFactor::fc_sample_uv(){ - compute_P(); - - //Rcpp::Rcout << "u " << arma::size(Pchol) << endl; - U = arma::randn(n, q) * Pchol.t(); - //Rcpp::Rcout << "v " << arma::size(Delta) << endl; - V = (*Y) + U * Pi * Lambda.t() * arma::diagmat(1.0/Delta); -} - -/* -void SparsePrecisionFactor::fc_sample_Lambda(){ - arma::mat rn = arma::randn(d, q); - arma::mat VtV = V.t() * V; - - // update Lambda by column (parallel OK) -#ifdef _OPENMP -#pragma omp parallel for -#endif - for(unsigned int i=0; i(0, d-1); - arma::rowvec vssq = arma::sum(V%V, 0); - // update Lambda by row - for(unsigned int i=0; i(0, d-1); - arma::rowvec vssq = arma::sum(V%V, 0); - arma::mat UVL = U - V*Lambda; - - // update Lambda by row - for(unsigned int i=0; in_rows; - d = Y->n_cols; - Iq = arma::eye(q,q); - - S_Lambda = arma::ones(d, q); - - } -}; diff --git a/src/spf_standalone.cpp b/src/spf_standalone.cpp deleted file mode 100644 index 4a8202e..0000000 --- a/src/spf_standalone.cpp +++ /dev/null @@ -1,59 +0,0 @@ -#include "spf.h" -using namespace std; - - -//[[Rcpp::export]] -Rcpp::List run_spf_model(arma::mat& Y, - unsigned int n_factors, - double delta_gamma_shape, - double delta_gamma_rate, - double dl_dirichlet_a, - const arma::mat& Lambda_start, - const arma::vec& Delta_start, - unsigned int mcmc=1000, - int print_every=1000, - bool seq_lambda=false){ - - SparsePrecisionFactor spf(&Y, n_factors, - delta_gamma_shape, delta_gamma_rate, dl_dirichlet_a); - - arma::cube Lambda_storage = arma::zeros(spf.d, spf.q, mcmc); - arma::mat Delta_storage = arma::zeros(spf.d, mcmc); - arma::cube Lambda_dlvar_storage = arma::zeros(spf.d, spf.q, mcmc); - - spf.Lambda_start(Lambda_start); - spf.Delta_start(Delta_start); - - for(unsigned int m=0; m0); - if(print_condition){ - print_condition = print_condition & (!(m % print_every)); - }; - if(print_condition){ - Rcpp::Rcout << "Iteration: " << m+1 << " of " << mcmc << endl; - } - } - - return Rcpp::List::create( - Rcpp::Named("Lambda") = Lambda_storage, - Rcpp::Named("Delta") = Delta_storage, - Rcpp::Named("Lambda_dlvar") = Lambda_dlvar_storage - ); -} \ No newline at end of file diff --git a/src/spiox.cpp b/src/spiox.cpp index 53e33ee..e3787ef 100644 --- a/src/spiox.cpp +++ b/src/spiox.cpp @@ -2,7 +2,7 @@ #include "interrupt.h" //[[Rcpp::export]] -Rcpp::List spiox_wishart(const arma::mat& Y, +Rcpp::List spiox_response(const arma::mat& Y, const arma::mat& X, const arma::mat& coords, @@ -15,7 +15,7 @@ Rcpp::List spiox_wishart(const arma::mat& Y, int mcmc = 1000, int print_every = 100, - bool matern = true, + int matern = 1, bool sample_iwish = true, bool sample_mvr = true, bool sample_theta_gibbs = true, @@ -131,7 +131,7 @@ Rcpp::List spiox_latent(const arma::mat& Y, int mcmc=1000, int print_every=100, - bool matern = true, + int matern = 1, bool sample_iwish=true, bool sample_mvr=true, bool sample_theta_gibbs=true, @@ -141,7 +141,7 @@ Rcpp::List spiox_latent(const arma::mat& Y, if(sampling==0){ - Rcpp::stop("Run the GP-IOX response model via spiox_wishart()"); + Rcpp::stop("Run the GP-IOX response model via spiox_response()"); } Rcpp::Rcout << "GP-IOX latent model, "; diff --git a/src/spiox.h b/src/spiox.h index b984fe6..e5d8861 100644 --- a/src/spiox.h +++ b/src/spiox.h @@ -1,8 +1,6 @@ #include "omp_import.h" #include "daggp.h" #include "ramadapt.h" -#include "spf.h" -#include "cholesky_lrupd.h" using namespace std; @@ -22,7 +20,6 @@ class SpIOX { // matrix of predictors dim n, p arma::mat X; - // metadata unsigned int n, q, p; int num_threads; @@ -32,10 +29,7 @@ class SpIOX { arma::mat B; arma::mat YXB; arma::mat B_Var; // prior variance on B, element by element - double B_a_dl; // dirichlet-laplace parameter for vec(B) - // SPF for sparse latent precision - SparsePrecisionFactor spf; arma::mat S, Si, Sigma, Q; // S^T * S = Sigma = Q^-1 = (Lambda*Lambda^T + Delta)^-1 = (Si * Si^T)^-1 // RadGP for spatial dependence @@ -47,13 +41,13 @@ class SpIOX { arma::mat V; // -------------- utilities - bool matern; + int matern; void sample_B(); // - void sample_Q_spf(); void sample_Sigma_wishart(); void compute_V(); // whitened void sample_theta_discr(); // gibbs for each outcome choosing from options void upd_theta_metrop(); + void upd_theta_metrop_conditional(); // n_options == q void init_theta_adapt(); // latent model @@ -68,15 +62,23 @@ class SpIOX { // bool phi_sampling, sigmasq_sampling, nu_sampling, tausq_sampling; + // adaptive metropolis to update theta atoms int theta_mcmc_counter; arma::uvec which_theta_elem; arma::mat theta_unif_bounds; - arma::mat theta_metrop_sd; + //arma::mat theta_metrop_sd; RAMAdapt theta_adapt; bool theta_adapt_active; // -------- + // adaptive metropolis (conditional update) to update theta atoms + // assume shared covariance functions and unknown parameters across variables + arma::mat c_theta_unif_bounds; + std::vector c_theta_adapt; + // -------- + + // -------------- run 1 gibbs iteration based on current values void gibbs(int it, int sample_sigma, bool sample_mvr, bool sample_theta_gibbs, bool upd_theta_opts); double logdens_eval(); @@ -94,7 +96,7 @@ class SpIOX { const arma::mat& Sigma_start, const arma::mat& mvreg_B_start, - bool use_matern, + int cov_model_matern, int num_threads_in) { num_threads = num_threads_in; @@ -129,7 +131,8 @@ class SpIOX { nu_sampling = arma::var(theta_options.row(2)) != 0; tausq_sampling = arma::var(theta_options.row(3)) != 0; - matern = use_matern; + matern = cov_model_matern; + //Rcpp::Rcout << "Covariance choice: " << matern << endl; for(unsigned int i=0; i(0); arma::uvec oneuv = arma::ones(1); + + if(q == 1) { + nu_sampling = true; + } + if(phi_sampling){ which_theta_elem = arma::join_vert(which_theta_elem, 0*oneuv); } @@ -176,7 +184,7 @@ inline void SpIOX::init_theta_adapt(){ if(tausq_sampling){ which_theta_elem = arma::join_vert(which_theta_elem, 3*oneuv); } - + int n_theta_par = n_options * which_theta_elem.n_elem; arma::mat bounds_all = arma::zeros(4, 2); // make bounds for all, then subset @@ -197,10 +205,22 @@ inline void SpIOX::init_theta_adapt(){ theta_unif_bounds = arma::join_vert(theta_unif_bounds, bounds_all); } - theta_metrop_sd = 0.05 * arma::eye(n_theta_par, n_theta_par); + arma::mat theta_metrop_sd = 0.05 * arma::eye(n_theta_par, n_theta_par); theta_adapt = RAMAdapt(n_theta_par, theta_metrop_sd, 0.24); theta_adapt_active = true; - // --- + + if(n_options == q){ + // conditional update + c_theta_unif_bounds = bounds_all; + int c_theta_par = which_theta_elem.n_elem; + arma::mat c_theta_metrop_sd = 0.05 * arma::eye(c_theta_par, c_theta_par); + c_theta_adapt = std::vector(n_options); + for(int j=0; j(1); + for(int j=0; j0){ + V_alt.col(j) = daggp_options_alt.at(spmap(j)).H_times_A(W.col(j));// * (Y.col(j) - X * B.col(j)); + } else { + V_alt.col(j) = daggp_options_alt.at(spmap(j)).H_times_A(YXB.col(j));// * (Y.col(j) - X * B.col(j)); + } + + double c_daggp_logdet = daggp_options.at(spmap(j)).precision_logdeterminant; + double c_daggp_alt_logdet = daggp_options_alt.at(spmap(j)).precision_logdeterminant; + + arma::vec Vjc = arma::zeros(n); + arma::vec Vjc_alt = arma::zeros(n); + for(int jc=0; jc 0){ - if(sample_sigma == 1){ - tstart = std::chrono::steady_clock::now(); - sample_Q_spf(); - timings(2) += time_count(tstart); - } else { - tstart = std::chrono::steady_clock::now(); - sample_Sigma_wishart(); - timings(2) += time_count(tstart); - } + tstart = std::chrono::steady_clock::now(); + sample_Sigma_wishart(); + timings(2) += time_count(tstart); } if(sample_mvr){ @@ -706,7 +783,11 @@ inline void SpIOX::gibbs(int it, int sample_sigma, bool sample_mvr, bool sample_ // update atoms for theta tstart = std::chrono::steady_clock::now(); if(upd_theta_opts){ - upd_theta_metrop(); + if( !sample_theta_gibbs & (q>3) & (n_options == q) ){ + upd_theta_metrop_conditional(); + } else { + upd_theta_metrop(); + } } timings(4) += time_count(tstart); diff --git a/src/spiox_predict.cpp b/src/spiox_predict.cpp index ec27d6a..2fe7d5b 100644 --- a/src/spiox_predict.cpp +++ b/src/spiox_predict.cpp @@ -14,7 +14,7 @@ Rcpp::List spiox_predict( const arma::cube& B, const arma::cube& Sigma, const arma::cube& theta, - bool matern = true, + int matern = 1, int num_threads = 1 ){ int q = Y.n_cols; @@ -51,7 +51,7 @@ Rcpp::List spiox_predict( for(int m=0; m::from( CC(j) - CPt(j).t() * ht )) ); // abs for numerical zeros - Y_out(i, j) = xb_new(i, j) + + Y_out(j) = xb_new(j) + arma::conv_to::from(ht.t() * yxb_old(px, outjx)); // post mean only Dj(j) = sqrtR; } - Y_out.row(i) += arma::trans( arma::diagmat(Dj) * Sigmauchol.t() * rndnorm_m.col(i) ); // pred uncert - Y_out_mcmc.subcube(i, 0, m, i, q-1, m) = arma::trans(Y_out.row(i)); + Y_out += arma::diagmat(Dj) * Sigmauchol.t() * rndnorm_m.col(i) ; // pred uncert + Y_out_mcmc.subcube(i, 0, m, i, q-1, m) = Y_out; theta_current = theta_m; } @@ -102,6 +102,158 @@ Rcpp::List spiox_predict( ); } + +//[[Rcpp::export]] +Rcpp::List spiox_predict_part( + const arma::mat& Y_new, + const arma::mat& X_new, + const arma::mat& coords_new, + + const arma::mat& Y, + const arma::mat& X, + const arma::mat& coords, + + const arma::field& dag, + + const arma::cube& B, + const arma::cube& Sigma, + const arma::cube& theta, + int matern = 1, + int num_threads = 1 +){ + // Y_new is a matrix with as many rows as X_new and coords_new + // and q columns. NA in Y_new are to be predicted + int q = Y.n_cols; + int p = X.n_cols; + int mcmc = B.n_slices; + + arma::uvec oneuv = arma::ones(1); + + int bessel_ws_inc = MAT_NU_MAX;//see bessel_k.c for working space needs + double * bessel_ws = (double *) R_alloc(num_threads*bessel_ws_inc, sizeof(double)); + + int ntrain = coords.n_rows; + int ntest = coords_new.n_rows; + + arma::mat cxall = arma::join_vert(coords, coords_new); + + arma::cube Y_out_mcmc = arma::zeros(ntest, q, mcmc); + arma::cube random_stdnormal = arma::randn(mcmc, q, ntest); + + // loop over test locations +#ifdef _OPENMP +#pragma omp parallel for num_threads(num_threads) +#endif + for(int i=0; i CC(q); + arma::field CPt(q); + arma::field PPi(q); + + for(int m=0; m 1e-15){ + CC(j) = Correlationf(cxall, ix, ix, theta_m.col(j), bessel_ws, matern, true); // 1 for matern + CPt(j) = Correlationf(cxall, px, ix, theta_m.col(j), bessel_ws, matern, false); + PPi(j) = arma::inv_sympd( Correlationf(cxall, px, px, theta_m.col(j), bessel_ws, matern, true) ); + } + } + + //Rcpp::Rcout << "building obvs" << endl; + for(int jx=0; jx::from( + CC(j) - CPt(j).t() * ht )) ); // abs for numerical zeros + Yo_resid(jx) = Y_new(i, j) - xb_new(j) - + arma::conv_to::from(ht.t() * yxb_old(px, outjx)); // post mean only + Djo(jx) = sqrtR; + } + + //Rcpp::Rcout << "building misss" << endl; + for(int jx=0; jx::from( + CC(j) - CPt(j).t() * ht )) ); // abs for numerical zeros + Ym_spatl(jx) = xb_new(j) + + arma::conv_to::from(ht.t() * yxb_old(px, outjx)); // post mean only + Djm(jx) = sqrtR; + } + + //Rcpp::Rcout << "finally" << endl; + + arma::mat S_mo = Sigmam(missing, observd); + arma::mat S_o_inv = arma::inv_sympd(Sigmam(observd, observd)); + arma::mat S_m = Sigmam(missing, missing); + + arma::mat Hmo = arma::diagmat(Djm) * S_mo * S_o_inv * arma::diagmat(1.0/Djo); + arma::mat Rmo = arma::diagmat(Djm) * (S_m - S_mo * S_o_inv * S_mo.t()) * arma::diagmat(Djm); + + arma::vec missing_post_mean = Hmo * Yo_resid; + arma::mat missing_post_cholvar = arma::chol(arma::symmatu(Rmo), "lower"); + + //Rcpp::Rcout << "finally 3" << endl; + arma::mat rnorm_all = random_stdnormal.subcube(m, 0, i, m, q-1, i) ; + arma::vec rnorm_miss = arma::trans( rnorm_all.cols(missing) ); + arma::vec Ym_pred = Ym_spatl + missing_post_mean + missing_post_cholvar * rnorm_miss; + + for(int jx=0; jx% lapply(\(i) (exp(-philist[i] * as.matrix(dist(cx))^explist[i]) ) ) -Llist <- Clist %>% lapply(\(C) t(chol(C))) - -#Q <- matrix(c(2,1,0,1,2,1,0,1,2), ncol=3) -#Sigma <- solve(Q) -Sigma <- rWishart(1, 10, diag(k))[,,1] %>% cov2cor() -#Sigma <- matrix(c(2,1.8,1.8,2), ncol=2) -#Sigma <- diag(k) -A <- t(chol(Sigma)) - -I <- \(x) diag(x) - -Lblocks <- Matrix::bdiag(Llist) - -Cblocks <- Lblocks %*% (Sigma %x% I(nr)) %*% t(Lblocks) - -#Ciblocks <- solve(t(Lblocks)) %*% (solve(Sigma) %x% I(nr)) %*% solve(Lblocks) -#C_LMC <- (A %x% I(nr)) %*% Matrix::bdiag(Clist) %*% (t(A) %x% I(nr)) -#Ci_LMC <- solve(C_LMC) - -v <- rnorm(nr*k) -y <- Lblocks %*% (A %x% I(nr)) %*% v -w <- (A %x% I(nr)) %*% Lblocks %*% v - - -Y <- matrix(y, ncol=k) -W <- matrix(w, ncol=k) - -df <- data.frame(coords, y=W) %>% - pivot_longer(cols=c(-Var1, -Var2)) - -ggplot(df, aes(Var1, Var2, fill=value)) + - geom_raster() + - scale_fill_viridis_c() + - facet_grid(~name) - - -newc <- c(0.124, 0.531) -Cnlist <- philist %>% lapply(\(phi) (exp(-phi * as.matrix(dist(rbind(cx,newc)))^1) ) ) -Lnlist <- Cnlist %>% lapply(\(C) t(chol(C))) -Lnblocks <- Matrix::bdiag(Lnlist) - - -S <- solve(Q) -B <- t(chol(S)) - -C_big <- Lnblocks %*% (S %x% I(nr+1)) %*% t(Lnblocks) -#C_big <- (B %x% I(nr+1)) %*% Matrix::bdiag(Cnlist) %*% (t(B) %x% I(nr+1)) -#C_big <- S %x% Cnlist[[1]] - -outix <- c(26, 52, 78) - -H <- C_big[outix, -outix] %*% solve(C_big[-outix, -outix]) -R <- C_big[outix, outix] - H %*% C_big[-outix, outix] - -H %*% as.vector(y) - -ych <- y -#ych[(nr+1):(k*nr)] <- ych[(nr+1):(k*nr)]*2 -#H %*% as.vector(ych) - - - - -Lnlist[[1]] -solve(Llist[[1]], Cnlist[[1]][1:nr, nr+1]) - - -Hcheck <- Matrix::bdiag(Lnlist %>% lapply(\(L) head(tail(L, 1), c(NA,-1)))) %*% - (Sigma %x% I(nr)) %*% - #Matrix::bdiag(Lnlist %>% lapply(\(L) t(head(L, c(nr, nr))))) %*% Ciblocks - Matrix::bdiag(Llist %>% lapply(\(L) t(L))) %*% Ciblocks - -Hcheck %*% as.vector(y) - - -Cnewold <- matrix(tail(Cnlist[[1]],1)[1:nr], ncol=1) - -- Cnewold %*% solve(Clist[[1]]) * -as.numeric(sqrt(1 - Cnewold %*% solve(Clist[[1]]) %*% t(Cnewold))) - - - - -outH <- function(C){ - return( - -C[nr+1,1:nr] %*% solve(C[1:nr,1:nr]) * (1/sqrt(C[nr+1, nr+1] - C[nr+1,1:nr] %*% solve(C[1:nr,1:nr]) %*% C[1:nr, nr+1])[1,1]) - ) -} - -C <- Cnlist[[1]] -ImH <- c(-C[nr+1,1:nr] %*% solve(C[1:nr,1:nr]), 1) -Rihalf <- 1/sqrt(C[nr+1, nr+1] - C[nr+1,1:nr] %*% solve(C[1:nr,1:nr]) %*% C[1:nr, nr+1])[1,1] - - -L1 <- outH(Cnlist[[1]]) -L2 <- outH(Cnlist[[2]]) -L3 <- outH(Cnlist[[3]]) -Matrix::bdiag(L1, L2, L3) - - - - - -