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 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 linear 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.

These equivalent modes make the regression parameters difficult to identify and can lead different initial values to different posterior regions. The repeated extrema below illustrate the source of the ambiguity.

2 Link-adjusted von Mises (LAvM) distribution

The link-adjusted circular construction takes a difference on the linear predictor scale and then maps the residual 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 transformed location is \(g(\eta)=2\arctan(\eta)\), which is also the median on the chosen interval \((-\pi,\pi)\). The Jacobian generally shifts the mode away from this location. When \(\eta=0\), the LAvM density reduces to a vM density centered at zero.

The inverse link is singular at \(|x|=\pi\). The regression examples therefore focus on concentrated data with an empty or nearly empty arc that can be placed at the boundary. Check a rose diagram before fitting and, if necessary, rotate the observations so that their main mass lies near zero. Rotate fitted locations and circular summaries back afterwards. The density is defined more broadly, but a response spread around the full circle requires particular care in choosing and assessing the link and its boundary.

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 transformed location through the inverse tangent link and also changes the 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\), rather than a common vM-scale concentration for every observation when \(\eta_i\ne0\).

2.1 Simulation utilities

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 density below uses a scaled Bessel function, with a large-concentration approximation. The CDF is approximated on a grid and interpolated on the log-odds scale; quantiles and random draws use a monotone spline approximation to its inverse.

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)

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
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\), matching the intended setting for these LAvM regression examples.

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 base model of no circular spread: a point mass at the transformed location, obtained as \(\kappa\rightarrow\infty\). Departure from that base is measured using the mean resultant length of the reference vM kernel,

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

which ranges from \(0\) for a circular-uniform kernel to \(1\) for a point mass. The implementation sets the exponential rate on the distance \(d(\kappa)=\sqrt{1-\rho(\kappa)}\) to \(\lambda=-\log(1-\alpha)/\sqrt{1-U}\), using \(U=\alpha=0.5\) here. For nonzero \(\eta\), \(\rho(\kappa)\) is a reference-kernel summary rather than the resultant length of the full LAvM density.

The PC-prior tutorial derives this prior. The function PC_vm_inf(kappa, lambda) supplies the log-density kernel for its positive-concentration component, with the rate lambda determined by \((U,\alpha)\).

library(cmdstanr)

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

functions {
  real PC_vm_inf(real kappa, real lambda) {
    real log_d_sq;
    real log_numerator;

    if (kappa < 10000) {
      real log_I0 = log_modified_bessel_first_kind(0, kappa);
      real log_Q = log_modified_bessel_first_kind(1, kappa) - log_I0;
      // I2 / I0 = 1 - 2 * Q / kappa, where Q = I1 / I0.
      real numerator = -expm1(2 * log_Q) - exp(log_Q) / kappa;
      log_d_sq = log1m_exp(log_Q);
      log_numerator = log(numerator);
    } else {
      real r = inv(kappa);
      log_d_sq = -log(2.0) - log(kappa) + log1p(0.25 * r);
      log_numerator = -log(2.0) - 2 * log(kappa) + log1p(0.5 * r);
    }

    return log(lambda) + log_numerator - log(2.0) - 0.5 * log_d_sq
           - lambda * exp(0.5 * log_d_sq);
  }

  vector lavm_log_densities(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);
    }

    vector[N] half_x = 0.5 * x;
    vector[N] eta_sq = square(eta);
    vector[N] A = cos(2 * atan(tan(half_x) - eta));
    vector[N] D_minus_1 = eta_sq - eta .* sin(x)
                          - eta_sq .* square(sin(half_x));
    return kappa * A - (log(2 * pi()) + log_I0_val) - log1p(D_minus_1);
  }

  real lavm_lpdf(vector x, vector eta, real kappa) {
    return sum(lavm_log_densities(x, eta, kappa));
  }

}

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 {
  real lambda = -log1m(alpha) / sqrt(1.0 - U);
}

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

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

model {
  beta0 ~ normal(0, 5);
  beta1 ~ normal(0, 5);

  target += PC_vm_inf(kappa, lambda);

  y ~ lavm(eta, kappa);
}

