Tutorial: Circular Distributions for Regression Models

This tutorial develops the link-adjusted von Mises (LAvM) regression model in Stan. It begins with the ambiguity of ordinary circular regression, introduces the LAvM data model, and then moves from fixed effects to latent temporal structure and a wind-direction application. The presentation assumes familiarity with a basic Stan workflow and focuses on the circular model itself.

After completing the tutorial, you will be able to:

• formulate an LAvM regression model through a real-valued linear predictor;
• implement fixed-effect and latent-process versions of the model in Stan; and
• interpret regression effects and concentration on the appropriate scales.

1 Preliminary

The model follows Ye, Van Niekerk, and Rue (2026), A Bayesian regression framework for circular models with INLA.

1.1 Von Mises distribution

The von Mises (vM) distribution is a standard model for a unimodal circular variable. Its probability density function is

\[ p_{\mathcal{VM}}\left(x \mid \mu, \kappa\right) = \frac{1}{2\pi \mathcal{I}_{0} \left(\kappa\right)} \exp \left\{ \kappa \cos \left( x - \mu \right) \right\}, \quad x \in \left[-\pi,\pi\right), \quad \kappa \in \left[0,\infty\right), \]

where \(\mu\) is the mean direction, \(\kappa\) is the concentration, and \(\mathcal I_0(\cdot)\) is the modified Bessel function of the first kind of order zero. Angles may equivalently be represented on \([0,2\pi)\). At \(\kappa=0\) the distribution is circular uniform; as \(\kappa\) increases, progressively more mass is placed near \(\mu\). Thus, \(\kappa\) behaves more like a precision than a variance parameter.

1.2 Inferential ambiguity

For a single circular variable, the signed angular residual can be represented uniquely on a chosen principal interval. Regression is more delicate because the mean direction varies by observation, for example \(\mu_i=g(\eta_i)\) for a real-valued predictor \(\eta_i\). Across observations, the unwrapped differences \(x_i-\mu_i\) can span \([-2\pi,2\pi)\), while the cosine term remains periodic. Distinct coefficient values can therefore give the same cosine values and, for some data configurations, equivalent modes of the likelihood.

This is an identifiability problem for the regression parameters rather than merely a choice of angular notation. Equivalent modes can lead different initial values to different solutions and make posterior exploration unnecessarily difficult. The repeated extrema below illustrate the source of the ambiguity.

2 Link-adjusted von Mises (LAvM) distribution

The link-adjusted circular construction resolves this ambiguity by taking the residual in linear space and only then mapping it back to the circle. With the inverse tangent link

\[ g(\eta)=2\arctan(\eta), \qquad g^{-1}(x)=\tan(x/2), \]

the circular residual entering the vM kernel is \(g\{g^{-1}(x)-\eta\}\). This residual has one principal representation, while the Jacobian of the transformation preserves a valid density. The resulting link-adjusted von Mises density is

\[ p_{\operatorname{LAvM}}\left(x \mid \eta, \kappa\right) = \frac{\exp\left\{ \kappa \cos\left( 2 \arctan \left( \tan \left( \frac{x}{2} \right) - \eta \right) \right) \right\}}{2\pi I_{0}\left(\kappa\right)\left(1 + \eta^{2} - \eta \sin\left(x\right) - \eta^{2}\sin^{2}\left(\frac{x}{2}\right)\right)},\quad x\in\left(-\pi,\pi\right), \quad \eta\in\mathbb{R}, \quad \kappa \in \left[0,\infty\right) \]

The mode is \(g(\eta)=2\arctan(\eta)\). When \(\eta=0\), the LAvM density reduces to a vM density centered at zero. The denominator is the Jacobian adjustment, so it is an essential part of the likelihood rather than a normalizing correction that can be omitted.

The inverse link is singular at \(|x|=\pi\), so the LAvM construction is intended for reasonably concentrated data with an empty or nearly empty arc that can be placed at the boundary. A rose diagram should therefore be checked before fitting. If necessary, rotate the observations so that their main mass lies near zero, then rotate fitted directions and circular summaries back afterwards. Data spread almost uniformly around the full circle are not suitable for this construction.

library(ggplot2)
library(patchwork)
library(INLAcircular)

col1 <- "#77AFA9"; col1_t <- adjustcolor(col1, alpha.f = 0.6)
col2 <- "#DD9871"; col2_t <- adjustcolor(col2, alpha.f = 0.6)
col3 <- "#839BB2"; col3_t <- adjustcolor(col3, alpha.f = 0.6)
col4 <- "#C08A98"; col4_t <- adjustcolor(col4, alpha.f = 0.6)
col5 <- "#E5C48A"; col5_t <- adjustcolor(col5, alpha.f = 0.6)

For fixed \(\kappa\), changing \(\eta\) moves the mode through the inverse tangent link and also changes the local shape through the Jacobian. Increasing \(\kappa\) concentrates all three curves more strongly. This distinction matters in the later examples: \(\kappa\) is the baseline concentration at \(\eta=0\), not a common vM-scale concentration for every observation when \(\eta_i\ne0\).

2.1 Simulation utilities

# Install INLAcircular
if (!requireNamespace("remotes", quietly = TRUE)) {
  install.packages("remotes")
}
remotes::install_github(
  "XiangYEstats/INLAcircular",
  dependencies = TRUE
)

library(INLAcircular)

The examples require vM density, distribution, quantile, and simulation operations. The following self-contained functions implement them through a monotone spline approximation to the CDF, so the simulations can also be followed without relying on package internals.

dvm <- function(x, mu, kappa, log = FALSE) {
  n <- max(length(x), length(mu), length(kappa))
  log.bessel.scaled <- numeric(n)
  kappa_long <- rep(kappa, length.out = n)
  large_kappa <- kappa_long > 1e5
  if (any(large_kappa)) {
    k_large <- kappa_long[large_kappa]
    log.bessel.scaled[large_kappa] <- -0.5 * log(2 * pi * k_large) + 1 / (8 * k_large)
  }
  if (any(!large_kappa)) {
    k_small <- kappa_long[!large_kappa]
    log.bessel.scaled[!large_kappa] <- log(besselI(k_small, 0, expon.scaled = TRUE))
  }
  log.dens <- -2 * kappa * sin((x - mu) / 2)^2 - log(2 * pi) - log.bessel.scaled
  if (log) {
    return(log.dens)
  } else {
    return(exp(log.dens))
  }
}

vm.spline.grid <- function(kappa, len = 2048L) {
  s <- 1 / sqrt(kappa)
  if (s < 1/sqrt(15)) {
    x.grid <- seq(-5*s, 5*s, length.out = len)
  } else {
    x.grid <- seq(-pi, pi, length.out = len)
  }
  fx <- dvm(x.grid, mu = 0, kappa = kappa, log = FALSE)
  Fx <- cumsum(fx)
  Fx.norm <- (Fx - Fx[1]) / (Fx[length(Fx)] - Fx[1])
  keep <- (Fx.norm > 0) & (Fx.norm < 1)
  if (sum(keep) < 4) {
    return(data.frame(x = x.grid, logit.Fx = qlogis(punif(seq(0,1,length.out=length(x.grid))))))
  }
  return(data.frame(x = x.grid[keep], logit.Fx = qlogis(Fx.norm[keep])))
}

