Go back to the About page.
Let us set some global options for all code chunks in this document.
# Set seed for reproducibility
set.seed(1982)
# Set global options for all code chunks
knitr::opts_chunk$set(
# Disable messages printed by R code chunks
message = FALSE,
# Disable warnings printed by R code chunks
warning = FALSE,
# Show R code within code chunks in output
echo = TRUE,
# Include both R code and its results in output
include = TRUE,
# Evaluate R code chunks
eval = TRUE,
# Enable caching of R code chunks for faster rendering
cache = FALSE,
# Align figures in the center of the output
fig.align = "center",
# Enable retina display for high-resolution figures
retina = 2,
# Show errors in the output instead of stopping rendering
error = TRUE,
# Do not collapse code and output into a single block
collapse = FALSE
)
# Start the figure counter
fig_count <- 0
# Define the captioner function
captioner <- function(caption) {
fig_count <<- fig_count + 1
paste0("Figure ", fig_count, ": ", caption)
}
For an implementation of the true covariance matrix using the KL expansion, see here.
# Function to build an interval graph and create a mesh
gets_graph_interval <- function(n){
edge <- rbind(c(0,0),c(1,0))
edges = list(edge)
graph <- metric_graph$new(edges = edges)
graph$build_mesh(n = n)
return(graph)
}
# Plot the interval graph
gets_graph_interval(n = 333)$plot() +
ggtitle("Interval graph") +
theme_minimal() +
theme(text = element_text(family = "Palatino"))
Figure 1: Interval graph.
Press the Show button below to reveal the code.
# Matern covariance function
matern.covariance <- function(h, kappa, nu, sigma) {
if (nu == 1 / 2) {
C <- sigma^2 * exp(-kappa * abs(h))
} else {
C <- (sigma^2 / (2^(nu - 1) * gamma(nu))) *
((kappa * abs(h))^nu) * besselK(kappa * abs(h), nu)
}
C[h == 0] <- sigma^2
return(as.matrix(C))
}
# Folded.matern.covariance.1d
folded.matern.covariance.1d.local <- function(x, kappa, nu, sigma,
L = 1, N = 10,
boundary = c("neumann",
"dirichlet", "periodic")) {
boundary <- tolower(boundary[1])
if (!(boundary %in% c("neumann", "dirichlet", "periodic"))) {
stop("The possible boundary conditions are 'neumann',
'dirichlet' or 'periodic'!")
}
addi = t(outer(x, x, "+"))
diff = t(outer(x, x, "-"))
s1 <- sapply(-N:N, function(j) {
diff + 2 * j * L
})
s2 <- sapply(-N:N, function(j) {
addi + 2 * j * L
})
if (boundary == "neumann") {
C <- rowSums(matern.covariance(h = s1, kappa = kappa,
nu = nu, sigma = sigma) +
matern.covariance(h = s2, kappa = kappa,
nu = nu, sigma = sigma))
} else if (boundary == "dirichlet") {
C <- rowSums(matern.covariance(h = s1, kappa = kappa,
nu = nu, sigma = sigma) -
matern.covariance(h = s2, kappa = kappa,
nu = nu, sigma = sigma))
} else {
C <- rowSums(matern.covariance(h = s1,
kappa = kappa, nu = nu, sigma = sigma))
}
return(matrix(C, nrow = length(x)))
}
# Function to get the true covariance matrix
gets_true_cov_mat = function(graph, kappa, nu, sigma, N, boundary){
#h = c(graph$mesh$V[-2,1], graph$mesh$V[2,1])
h = graph$mesh$V[,1]
true_cov_mat = folded.matern.covariance.1d.local(x = h, kappa = kappa, nu = nu, sigma = sigma, N = N, boundary = boundary)
return(true_cov_mat)
}
# Parameters of the approximation
n = 333
type = "covariance"
type_rational_approximation = "chebfun"
rho = 1
m = 4
nu = 0.6
sigma = 1
N.folded = 10
boundary = "neumann"
# Get the graph and build the mesh
graph = gets_graph_interval(n = n)
kappa = sqrt(8*nu)/rho
tau = sqrt(gamma(nu) / (sigma^2 * kappa^(2*nu) * (4*pi)^(1/2) * gamma(nu + 1/2))) #sigma = 1, d = 1
alpha = nu + 1/2
# Get the true covariance
true_cov_mat = gets_true_cov_mat(graph = graph,
kappa = kappa,
nu = nu,
sigma = sigma,
N = N.folded,
boundary = boundary)
# Get the approximate covariance matrix
op = matern.operators(alpha = alpha,
kappa = kappa,
tau = tau,
m = m,
graph = graph,
type = type,
type_rational_approximation = type_rational_approximation)
appr_cov_mat = op$covariance_mesh()
# Plot the true covariance and the approximation
point <- c(1,0.5)
c_cov <- op$cov_function_mesh(matrix(point,1,2))
loc <- graph$coordinates(PtE = point)
m1 <- which.min((graph$mesh$V[,1]-loc[1])^2 + (graph$mesh$V[,2]-loc[2])^2)
p = graph$plot_function(true_cov_mat[m1,], plotly = TRUE, line_width = 3, edge_width = 3) # blue is the true one
graph$plot_function(c_cov, p = p, line_color = "red", plotly = TRUE, line_width = 3, edge_width = 3) %>%
config(mathjax = 'cdn') %>%
layout(font = list(family = "Palatino"),
scene = list(
annotations = list(
list(
x = 0, y = 0, z = 0,
text = TeX("v_1"),
textangle = 0, ax = 0, ay = 15,
font = list(color = "black", size = 16),
arrowcolor = "black", arrowsize = 1, arrowwidth = 0.1, arrowhead = 1),
list(
x = 0, y = 1, z = 0,
text = TeX("v_2"),
textangle = 0, ax = 0, ay = 15,
font = list(color = "black", size = 16),
arrowcolor = "black", arrowsize = 1, arrowwidth = 0.1, arrowhead = 1),
list(
x = 0, y = 0.5, z = 0,
text = TeX("e_1"),
textangle = 0, ax = 0, ay = 15,
font = list(color = "black", size = 16),
arrowcolor = "white", arrowsize = 1, arrowwidth = 0.1, arrowhead = 1),
list(
x = 0, y = 0.6, z = 1.25,
text = TeX("\\bullet\\mbox{ exact}"),
textangle = 0, ax = 60, ay = 0,
font = list(color = "blue", size = 16),
arrowcolor = "white", arrowsize = 1, arrowwidth = 0.1, arrowhead = 1),
list(
x = 0, y = 0.745, z = 1.20,
text = TeX("\\bullet\\mbox{ approximated}"),
textangle = 0, ax = 60, ay = 0,
font = list(color = "red", size = 16),
arrowcolor = "white", arrowsize = 1, arrowwidth = 0.1, arrowhead = 1))))
Figure 2: Exact and approximated covariance function on the interval graph.
For an implementation of the true covariance matrix using the KL expansion, see here.
# Function to build a circle graph and create a mesh
gets_graph_circle <- function(n){
r = 1/(pi)
theta <- seq(from=-pi,to=pi,length.out = 100)
edge <- cbind(1+r+r*cos(theta),r*sin(theta))
edges = list(edge)
graph <- metric_graph$new(edges = edges)
graph$build_mesh(n = n)
return(graph)
}
# Plot the circle graph
gets_graph_circle(n = 666)$plot() +
ggtitle("Circle graph") +
theme_minimal() +
theme(text = element_text(family = "Palatino"))
Figure 3: Circle graph.
Press the Show button below to reveal the code.
# Matern covariance function
matern.covariance <- function(h, kappa, nu, sigma) {
if (nu == 1 / 2) {
C <- sigma^2 * exp(-kappa * abs(h))
} else {
C <- (sigma^2 / (2^(nu - 1) * gamma(nu))) *
((kappa * abs(h))^nu) * besselK(kappa * abs(h), nu)
}
C[h == 0] <- sigma^2
return(as.matrix(C))
}
# Folded.matern.covariance.1d
folded.matern.covariance.1d.local <- function(x, kappa, nu, sigma, L = 1, N = 10, boundary = c("neumann",
"dirichlet", "periodic")) {
boundary <- tolower(boundary[1])
if (!(boundary %in% c("neumann", "dirichlet", "periodic"))) {
stop("The possible boundary conditions are 'neumann',
'dirichlet' or 'periodic'!")
}
addi = t(outer(x, x, "+"))
diff = t(outer(x, x, "-"))
s1 <- sapply(-N:N, function(j) { # s1 is a matrix of size length(h)x(2N+1)
diff + 2 * j * L
})
s2 <- sapply(-N:N, function(j) {
addi + 2 * j * L
})
if (boundary == "neumann") {
C <- rowSums(matern.covariance(h = s1, kappa = kappa,
nu = nu, sigma = sigma) +
matern.covariance(h = s2, kappa = kappa,
nu = nu, sigma = sigma))
} else if (boundary == "dirichlet") {
C <- rowSums(matern.covariance(h = s1, kappa = kappa,
nu = nu, sigma = sigma) -
matern.covariance(h = s2, kappa = kappa,
nu = nu, sigma = sigma))
} else {
C <- rowSums(matern.covariance(h = s1,
kappa = kappa, nu = nu, sigma = sigma))
}
return(matrix(C, nrow = length(x)))
}
# Function to get the true covariance matrix
gets_true_cov_mat = function(graph, kappa, nu, sigma, N, boundary){
h = c(0,graph$get_edge_lengths()[1]*graph$mesh$PtE[,2])
true_cov_mat = folded.matern.covariance.1d.local(x = h, kappa = kappa, nu = nu, sigma = sigma, N = N, boundary = boundary)
return(true_cov_mat)
}
# Parameters of the approximation
n = 666
type = "covariance"
type_rational_approximation = "chebfun"
rho = 1
m = 4
nu = 0.6
sigma = 1
N.folded = 10
boundary = "periodic"
# Get the graph and build the mesh
graph = gets_graph_circle(n = n)
kappa = sqrt(8*nu)/rho
tau = sqrt(gamma(nu) / (sigma^2 * kappa^(2*nu) * (4*pi)^(1/2) * gamma(nu + 1/2))) #sigma = 1, d = 1
alpha = nu + 1/2
# Get the "true" covariance matrix
true_cov_mat = gets_true_cov_mat(graph = graph,
kappa = kappa,
nu = nu,
sigma = sigma,
N = N.folded,
boundary = boundary)
# Get the approximate covariance matrix
op = matern.operators(alpha = alpha,
kappa = kappa,
tau = tau,
m = m,
graph = graph,
type = type,
type_rational_approximation = type_rational_approximation)
appr_cov_mat = op$covariance_mesh()
# Plot the true covariance and the approximation
point <- c(1,0)
c_cov <- op$cov_function_mesh(matrix(point,1,2))
loc <- graph$coordinates(PtE = point)
m1 <- which.min((graph$mesh$V[,1]-loc[1])^2 + (graph$mesh$V[,2]-loc[2])^2)
p = graph$plot_function(true_cov_mat[m1,], plotly = TRUE, line_width = 3, edge_width = 3)
graph$plot_function(c_cov,p = p, line_color = "red", plotly = TRUE,line_width = 3, edge_width = 3) %>%
config(mathjax = 'cdn') %>%
layout(font = list(family = "Palatino"),
scene = list(
annotations = list(
list(
x = -0.03, y = 1.03, z = 0,
text = TeX("v_1"),
textangle = 0, ax = 0, ay = 15,
font = list(color = "black", size = 16),
arrowcolor = "white", arrowsize = 1, arrowwidth = 0.1, arrowhead = 1),
list(
x = 0.26, y = 1.54, z = 0,
text = TeX("e_1"),
textangle = 0, ax = 0, ay = 15,
font = list(color = "black", size = 16),
arrowcolor = "white", arrowsize = 1, arrowwidth = 0.1, arrowhead = 1),
list(
x = -0.05, y = 1.085, z = 0.95,
text = TeX("\\bullet\\mbox{ exact}"),
textangle = 0, ax = 60, ay = 0,
font = list(color = "blue", size = 16),
arrowcolor = "white", arrowsize = 1, arrowwidth = 0.1, arrowhead = 1),
list(
x = -0.05, y = 1.2, z = 0.9,
text = TeX("\\bullet\\mbox{ approximated}"),
textangle = 0, ax = 60, ay = 0,
font = list(color = "red", size = 16),
arrowcolor = "white", arrowsize = 1, arrowwidth = 0.1, arrowhead = 1))))
Figure 4: Exact and approximated covariance function on the circle graph.
# Function to build a tadpole graph and create a mesh
gets_graph_tadpole <- function(h){
edge1 <- rbind(c(0,0),c(1,0))
theta <- seq(from=-pi,to=pi,length.out = 100)
edge2 <- cbind(1+1/pi+cos(theta)/pi,sin(theta)/pi)
edges = list(edge1, edge2)
graph <- metric_graph$new(edges = edges)
graph$build_mesh(h = h)
return(graph)
}
# Plot the tadpole graph
gets_graph_tadpole(h = 0.003)$plot() +
ggtitle("Tadpole graph") +
theme_minimal() +
theme(text = element_text(family = "Palatino"))
Figure 5: Tadpole graph.
Press the Show button below to reveal the code.
# Function to compute the eigenfunctions
tadpole.eig <- function(k,graph){
x1 <- c(0,graph$get_edge_lengths()[1]*graph$mesh$PtE[graph$mesh$PtE[,1]==1,2])
x2 <- c(0,graph$get_edge_lengths()[2]*graph$mesh$PtE[graph$mesh$PtE[,1]==2,2])
if(k==0){
f.e1 <- rep(1,length(x1))
f.e2 <- rep(1,length(x2))
f1 = c(f.e1[1],f.e2[1],f.e1[-1], f.e2[-1])
f = list(phi=f1/sqrt(3))
} else {
f.e1 <- -2*sin(pi*k*1/2)*cos(pi*k*x1/2)
f.e2 <- sin(pi*k*x2/2)
f1 = c(f.e1[1],f.e2[1],f.e1[-1], f.e2[-1])
if((k %% 2)==1){
f = list(phi=f1/sqrt(3))
} else {
f.e1 <- (-1)^{k/2}*cos(pi*k*x1/2)
f.e2 <- cos(pi*k*x2/2)
f2 = c(f.e1[1],f.e2[1],f.e1[-1],f.e2[-1])
f <- list(phi=f1,psi=f2/sqrt(3/2))
}
}
return(f)
}
# Function to get the true covariance matrix
gets_true_cov_mat <- function(graph, kappa, tau, alpha, n.overkill){
Sigma.kl <- matrix(0,nrow = dim(graph$mesh$V)[1],ncol = dim(graph$mesh$V)[1])
for(i in 0:n.overkill){
phi <- tadpole.eig(i,graph)$phi
Sigma.kl <- Sigma.kl + (1/(kappa^2 + (i*pi/2)^2)^(alpha))*phi%*%t(phi)
if(i>0 && (i %% 2)==0){
psi <- tadpole.eig(i,graph)$psi
Sigma.kl <- Sigma.kl + (1/(kappa^2 + (i*pi/2)^2)^(alpha))*psi%*%t(psi)
}
}
Sigma.kl <- Sigma.kl/tau^2
return(Sigma.kl)
}
# Parameters of the approximation
h = 0.003
type = "covariance"
type_rational_approximation = "chebfun"
rho = 1
m = 4
nu = 0.6
# Get the graph and build the mesh
graph = gets_graph_tadpole(h = h)
n.overkill = 10000
kappa = sqrt(8*nu)/rho
tau = sqrt(gamma(nu) / (1^2 * kappa^(2*nu) * (4*pi)^(1/2) * gamma(nu + 1/2))) #sigma = 1, d = 1
alpha = nu + 1/2
# Get the true covariance
true_cov_mat = gets_true_cov_mat(graph = graph,
kappa = kappa,
tau = tau,
alpha = alpha,
n.overkill = n.overkill)
# Get the approximate covariance matrix
op = matern.operators(alpha = alpha,
kappa = kappa,
tau = tau,
m = m,
graph = graph,
type = type,
type_rational_approximation = type_rational_approximation)
appr_cov_mat = op$covariance_mesh()
# Plot the true covariance and the approximation
point <- c(1,1)
c_cov <- op$cov_function_mesh(matrix(point,1,2))
loc <- graph$coordinates(PtE = point)
m1 <- which.min((graph$mesh$V[,1]-loc[1])^2 + (graph$mesh$V[,2]-loc[2])^2)
p <- graph$plot_function(true_cov_mat[m1,], plotly = TRUE, line_width = 3, edge_width = 3) # blue is the true one
graph$plot_function(c_cov, p = p, line_color = "red", plotly = TRUE,line_width = 3, edge_width = 3) %>%
config(mathjax = 'cdn') %>%
layout(font = list(family = "Palatino"),
scene = list(
annotations = list(
list(
x = 0, y = 0, z = 0,
text = TeX("v_1"),
textangle = 0, ax = 0, ay = 15,
font = list(color = "black", size = 16),
arrowcolor = "black", arrowsize = 1, arrowwidth = 0.1, arrowhead = 1),
list(
x = -0.03, y = 1.03, z = 0,
text = TeX("v_2"),
textangle = 0, ax = 0, ay = 15,
font = list(color = "black", size = 16),
arrowcolor = "white", arrowsize = 1, arrowwidth = 0.1, arrowhead = 1),
list(
x = 0, y = 0.5, z = 0,
text = TeX("e_1"),
textangle = 0, ax = 0, ay = 15,
font = list(color = "black", size = 16),
arrowcolor = "white", arrowsize = 1, arrowwidth = 0.1, arrowhead = 0),
list(
x = 0.26, y = 1.54, z = 0,
text = TeX("e_2"),
textangle = 0, ax = 0, ay = 15,
font = list(color = "black", size = 16),
arrowcolor = "white", arrowsize = 1, arrowwidth = 0.1, arrowhead = 1),
list(
x = -0.08, y = 1.135, z = 0.65,
text = TeX("\\bullet\\mbox{ exact}"),
textangle = 0, ax = 60, ay = 0,
font = list(color = "blue", size = 16),
arrowcolor = "white", arrowsize = 1, arrowwidth = 0.1, arrowhead = 1),
list(
x = -0.08, y = 1.25, z = 0.6,
text = TeX("\\bullet\\mbox{ approximated}"),
textangle = 0, ax = 60, ay = 0,
font = list(color = "red", size = 16),
arrowcolor = "white", arrowsize = 1, arrowwidth = 0.1, arrowhead = 1))))
Figure 6: Exact and approximated covariance function on the tadpole graph.
We used R version 4.4.1 (R Core Team 2024a) and the following R packages: cowplot v. 1.1.3 (Wilke 2024), ggmap v. 4.0.0.900 (Kahle and Wickham 2013), ggpubr v. 0.6.0 (Kassambara 2023), ggtext v. 0.1.2 (Wilke and Wiernik 2022), grid v. 4.4.1 (R Core Team 2024b), here v. 1.0.1 (Müller 2020), htmltools v. 0.5.8.1 (Cheng et al. 2024), INLA v. 24.12.11 (Rue, Martino, and Chopin 2009; Lindgren, Rue, and Lindström 2011; Martins et al. 2013; Lindgren and Rue 2015; De Coninck et al. 2016; Rue et al. 2017; Verbosio et al. 2017; Bakka et al. 2018; Kourounis, Fuchs, and Schenk 2018), inlabru v. 2.12.0.9002 (Yuan et al. 2017; Bachl et al. 2019), knitr v. 1.48 (Xie 2014, 2015, 2024), latex2exp v. 0.9.6 (Meschiari 2022), Matrix v. 1.6.5 (Bates, Maechler, and Jagan 2024), MetricGraph v. 1.4.0.9000 (Bolin, Simas, and Wallin 2023b, 2023a, 2023c, 2024; Bolin et al. 2024), OpenStreetMap v. 0.4.0 (Fellows and JMapViewer library by Jan Peter Stotz 2023), osmdata v. 0.2.5 (Mark Padgham et al. 2017), patchwork v. 1.2.0 (Pedersen 2024), plotly v. 4.10.4 (Sievert 2020), plotrix v. 3.8.4 (J 2006), reshape2 v. 1.4.4 (Wickham 2007), rmarkdown v. 2.28 (Xie, Allaire, and Grolemund 2018; Xie, Dervieux, and Riederer 2020; Allaire et al. 2024), rSPDE v. 2.4.0.9000 (Bolin and Kirchner 2020; Bolin and Simas 2023; Bolin, Simas, and Xiong 2024), scales v. 1.3.0 (Wickham, Pedersen, and Seidel 2023), sf v. 1.0.19 (E. Pebesma 2018; E. Pebesma and Bivand 2023), sp v. 2.1.4 (E. J. Pebesma and Bivand 2005; Bivand, Pebesma, and Gomez-Rubio 2013), tidyverse v. 2.0.0 (Wickham et al. 2019), viridis v. 0.6.4 (Garnier et al. 2023), xaringanExtra v. 0.8.0 (Aden-Buie and Warkentin 2024).