mod <- cmdstan_model(sim1_stan_file)

lavm_lpdf combines the von Mises kernel with the Jacobian of the link adjustment. lavm_log_densities gives the observation-level contributions used for predictive assessment.

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.0118 0.0118  0.0875  0.126 1.00     4054.    3015.
## 2 beta1     1.26   1.26  0.113  0.107   1.08    1.45  1.00     3462.    2952.
## 3 kappa    19.9   19.8   2.80   2.87   15.5    24.6   1.000    3607.    2848.

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)

The posterior marginals show parameter uncertainty; dashed lines mark the generating values.

plot_density <- function(df, var_name, title_expr, true_value = NULL,
                         fill = col1_t, alpha = 0.8, y_label = "") {
  p <- ggplot(df, aes(x = .data[[var_name]])) +
    geom_density(fill = fill, alpha = alpha, color = "white") +
    labs(title = title_expr, x = "", y = y_label) +
    theme_minimal() +
    theme(plot.title = element_text(face = "bold", hjust = 0.5))
  if (!is.null(true_value)) {
    p <- p + geom_vline(xintercept = true_value, linetype = "dashed",
                        color = "black", linewidth = 0.6)
  }
  p
}

draws <- fit$draws(variables = c("beta0", "beta1", "kappa"), format = "df")

p_beta0 <- plot_density(draws, "beta0", expression(beta[0]), 0.1, col1,
                        y_label = "Posterior Density")
p_beta1 <- plot_density(draws, "beta1", expression(beta[1]), 1.2, col2)
p_kappa <- plot_density(draws, "kappa", expression(kappa), 20, col3)

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

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\) follows a second-order random walk (RW2), which penalizes second differences to encourage a smooth trajectory.

The simulated path uses a slightly regularized RW2 precision, with its constant and linear components removed and its standard deviation set to \(0.4\). The fitted path has a zero-sum constraint; its linear component is governed by the initial-value priors below.

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)

w <- as.numeric(scale(w_proj)) * sigma_w

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.

The priors are \(N(0,5^2)\) for the fixed effects and a PC prior with \(U=\alpha=0.5\) for concentration. The first two latent values have \(N(0,10^2)\) prior terms. The positive RW2 innovation scale sigma has a half-normal prior with scale 1 and controls the variability of second differences. sigma_w = sd(w) measures the empirical spread of the path.

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

functions {
  real PC_vm_inf(real kappa, real lambda) {
    real log_d_sq;
    real log_numerator;

    if (kappa < 10000) {
      real log_I0 = log_modified_bessel_first_kind(0, kappa);
      real log_Q = log_modified_bessel_first_kind(1, kappa) - log_I0;
      // I2 / I0 = 1 - 2 * Q / kappa, where Q = I1 / I0.
      real numerator = -expm1(2 * log_Q) - exp(log_Q) / kappa;
      log_d_sq = log1m_exp(log_Q);
      log_numerator = log(numerator);
    } else {
      real r = inv(kappa);
      log_d_sq = -log(2.0) - log(kappa) + log1p(0.25 * r);
      log_numerator = -log(2.0) - 2 * log(kappa) + log1p(0.5 * r);
    }

    return log(lambda) + log_numerator - log(2.0) - 0.5 * log_d_sq
           - lambda * exp(0.5 * log_d_sq);
  }

  vector lavm_log_densities(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);
    }

    vector[N] half_x = 0.5 * x;
    vector[N] eta_sq = square(eta);
    vector[N] A = cos(2 * atan(tan(half_x) - eta));
    vector[N] D_minus_1 = eta_sq - eta .* sin(x)
                          - eta_sq .* square(sin(half_x));
    return kappa * A - (log(2 * pi()) + log_I0_val) - log1p(D_minus_1);
  }

  real lavm_lpdf(vector x, vector eta, real kappa) {
    return sum(lavm_log_densities(x, eta, kappa));
  }


  array[] real lavm_rng(vector eta, real kappa) {
    int N = num_elements(eta);
    int G = 400;
    real step_size = 2 * pi() / G;
    real start_val = -pi() + 0.5 * step_size;
    vector[G] x = start_val + linspaced_vector(G, 0, G - 1) * step_size;
    vector[G] tan_half_x = tan(0.5 * x);
    vector[G] sin_x = sin(x);
    vector[G] sin_half_x_sq = square(sin(0.5 * x));
    array[N] real draws;

    for (n in 1:N) {
      real eta_sq = square(eta[n]);
      vector[G] A = cos(2 * atan(tan_half_x - eta[n]));
      vector[G] D_minus_1 = eta_sq - eta[n] * sin_x
                            - eta_sq * sin_half_x_sq;
      vector[G] log_probs = kappa * A - log1p(D_minus_1);
      draws[n] = x[categorical_logit_rng(log_probs)];
    }
    return draws;
  }
}

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