pvm <- function(q, mu, kappa, strategy = "circular", log = FALSE, len = 2048L) {
  grids <- vm.spline.grid(kappa = kappa, len = len)
  L.spline <- splinefun(grids$x, grids$logit.Fx)
  if (strategy == "linear") {
    logit_val <- L.spline(q - mu)
  } else if (strategy == "circular") {
        z <- (q - mu + pi) %% (2 * pi) - pi
    logit_val <- L.spline(z)
  } else {
    stop("strategy must be 'circular' or 'linear'")
  }
  if (log) {
    return(plogis(logit_val, log.p = TRUE))
  } else {
    return(plogis(logit_val))
  }
}

qvm <- function(p, mu, kappa, len = 2048L) {
  n_p <- length(p)
  n_mu <- length(mu)
  if (n_mu == 1) {
    mu_vec <- rep(mu, n_p)
  } else if (n_mu == n_p) {
    mu_vec <- mu
  } else {
    stop("Error: length(mu) must be 1 or equal to length(p).")
  }
  grids <- vm.spline.grid(kappa = kappa, len = len)
  inv.spline <- splinefun(grids$logit.Fx, grids$x, method = "monoH.FC")
  p_bounded <- pmin(pmax(p, 0), 1)
  logit_p <- qlogis(p_bounded)
  min_logit <- min(grids$logit.Fx)
  max_logit <- max(grids$logit.Fx)
  logit_p_clamped <- pmin(pmax(logit_p, min_logit), max_logit)
  z <- inv.spline(logit_p_clamped)
  q_val <- z + mu_vec
  q_val[p <= 0] <- mu_vec[p <= 0] - pi
  q_val[p >= 1] <- mu_vec[p >= 1] + pi
  q_val[is.na(p)] <- NaN
  return(q_val)
}

rvm <- function(n, mu, kappa, len = 2048L) {
  if (length(mu) == 1) {
    mu <- rep(mu, n)
  } else if (length(mu) != n) {
    stop("length(mu) must be 1 or n")
  }
  if (length(kappa) == 1) {
    kappa <- rep(kappa, n)
  } else if (length(kappa) != n) {
    stop("length(kappa) must be 1 or n")
  }
  unique_k <- unique(kappa)
  unique_mu <- unique(mu)
  if (length(unique_k) == 1 && length(unique_mu) == 1) {
    u <- runif(n)
    return(qvm(p = u, mu = unique_mu, kappa = unique_k, len = len))
  }
  u <- runif(n)
  res <- numeric(n)
  for (k in unique_k) {
    idx <- which(kappa == k)
    res[idx] <- qvm(p = u[idx], mu = mu[idx], kappa = k, len = len)
  }
  return(res)
}

The simulation function accepts an observation-specific mean direction through vector-valued mu.

2.2 Simulation example

The first simulation isolates the fixed-effects part of the model. Let

\[ \eta_i=\beta_0+\beta_1z_i, \qquad \mu_i=g(\eta_i)=2\arctan(\eta_i), \qquad y_i\sim\operatorname{vM}(\mu_i,20). \]

We use \(n=100\), \(\beta_0=0.1\), \(\beta_1=1.2\), and a narrowly distributed covariate \(z_i\sim N(0,0.1^2)\). The data-generating distribution is the standard vM rather than LAvM. Because the predictors remain close to zero, this is a transparent first check that LAvM regression can recover the mean relationship and concentration in a setting where the link adjustment is modest.

library(circular)

# inverse tangent link for circular response
lfun <- function(x) {2 * atan(x)}
n <- 100
set.seed(1314)
z <- rnorm(n, mean = 0, sd = 0.1)
beta0 <- 0.1; beta1 <- 1.2
# eta = beta0 + beta1*z
y <- rvm(n, mu = lfun(beta0 + beta1*z), kappa = 20)

par(mfrow = c(1, 2), pty = "s", mar = c(2, 0.5, 2, 0.5))
hist(y, breaks = 10, col = col1_t, border = col1, main="", xlab = "", ylab = "")
rose.diag(y, bins = 20, col = col1_t, border = col1, prop = 1.5)

The ordinary histogram treats an angle as a linear number and can split a cluster at the chosen boundary, whereas the rose diagram respects the circular geometry. Here both displays show a concentrated arc away from \(\pm\pi\), which is the setting required by the LAvM construction.

We fit

\[ y_i\mid\eta_i,\kappa\sim\operatorname{LAvM}(\eta_i,\kappa), \qquad \eta_i=\beta_0+\beta_1z_i. \]

We use \(N(0,5^2)\) priors for the regression coefficients. For \(\kappa\), a penalized-complexity (PC) prior shrinks toward the simpler base model of no circular spread: a point mass at the mean direction, obtained as \(\kappa\rightarrow\infty\). Departure from that base is measured using the mean resultant length

\[ \rho(\kappa)=\frac{\mathcal I_1(\kappa)}{\mathcal I_0(\kappa)}, \]

which ranges from \(0\) for a circular-uniform distribution to \(1\) for a point mass. The prior is calibrated by the interpretable statement \(\Pr\{\rho(\kappa)>U\}=\alpha\); here \(U=0.5\) and \(\alpha=0.5\) give a weak default calibration.

library(cmdstanr)

sim1_stan_file <- file.path("stan", "lavm_simulation_1.stan")

functions {
  real lavm_lpdf(vector x, vector eta, real kappa) {
    int N = num_elements(x);
    real log_lik = 0;
    real log_I0_val;
    if (kappa < 10000) {
      log_I0_val = log_modified_bessel_first_kind(0, kappa);
    } else {
      log_I0_val = kappa - 0.5 * log(2 * pi() * kappa) + 1.0 / (8.0 * kappa);
    }
    real log_bessel = log(2 * pi()) + log_I0_val;
    for (n in 1:N) {
      real tan_half_x = tan(x[n] / 2.0);
      real A = cos(2.0 * atan(tan_half_x - eta[n]));
      real D = 1.0 + square(eta[n]) - eta[n] * sin(x[n]) - square(eta[n]) * square(sin(x[n] / 2.0));
      log_lik += (kappa * A) - log_bessel - log(D);
    }
    return log_lik;
  }
}

data {
  int<lower=1> N;
  vector[N] y;
  vector[N] z;
  real<lower=0, upper=1> U;
  real<lower=0, upper=1> alpha;
}

transformed data {
  // Compute lambda based on mean resultant length
  real lambda = -log(1.0 - alpha) / sqrt(1.0 - U);
}

parameters {
  real beta0;
  real beta1;
  real<lower=0> kappa;
}

transformed parameters {
  vector[N] eta = beta0 + beta1 * z;
}

