## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)


## ----eval=TRUE,  warning=FALSE, message=FALSE---------------------------------
library(BMEmapping)
library(ggplot2)
library(sf)
library(gstat)
library(dplyr)
library(tidyr)
library(scales)
library(knitr)
library(gridExtra)

data("utsnowload")


## ----eval=TRUE----------------------------------------------------------------
# Load the sample spatial dataset
data("utsnowload")

# Display baseline distributions and structural summaries
summary(utsnowload)


## ----eval=FALSE---------------------------------------------------------------
# ?utsnowload
# 

## ----eval=TRUE----------------------------------------------------------------
# hard data locations
ch <- utsnowload[1:30, c("longitude", "latitude")]

# hard data values
zh <- utsnowload[1:30, c("hard")]

# soft data locations
cs <- utsnowload[68:167, c("longitude", "latitude")]

# lower and upper bounds of soft data (intervals)
a <- utsnowload[68:167, c("lower")]
b <- utsnowload[68:167, c("upper")]


## ----eval=TRUE----------------------------------------------------------------
data_object <- bme_map(ch, cs, zh, a, b)


## ----eval=TRUE, fig.width = 6, fig.height = 6, fig.align='center'-------------
plot(data_object)


## ----eval=TRUE----------------------------------------------------------------
xk <- utsnowload[201:205, c("longitude", "latitude")]
xk


## ----eval=TRUE----------------------------------------------------------------
df <- data.frame(rbind(ch, cs), z = c(zh, (a + b) / 2))
sf_data <- sf::st_as_sf(df, coords = c("longitude", "latitude"))

vg <- variogram(z ~ 1, data = sf_data)
vg_model <- fit.variogram(vg, model = vgm(c("Exp", "Sph")))
vg_model


## ----eval=TRUE, fig.width = 6, fig.height = 5, fig.align='center'-------------
# Extract spatial parameters
model  <- as.character(vg_model[2, 1])
nugget <- vg_model[1, 2]
sill   <- vg_model[2, 2]
range  <- vg_model[2, 3]

# default zk_range
p_1 <- prob_zk(xk[1,], data_object, model, nugget, sill, range)
p_2 <- prob_zk(xk[2,], data_object, model, nugget, sill, range)
p_3 <- prob_zk(xk[3,], data_object, model, nugget, sill, range)
p_4 <- prob_zk(xk[4,], data_object, model, nugget, sill, range)
p_df <- cbind.data.frame(p_1, p_2[, 2], p_3[, 2], p_4[, 2])
names(p_df) <- c("zk_i", "p1", "p2", "p3", "p4")

# Function to generate ggplot for a given column name
plot_prob_curve <- function(df, pi_col) {
  ggplot(df, aes(x = zk_i, y = .data[[pi_col]])) +
    geom_line(color = "darkblue", linewidth = 0.7) +
    labs(x = "Spatial Field Value (z)", y = "Posterior Density f(z)") +
    theme_minimal() +
    theme(
      panel.background = ggplot2::element_rect(fill = "white", color = "black")
    )
}

# Generate individual plots
p1 <- plot_prob_curve(p_df, "p1")
p2 <- plot_prob_curve(p_df, "p2")
p3 <- plot_prob_curve(p_df, "p3")
p4 <- plot_prob_curve(p_df, "p4")

# Arrange in a 2x2 grid
grid.arrange(p1, p2, p3, p4, ncol = 2)


## ----eval=TRUE, fig.width = 6, fig.height = 5, fig.align='center'-------------
# updated zk_range: [-2, 2]
q_1 <- prob_zk(xk[1,], data_object, model, nugget, sill, range, zk_range = c(-2, 2))
q_2 <- prob_zk(xk[2,], data_object, model, nugget, sill, range, zk_range = c(-2, 2))
q_3 <- prob_zk(xk[3,], data_object, model, nugget, sill, range, zk_range = c(-2, 2))
q_4 <- prob_zk(xk[4,], data_object, model, nugget, sill, range, zk_range = c(-2, 2))
q_df <- cbind.data.frame(q_1, q_2[, 2], q_3[, 2], q_4[, 2])
names(q_df) <- c("zk_i", "q1", "q2", "q3", "q4")

# Generate individual plots
q1 <- plot_prob_curve(q_df, "q1")
q2 <- plot_prob_curve(q_df, "q2")
q3 <- plot_prob_curve(q_df, "q3")
q4 <- plot_prob_curve(q_df, "q4")

grid.arrange(q1, q2, q3, q4, ncol = 2)


## ----eval=TRUE----------------------------------------------------------------
CBME_mode <- bme_predict(xk, data_object, model, nugget, sill, range,
                         n = 100, zk_range = c(-2, 2), type = "mode")
head(CBME_mode)


## ----eval=TRUE----------------------------------------------------------------
CBME_mean <- bme_predict(xk, data_object, model, nugget, sill, range,
                         n = 100, zk_range = c(-2, 2), type = "mean")
head(CBME_mean)


## ----eval=TRUE----------------------------------------------------------------
CBME_median <- bme_predict(xk, data_object, model, nugget, sill, range,
                           n = 100, zk_range = c(-2, 2), type = "median")
head(CBME_median)


## ----eval=TRUE----------------------------------------------------------------
CBME_interval <- bme_predict_ci(xk, data_object, model, nugget, sill, range,
                                n = 100, zk_range = c(-2, 2), level = 0.90)
head(CBME_interval)


## ----eval=FALSE---------------------------------------------------------------
# # Evaluate localized posterior density curve using nq = 8 quantile slices
# q_1 <- q_prob_zk(xk[1,], data_object, nq = 8)
# 

## ----eval=FALSE---------------------------------------------------------------
# # Compute posterior mode point predictions
# QBME_mode   <- q_bme_predict(xk, data_object, n = 100, type = "mode", nq = 8)
# 
# # Compute posterior mean point predictions
# QBME_mean   <- q_bme_predict(xk, data_object, n = 100, type = "mean", nq = 8)
# 
# # Compute posterior median point predictions
# QBME_median <- q_bme_predict(xk, data_object, n = 100, type = "median", nq = 8)
# 
# # Integrated 90% probabilistic credible interval bounds
# QBME_int    <- q_bme_predict_ci(xk, data_object, n = 100, level = 0.90, nq = 8)
# 

## ----eval=TRUE----------------------------------------------------------------
# Execute 5-fold cross-validation across the hard data locations
CBME_cv <- bme_cv(data_object, model, nugget, sill, range, n = 100, 
                  zk_range = c(-2, 2), type = "mean", k = 5)
CBME_cv


## ----eval=TRUE----------------------------------------------------------------
# Extract comprehensive cross-validation error metrics
summary(CBME_cv)


## ----eval=TRUE, fig.width = 6, fig.height = 5, fig.align='center'-------------
# Render diagnostic residual charts
plot(CBME_cv)


## ----eval=FALSE---------------------------------------------------------------
# # Execute Leave-One-Out Cross-Validation (LOOCV) under the QBME framework
# QBME_cv <- q_bme_cv(data_object, n = 100, nq = 8, type = "mean", k = nrow(ch))
# QBME_cv
# 

## ----eval=FALSE---------------------------------------------------------------
# # Extract diagnostic summary metrics for the QBME validation routine
# summary(QBME_cv)
# 

## ----eval=FALSE---------------------------------------------------------------
# # Render diagnostic residual plots to assess QBME prediction adequacy
# plot(QBME_cv)
# 

## ----eval=TRUE, fig.width = 6, fig.height = 5, fig.align='center'-------------
# Map the continuous spatial prediction surface across the target geographic area
plot(CBME_mean)