Check integrated kernel computation
Adam Howes
Abstract
Task We check that the Stan code to calculate the integrated kernel is correct.
Performing the check
sf <- mw
L <- 10
type = "random"
n <- nrow(sf)
samples <- sf::st_sample(sf, type = type, exact = TRUE, size = rep(L, n))
S <- sf::st_distance(samples, samples)
sample_index <- sf::st_intersects(sf, samples)
sample_lengths <- lengths(sample_index)
start_index <- sapply(sample_index, function(x) x[1])First, calculate the integrated kernel directly using R only:
D <- sf::st_distance(samples, samples)
kD <- matern(D, l = 1)
K <- matrix(nrow = n, ncol = n)
for(i in 1:n) {
K[i, i] <- mean(kD[sample_index[[i]], sample_index[[i]]])
}
for(i in 1:(n - 1)) {
for(j in (i + 1):n) {
K[i, j] <- mean(kD[sample_index[[i]], sample_index[[j]]]) #' Fill the upper triangle
K[j, i] <- K[i, j] #' Fill the lower triangle
}
}
cov_r <- KNow, calculate it in Stan via two different methods. We use
rstan::expose_stan_functions to put the C++ functions into
the R local workspace to use:
rstan::expose_stan_functions("cov_sample_average.stan")Newer method (using cov_sample_average):
cov_stan <- cov_sample_average(
S = S,
l = 1,
n = n,
start_index = start_index,
sample_lengths = sample_lengths,
total_samples = sum(sample_lengths)
)Old method (using cov_sample_average_old):
cov_stan_old <- cov_sample_average_old(
n = n,
L = 10,
l = 1,
S = S
)Comparison (these statements should all be TRUE):
all.equal(cov_r, cov_stan)## [1] TRUEall.equal(cov_r, cov_stan_old)## [1] TRUEall.equal(cov_stan, cov_stan_old)## [1] TRUEAnd here are the images (it’s hard to tell from an image if something is identical, but anyway):
cowplot::plot_grid(
plot_matrix(cov_r) + labs(title = "R"),
plot_matrix(cov_stan_old) + labs(title = "Stan old"),
plot_matrix(cov_stan) + labs(title = "Stan new"),
ncol = 2
)
cov_sample_average.stan
// cov_sample_average.stan
functions {
matrix cov_matern32(matrix D, real l) {
int n = rows(D);
matrix[n, n] K;
real norm_K;
real sqrt3;
sqrt3 = sqrt(3.0);
for(i in 1:(n - 1)){
// Diagonal entries (apart from the last)
K[i, i] = 1;
for(j in (i + 1):n){
// Off-diagonal entries
norm_K = D[i, j] / l;
K[i, j] = (1 + sqrt3 * norm_K) * exp(-sqrt3 * norm_K); // Fill lower triangle
K[j, i] = K[i, j]; // Fill upper triangle
}
}
K[n, n] = 1;
return(K);
}
matrix cov_sample_average(matrix S, real l, int n, int[] start_index, int[] sample_lengths, int total_samples) {
matrix[n, n] K;
matrix[total_samples, total_samples] kS = cov_matern32(S, l);
for(i in 1:(n - 1)) {
// Diagonal entries (apart from the last)
int start_i = start_index[i];
int length_i = sample_lengths[i];
K[i, i] = mean(block(kS, start_i, start_i, length_i, length_i));
for(j in (i + 1):n) {
// Off-diagonal entries
K[i, j] = mean(block(kS, start_i, start_index[j], length_i, sample_lengths[j]));
K[j, i] = K[i, j];
}
}
K[n, n] = mean(block(kS, start_index[n], start_index[n], sample_lengths[n], sample_lengths[n]));
return(K);
}
matrix cov_sample_average_old(int n, int L, real l, matrix S) {
matrix[n, n] K;
matrix[L * n, L * n] kS = cov_matern32(S, l);
// Diagonal entries
for(i in 1:n){
K[i, i] = mean(kS[(L * (i - 1) + 1):(i * L), (L * (i - 1) + 1):(i * L)]);
}
// Off-diagonal entries
for(i in 1:(n - 1)) {
for(j in (i + 1):n) {
K[i, j] = mean(kS[(L * (i - 1) + 1):(i * L), (L * (j - 1) + 1):(j * L)]);
K[j, i] = K[i, j];
}
}
return(K);
}
}Original computing environment
sessionInfo()## R version 4.2.0 (2022-04-22)
## Platform: x86_64-apple-darwin17.0 (64-bit)
## Running under: macOS Monterey 12.6
##
## Matrix products: default
## LAPACK: /Library/Frameworks/R.framework/Versions/4.2/Resources/lib/libRlapack.dylib
##
## locale:
## [1] en_GB.UTF-8/en_GB.UTF-8/en_GB.UTF-8/C/en_GB.UTF-8/en_GB.UTF-8
##
## attached base packages:
## [1] stats4 stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] rmarkdown_2.18 bayesplot_1.9.0 rstan_2.21.5 StanHeaders_2.21.0-7 multi.utils_0.1.0
## [6] purrr_0.3.5 sn_2.0.2 tidyr_1.2.1 gt_0.8.0 magrittr_2.0.3
## [11] cubature_2.0.4.4 tikzDevice_0.12.3.1 viridis_0.6.2 viridisLite_0.4.1 reshape2_1.4.4
## [16] dplyr_1.0.10 ggplot2_3.4.0 bsae_0.2.7 sf_1.0-9
##
## loaded via a namespace (and not attached):
## [1] matrixStats_0.62.0 numDeriv_2016.8-1.1 tools_4.2.0 bslib_0.4.1 utf8_1.2.2
## [6] R6_2.5.1 KernSmooth_2.23-20 DBI_1.1.3 colorspace_2.0-3 withr_2.5.0
## [11] sp_1.5-1 tidyselect_1.2.0 gridExtra_2.3 prettyunits_1.1.1 mnormt_2.0.2
## [16] processx_3.5.3 compiler_4.2.0 cli_3.4.1 labeling_0.4.2 sass_0.4.4
## [21] scales_1.2.1 classInt_0.4-8 ggridges_0.5.3 callr_3.7.0 proxy_0.4-27
## [26] askpass_1.1 stringr_1.5.0 digest_0.6.30 pkgconfig_2.0.3 htmltools_0.5.3
## [31] fastmap_1.1.0 highr_0.9 rlang_1.0.6 rstudioapi_0.13 jquerylib_0.1.4
## [36] generics_0.1.3 farver_2.1.1 jsonlite_1.8.4 inline_0.3.19 loo_2.5.1
## [41] Matrix_1.5-1 Rcpp_1.0.9 munsell_0.5.0 fansi_1.0.3 lifecycle_1.0.3
## [46] stringi_1.7.8 yaml_2.3.6 pkgbuild_1.3.1 plyr_1.8.8 grid_4.2.0
## [51] parallel_4.2.0 crayon_1.5.2 lattice_0.20-45 cowplot_1.1.1 INLA_22.05.07
## [56] splines_4.2.0 BH_1.78.0-0 tmvnsim_1.0-2 knitr_1.41 ps_1.7.0
## [61] pillar_1.8.1 uuid_1.1-0 codetools_0.2-18 glue_1.6.2 evaluate_0.18
## [66] RcppParallel_5.1.5 vctrs_0.5.1 ids_1.0.1 MatrixModels_0.5-0 gtable_0.3.1
## [71] openssl_2.0.5 assertthat_0.2.1 cachem_1.0.6 xfun_0.35 RcppEigen_0.3.3.9.3
## [76] e1071_1.7-12 filehash_2.4-3 class_7.3-20 tibble_3.1.8 orderly_1.4.3
## [81] units_0.8-0 ellipsis_0.3.2