model {
  // Linear Priors
  beta0 ~ normal(0, 5);
  beta1 ~ normal(0, 5);
  
  // PC prior for kappa (Point Mass Base Model)
  real numerator;  
  real d_sq; 
  
  if (kappa < 10000) {
    real log_I0 = log_modified_bessel_first_kind(0, kappa);
    real log_I1 = log_modified_bessel_first_kind(1, kappa);
    real log_I2 = log_modified_bessel_first_kind(2, kappa);
    real I1_over_I0 = exp(log_I1 - log_I0);
    real I2_over_I0 = exp(log_I2 - log_I0);
    d_sq = 1.0 - I1_over_I0;
    numerator = 0.5 * (1.0 + I2_over_I0) - square(I1_over_I0);
  } else {
    d_sq = 0.5 / kappa + 0.125 / square(kappa);
    numerator = 0.5 / square(kappa) + 0.25 / (kappa^3);
  }
  
  if (d_sq < 1e-12) d_sq = 1e-12;
  if (numerator < 1e-12) numerator = 1e-12;
  
  target += log(lambda) + log(numerator) - log(2.0) - 0.5 * log(d_sq) - lambda * sqrt(d_sq);
  
  // Data model
  y ~ lavm(eta, kappa);
}

mod <- cmdstan_model(sim1_stan_file)

The custom density is the logarithm of the LAvM expression introduced above: its cosine term supplies the vM kernel and -log(D) supplies the Jacobian adjustment. The transformed-data block converts \((U,\alpha)\) to the PC-prior rate, and an asymptotic expansion for \(\log I_0(\kappa)\) avoids overflow at large \(\kappa\).

fit <- mod$sample(
  data = list(N = n, y = y, z = z, U = 0.5, alpha = 0.5),
  seed = 131401,
  refresh = 0
)

fit$summary(variables = c("beta0", "beta1", "kappa"))

## # A tibble: 3 × 10
##   variable   mean median     sd    mad      q5    q95  rhat ess_bulk ess_tail
##   <chr>     <dbl>  <dbl>  <dbl>  <dbl>   <dbl>  <dbl> <dbl>    <dbl>    <dbl>
## 1 beta0     0.107  0.107 0.0115 0.0115  0.0884  0.126  1.00    3454.    2750.
## 2 beta1     1.26   1.26  0.113  0.114   1.08    1.45   1.00    2964.    2901.
## 3 kappa    19.9   19.7   2.85   2.88   15.4    24.8    1.00    3935.    3063.

The posterior means are approximately \(0.107\) for \(\beta_0\), \(1.26\) for \(\beta_1\), and \(19.9\) for \(\kappa\), close to the generating values \(0.1\), \(1.2\), and \(20\).

library(bayesplot)

# 1. Extract draws directly into a data frame
draws <- fit$draws(variables = c("beta0", "beta1", "kappa"), format = "df")

# 2. Build individual density plots
p_beta0 <- ggplot(draws, aes(x = beta0)) +
  geom_density(fill = col1, alpha = 0.8, color = "white") +
  geom_vline(xintercept = 0.1, linetype = "dashed", color = "black", linewidth = 0.6) +
  labs(title = expression(beta[0]), x = "", y = "Posterior Density") +
  theme_minimal() +
  theme(plot.title = element_text(face = "bold", hjust = 0.5))

p_beta1 <- ggplot(draws, aes(x = beta1)) +
  geom_density(fill = col2, alpha = 0.8, color = "white") +
  geom_vline(xintercept = 1.2, linetype = "dashed", color = "black", linewidth = 0.6) +
  labs(title = expression(beta[1]), x = "", y = "") +
  theme_minimal() +
  theme(plot.title = element_text(face = "bold", hjust = 0.5))

p_kappa <- ggplot(draws, aes(x = kappa)) +
  geom_density(fill = col3, alpha = 0.8, color = "white") +
  geom_vline(xintercept = 20, linetype = "dashed", color = "black", linewidth = 0.6) +
  labs(title = expression(kappa), x = "", y = "") +
  theme_minimal() +
  theme(plot.title = element_text(face = "bold", hjust = 0.5))

# 3. Combine into a 1x3 layout using patchwork
(p_beta0 | p_beta1 | p_kappa) + 
  plot_annotation(
    title = "Posterior Distributions",
    theme = theme(plot.title = element_text(face = "bold", size = 15, hjust = 0.5))
  )

The marginal posteriors are centred near the generating values, and the slope is clearly separated from zero.

2.3 Simulation example 2

The second simulation adds a smooth latent component. The data-generating model is

\[ y_i \sim \operatorname{vM}\left(g(\eta_i),100\right), \qquad \eta_i=\beta_0+\beta_1z_i+w_i, \qquad z_i\sim N(0,1), \]

where \(g(\eta)=2\arctan(\eta)\) and \(\mathbf w\) is a second-order random walk (RW2). An RW2 penalizes changes in the first difference, making it suitable for a smoothly evolving trajectory rather than independent observation-level noise. Its precision matrix is intrinsically rank deficient because constant and linear trends lie in its null space. Level and linear-trend constraints are therefore needed to keep those components from being confounded with the fixed effects.

The simulation constructs the second-difference precision \(\mathbf Q=\mathbf D^\mathsf{T}\mathbf D\), projects the draw away from the constant and linear null-space components, and scales the resulting path to standard deviation \(\sigma_w=0.4\). The fitted model imposes the analogous constraints on the latent trajectory.

n <- 100
sigma_w <- 0.4

# Construct the second-difference matrix for an RW2 process
D <- Matrix::sparseMatrix(
  i = rep(seq_len(n - 2), each = 3),
  j = as.vector(rbind(seq_len(n - 2), 2:(n - 1), 3:n)),
  x = rep(c(1, -2, 1), times = n - 2),
  dims = c(n - 2, n)
)
Q <- Matrix::crossprod(D)
Q_eps <- as.matrix(Q + Matrix::Diagonal(n, x = 1e-5))
R_Q <- chol(Q_eps)

set.seed(520)
z_innov <- rnorm(n)
w_unconstrained <- as.numeric(backsolve(R_Q, z_innov))
A <- matrix(c(rep(1, n), 1:n), nrow = 2, byrow = TRUE)
QA_t <- solve(Q_eps, t(A)) 
w_proj <- w_unconstrained - QA_t %*% solve(A %*% QA_t) %*% (A %*% w_unconstrained)

# Standardize to sd=1 before multiplying by sigma_w
w <- as.numeric(scale(w_proj)) * sigma_w

# Simulate covariate and response
z <- rnorm(n, mean = 0, sd = 1)
beta0 <- 0.1; beta1 <- 0.2
kappa_true <- 100

eta <- beta0 + beta1 * z + w
y <- rvm(n, mu = lfun(eta), kappa = kappa_true)

par(mfrow = c(1, 3), pty = "s", mar = c(2, 2, 2, 0.5))
plot(1:n, w, type = "l", col = "black", lwd = 2, main = "Latent RW2 Path (w)", xlab = "", ylab = "")
hist(y, breaks = 10, col = col1_t, border = col1, main="y", xlab = "", ylab = "")
rose.diag(y, bins = 20, col = col1_t, border = col1, prop = 1.5, main = "x")

The first panel displays the latent path that the model must recover; the other two show the response on linear and circular scales. Together they help distinguish smooth movement in the mean direction from unexplained angular dispersion and verify that the responses remain away from the branch boundary.

For \(\kappa\), we again use the PC prior whose base model is a point mass, corresponding to \(\kappa\rightarrow\infty\):