transformed data {
  real lambda = -log1m(alpha) / sqrt(1.0 - U);
}

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

  sum_to_zero_vector[N] w;
  real<lower=0> sigma;    // RW2 second-difference scale
}

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

model {
  target += PC_vm_inf(kappa, lambda);

  beta0 ~ normal(0, 5);
  beta1 ~ normal(0, 5);

  sigma ~ std_normal();

  w[1:2] ~ normal(0, 10);

  w[3:N] ~ normal(2 * w[2:(N - 1)] - w[1:(N - 2)], sigma);

  y ~ lavm(eta, kappa);
}

generated quantities {
  vector[N] log_lik = lavm_log_densities(y, eta, kappa);
  array[N] real y_rep = lavm_rng(eta, kappa);

  real sigma_w = sd(w);
  real kappa_eff;

  { // Local scope to prevent saving intermediate variables to MCMC draws

    // Gaussian approximation to the effective concentration
    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);

  }
}

mod <- cmdstan_model(sim2_stan_file)

stan_data <- list(
  N = n,
  y = as.numeric(y),
  z = z,
  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", "sigma_w"))

## # A tibble: 6 × 10
##   variable      mean   median       sd      mad       q5      q95  rhat ess_bulk
##   <chr>        <dbl>    <dbl>    <dbl>    <dbl>    <dbl>    <dbl> <dbl>    <dbl>
## 1 beta0      0.104    0.104    0.00639  0.00645  0.0936   1.15e-1 1.000    4800.
## 2 beta1      0.194    0.194    0.00773  0.00776  0.181    2.06e-1 1.00     4616.
## 3 kappa     66.2     65.8      9.96     9.76    50.8      8.39e+1 1.00     4720.
## 4 kappa_eff 98.2     97.5     14.8     14.4     75.3      1.25e+2 1.00     4751.
## 5 sigma      0.00573  0.00541  0.00172  0.00140  0.00357  8.84e-3 1.03      167.
## 6 sigma_w    0.403    0.403    0.00647  0.00635  0.392    4.13e-1 1.00     2898.
## # ℹ 1 more variable: ess_tail <dbl>

The posterior means are approximately \(0.104\) for \(\beta_0\), \(0.194\) for \(\beta_1\), \(98.2\) for the average effective concentration, and \(0.403\) 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.2\); its different scale is explained below.

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

p_beta0 <- plot_density(draws, "beta0", expression(beta[0]), beta0, col1_t,
                        y_label = "Density")
p_beta1 <- plot_density(draws, "beta1", expression(beta[1]), beta1, col2_t)
p_sigma_w <- plot_density(draws, "sigma_w", expression(sigma[w]), sigma_w, col4_t)
p_kappa <- plot_density(draws, "kappa", expression(kappa), kappa_true, col3_t)
p_kappa_eff <- plot_density(draws, "kappa_eff", expression(kappa[eff]), kappa_true,
                            col5_t, y_label = "Density")

(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))
  )

The LAvM parameter \(\kappa\) describes concentration at \(\eta=0\). For comparison with the generating von Mises distribution, we report a curvature-based effective concentration evaluated at the transformed location \(x_i=2\arctan(\eta_i)\):

\[ \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 \]

