TMB fitTMB::compile("naomi_simple.cpp")## [1] 0[1] 0dyn.load(TMB::dynlib("naomi_simple"))
# tmb <- readRDS("depends/tmb.rds")
@@ -2939,35 +2961,35 @@ 2 Point concentration proportiona
)
base_grid
## Grid for Numerical Integration
-## # dimensions: 24 # gridpoints: 2 used memory:832 bytes
-##
-## nD-Construction principle: 'product'
-## individual quadrature rule for each dimension
-## dim = 1 --> type: GHe level: 1 initial domain: -Inf Inf
-## dim = 2 --> type: GHe level: 1 initial domain: -Inf Inf
-## dim = 3 --> type: GHe level: 1 initial domain: -Inf Inf
-## dim = 4 --> type: GHe level: 1 initial domain: -Inf Inf
-## dim = 5 --> type: GHe level: 1 initial domain: -Inf Inf
-## dim = 6 --> type: GHe level: 1 initial domain: -Inf Inf
-## dim = 7 --> type: GHe level: 1 initial domain: -Inf Inf
-## dim = 8 --> type: GHe level: 1 initial domain: -Inf Inf
-## dim = 9 --> type: GHe level: 1 initial domain: -Inf Inf
-## dim = 10 --> type: GHe level: 1 initial domain: -Inf Inf
-## dim = 11 --> type: GHe level: 1 initial domain: -Inf Inf
-## dim = 12 --> type: GHe level: 1 initial domain: -Inf Inf
-## dim = 13 --> type: GHe level: 1 initial domain: -Inf Inf
-## dim = 14 --> type: GHe level: 1 initial domain: -Inf Inf
-## dim = 15 --> type: GHe level: 1 initial domain: -Inf Inf
-## dim = 16 --> type: GHe level: 2 initial domain: -Inf Inf
-## dim = 17 --> type: GHe level: 1 initial domain: -Inf Inf
-## dim = 18 --> type: GHe level: 1 initial domain: -Inf Inf
-## dim = 19 --> type: GHe level: 1 initial domain: -Inf Inf
-## dim = 20 --> type: GHe level: 1 initial domain: -Inf Inf
-## dim = 21 --> type: GHe level: 1 initial domain: -Inf Inf
-## dim = 22 --> type: GHe level: 1 initial domain: -Inf Inf
-## dim = 23 --> type: GHe level: 1 initial domain: -Inf Inf
-## dim = 24 --> type: GHe level: 1 initial domain: -Inf InfGrid for Numerical Integration
+ # dimensions: 24 # gridpoints: 2 used memory:832 bytes
+
+nD-Construction principle: 'product'
+ individual quadrature rule for each dimension
+ dim = 1 --> type: GHe level: 1 initial domain: -Inf Inf
+ dim = 2 --> type: GHe level: 1 initial domain: -Inf Inf
+ dim = 3 --> type: GHe level: 1 initial domain: -Inf Inf
+ dim = 4 --> type: GHe level: 1 initial domain: -Inf Inf
+ dim = 5 --> type: GHe level: 1 initial domain: -Inf Inf
+ dim = 6 --> type: GHe level: 1 initial domain: -Inf Inf
+ dim = 7 --> type: GHe level: 1 initial domain: -Inf Inf
+ dim = 8 --> type: GHe level: 1 initial domain: -Inf Inf
+ dim = 9 --> type: GHe level: 1 initial domain: -Inf Inf
+ dim = 10 --> type: GHe level: 1 initial domain: -Inf Inf
+ dim = 11 --> type: GHe level: 1 initial domain: -Inf Inf
+ dim = 12 --> type: GHe level: 1 initial domain: -Inf Inf
+ dim = 13 --> type: GHe level: 1 initial domain: -Inf Inf
+ dim = 14 --> type: GHe level: 1 initial domain: -Inf Inf
+ dim = 15 --> type: GHe level: 1 initial domain: -Inf Inf
+ dim = 16 --> type: GHe level: 2 initial domain: -Inf Inf
+ dim = 17 --> type: GHe level: 1 initial domain: -Inf Inf
+ dim = 18 --> type: GHe level: 1 initial domain: -Inf Inf
+ dim = 19 --> type: GHe level: 1 initial domain: -Inf Inf
+ dim = 20 --> type: GHe level: 1 initial domain: -Inf Inf
+ dim = 21 --> type: GHe level: 1 initial domain: -Inf Inf
+ dim = 22 --> type: GHe level: 1 initial domain: -Inf Inf
+ dim = 23 --> type: GHe level: 1 initial domain: -Inf Inf
+ dim = 24 --> type: GHe level: 1 initial domain: -Inf Inf The total number of points for this grid is 2. Let’s try fitting this, and see how long it takes. First prepare the data:
@@ -3030,13 +3052,13 @@Now let’s fit the model with aghq, and time how long it takes:
start <- Sys.time()
quad <- fit_aghq(tmb_inputs, k = 1, basegrid = base_grid)## Warning: '.T2Cmat' is deprecated.
-## Use '.T2CR' instead.
-## See help("Deprecated") and help("Matrix-deprecated").Warning: '.T2Cmat' is deprecated.
+Use '.T2CR' instead.
+See help("Deprecated") and help("Matrix-deprecated").end <- Sys.time()
end - start## Time difference of 1.329712 minsTime difference of 58.11896 secsVerifying the number of grid points used:
stopifnot(nrow(quad$modesandhessians) == 2)How are these grid points spaced across the 24 dimensions?1
@@ -3048,37 +3070,37 @@
Interesting to note that for the Cholesky adapted grid with levels = c(1, 1, ..., 2, 1, ..., 1) there is variation in all of the dimensions greater than or equal to the 16th.
This is related to the Cholesky adaption being based on a lower triangular matrix.
To verify this, we can look at the differences between the two grid points.
For many of the dimensions the difference is too small to see on the plot above.
t(unname(diff(as.matrix(quad$modesandhessians[, 1:24]))))## [,1]
-## [1,] 0.000000e+00
-## [2,] 0.000000e+00
-## [3,] 0.000000e+00
-## [4,] 0.000000e+00
-## [5,] 0.000000e+00
-## [6,] 0.000000e+00
-## [7,] 0.000000e+00
-## [8,] 0.000000e+00
-## [9,] 0.000000e+00
-## [10,] 0.000000e+00
-## [11,] 0.000000e+00
-## [12,] 0.000000e+00
-## [13,] 0.000000e+00
-## [14,] 0.000000e+00
-## [15,] 0.000000e+00
-## [16,] 3.814905e+00
-## [17,] 1.606655e+00
-## [18,] 1.518063e-03
-## [19,] -3.494600e-03
-## [20,] -7.841190e-06
-## [21,] 4.996001e-05
-## [22,] -2.818483e-04
-## [23,] -6.447179e-03
-## [24,] -1.368464e-04 [,1]
+ [1,] 0.000000e+00
+ [2,] 0.000000e+00
+ [3,] 0.000000e+00
+ [4,] 0.000000e+00
+ [5,] 0.000000e+00
+ [6,] 0.000000e+00
+ [7,] 0.000000e+00
+ [8,] 0.000000e+00
+ [9,] 0.000000e+00
+[10,] 0.000000e+00
+[11,] 0.000000e+00
+[12,] 0.000000e+00
+[13,] 0.000000e+00
+[14,] 0.000000e+00
+[15,] 0.000000e+00
+[16,] 3.814905e+00
+[17,] 1.606655e+00
+[18,] 1.518063e-03
+[19,] -3.494600e-03
+[20,] -7.841190e-06
+[21,] 4.996001e-05
+[22,] -2.818483e-04
+[23,] -6.447179e-03
+[24,] -1.368464e-04Let’s also verify that adaption is happening as we expect it to be:
base_grid_copy <- rlang::duplicate(base_grid)
@@ -3091,31 +3113,31 @@ 2 Point concentration proportiona
Notice that there is a difference for all of the 24 dimenions when dec.type = 1 (spectral) is used:
mvQuad::rescale(base_grid_copy, m = quad$optresults$mode, C = Matrix::forceSymmetric(solve(quad$optresults$hessian)), dec.type = 1)
t(diff(mvQuad::getNodes(base_grid_copy)))
-## [,1]
-## [1,] -4.562110e-03
-## [2,] -4.788531e-03
-## [3,] -5.313556e-03
-## [4,] 3.876109e-01
-## [5,] 6.925962e-04
-## [6,] -1.955305e-03
-## [7,] 5.923101e-03
-## [8,] -1.301535e-02
-## [9,] 0.000000e+00
-## [10,] 8.809706e-03
-## [11,] -9.468319e-02
-## [12,] 2.367271e-02
-## [13,] -3.194831e-01
-## [14,] -1.485009e-02
-## [15,] 3.214202e-02
-## [16,] -7.588094e-02
-## [17,] 1.752933e-01
-## [18,] -1.072913e-02
-## [19,] 2.496620e-04
-## [20,] 3.893256e-05
-## [21,] -3.087050e-04
-## [22,] 4.770154e-03
-## [23,] 4.556001e-01
-## [24,] -3.400338e-03
+ [,1]
+ [1,] -4.562110e-03
+ [2,] -4.788531e-03
+ [3,] -5.313556e-03
+ [4,] 3.876109e-01
+ [5,] 6.925962e-04
+ [6,] -1.955305e-03
+ [7,] 5.923101e-03
+ [8,] -1.301535e-02
+ [9,] 0.000000e+00
+[10,] 8.809706e-03
+[11,] -9.468319e-02
+[12,] 2.367271e-02
+[13,] -3.194831e-01
+[14,] -1.485009e-02
+[15,] 3.214202e-02
+[16,] -7.588094e-02
+[17,] 1.752933e-01
+[18,] -1.072913e-02
+[19,] 2.496620e-04
+[20,] 3.893256e-05
+[21,] -3.087050e-04
+[22,] 4.770154e-03
+[23,] 4.556001e-01
+[24,] -3.400338e-03
This model fit in a very reasonable amount of time.
More generally, it would be of interest to know more about the relationship between the number of points in the grid, and the length of time it takes run fit_aghq.
Are there any theoretical considerations statements that can be made e.g. scaling \(\mathcal{O}(n)\) where \(n\) is the number of grid points?
@@ -3138,7 +3160,7 @@ 3 Eigendecomposition
E <- eigenC$vectorsSuch that C can be obtained by \(E \Lambda E^\top\), or in code E %*% diag(lambda) %*% t(E):
max(C - (E %*% Lambda %*% t(E))) < 10E-12## [1] TRUE[1] TRUEThe relative contributions of each principle component are given by \(\lambda_i / \sum_i \lambda_i\):
plot_total_variation(eigenC, label_x = 20)This grid has \(3^8 = 6561\) nodes:
(n_nodes <- nrow(mvQuad::getNodes(gg_s)))## [1] 6561[1] 6561How does the reconstruction of the covariance matrix with 8 components look?
E_s <- E[, 1:s]
Lambda_s <- Lambda[1:s, 1:s]
@@ -3195,16 +3217,16 @@ 3.2 Scoping the application to Na
gg3 <- pca_rescale(m = m, C = C, s = 3, k = 3)
nrow(mvQuad::getNodes(gg3))
## [1] 27[1] 27gg5 <- pca_rescale(m = m, C = C, s = 5, k = 3)
nrow(mvQuad::getNodes(gg5))## [1] 243[1] 243gg7 <- pca_rescale(m = m, C = C, s = 7, k = 3)
nrow(mvQuad::getNodes(gg7))## [1] 2187[1] 2187gg9 <- pca_rescale(m = m, C = C, s = 9, k = 3)
nrow(mvQuad::getNodes(gg9))## [1] 19683[1] 19683
Want to use a particular coordinate system, but still the covariance matrix scaling as in adaptive quadrature.
+m <- c(1, 1.5)
+C <- matrix(c(2, 1, 1, 1), ncol = 2)
+
+obj <- function(theta) {
+ mvtnorm::dmvnorm(theta, mean = m, sigma = C, log = TRUE)
+}
+
+ff <- list(
+ fn = obj,
+ gr = function(theta) numDeriv::grad(obj, theta),
+ he = function(theta) numDeriv::hessian(obj, theta)
+)
+
+grid <- expand.grid(
+ theta1 = seq(-6, 6, length.out = 600),
+ theta2 = seq(-6, 6, length.out = 600)
+)
+
+ground_truth <- cbind(grid, pdf = exp(obj(grid)))
+
+plot <- ggplot(ground_truth, aes(x = theta1, y = theta2, z = pdf)) +
+ geom_contour(col = multi.utils::cbpalette()[1]) +
+ coord_fixed(xlim = c(-6, 6), ylim = c(-6, 6), ratio = 1) +
+ labs(x = "", y = "") +
+ theme_minimal()
+
+grid_5 <- mvQuad::createNIGrid(dim = 2, type = "GHe", level = 5)
+mvQuad::rescale(grid_5, m = m, C = C, dec.type = 1)
+
+pca_grid_5 <- pca_rescale(m = m, C = C, s = 1, k = 5)
+
+plot_points <- function(gg) {
+ plot +
+ geom_point(
+ data = mvQuad::getNodes(gg) %>%
+ as.data.frame() %>%
+ mutate(weights = mvQuad::getWeights(gg)),
+ aes(x = V1, y = V2, size = weights),
+ alpha = 0.8,
+ col = multi.utils::cbpalette()[2],
+ inherit.aes = FALSE
+ ) +
+ scale_size_continuous(range = c(1, 2)) +
+ labs(x = "", y = "", size = "Weight") +
+ theme_minimal()
+}
+
+plot_points(grid_5) + plot_points(pca_grid_5) + plot_annotation(tag_levels = "A")
Try to create a grid similar to A above but using the second principal component rather than the first:
+eigenC <- eigen(C)
+E <- eigenC$vectors
+lambda <- eigenC$values
+
+gg_s <- mvQuad::createNIGrid(dim = 1, type = "GHe", level = 5)
+nodes_out <- t(E[, 2] %*% diag(sqrt(lambda[2]), ncol = 1) %*% t(mvQuad::getNodes(gg_s)))
+for(j in 1:2) nodes_out[, j] <- nodes_out[, j] + m[j] # Adaption
+weights_out <- mvQuad::getWeights(gg_s) * as.numeric(mvQuad::getWeights(mvQuad::createNIGrid(dim = 1, type = "GHe", level = 1)))
+weights_out <- det(chol(C)) * weights_out # Adaption
+
+custom_grid_5 <- mvQuad::createNIGrid(dim = 2, type = "GHe", level = 1)
+custom_grid_5$level <- rep(NA, times = 2)
+custom_grid_5$ndConstruction <- "custom"
+custom_grid_5$nodes <- nodes_out
+custom_grid_5$weights <- weights_out
+
+plot_points(grid_5) / plot_points(pca_grid_5) / plot_points(custom_grid_5) + plot_annotation(tag_levels = "A")
ggsave("pca-demo.png", h = 7, w = 6.25)We know that \(C = E \Lambda E^\top = \sum \lambda_i \mathbf{e}_i \mathbf{e}_i^\top\)
+lambda[1] * E[, 1] %*% t(E[, 1]) + lambda[2] * E[, 2] %*% t(E[, 2]) [,1] [,2]
+[1,] 2 1
+[2,] 1 1We also know that any vector can be written in terms of the eigendecomposition.
+For example, the vector c(1, 1) can be written as:
a <- c(1, 1) %*% E
+
+a[1] * E[, 1] + a[2] * E[, 2][1] 1 1Remove any nodes which have small enough weight.
+See Jäckel (2005).
+R45 <- matrix(
+ c(cos(pi / 2), sin(pi / 2), -sin(pi / 2), cos(pi / 2)), ncol = 2, nrow = 2
+)#' TODO!Remove any nodes which have small enough weight (Jäckel 2005).
+#' TODO!sessionInfo()R version 4.2.0 (2022-04-22)
+Platform: x86_64-apple-darwin17.0 (64-bit)
+Running under: macOS 13.3.1
+
+Matrix products: default
+LAPACK: /Library/Frameworks/R.framework/Versions/4.2/Resources/lib/libRlapack.dylib
+
+locale:
+[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
+
+attached base packages:
+[1] stats graphics grDevices utils datasets methods base
+
+other attached packages:
+ [1] multi.utils_0.1.0 mvQuad_1.0-6 patchwork_1.1.2 naomi_2.8.5 sf_1.0-9
+ [6] aghq_0.4.1 tmbstan_1.0.4 rstan_2.21.5 StanHeaders_2.21.0-7 TMB_1.9.2
+[11] Matrix_1.5-4.1 stringr_1.5.0 purrr_1.0.1 tidyr_1.3.0 readr_2.1.3
+[16] ggplot2_3.4.0 forcats_0.5.2 dplyr_1.0.10
+
+loaded via a namespace (and not attached):
+ [1] colorspace_2.0-3 deldir_1.0-6 ellipsis_0.3.2 class_7.3-20 eppasm_0.7.1 rstudioapi_0.14
+ [7] proxy_0.4-27 farver_2.1.1 bit64_4.0.5 fansi_1.0.4 mvtnorm_1.1-3 xml2_1.3.3
+ [13] splines_4.2.0 codetools_0.2-18 cachem_1.0.6 knitr_1.41 jsonlite_1.8.4 data.tree_1.0.0
+ [19] httr_1.4.4 compiler_4.2.0 assertthat_0.2.1 fastmap_1.1.0 cli_3.6.1 s2_1.1.1
+ [25] htmltools_0.5.3 prettyunits_1.1.1 tools_4.2.0 gtable_0.3.1 glue_1.6.2 reshape2_1.4.4
+ [31] wk_0.7.0 V8_4.2.2 fastmatch_1.1-3 Rcpp_1.0.10 jquerylib_0.1.4 vctrs_0.6.1
+ [37] spdep_1.2-7 svglite_2.1.0 xfun_0.37 ps_1.7.3 brio_1.1.3 rvest_1.0.2
+ [43] lifecycle_1.0.3 statmod_1.4.36 orderly_1.4.3 ids_1.0.1 zoo_1.8-11 scales_1.2.1
+ [49] vroom_1.6.0 ragg_1.2.2 hms_1.1.2 parallel_4.2.0 inline_0.3.19 yaml_2.3.7
+ [55] curl_5.0.0 gridExtra_2.3 sass_0.4.4 loo_2.5.1 stringi_1.7.8 highr_0.9
+ [61] e1071_1.7-12 boot_1.3-28 pkgbuild_1.3.1 zip_2.2.2 spData_2.2.1 rlang_1.1.0
+ [67] pkgconfig_2.0.3 systemfonts_1.0.4 inf.utils_0.1.0 matrixStats_0.62.0 evaluate_0.20 lattice_0.20-45
+ [73] labeling_0.4.2 cowplot_1.1.1 bit_4.0.5 processx_3.8.0 tidyselect_1.2.0 traduire_0.0.6
+ [79] bookdown_0.26 plyr_1.8.8 magrittr_2.0.3 R6_2.5.1 generics_0.1.3 DBI_1.1.3
+ [85] pillar_1.9.0 withr_2.5.0 units_0.8-0 sp_1.5-1 tibble_3.2.1 crayon_1.5.2
+ [91] uuid_1.1-0 first90_1.5.3 KernSmooth_2.23-20 utf8_1.2.3 tzdb_0.3.0 rmarkdown_2.18
+ [97] grid_4.2.0 isoband_0.2.6 data.table_1.14.6 callr_3.7.3 webshot_0.5.3 digest_0.6.31
+[103] classInt_0.4-8 numDeriv_2016.8-1.1 textshaping_0.3.6 openssl_2.0.5 RcppParallel_5.1.5 stats4_4.2.0
+[109] munsell_0.5.0 kableExtra_1.3.4 viridisLite_0.4.1 bslib_0.4.1 askpass_1.1