\[ p\left(\kappa\right) = \lambda \exp\left\{ - \lambda \sqrt{1 - \frac{\mathcal{I}_{1}(\kappa)}{\mathcal{I}_{0}(\kappa)}} \right\} \frac{ \frac{\mathcal{I}_{0}(\kappa) + \mathcal{I}_{2}(\kappa)}{2\mathcal{I}_{0}(\kappa)} - \frac{\mathcal{I}_{1}(\kappa)^{2}}{\mathcal{I}_{0}(\kappa)^{2}} }{2\sqrt{1 - \frac{\mathcal{I}_{1}(\kappa)}{\mathcal{I}_{0}(\kappa)}}}, \]

The rate \(\lambda\) is determined by the same probability statement \(\Pr\{\rho(\kappa)>U\}=\alpha\) used in the first example. We again set \(U=\alpha=0.5\). The fixed effects have \(N(0,5^2)\) priors, and the positive RW2 innovation scale has a half-normal prior.

sim2_stan_file <- file.path("stan", "lavm_simulation_2.stan")

functions {
  // 1. LAvM log density
  real lavm_lpdf(vector x, vector eta, real kappa) {
    int N = num_elements(x);
    real log_lik = 0;
    
    real log_I0_val;
    if (kappa < 10000) {
      log_I0_val = log_modified_bessel_first_kind(0, kappa);
    } else {
      log_I0_val = kappa - 0.5 * log(2 * pi() * kappa) + 1.0 / (8.0 * kappa);
    }
    real log_bessel = log(2 * pi()) + log_I0_val;
    
    for (n in 1:N) {
      real tan_half_x = tan(x[n] / 2.0);
      real A = cos(2.0 * atan(tan_half_x - eta[n]));
      real D = 1.0 + square(eta[n]) - eta[n] * sin(x[n]) - square(eta[n]) * square(sin(x[n] / 2.0));
      log_lik += (kappa * A) - log_bessel - log(D);
    }
    return log_lik;
  }
  
  // 2. Custom Random Number Generator for Posterior Predictive Checks
  real lavm_rng(real eta, real kappa) {
    int G = 400; // Grid resolution for numerical sampling
    vector[G] log_probs; 
    real step_size = (2.0 * pi()) / G;
    real start_val = -pi() + step_size / 2.0;
    
    for (g in 1:G) {
      real x = start_val + (g - 1) * step_size;
      real tan_half_x = tan(x / 2.0);
      real A = cos(2.0 * atan(tan_half_x - eta));
      real D = 1.0 + square(eta) - eta * sin(x) - square(eta) * square(sin(x / 2.0));
      
      // Calculate securely in log space to prevent overflow
      log_probs[g] = (kappa * A) - log(D);
    }
    
    int idx = categorical_rng(softmax(log_probs));
    return start_val + (idx - 1) * step_size;
  }
}

data {
  int<lower=1> N;
  vector[N] y;
  vector[N] z;
  vector[N] seq_idx;
  real<lower=0, upper=1> U;
  real<lower=0, upper=1> alpha;
}

transformed data {
  // Hyperparameter for the Penalized Complexity (PC) prior
  real lambda = -log(1.0 - alpha) / sqrt(1.0 - U);
}

parameters {
  real beta0;
  real beta1;
  real<lower=0> kappa; // Sampled directly on the natural scale
  
  vector[N] w;               // The latent Gaussian process
  real<lower=0> sigma_wi;    // SD of the individual step-to-step innovations (roughness)
}

transformed parameters {
  vector[N] eta = beta0 + beta1 * z + w;
}

model {
  // PC Prior Implementation for kappa
  real numerator;
  real d_sq;
  
  if (kappa < 10000) {
    real log_I0 = log_modified_bessel_first_kind(0, kappa);
    real log_I1 = log_modified_bessel_first_kind(1, kappa);
    real log_I2 = log_modified_bessel_first_kind(2, kappa);
    
    real I1_over_I0 = exp(log_I1 - log_I0);
    real I2_over_I0 = exp(log_I2 - log_I0);
    
    d_sq = 1.0 - I1_over_I0;
    numerator = 0.5 * (1.0 + I2_over_I0) - square(I1_over_I0);
  } else {
    // Asymptotic expansion for numerical stability at large kappa
    d_sq = 0.5 / kappa + 0.125 / square(kappa);
    numerator = 0.5 / square(kappa) + 0.25 / (kappa^3);
  }
  
  if (d_sq < 1e-12) d_sq = 1e-12;
  if (numerator < 1e-12) numerator = 1e-12;
  
  target += log(lambda) + log(numerator) - log(2.0) - 0.5 * log(d_sq) - lambda * sqrt(d_sq);
  
  // Fixed Effects Priors
  beta0 ~ normal(0, 5);
  beta1 ~ normal(0, 5);
  
  // Variance Prior
  sigma_wi ~ std_normal(); 
  
  // --- Latent Field (RW2) Generative Story ---
  w[1] ~ normal(0, 10);
  w[2] ~ normal(0, 10);
  
  for (i in 3:N) {
    w[i] ~ normal(2.0 * w[i-1] - w[i-2], sigma_wi);
  }
  
  // Soft constraints to preserve identifiability
  sum(w) ~ normal(0, 0.001 * N);
  dot_product(w, seq_idx) ~ normal(0, 0.001 * N);
  
  // Data model
  y ~ lavm(eta, kappa);
}

generated quantities {
  vector[N] log_lik;
  array[N] real y_rep;
  
  // Derived Global Parameters
  real sigma_w = sd(w); 
  real kappa_eff; 
  
  { // Local scope to prevent saving intermediate variables to MCMC draws
  
    // 1. Gaussian-approximation effective concentration of the LAvM distribution
    vector[N] eta_sq = square(eta);
    vector[N] term1 = kappa * square(1.0 + eta_sq);
    vector[N] term2 = 0.5 * (eta_sq .* (1.0 + eta_sq));
    kappa_eff = mean(term1 + term2);
  
    // 2. PPC and LOO Calculations
    real log_I0_gen = kappa < 10000 ? log_modified_bessel_first_kind(0, kappa) 
                                    : (kappa - 0.5 * log(2 * pi() * kappa) + 1.0 / (8.0 * kappa));
                                    
    real log_bessel = log(2.0 * pi()) + log_I0_gen;
    
    for (n in 1:N) {
      // Pointwise log-likelihood for LOO-CV
      real tan_half_x = tan(y[n] / 2.0);
      real A = cos(2.0 * atan(tan_half_x - eta[n]));
      real D = 1.0 + square(eta[n]) - eta[n] * sin(y[n]) - square(eta[n]) * square(sin(y[n] / 2.0));
      log_lik[n] = (kappa * A) - log_bessel - log(D);
      
      // PPC Generation via grid sampler
      y_rep[n] = lavm_rng(eta[n], kappa);
    }
  }
}

mod <- cmdstan_model(sim2_stan_file)

Relative to the first Stan program, this model estimates the constrained RW2 trajectory and its innovation scale. The two soft constraints keep the latent level and trend separate from \(\beta_0\) and \(\beta_1\). Its generated quantities provide pointwise log_lik, replicated angles y_rep, and the draw-specific average effective concentration kappa_eff used below.

stan_data <- list(
  N = n,
  y = as.numeric(y),
  z = z,
  seq_idx = 1:n,
  U = 0.5,
  alpha = 0.5
)

fit <- mod$sample(
  data = stan_data,
  seed = 52002,
  refresh = 0
)