The leading term is \(\kappa(1+\eta_i^2)^2\) at high concentration; the additional term comes from the Jacobian. Because the transformed location is generally not the mode, this quantity is a local comparison, not an exact vM concentration or inverse circular variance. 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.

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
)

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)

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))

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))

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))

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, reliability of the PSIS-LOO approximation through Pareto-\(k\), and predictive calibration through LOO-PIT.

2.4 Wind application example

The final example uses January hourly climate normals from John F. Kennedy International Airport. These are summaries for each calendar hour, including mean wind-vector direction and temperature, rather than observations from one January; see the NOAA hourly normals documentation. The prepared data express direction in radians and attach a calendar date for plotting. The accompanying article develops a joint model for wind direction and speed; here we fit the circular direction component.

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

Rotation places the sample circular mean at zero and wraps the observations to \([-\pi,\pi)\), moving the occupied arc away from the LAvM singularity at the endpoints. Use rotation to translate fitted or predicted directions 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 = 0, direction = 1) +
  scale_x_continuous(limits = c(0, 2 * pi), 
                     breaks = c(0, pi / 2, pi, 3 * pi / 2), 
                     labels = c("N", "E", "S", "W")) +
  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("0", "π/2", "−π", "−π/2")) +
  labs(x = "direction relative\nto circular mean", 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 normal", 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 show the angular concentration before and after rotation. The temperature histogram shows the scale and support of the covariate, while the calendar-hour trace shows smooth variation that motivates a temporal latent component. These patterns describe climate normals, not the progression of a particular weather event. 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} &= \widetilde{\mathbf w}-\overline{\widetilde w}\,\mathbf 1, \qquad \mathbf{w}_{2} \sim N\left( 0, \mathbf{Q}_{\text{cyclic 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 at zero temperature, \(\alpha\) is the temperature coefficient, and \(a_1\) and \(a_2\) scale the calendar-time 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 transformed location 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 \(U=0.5\) and \(\alpha=0.99\) in the distance-rate formula above. The exponential priors on the positive loadings \(a_1\) and \(a_2\) have median \(0.5\) and regularize the sizes of the two latent contributions.

The calendar-time path follows an AR(2) model before centering: \[ \widetilde w_i = \phi_1\widetilde w_{i-1}+\phi_2\widetilde w_{i-2}+\varepsilon_i, \qquad \varepsilon_i\sim N(0,1), \quad i\geq3. \] Its first two states have independent \(N(0,5^2)\) priors, and \(\mathcal S_{\mathrm{AR2}}\) restricts the coefficients to \(-1<\phi_2<1\) and \(\phi_2-1<\phi_1<1-\phi_2\). Thus \(a_1\) scales the innovations; the initial states are specified separately from the stationary distribution. Subtracting the path mean separates its level from the intercept.

The hourly effects have a cyclic RW2 prior on the zero-sum subspace, \[ p(\mathbf w_2)\propto\exp\left\{-\tfrac12\sum_h(w_{2,h+1}-2w_{2,h}+w_{2,h-1})^2\right\}, \] with hour indices wrapping across midnight. The inverse precision in the model above is restricted to this subspace.

The implementation integrates Gaussian coefficients and hourly effects analytically, then recovers their conditional posterior draws for summaries.

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

functions {
  // Observation likelihood plus the residual transformation Jacobian.
  // Combining the kappa terms here avoids subtracting nearly equal large values.
  real lavm_scaled_residual_lpdf(vector x_residual, data vector x, real log_kappa) {
    int N = num_elements(x);
    real kappa = exp(log_kappa);
    real log_I0_scaled;
    vector[N] residual_sq = square(x_residual);
    vector[N] delta_sq = exp(-log_kappa) * residual_sq;
    if (kappa < 10000) {
      log_I0_scaled = log_modified_bessel_first_kind(0, kappa) - kappa
                     + 0.5 * log_kappa;
    } else {
      log_I0_scaled = -0.5 * log(2 * pi()) + 0.125 * exp(-log_kappa);
    }
    return -2 * sum(residual_sq ./ (1 + delta_sq))
           - N * (log(2 * pi()) + log_I0_scaled)
           - 2 * sum(log(abs(cos(0.5 * x)))) - sum(log1p(delta_sq));
  }

  real PC_vm_inf(real kappa, real lambda) {
    real log_d_sq;
    real log_numerator;

    if (kappa < 10000) {
      real log_I0 = log_modified_bessel_first_kind(0, kappa);
      real log_Q = log_modified_bessel_first_kind(1, kappa) - log_I0;
      // I2 / I0 = 1 - 2 * Q / kappa, where Q = I1 / I0.
      real numerator = -expm1(2 * log_Q) - exp(log_Q) / kappa;
      log_d_sq = log1m_exp(log_Q);
      log_numerator = log(numerator);
    } else {
      real r = inv(kappa);
      log_d_sq = -log(2.0) - log(kappa) + log1p(0.25 * r);
      log_numerator = -log(2.0) - 2 * log(kappa) + log1p(0.5 * r);
    }

    return log(lambda) + log_numerator - log(2.0) - 0.5 * log_d_sq
           - lambda * exp(0.5 * log_d_sq);
  }

  array[] real lavm_rng(vector eta, real kappa) {
    int N = num_elements(eta);
    int G = 400;
    real step_size = 2 * pi() / G;
    real start_val = -pi() + 0.5 * step_size;
    vector[G] x = start_val + linspaced_vector(G, 0, G - 1) * step_size;
    vector[G] tan_half_x = tan(0.5 * x);
    vector[G] sin_x = sin(x);
    vector[G] sin_half_x_sq = square(sin(0.5 * x));
    array[N] real draws;

    for (n in 1:N) {
      real eta_sq = square(eta[n]);
      vector[G] A = cos(2 * atan(tan_half_x - eta[n]));
      vector[G] D_minus_1 = eta_sq - eta[n] * sin_x
                            - eta_sq * sin_half_x_sq;
      vector[G] log_probs = kappa * A - log1p(D_minus_1);
      draws[n] = x[categorical_logit_rng(log_probs)];
    }
    return draws;
  }
  matrix cyclic_rw2_basis(int H) {
    matrix[H, H - 1] basis = rep_matrix(0, H, H - 1);
    for (j in 1:(H - 1)) {
      basis[1:j, j] = rep_vector(inv_sqrt(j * (j + 1.0)), j);
      basis[j + 1, j] = -sqrt(j / (j + 1.0));
    }
    matrix[H, H - 1] d2 = append_row(basis[2:H], basis[1]) - 2 * basis
                           + append_row(basis[H], basis[1:(H - 1)]);
    matrix[H - 1, H - 1] L = cholesky_decompose(crossprod(d2));
    return mdivide_left_tri_low(L, basis')';
  }

  array[] matrix ar2_static_grams(matrix Xc) {
    int N = rows(Xc);
    int P = cols(Xc);
    array[7] matrix[P, P] grams;
    matrix[P, P] g01 = Xc[3:N]' * Xc[2:(N - 1)];
    matrix[P, P] g02 = Xc[3:N]' * Xc[1:(N - 2)];
    matrix[P, P] g12 = Xc[2:(N - 1)]' * Xc[1:(N - 2)];
    grams[1] = crossprod(Xc[1:2]) / 25;
    grams[2] = crossprod(Xc[3:N]);
    grams[3] = crossprod(Xc[2:(N - 1)]);
    grams[4] = crossprod(Xc[1:(N - 2)]);
    grams[5] = g01 + g01';
    grams[6] = g02 + g02';
    grams[7] = g12 + g12';
    return grams;
  }

  matrix ar2_static_sums(matrix Xc) {
    int N = rows(Xc);
    int P = cols(Xc);
    matrix[4, P] sums;
    sums[1] = (Xc[1] + Xc[2]) / 25;
    sums[2] = rep_row_vector(1, N - 2) * Xc[3:N];
    sums[3] = rep_row_vector(1, N - 2) * Xc[2:(N - 1)];
    sums[4] = rep_row_vector(1, N - 2) * Xc[1:(N - 2)];
    return sums;
  }

  vector ar2_static_standardize(vector path, real phi1, real phi2, real sigma) {
    int N = num_elements(path);
    vector[N] innovations;
    innovations[1:2] = path[1:2] / 5;
    innovations[3:N] = (path[3:N] - phi1 * path[2:(N - 1)]
                        - phi2 * path[1:(N - 2)]) / sigma;
    return innovations;
  }

  vector ar2_static_cross(data matrix Xc, vector innovations,
                          real phi1, real phi2, real sigma) {
    int N = rows(Xc);
    return Xc[1:2]' * (innovations[1:2] / 5)
           + (Xc[3:N]' * innovations[3:N]
              - phi1 * (Xc[2:(N - 1)]' * innovations[3:N])
              - phi2 * (Xc[1:(N - 2)]' * innovations[3:N])) / sigma;
  }

  tuple(real, vector, matrix, vector) gaussian_regression_distribution(
      vector predictor, data matrix Xc, data row_vector mean_X,
      data array[] matrix grams, data matrix sums, vector scales,
      real loading, real phi1, real phi2, real sigma) {
    int N = num_elements(predictor);
    int P = cols(Xc);
    real mean_predictor = mean(predictor);
    real ar_constant = 1 - phi1 - phi2;
    real precision_mean = 2.0 / 25 + (N - 2) * square(ar_constant / sigma);
    vector[P] mean_cross = (sums[1] + ar_constant / square(sigma)
                              * (sums[2] - phi1 * sums[3] - phi2 * sums[4]))';
    matrix[P, P] gram_raw = grams[1]
      + (grams[2] + square(phi1) * grams[3] + square(phi2) * grams[4]
         - phi1 * grams[5] - phi2 * grams[6] + phi1 * phi2 * grams[7])
        / square(sigma);
    matrix[P, P] gram_projected = gram_raw - mean_cross * mean_cross' / precision_mean;
    vector[P] scaled_columns = scales / loading;
    row_vector[P] mean_design = mean_X .* scales';
    matrix[P, P] precision = quad_form_diag(gram_projected, scaled_columns)
                            + mean_design' * mean_design
                            + diag_matrix(rep_vector(1, P));
    matrix[P, P] L = cholesky_decompose(precision);
    vector[N] r = ar2_static_standardize((predictor - mean_predictor) / loading,
                                          phi1, phi2, sigma);
    real mean_response_cross = sum(r[1:2]) / 5 + ar_constant / sigma * sum(r[3:N]);
    vector[P] rhs = mean_design' * mean_predictor
      + scaled_columns .* (ar2_static_cross(Xc, r, phi1, phi2, sigma)
                           - mean_cross * (mean_response_cross / precision_mean));
    vector[P] conditional_mean = mdivide_right_tri_low(mdivide_left_tri_low(L, rhs)', L)';
    vector[N] residual = ar2_static_standardize(
      (predictor - mean_predictor - Xc * (scales .* conditional_mean)) / loading,
      phi1, phi2, sigma);
    real residual_mean = (sum(residual[1:2]) / 5
                         + ar_constant / sigma * sum(residual[3:N])) / precision_mean;
    residual[1:2] -= rep_vector(residual_mean / 5, 2);
    residual[3:N] -= rep_vector(residual_mean * ar_constant / sigma, N - 2);
    real residual_ss = dot_self(residual)
                      + square(mean_predictor - mean_design * conditional_mean)
                      + dot_self(conditional_mean);
    real log_density = -0.5 * N * log(2 * pi()) - log(25)
                      - (N - 2) * log(sigma) - (N - 1) * log(loading)
                      - 0.5 * (log(precision_mean) + residual_ss)
                      - sum(log(diagonal(L)));
    vector[N] score = rep_vector(0, N);
    score[1:2] = residual[1:2] / 5;
    score[3:N] += residual[3:N] / sigma;
    score[2:(N - 1)] -= phi1 * residual[3:N] / sigma;
    score[1:(N - 2)] -= phi2 * residual[3:N] / sigma;
    score = -(score - mean(score)) / loading
      - rep_vector((mean_predictor - mean_design * conditional_mean) / N, N);
    return (log_density, conditional_mean, L, score);
  }

  tuple(vector, real) gaussian_predictor_coordinates(
      vector coordinates, data vector location, data matrix Xc,
      data row_vector mean_X, data array[] matrix grams, data matrix sums,
      vector scales, real loading, real phi1, real phi2, real log_kappa) {
    int N = num_elements(location);
    vector[N] score = gaussian_regression_distribution(
      location, Xc, mean_X, grams, sums, scales, loading, phi1, phi2, 1).4;
    vector[N] precision_diag = rep_vector(0, N);
    precision_diag[1:2] = rep_vector(1.0 / 25, 2);
    precision_diag[3:N] += rep_vector(1, N - 2);
    precision_diag[2:(N - 1)] += rep_vector(square(phi1), N - 2);
    precision_diag[1:(N - 2)] += rep_vector(square(phi2), N - 2);
    precision_diag /= square(loading);
    vector[N] denominator = 4 + exp(-log_kappa) * (precision_diag + 2);
    vector[N] residual = coordinates ./ sqrt(denominator)
      + exp(-0.5 * log_kappa) * (score ./ denominator);
    return (residual, -0.5 * sum(log(denominator)));
  }

  real gaussian_regression_lpdf(vector predictor, data matrix Xc, data row_vector mean_X,
                       data array[] matrix grams, data matrix sums, vector scales,
                       real loading, real phi1, real phi2, real sigma) {
    return gaussian_regression_distribution(predictor, Xc, mean_X, grams, sums, scales,
                                   loading, phi1, phi2, sigma).1;
  }

  vector gaussian_regression_coefficients_rng(vector predictor, data matrix Xc,
      data row_vector mean_X, data array[] matrix grams, data matrix sums,
      vector scales, real loading, real phi1, real phi2, real sigma) {
    int P = cols(Xc);
    tuple(real, vector[P], matrix[P, P], vector[num_elements(predictor)]) conditional = gaussian_regression_distribution(
      predictor, Xc, mean_X, grams, sums, scales, loading, phi1, phi2, sigma);
    vector[P] noise = to_vector(normal_rng(rep_vector(0, P), 1));
    return conditional.2 + mdivide_right_tri_low(noise', conditional.3)';
  }
}

data {
  int<lower=3> 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 = -log1m(alpha_pc) / sqrt(1 - U);
  real lambda_tau = -log(0.5) / 0.5;
  vector[N] x_location = tan(0.5 * x);
  matrix[24, 23] hour_basis = cyclic_rw2_basis(24);
  matrix[N, 23] hour_design = hour_basis[hour_idx];
  matrix[N, 24] fixed_design = append_col(z, hour_design);
  row_vector[24] mean_X = rep_row_vector(1.0 / N, N) * fixed_design;
  matrix[N, 24] Xc = fixed_design - rep_matrix(mean_X, N);
  array[7] matrix[24, 24] static_grams = ar2_static_grams(Xc);
  matrix[4, 24] static_sums = ar2_static_sums(Xc);
}

parameters {
  vector[N] x_residual;
  real<lower=0> a1;
  real<lower=0> a2;
  real log_kappa;
  real<lower=-1, upper=1> phi2;
  real<lower=phi2-1, upper=1-phi2> phi1;
}

transformed parameters {
  real<lower=0> kappa = exp(log_kappa);
  vector[N] eta_x;
  vector[N] effective_residual;
  real residual_log_jacobian;
  {
    vector[24] scales = rep_vector(a2, 24);
    scales[1] = 1;
    tuple(vector[N], real) coordinates = gaussian_predictor_coordinates(
      x_residual, x_location, Xc, mean_X, static_grams, static_sums,
      scales, a1, phi1, phi2, log_kappa);
    effective_residual = coordinates.1;
    residual_log_jacobian = coordinates.2;
    eta_x = x_location + exp(-0.5 * log_kappa) * effective_residual;
  }
}

model {
  target += log_kappa;
  target += PC_vm_inf(kappa, lambda);
  a1 ~ exponential(lambda_tau);
  a2 ~ exponential(lambda_tau);
  vector[24] scales = rep_vector(a2, 24);
  scales[1] = 1;
  target += gaussian_regression_lpdf(eta_x | Xc, mean_X, static_grams, static_sums,
                          scales, a1, phi1, phi2, 1);
  target += lavm_scaled_residual_lpdf(effective_residual | x, log_kappa);
  target += residual_log_jacobian;
}

generated quantities {
  real alpha;
  real a0;
  vector[N] w;
  vector[24] w2;
  {
    vector[24] scales = rep_vector(a2, 24);
    scales[1] = 1;
    vector[24] direction_coef = gaussian_regression_coefficients_rng(
      eta_x, Xc, mean_X, static_grams, static_sums, scales, a1, phi1, phi2, 1);
    vector[N] fitted_effect = fixed_design * (scales .* direction_coef);
    alpha = direction_coef[1];
    w2 = hour_basis * direction_coef[2:24];
    a0 = mean(eta_x - fitted_effect);
    w = (eta_x - a0 - fitted_effect) / a1;
  }
  array[N] real x_rep = lavm_rng(eta_x, kappa);
}

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

Scaled predictor coordinates keep the circular likelihood evaluation stable at large concentration. Posterior predictive angles use a 400-point grid approximation to the continuous LAvM distribution.

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(
    x_residual = rnorm(n, sd = 0.1),
    a1 = runif(1, min = 0.05, max = 0.2),
    a2 = runif(1, min = 0.05, max = 0.2),
    log_kappa = log(runif(1, min = 3, max = 7)),
    phi2 = runif(1, min = -0.1, max = 0.1),
    phi1 = runif(1, min = 0.3, max = 0.7)
  )
}

output_dir <- file.path("results", "wind_csv")
dir.create(output_dir, recursive = TRUE, showWarnings = FALSE)

set.seed(1314)
fit <- mod$sample(
  data = stan_data,
  output_dir = output_dir,
  init = init_fun,
  seed = 1314,
  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.412   0.414  0.100    0.101     0.242    0.571    1.00    8306.
## 2 a1         0.0103  0.0104 0.000591 0.000395  0.00903  0.0110   1.03     246.
## 3 a2         0.0223  0.0219 0.00379  0.00356   0.0171   0.0293   1.01     812.
## 4 alpha     -0.0125 -0.0126 0.00303  0.00305  -0.0173  -0.00734  1.00    8292.
## 5 log_kappa 12.2    11.8    2.25     1.82      9.49    16.5      1.02     252.
## 6 phi1       0.746   0.736  0.0642   0.0471    0.667    0.872    1.02     196.
## 7 phi2       0.197   0.207  0.0640   0.0480    0.0762   0.277    1.02     187.
## # ℹ 1 more variable: ess_tail <dbl>

These examples illustrate circular regression; they are not presented as the best models for these data. Practical analyses may require different model choices, more warmup and sampling, or further tuning. Detailed long-term time-series modelling is beyond this tutorial’s scope.

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

p_a0 <- plot_density(draws, "a0", expression(a[0]), alpha = 0.6) + labs(y = "Density")
p_a1 <- plot_density(draws, "a1", expression(a[1]), alpha = 0.6)
p_a2 <- plot_density(draws, "a2", expression(a[2]), alpha = 0.6)
p_alpha <- plot_density(draws, "alpha", expression(alpha), alpha = 0.6) +
  scale_x_continuous(n.breaks = 3, guide = guide_axis(check.overlap = TRUE))
p_log_kappa <- plot_density(draws, "log_kappa", expression(log(kappa)), alpha = 0.6) + labs(y = "Density")
p_phi1 <- plot_density(draws, "phi1", expression(phi[1]), alpha = 0.6)
p_phi2 <- plot_density(draws, "phi2", expression(phi[2]), alpha = 0.6)

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.

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
w2_summary$hour <- w2_summary$idx - 1

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 = hour, 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 = c(0, 6, 12, 18, 23)) +
  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 centered AR(2) trajectory across January calendar hours; the right panel shows the 24 hour-of-day states, which represent a recurring diurnal pattern in the normals. 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.

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.

3 Conclusion

Link-adjusted circular distributions form a regression family particularly suited to reasonably concentrated circular data. The LAvM examples show how a linear predictor combines covariates and latent effects, with the transformed location and concentration interpreted together.