---
title: "Sampling from Gaussian CAR and ICAR Models"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Sampling from Gaussian CAR and ICAR Models}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  message = FALSE,
  warning = FALSE,
  fig.width = 6,
  fig.height = 4,
  dpi = 72,
  fig.retina = 1
)
```

## Purpose

This vignette demonstrates how to sample Gaussian random effects defined on a graph using **trafficCAR**.
The emphasis is on models specified through a **sparse precision matrix**:

\[
x \sim \mathrm{N}(Q^{-1} b,\; Q^{-1}),
\]

where \(Q\) is a sparse precision matrix and \(b\) is an optional linear term.
Both **proper CAR** and **intrinsic CAR (ICAR)** specifications are illustrated.

```{r}
library(trafficCAR)
library(Matrix)
```

## CAR and ICAR precision matrices

Let \(A\) be a symmetric adjacency matrix with zero diagonal, and let \(D = \mathrm{diag}(A\mathbf{1})\).
A common CAR family is

\[
Q = \tau (D - \rho A),
\]

with \(\tau > 0\). For a **proper CAR**, choose \(\rho\) so that \(D - \rho A\) is positive definite.
For an **ICAR**, the conventional choice is \(\rho = 1\), which yields a singular precision matrix.

The function `car_precision()` constructs \(Q\) from \(A\).

## A small graph example

Construct a simple path graph with 3 nodes.

```{r}
A_small <- matrix(
  c(0, 1, 0,
    1, 0, 1,
    0, 1, 0),
  nrow = 3,
  byrow = TRUE
)
A_small
```

### Proper CAR

```{r}
Q_car <- car_precision(A_small, type = "proper", rho = 0.5, tau = 1)
Q_car
```

### ICAR

```{r}
Q_icar <- car_precision(A_small, type = "icar", tau = 1)
Q_icar
```

## Sampling from a sparse precision matrix

Sampling is performed by `rmvnorm_prec()`, which draws from \(\mathcal{N}(Q^{-1} b, Q^{-1})\) using sparse factorizations.
A dense covariance matrix is not formed.

### Zero-mean sampling

```{r}
x0 <- trafficCAR:::rmvnorm_prec(Q_car)
x0
```

### Nonzero mean (linear term)

```{r}
b <- c(1, 0, -1)
x1 <- trafficCAR:::rmvnorm_prec(Q_car, b = b)
x1
```

## Notes on ICAR sampling

For ICAR models, \(Q\) is singular and the implied Gaussian distribution is improper.
Samples are identifiable only up to an additive constant. A common post-processing step is to center a draw.

```{r}
# ICAR precision matrices are singular, so direct Cholesky-based sampling may fail.
# A practical approach for simulation is to add a small ridge and then center.
# This yields an approximate draw on the ICAR subspace.

n_icar <- nrow(Q_icar)
eps <- 1e-6
Q_icar_ridge <- Q_icar + Matrix::Diagonal(n_icar, x = eps)

x_icar <- trafficCAR:::rmvnorm_prec(Q_icar_ridge)
x_icar_centered <- x_icar - mean(x_icar)
x_icar_centered
```

## A road-like graph example

This section illustrates CAR/ICAR sampling on a graph resembling a small road network.
If `sf` is available, we build an adjacency matrix from the bundled `roads_small` dataset.
Otherwise, we fall back to a small synthetic graph.

```{r}
A_road <- NULL

has_sf <- requireNamespace("sf", quietly = TRUE)

if (has_sf) {
  data("roads_small", package = "trafficCAR")
  roads <- roads_small

  segments <- roads_to_segments(
    roads,
    crs_m = 3857,
    split_at_intersections = TRUE
  )

  if (nrow(segments) > 200) {
    segments <- segments[seq_len(200), ]
  }

  adjacency <- build_adjacency(segments, crs_m = 3857)

  # Drop isolated segments to keep CAR models well-defined.
  if (any(adjacency$isolates)) {
    segments <- segments[!adjacency$isolates, ]
    adjacency <- build_adjacency(segments, crs_m = 3857)
  }

  A_road <- adjacency$A
}
```

### Option 2: Synthetic road-like graph (fallback)

The following graph mimics a small road network with a T-junction.

```{r}
if (is.null(A_road)) {
  # nodes: 1--2--3 and 2--4 (T-junction at node 2)
  edges <- matrix(c(
    1, 2,
    2, 3,
    2, 4
  ), ncol = 2, byrow = TRUE)

  if (requireNamespace("igraph", quietly = TRUE)) {
    g <- igraph::graph_from_edgelist(edges, directed = FALSE)
    A_road <- igraph::as_adjacency_matrix(g, sparse = TRUE, attr = NULL)
  } else {
    # minimal fallback without igraph (dense, small only)
    n <- max(edges)
    A_road <- matrix(0, n, n)
    for (k in seq_len(nrow(edges))) {
      i <- edges[k, 1]; j <- edges[k, 2]
      A_road[i, j] <- 1
      A_road[j, i] <- 1
    }
    A_road <- Matrix(A_road, sparse = TRUE)
  }
}

c(n = nrow(A_road), nnz = Matrix::nnzero(A_road))
```

### CAR/ICAR sampling on the road graph

```{r}
# proper CAR (users should choose rho consistent with their graph and model)
Q_car_road <- car_precision(A_road, type = "proper", rho = 0.5, tau = 1)

# ICAR
Q_icar_road <- car_precision(A_road, type = "icar", tau = 1)

x_car_road <- trafficCAR:::rmvnorm_prec(Q_car_road)

n_road <- nrow(Q_icar_road)
eps <- 1e-6
Q_icar_road_ridge <- Q_icar_road + Matrix::Diagonal(n_road, x = eps)
x_icar_road <- trafficCAR:::rmvnorm_prec(Q_icar_road_ridge)
x_icar_road <- x_icar_road - mean(x_icar_road)

summary(x_car_road)
summary(x_icar_road)
```

## Summary

- Proper CAR models yield positive definite precision matrices and support direct Gaussian sampling.
- ICAR models yield singular precisions; sampled fields are defined up to an additive constant and are typically centered or constrained.
- Sampling from sparse precision matrices avoids dense covariance calculations and scales well to larger graphs.