fit$summary(variables = c("beta0", "beta1", "kappa", "kappa_eff", "sigma_w"))

## # A tibble: 5 × 10
##   variable   mean median      sd     mad      q5     q95  rhat ess_bulk ess_tail
##   <chr>     <dbl>  <dbl>   <dbl>   <dbl>   <dbl>   <dbl> <dbl>    <dbl>    <dbl>
## 1 beta0     0.104  0.104 6.50e-3 6.28e-3  0.0936   0.115  1.00    5004.    2191.
## 2 beta1     0.194  0.194 7.81e-3 7.60e-3  0.182    0.207  1.00    3804.    2888.
## 3 kappa    66.6   66.1   1.03e+1 1.04e+1 50.7     84.1    1.00    4093.    2885.
## 4 kappa_e… 98.9   98.2   1.54e+1 1.54e+1 74.7    125.     1.00    4079.    2808.
## 5 sigma_w   0.402  0.403 6.45e-3 6.26e-3  0.392    0.413  1.00    2272.    2656.

The posterior means are approximately \(0.104\) for \(\beta_0\), \(0.194\) for \(\beta_1\), \(98.9\) for the average effective concentration, and \(0.402\) for the empirical standard deviation of \(\mathbf w\), close to the generating values \(0.1\), \(0.2\), \(100\), and \(0.4\). The raw LAvM baseline concentration has posterior mean \(66.6\); its different scale is explained below.

# Extract parameters
draws <- fit$draws(variables = c("beta0", "beta1", "kappa", "kappa_eff", "sigma_w"), format = "df")

p_beta0 <- ggplot(draws, aes(x = beta0)) +
  geom_density(fill = col1_t, alpha = 0.8, color = "white") +
  geom_vline(xintercept = beta0, linetype = "dashed", color = "black", linewidth = 0.6) +
  labs(title = expression(beta[0]), x = "", y = "Density") +
  theme_minimal() + theme(plot.title = element_text(face = "bold", hjust = 0.5))

p_beta1 <- ggplot(draws, aes(x = beta1)) +
  geom_density(fill = col2_t, alpha = 0.8, color = "white") +
  geom_vline(xintercept = beta1, linetype = "dashed", color = "black", linewidth = 0.6) +
  labs(title = expression(beta[1]), x = "", y = "") +
  theme_minimal() + theme(plot.title = element_text(face = "bold", hjust = 0.5))

p_sigma_w <- ggplot(draws, aes(x = sigma_w)) +
  geom_density(fill = col4_t, alpha = 0.8, color = "white") +
  geom_vline(xintercept = sigma_w, linetype = "dashed", color = "black", linewidth = 0.6) +
  labs(title = expression(sigma[w]), x = "", y = "") +
  theme_minimal() + theme(plot.title = element_text(face = "bold", hjust = 0.5))

p_kappa <- ggplot(draws, aes(x = kappa)) +
  geom_density(fill = col3_t, alpha = 0.8, color = "white") +
  geom_vline(xintercept = kappa_true, linetype = "dashed", color = "black", linewidth = 0.6) +
  labs(title = expression(kappa), x = "", y = "") +
  theme_minimal() + theme(plot.title = element_text(face = "bold", hjust = 0.5))

p_kappa_eff <- ggplot(draws, aes(x = kappa_eff)) +
  geom_density(fill = col5_t, alpha = 0.8, color = "white") +
  geom_vline(xintercept = kappa_true, linetype = "dashed", color = "black", linewidth = 0.6) +
  labs(title = expression(kappa[eff]), x = "", y = "Density") +
  theme_minimal() + theme(plot.title = element_text(face = "bold", hjust = 0.5))


(p_beta0 | p_beta1 | p_sigma_w) / (p_kappa | p_kappa_eff | plot_spacer()) + 
  plot_annotation(
    title = "Posterior Distributions",
    theme = theme(plot.title = element_text(face = "bold", size = 15, hjust = 0.5))
  )

Raw \(\kappa\) is the LAvM baseline at \(\eta=0\), whereas \(\kappa_{\mathrm{eff}}\) is on the vM-comparable scale. The posterior of \(\kappa\) should not be expected to recover the generating vM concentration directly because the data were generated from a vM distribution but fitted with an LAvM distribution. Under a local Gaussian approximation, observation \(i\) has effective concentration

\[ \kappa_{\mathrm{eff},i} =\kappa(1+\eta_i^2)^2+\frac{1}{2}\eta_i^2(1+\eta_i^2) \approx\kappa(1+\eta_i^2)^2 \]

for large \(\kappa\). The reported value averages \(\kappa_{\mathrm{eff},i}\) within each posterior draw, preserving its dependence on \(\kappa\) and \(\boldsymbol\eta\).

library(loo)

if (packageVersion("bayesplot") < package_version("1.16.0")) {
  stop("The correlated LOO-PIT ECDF requires bayesplot 1.16.0 or later.")
}

The LOO-PIT ECDF below uses the dependence-aware correlated method of Tesso and Vehtari (2026), following the Bayesian Workflow example. This method requires bayesplot 1.16.0 or later.

# Latent field w
w_summary <- fit$summary("w")
w_df <- data.frame(
  Index = 1:n,
  True_w = w,
  Est_w = w_summary$mean,
  Lower = w_summary$q5,
  Upper = w_summary$q95
)

# LOO and y_rep
log_lik_matrix <- fit$draws("log_lik", format = "matrix")
loo_res <- loo(log_lik_matrix, save_psis = TRUE)
yrep <- fit$draws("y_rep", format = "matrix")

wrap_to_pi <- function(x) { (x + pi) %% (2 * pi) - pi }
y_shifted <- wrap_to_pi(as.numeric(y))
yrep_shifted <- wrap_to_pi(yrep)

# Top-Left: Latent Process
p_latent <- ggplot(w_df, aes(x = Index)) +
  geom_ribbon(aes(ymin = Lower, ymax = Upper), fill = col1, alpha = 0.3) +
  geom_line(aes(y = Est_w, color = "Posterior Mean"), linewidth = 1) +
  geom_line(aes(y = True_w, color = "True Latent Field"), linewidth = 0.8, linetype = "dashed") +
  scale_color_manual(values = c("Posterior Mean" = col1, "True Latent Field" = "black")) +
  labs(title = "Latent Process: True vs. Estimated", x = "Index", y = expression(w), color = "") +
  theme_minimal() +
  theme(plot.title = element_text(face = "bold", size = 13, hjust = 0.5))

# Top-Right: Posterior Predictive Check
set.seed(520)
yrep_sample <- yrep_shifted[sample(nrow(yrep_shifted), 100), ]
p_ppc <- ppc_dens_overlay(y = y_shifted, yrep = yrep_sample) +
  scale_x_continuous(limits = c(-pi, pi), breaks = c(-pi, 0, pi), 
                     labels = c(expression(-pi), "0", expression(pi))) +
  labs(title = "Posterior Predictive Check", x = "Circular Response", y = "Density") +
  theme_minimal() +
  theme(legend.position = "none", plot.title = element_text(face = "bold", size = 13, hjust = 0.5))

