Performing the check

type <- "random"
L <- 200
grid <- readRDS("depends/grid.rds")
n <- nrow(grid)
samples <- sf::st_sample(grid, type = type, exact = TRUE, size = rep(L, n))
length(samples) # n * L
## [1] 7200
sf::st_crs(samples) <- NA

Add an id column:

samples <- st_sf(id = rep(seq_along(1:L), n), geom = samples)

sample_sets <- lapply(seq(10, L, by = 10), function(i) subset(samples, id <= i))

Check that the samples look OK:

cowplot::plot_grid(
  plot_samples(sample_sets[[1]]),
  plot_samples(sample_sets[[10]]),
  plot_samples(sample_sets[[20]]),
  ncol = 3
)

Calculate the Gram matrices:

covs <- lapply(sample_sets, adapted_integrated_covariance)

cowplot::plot_grid(
  plot_matrix(covs[[1]]),
  plot_matrix(covs[[10]]),
  plot_matrix(covs[[20]]),
  ncol = 3
)

matrix_comparison <- outer(covs, covs, Vectorize(frobenius))
plot_matrix(matrix_comparison)

Perform an experiment with 10, 20 and 50 points versus ground truth 300 points:

L_settings <- list(10, 20, 50)
gt_samples <- sf::st_sample(grid, type = type, exact = TRUE, size = rep(300, n))
gt_cov <- adapted_integrated_covariance(gt_samples)
result <- lapply(L_settings, int_convergence_experiment)

names(result) <- c("L_10", "L_20", "L_50")

result <- list_to_df(result) %>%
  pivot_longer(
    cols = 1:3,
    names_to = "L",
    names_prefix= "L_",
    values_to = "value"
  ) %>%
  mutate(L = as.numeric(L))

ggplot(result, aes(x = value)) +
  geom_histogram() +
  facet_grid(~L) +
  labs(x = "Frobenius norm", y = "") +
  theme_minimal()
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

functions.R

adapted_integrated_covariance <- function(samples) {
  sample_index <- sf::st_intersects(grid, samples)
  D <- sf::st_distance(samples, samples)
  kD <- matern(D, l = 1)
  K <- matrix(nrow = n, ncol = n)
  #' Diagonal entries
  for(i in 1:(n - 1)) {
    K[i, i] <- mean(kD[sample_index[[i]], sample_index[[i]]])
    for(j in (i + 1):n) {
      #' Off-diagonal entries
      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
    }
  }
  K[n, n] <- mean(kD[sample_index[[n]], sample_index[[n]]])
  return(K)
}

plot_samples <- function(samples){
  ggplot(grid) +
    geom_sf(fill = "lightgrey") +
    geom_sf(data = samples, alpha = 0.5, shape = 4) +
    labs(x = "Longitude", y = "Latitude") +
    theme_minimal() +
    labs(subtitle = paste0(length(samples$geom), " samples"), fill = "") +
    theme(axis.title = element_blank(),
          axis.text = element_blank(),
          axis.ticks = element_blank())
}

#' Simple metric between matrices of identical dimension
frobenius <- function(M1, M2) {
  diff <- M1 - M2
  sqrt(sum(diff^2))
}

int_convergence_experiment <- function(L, n_sim = 100) {
  sample_sets <- lapply(1:n_sim, function(i) sf::st_sample(grid, type = type, exact = TRUE, size = rep(L, n)))
  covs <- lapply(sample_sets, adapted_integrated_covariance)
  frobs <- sapply(covs, function(cov) frobenius(gt_cov, cov))
  return(frobs)
}

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