# Bottom-Left: Pareto K Diagnostics
k_df <- data.frame(Observation = 1:n, k = loo_res$diagnostics$pareto_k)
p_k <- ggplot(k_df, aes(x = Observation, y = k)) +
  geom_point(color = col3, size = 2, alpha = 0.8) +
  geom_hline(yintercept = 0.7, linetype = "dashed", color = col2, linewidth = 0.8) +
  geom_hline(yintercept = 0, color = "black", linewidth = 0.3) +
  labs(title = "PSIS Diagnostics", subtitle = "Pareto k values (Threshold = 0.7)", 
       x = "Data Point Index", y = expression(Pareto~k)) +
  theme_minimal() +
  theme(plot.title = element_text(face = "bold", size = 13, hjust = 0.5),
        plot.subtitle = element_text(hjust = 0.5))

# Bottom-Right: LOO-PIT ECDF
p_pit <- bayesplot::ppc_loo_pit_ecdf(
  y = y_shifted,
  yrep = yrep_shifted,
  psis_object = loo_res$psis_object,
  method = "correlated"
) +
  labs(title = "LOO-PIT ECDF Calibration") +
  theme_minimal() +
  theme(legend.position = "none", plot.title = element_text(face = "bold", size = 13, hjust = 0.5))

(p_latent | p_ppc) / (p_k | p_pit) + 
  plot_annotation(
    title = "Model Diagnostics and Validation",
    theme = theme(plot.title = element_text(face = "bold", size = 16, hjust = 0.5))
  ) +
  plot_layout(guides = 'collect') & theme(legend.position = 'bottom')

The grid checks four aspects of the toy fit: recovery of the known latent path, agreement between observed and replicated angular distributions, influential observations through Pareto-\(k\), and predictive calibration through LOO-PIT. These are complementary summaries of this example rather than a detailed model-comparison exercise.

2.4 Wind application example

The final example uses the January subset of hourly wind data from John F. Kennedy International Airport. The data include wind direction in radians, temperature, time, and hour of day. The accompanying article develops a joint model for wind direction and speed; to keep the example focused on LAvM regression, we fit only the circular wind-direction component here.

library(dplyr)

Download the prepared wind_data.rds file and save it in the same folder as this tutorial before running the following code.

wind_data <- readRDS("wind_data.rds")
dat <- wind_data %>% filter(month %in% c(1))
n <- nrow(dat)

x <- as.numeric(dat$HLY.WIND.VCTDIR)
z <- dat$HLY.TEMP.NORMAL
hour_idx <- dat$hour + 1 

# Rotate around the circular mean and wrap to (-pi, pi]
rotation <- as.numeric(circular::mean.circular(circular::circular(x)))
x_centered <- (x - rotation + pi) %% (2 * pi) - pi

This rotation is a modeling step, not a cosmetic change. It places the sample circular mean at zero and wraps the observations to \([-\pi,\pi)\), moving the occupied part of the circle away from the LAvM singularity at the endpoints. The value rotation is retained so fitted or predicted directions can be translated back to the original compass coordinates.

plot_dat <- data.frame(
  direction = x,
  direction_rotated = x_centered %% (2 * pi),
  temp = z
)

p1 <- ggplot(plot_dat, aes(x = direction)) +
  geom_histogram(breaks = seq(0, 2 * pi, length.out = 37), fill = "#2A9D8F", color = "white", alpha = 0.9) +
  coord_polar(theta = "x", start = -pi / 2, direction = -1) +
  scale_x_continuous(limits = c(0, 2 * pi), 
                     breaks = c(0, pi / 2, pi, 3 * pi / 2), 
                     labels = c("E", "N", "W", "S")) +
  labs(x = "wind direction", y = NULL) +
  theme_minimal(base_size = 14) + 
  theme(axis.text.y = element_blank(), axis.ticks.y = element_blank(), aspect.ratio = 1)

p2 <- ggplot(plot_dat, aes(x = direction_rotated)) +
  geom_histogram(breaks = seq(0, 2 * pi, length.out = 37), fill = "#2A9D8F", color = "white", alpha = 0.9) +
  coord_polar(theta = "x", start = -pi / 2, direction = -1) +
  scale_x_continuous(limits = c(0, 2 * pi), 
                     breaks = c(0, pi / 2, pi, 3 * pi / 2), 
                     labels = c("E", "N", "W", "S")) +
  labs(x = "direction (rotated)", y = NULL) +
  theme_minimal(base_size = 14) + 
  theme(axis.text.y = element_blank(), axis.ticks.y = element_blank(), aspect.ratio = 1)

p3 <- ggplot(plot_dat, aes(x = temp)) +
  geom_histogram(bins = 30, fill = "#F4A261", color = "white", alpha = 0.9) +
  scale_x_continuous(n.breaks = 4, guide = guide_axis(check.overlap = TRUE)) +
  labs(x = "temperature (℃)", y = NULL) +
  theme_minimal(base_size = 14)

trace_dat <- data.frame(time = dat$datetime, direction = x)

p4 <- ggplot(trace_dat, aes(x = time, y = direction)) +
  geom_line(color = "#2A9D8F", alpha = 0.9, linewidth = 0.8) + 
  scale_y_continuous(limits = c(pi, 2 * pi), breaks = c(pi, 3 * pi / 2, 2 * pi), labels = c("S", "W", "N")) +
  scale_x_datetime(breaks = scales::breaks_pretty(n = 3), guide = guide_axis(check.overlap = TRUE)) +
  labs(x = NULL, y = "wind direction") + 
  theme_minimal(base_size = 14)

(p1 | p3) / (p2 | p4)

The two polar histograms verify the angular concentration before and after rotation. The temperature histogram shows the scale and support of the linear covariate, while the time trace reveals persistence that an independent-response model would miss. Rotation changes the coordinate origin but not the ordering of directions around the circle.

For centered wind direction \(x_i\), the fitted temporal model is

\[ \begin{aligned} x_i\mid w_i,w_{2,h(i)},z_i &\sim \operatorname{LAvM}\left(\eta_i,\kappa\right), \\ \eta_i &= a_0+a_1w_i+a_2w_{2,h(i)}+\alpha z_i, \\ \mathbf{w} &\sim N\left( 0, \mathbf{Q}_{\text{AR2}}^{-1} \right), \mathbf{w}_{2} \sim N\left( 0, \mathbf{Q}_{\text{RW2}}^{-1} \right), \\ a_{0}, \alpha &\sim N\left(0, 1\right), \quad \kappa \sim \mathcal{PC}_{\kappa}\left(0.5,0.99\right), \\ a_{1}, a_{2} &\sim \text{Exp}\left(-\log(0.5)/0.5\right), \\ (\phi_{1}, \phi_{2}) &\sim \operatorname{Uniform}(\mathcal{S}_{\text{AR2}}). \end{aligned} \]

Here \(a_0\) is the baseline linear predictor, \(\alpha\) is the temperature coefficient, and \(a_1\) and \(a_2\) scale the long-term and hour-of-day latent effects. Their positivity constraints remove a sign ambiguity between each latent process and its loading. On the centered circular scale, the conditional mode is \(2\arctan(\eta_i)\); on the original compass scale, add rotation and wrap the result modulo \(2\pi\).

The PC prior for \(\kappa\) uses the same mean-resultant-length calibration as before, now with \(\Pr\{\rho(\kappa)>0.5\}=0.99\). The exponential priors on \(a_1\) and \(a_2\) have median \(0.5\), regularizing the sizes of the two latent contributions while retaining their positive orientation.

The stationary coefficient region is the triangle

\[\mathcal{S}_{\text{AR2}} = \left\{ (\phi_{1}, \phi_{2}) \in \mathbb{R}^{2} \;\middle|\; -1 < \phi_{2} < 1 \quad \text{and} \quad \phi_{2} - 1 < \phi_{1} < 1 - \phi_{2} \right\}.\]

Within this region, the long-term latent process follows an autoregression of order two,

\[ w_i = \phi_1w_{i-1}+\phi_2w_{i-2}+\varepsilon_{w,i}, \qquad \varepsilon_{w,i}\sim N(0,1). \]

The 24-dimensional vector \(\mathbf w_2\) is an RW2 indexed by hour of day, with soft sum-to-zero and linear-trend constraints. The implementation is non-cyclic: it smooths adjacent hours but does not explicitly connect hour 24 to hour 1. A cyclic RW2 would be more appropriate if continuity across midnight were a central scientific requirement.

wind_stan_file <- file.path("stan", "lavm_wind_application.stan")

functions {
  // 1. Vectorized LAvM log density
  real lavm_lpdf(vector x, vector eta, real kappa) {
    int N = num_elements(x);
    real log_I0_val;
    
    if (kappa < 10000) {
      log_I0_val = log_modified_bessel_first_kind(0, kappa);
    } else {
      log_I0_val = kappa - 0.5 * log(2 * pi() * kappa) + 1.0 / (8.0 * kappa);
    }
    real log_bessel = log(2 * pi()) + log_I0_val;
    
    // Vectorized calculations
    vector[N] tan_half_x = tan(x / 2.0);
    vector[N] sin_x = sin(x);
    vector[N] sin_half_x_sq = square(sin(x / 2.0));
    
    vector[N] A = cos(2.0 * atan(tan_half_x - eta));
    vector[N] D = 1.0 + square(eta) - eta .* sin_x - square(eta) .* sin_half_x_sq;
    
    return sum((kappa * A) - log_bessel - log(D));
  }
  
  // 2. LAvM RNG
  real lavm_rng(real eta, real kappa) {
    int G = 400; 
    vector[G] log_probs; 
    real step_size = (2.0 * pi()) / G;
    real start_val = -pi() + step_size / 2.0;
    
    for (g in 1:G) {
      real x_val = start_val + (g - 1) * step_size;
      real tan_half_x = tan(x_val / 2.0);
      real A = cos(2.0 * atan(tan_half_x - eta));
      real D = 1.0 + square(eta) - eta * sin(x_val) - square(eta) * square(sin(x_val / 2.0));
      log_probs[g] = (kappa * A) - log(D);
    }
    int idx = categorical_rng(softmax(log_probs));
    return start_val + (idx - 1) * step_size;
  }
}

data {
  int<lower=1> N;
  vector[N] x;
  vector[N] z;
  array[N] int<lower=1, upper=24> hour_idx; 
  
  real<lower=0, upper=1> U;
  real<lower=0, upper=1> alpha_pc;
}

transformed data {
  real lambda = -log(1.0 - alpha_pc) / sqrt(1.0 - U);
  real lambda_tau = -log(0.5) / 0.5;
  
  vector[24] seq_idx_24;
  for (h in 1:24) seq_idx_24[h] = h;
}

parameters {
  // Fixed Effects
  real a0;
  real<lower=0> a1; 
  real<lower=0> a2; 
  real alpha;
  
  // Dispersion
  real log_kappa; 
  
  // AR(2) Process (w)
  real<lower=-1, upper=1> phi2;
  real<lower=phi2-1, upper=1-phi2> phi1;
  vector[N] z_w; 
  
  // RW2 Diurnal Process (w2)
  vector[24] w2;
}

transformed parameters {
  real<lower=0> kappa = exp(log_kappa);
  
  vector[N] w_unc;
  w_unc[1] = z_w[1] * 5.0;
  w_unc[2] = z_w[2] * 5.0;
  
  // Sequential dependency prevents vectorization here
  for (i in 3:N) {
    w_unc[i] = phi1 * w_unc[i-1] + phi2 * w_unc[i-2] + z_w[i];
  }
  
  vector[N] w = w_unc - mean(w_unc);
  vector[N] eta_x = a0 + a1 * w + a2 * w2[hour_idx] + alpha * z;
}

model {
  target += log_kappa; 
  
  // PC Prior for kappa
  real numerator;
  real d_sq;
  if (kappa < 10000) {
    real log_I0 = log_modified_bessel_first_kind(0, kappa);
    real log_I1 = log_modified_bessel_first_kind(1, kappa);
    real log_I2 = log_modified_bessel_first_kind(2, kappa);
    
    real I1_over_I0 = exp(log_I1 - log_I0);
    real I2_over_I0 = exp(log_I2 - log_I0);
    
    d_sq = 1.0 - I1_over_I0;
    numerator = 0.5 * (1.0 + I2_over_I0) - square(I1_over_I0);
  } else {
    d_sq = 0.5 / kappa + 0.125 / square(kappa);
    numerator = 0.5 / square(kappa) + 0.25 / (kappa^3);
  }
  if (d_sq < 1e-12) d_sq = 1e-12;
  if (numerator < 1e-12) numerator = 1e-12;
  target += log(lambda) + log(numerator) - log(2.0) - 0.5 * log(d_sq) - lambda * sqrt(d_sq);
  
  // Structural Priors (Optimized with std_normal)
  a0 ~ std_normal();
  alpha ~ std_normal();
  
  a1 ~ exponential(lambda_tau);
  a2 ~ exponential(lambda_tau);
  
  z_w ~ std_normal();
  
  // Vectorized RW2 Generative Model
  w2[1:2] ~ normal(0, 10);
  w2[3:24] ~ normal(2.0 * w2[2:23] - w2[1:22], 1.0); 
  sum(w2) ~ normal(0, 0.001 * 24); 
  dot_product(w2, seq_idx_24) ~ normal(0, 0.001 * 24);
  
  // Vectorized data model
  x ~ lavm(eta_x, kappa);
}

generated quantities {
  array[N] real x_rep;
  for (i in 1:N) {
    x_rep[i] = lavm_rng(eta_x[i], kappa);
  }
}

mod <- cmdstan_model(
  stan_file = wind_stan_file,
  cpp_options = list(
    STAN_CPP_OPTIMS = TRUE,
    STAN_NO_RANGE_CHECKS = TRUE
  )
)

The model samples log_kappa for numerical stability and adds the corresponding Jacobian term before applying the PC prior on the natural \(\kappa\) scale. The triangular parameter bounds keep the AR(2) coefficients in the stationary region, and centering \(\mathbf w\) separates its level from \(a_0\). The replicated angles x_rep are generated from the fitted LAvM density using a 400-point numerical grid.

stan_data <- list(
  N = n,
  x = x_centered, 
  z = z,
  hour_idx = hour_idx,
  U = 0.5,
  alpha_pc = 0.99
)

init_fun <- function() {
  list(
    a0 = rnorm(1, mean = 0, sd = 0.1),
    a1 = runif(1, min = 0.05, max = 0.2),
    a2 = runif(1, min = 0.05, max = 0.2),
    alpha = rnorm(1, mean = 0, sd = 0.1),
    log_kappa = log(runif(1, min = 3.0, max = 7.0)),
    
    phi2 = runif(1, min = -0.1, max = 0.1),
    phi1 = runif(1, min = 0.3, max = 0.7),
    z_w = rnorm(n, mean = 0, sd = 0.1),
    
    w2 = rnorm(24, mean = 0, sd = 0.1)
  )
}

set.seed(202601)
fit <- mod$sample(
  data = stan_data,
  init = init_fun,
  seed = 202601,
  refresh = 100,             
  iter_warmup = 500,       
  iter_sampling = 3000,     
  adapt_delta = 0.95,       
  chains = 4,
  parallel_chains = 4
)

structural_vars <- c("a0", "a1", "a2", "alpha", "log_kappa", "phi1", "phi2")
fit$summary(variables = structural_vars)

## # A tibble: 7 × 10
##   variable      mean   median       sd      mad       q5      q95  rhat ess_bulk
##   <chr>        <dbl>    <dbl>    <dbl>    <dbl>    <dbl>    <dbl> <dbl>    <dbl>
## 1 a0         0.647    0.646   0.0222   0.0213    0.611    0.685    1.11     36.0
## 2 a1         0.00672  0.00671 0.000826 0.000833  0.00537  0.00813  1.22     13.9
## 3 a2         0.0242   0.0235  0.00505  0.00379   0.0178   0.0320   1.22     16.9
## 4 alpha     -0.0196  -0.0196  0.000673 0.000644 -0.0208  -0.0185   1.11     36.0
## 5 log_kappa  8.82     8.81    0.161    0.159     8.57     9.11     1.21     14.4
## 6 phi1       1.19     1.19    0.119    0.120     0.991    1.39     1.26     11.6
## 7 phi2      -0.213   -0.212   0.117    0.118    -0.409   -0.0166   1.26     11.7
## # ℹ 1 more variable: ess_tail <dbl>

Warning: This wind fit has not converged adequately (\(\widehat R\approx1.11\)–\(1.26\) with low bulk effective sample sizes), so its summaries and figures are illustrative only.

fit$diagnostic_summary()

# --- Posterior marginals ---
draws <- fit$draws(variables = structural_vars, format = "df")

plot_density <- function(df, var_name, title_expr) {
  ggplot(df, aes(x = .data[[var_name]])) +
    geom_density(fill = col1_t, alpha = 0.6, color = "white") +
    labs(title = title_expr, x = "", y = "") +
    theme_minimal() + 
    theme(plot.title = element_text(face = "bold", hjust = 0.5))
}

p_a0 <- plot_density(draws, "a0", expression(a[0])) + labs(y = "Density")
p_a1 <- plot_density(draws, "a1", expression(a[1]))
p_a2 <- plot_density(draws, "a2", expression(a[2]))
p_alpha <- plot_density(draws, "alpha", expression(alpha))
p_log_kappa <- plot_density(draws, "log_kappa", expression(log(kappa))) + labs(y = "Density")
p_phi1 <- plot_density(draws, "phi1", expression(phi[1]))
p_phi2 <- plot_density(draws, "phi2", expression(phi[2]))

density_grid <- (p_a0 | p_a1 | p_a2 | p_alpha) /
                (p_log_kappa | p_phi1 | p_phi2 | plot_spacer()) + 
  plot_annotation(
    title = "Posterior Marginal Distributions",
    theme = theme(plot.title = element_text(face = "bold", size = 18, hjust = 0.5))
  )

density_grid

The panels show the baseline and temperature effects, the two latent loadings, log concentration, and AR(2) coefficients.

# --- Latent Processes ---
get_latent_summary <- function(fit, param_name, length) {
  draws_mat <- fit$draws(param_name, format = "matrix")
  data.frame(
    idx = 1:length,
    mean = apply(draws_mat, 2, mean),
    lower = apply(draws_mat, 2, quantile, probs = 0.025),
    upper = apply(draws_mat, 2, quantile, probs = 0.975)
  )
}

w_summary <- get_latent_summary(fit, "w", n)
w2_summary <- get_latent_summary(fit, "w2", 24)
w_summary$time <- dat$datetime

p_w <- ggplot(w_summary, aes(x = time, y = mean)) +
  geom_ribbon(aes(ymin = lower, ymax = upper), fill = col1_t, alpha = 0.3) +
  geom_line(color = col1, linewidth = 0.8) +
  labs(title = expression("AR(2) Process (" * w * ")"), x = "", y = "Value") +
  theme_minimal(base_size = 14) + theme(plot.title = element_text(face = "bold"))

p_w2 <- ggplot(w2_summary, aes(x = idx, y = mean)) +
  geom_ribbon(aes(ymin = lower, ymax = upper), fill = col3_t, alpha = 0.3) +
  geom_line(color = col3, linewidth = 0.8) +
  geom_point(color = col5, size = 2) +
  scale_x_continuous(breaks = seq(0, 24, by = 4)) +
  labs(title = expression("Diurnal RW2 (" * w[2] * ")"), x = "Hour of Day", y = "") +
  theme_minimal(base_size = 14) + theme(plot.title = element_text(face = "bold"))

(p_w | p_w2) +
  plot_annotation(
    title = "Latent Temporal Trajectories",
    subtitle = "Solid lines indicate posterior means; shaded regions indicate 95% credible intervals.",
    theme = theme(plot.title = element_text(face = "bold", size = 18, hjust = 0.5),
                  plot.subtitle = element_text(hjust = 0.5))
  )

The left panel shows the observation-indexed AR(2) trajectory, which represents persistence across the month; the right panel shows the 24 hour-of-day states, which represent a recurring diurnal pattern. Their contributions to the predictor are \(a_1w_i\) and \(a_2w_{2,h(i)}\), so the latent states and their loadings must be interpreted together.

# --- Posterior Predictive Check ---
x_rep <- fit$draws("x_rep", format = "matrix")
set.seed(123)
sample_idx <- sample(nrow(x_rep), 100)
x_rep_sample <- x_rep[sample_idx, ]

ppc_dens_overlay(x_centered, x_rep_sample) +
  scale_x_continuous(
    breaks = c(-pi/8, 0, pi/8),
    labels = c("-π/8", "0", "π/8")
  ) +
  coord_cartesian(xlim = c(-pi/8, pi/8)) +
  labs(
    title = "Posterior Predictive Check",
    subtitle = "LAvM Distribution",
    x = "Radians",
    y = "Density"
  ) +
  theme_minimal(base_size = 14) +
  theme(plot.title = element_text(face = "bold"), legend.position = "none")

The black curve is the centred observed density, and the coloured curves are posterior replications of the same marginal angular feature.

Across the three examples, the same LAvM likelihood carries a real-valued predictor into a circular response model without reintroducing the periodic regression ambiguity. The main modeling choices are to keep the observed arc away from the branch boundary, interpret \(2\arctan(\eta_i)\) as the conditional direction, distinguish baseline \(\kappa\) from observation-specific effective concentration, and constrain latent processes so that they remain identifiable from the fixed effects.