This tutorial shows how a circular variable can enter a regression through a joint circular–linear model rather than as a raw angular covariate. We begin with an identifiable fixed-effect model, add shared and response-specific AR(2) processes, and finish with a joint analysis of hourly wind direction and speed. Familiarity with a basic Stan workflow is assumed.

After completing the tutorial, you will be able to:

represent a circular covariate through its LAvM linear predictor;
construct an identifiable joint model for circular and linear responses; and
share structured temporal effects between responses and interpret their association.

1 Preliminary

The joint framework follows Xiang Ye, Janet Van Niekerk, and Håvard Rue (2026), A Bayesian regression framework for circular models with INLA.

1.1 Simulation utilities

The simulations can use either INLAcircular or the standalone von Mises functions below.

Install INLAcircular

options(
  repos = c(
    getOption("repos"), 
    INLA = "https://inla.r-inla-download.org/R/testing"
    )
  )
if (!requireNamespace("remotes", quietly = TRUE)) {
  install.packages("remotes")
}
remotes::install_github("XiangYEstats/INLAcircular")

The INLAcircular package provides the LAvM distribution functions used here.

Standalone von Mises functions

## Functions for von Mises distribution
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)
}

Both versions of rvm() accept an observation-specific mean direction.

1.2 Link-adjusted von Mises (LAvM) distribution

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

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)
col6 <- "#998B9F"; col6_t <- adjustcolor(col6, alpha.f = 0.6)
col7 <- "#A3B899"; col7_t <- adjustcolor(col7, alpha.f = 0.6)
col8 <- "#D4A373"; col8_t <- adjustcolor(col8, alpha.f = 0.6)

For the inverse-tangent link \(g(\eta)=2\arctan(\eta)\), the link-adjusted von Mises (LAvM) distribution has density \[ 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) \] where \(\eta\) is the linear predictor and \(I_0\) is the modified Bessel function of the first kind. The modal direction is \(g(\eta)\), while \(\kappa\) controls conditional concentration.

The link adjustment compares \(g^{-1}(x)\) with \(\eta\) on the real line and then maps the difference back to the circle. This keeps an unrestricted regression predictor compatible with a circular response, without treating angular and linear quantities as if they belonged to the same space. Notice that \(\eta\) is not itself an angle: regression coefficients act linearly on the \(\eta\) scale, whereas their effect on the modal direction is nonlinear and bounded after applying \(g\).

When \(\eta=0\), the LAvM density is the corresponding von Mises density centred at zero. Moving \(\eta\) shifts the mode through \(g(\eta)\) and also changes the local concentration; increasing \(\kappa\) sharpens the density in both panels. This predictor-dependent shape is why the fitted LAvM \(\kappa\) is not always directly comparable with a von Mises concentration.

2 Circular joint regression framework

Suppose \(y\) is Gaussian, \(z\) is a linear covariate, and \(x\) is a circular variable that we want to use as a covariate. Directly inserting the numeric representation of \(x\) into a linear predictor depends on the angular origin and branch cut. Angles just below \(2\pi\) and just above zero, for example, are close on the circle but far apart numerically, so we avoid \[ y_{i} \sim N(\beta_{0} + \beta_{1}z_{i} + \beta_{2}x_{i}, \sigma^{2}) \] An origin-invariant first-harmonic alternative is \[ y_{i} \sim N(\beta_{0} + \beta_{1}z_{i} + \beta_{2}\cos\left(x_{i}\right) + \beta_{3}\sin\left(x_{i}\right), \sigma^{2}) \] This sine–cosine representation is useful when a first-harmonic association is scientifically adequate. It is less convenient when the systematic structure underlying \(x\) contains temporal, spatial, or other random effects that should also inform \(y\). We then treat \(x\) as a second response: \[ \begin{align} y_{i} &\sim N(\beta_{0} + \beta_{1}z_{i} + b_{1}\eta^{x}_{i}, \sigma^{2}) \\ x_{i} &\sim \text{LAvM}(\eta^{x}_{i}, \kappa) \end{align} \] Here \(\eta_i^x\) is the systematic component of the circular response on the linear scale. The observed angle is a noisy circular measurement of that component, and \(b_1\) carries the component into the Gaussian mean. Because both models are fitted jointly, uncertainty about \(\eta_i^x\) is propagated to \(b_1\) and to predictions of \(y_i\).

The interpretation is therefore different from an ordinary angular slope. A positive \(b_1\) means that larger values of the linear predictor for \(x_i\) are associated with larger expected \(y_i\); it is not a constant change in \(y_i\) per radian of observed angle. Angular predictions are obtained only after applying \(g(\eta)=2\arctan(\eta)\).

2.1 Simulation example 1 (joint linear covariate)

We first consider an identifiable example with two independent linear covariates: \(z_i^x\) drives the circular predictor and \(z_i^y\) enters the linear response directly: \[ \begin{align} y_{i} &\sim N\left(\beta_{0} + \beta_{1}z_{i}^{y} + b_{1}\left(\alpha_{0} + \alpha_{1}z_{i}^{x}\right), \sigma^{2}\right), \\ x_{i} &\sim \text{LAvM}\left(\alpha_{0} + \alpha_{1}z_{i}^{x}, \kappa\right). \end{align} \] If the same covariate were used in both places, its direct contribution \(\beta_1z_i\) could not be separated from the mediated contribution \(b_1\alpha_1z_i\) without further constraints or prior information. Using independent \(z_i^x\) and \(z_i^y\) makes the teaching example deliberately identifiable: \(\alpha_1\) controls the circular predictor, \(\beta_1\) is the direct effect on \(y_i\), and \(b_1\) transfers information from the circular submodel to the Gaussian submodel.

library(circular)

We simulate \(n=100\) observations with high von Mises concentration. This example is designed for parameter recovery: the LAvM fit should recover the regression effects, although its raw \(\kappa\) needs a scale adjustment before it can be compared with the generating von Mises concentration.

lfun <- function(x) {2 * atan(x)}
n <- 100

beta0 <- 2; beta1 <- 1.2; b1 <- 0.7
alpha0 <- 0.1; alpha1 <- 0.5
kappa <- 200
sd_y <- 0.5

set.seed(1314)
z_x <- rnorm(n, mean = 0, sd = 1)
z_y <- rnorm(n, mean = 0, sd = 1)
eta_x <- alpha0 + alpha1 * z_x
# x ~ vM(g(alpha0 + alpha1*z_x), kappa)
# y ~ N(beta0 + beta1*z_y + b1*eta_x, sd_y^2)
x <- rvm(n, mu = lfun(eta_x), kappa = kappa)
y <- rnorm(n, mean = beta0 + beta1 * z_y + b1 * eta_x, sd = sd_y)

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

For the concentration parameter \(\kappa\), we use a penalised-complexity (PC) prior. A PC prior shrinks toward a simpler reference model and controls departures from it through an interpretable probability statement. Here the reference is the limiting point mass, corresponding to perfect concentration: \[ 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)}}}, \] where \(\mathcal{I}_j\) is the modified Bessel function of order \(j\). The rate \(\lambda\) is calibrated through \(\Pr(R>U)=\alpha\), where the mean resultant length \(R\in[0,1]\) measures circular concentration: values near zero indicate weak concentration and values near one indicate angles clustered around a common direction. We use \(U=0.5\) and \(\alpha=0.5\) for this illustration; in an application, these values should express plausible prior knowledge about concentration.

library(cmdstanr)

Stan evaluates the LAvM density directly. The model samples log_kappa for more stable geometry and adds log_kappa as the Jacobian because the PC prior is defined on the natural \(\kappa\) scale. Asymptotic Bessel expressions avoid overflow at large \(\kappa\). The lavm_rng function discretises the circle into 400 cells, which is sufficient for these illustrative predictive plots but should be checked if fine angular resolution matters.

stan_file <- file.path("stan", "joint_simulation_1.stan")

functions {
  real lavm_lpdf(vector x, vector eta, real kappa) {
    int N = num_elements(x);
    real log_lik = 0;
    
    // Stable log(I0) calculation with asymptotic safeguard
    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;
  }
  
  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 = 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));
      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] x;
  vector[N] z_x;
  vector[N] z_y;
  real<lower=0, upper=1> U;
  real<lower=0, upper=1> alpha;
}

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

parameters {
  real beta0;
  real beta1;
  real b1;
  
  real alpha0;
  real alpha1;
  
  // Sampled on unconstrained log scale for optimal HMC geometry
  real log_kappa; 
  real<lower=0> sd_y;  
}

transformed parameters {
  // Transformed to natural scale for likelihood and outputs
  real<lower=0> kappa = exp(log_kappa);
  
  vector[N] eta_x = alpha0 + alpha1 * z_x;
  vector[N] mu_y = beta0 + beta1 * z_y + b1 * eta_x;
}

model {
  // Jacobian adjustment for log_kappa
  target += log_kappa;
  
  // PC Prior Implementation (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);
  
  // Fixed Effects Priors
  beta0 ~ normal(0, 5);
  beta1 ~ normal(0, 5);
  b1 ~ normal(0, 5);
  alpha0 ~ normal(0, 5);
  alpha1 ~ normal(0, 5);
  sd_y ~ std_normal();
  
  // Joint model
  y ~ normal(mu_y, sd_y);
  x ~ lavm(eta_x, kappa);
}

generated quantities {
  vector[N] log_lik; 
  array[N] real y_rep;
  array[N] real x_rep;
  
  // Returned directly on natural scale
  real kappa_eff; 
  
  { 
    vector[N] eta_sq = square(eta_x);
    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);
  
    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) {
      real ll_y = normal_lpdf(y[n] | mu_y[n], sd_y);
      
      real tan_half_x = tan(x[n] / 2.0);
      real A = cos(2.0 * atan(tan_half_x - eta_x[n]));
      real D = 1.0 + square(eta_x[n]) - eta_x[n] * sin(x[n]) - square(eta_x[n]) * square(sin(x[n] / 2.0));
      real ll_x = (kappa * A) - log_bessel - log(D);
      
      log_lik[n] = ll_y + ll_x;
      
      y_rep[n] = normal_rng(mu_y[n], sd_y);
      x_rep[n] = lavm_rng(eta_x[n], kappa);
    }
  }
}

mod <- cmdstan_model(stan_file)

stan_data <- list(
  N = n,
  y = y,
  x = as.numeric(x),
  z_x = z_x,
  z_y = z_y,
  U = 0.5,
  alpha = 0.5
)

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

fit$summary(variables = c("beta0", "beta1", "b1", "alpha0", "alpha1", "sd_y", "kappa", "kappa_eff"))

## # A tibble: 8 × 10
##   variable     mean  median       sd      mad       q5     q95  rhat ess_bulk
##   <chr>       <dbl>   <dbl>    <dbl>    <dbl>    <dbl>   <dbl> <dbl>    <dbl>
## 1 beta0       1.99    1.99   0.0562   0.0556    1.90     2.08  1.00     5635.
## 2 beta1       1.19    1.19   0.0583   0.0568    1.09     1.28  1.000    6443.
## 3 b1          0.694   0.693  0.101    0.0995    0.528    0.862 1.00     5636.
## 4 alpha0      0.102   0.102  0.00446  0.00440   0.0951   0.110 1.00     6790.
## 5 alpha1      0.500   0.500  0.00432  0.00435   0.493    0.507 1.00     5548.
## 6 sd_y        0.540   0.537  0.0386   0.0379    0.481    0.606 1.00     6329.
## 7 kappa     132.    132.    18.4     18.0     103.     164.    1.000    5617.
## 8 kappa_eff 251.    251.    35.2     34.6     196.     312.    1.000    5639.
## # ℹ 1 more variable: ess_tail <dbl>

The posterior centres for the regression coefficients and \(\sigma_y\) are close to their generating values, with the expected uncertainty from a sample of 100 observations.

library(bayesplot)

draws <- fit$draws(variables = c("beta0", "beta1", "b1", "alpha0", "alpha1", "sd_y", "kappa", "kappa_eff"), 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_b1 <- ggplot(draws, aes(x = b1)) +
  geom_density(fill = col3_t, alpha = 0.8, color = "white") +
  geom_vline(xintercept = b1, linetype = "dashed", color = "black", linewidth = 0.6) +
  labs(title = expression(b[1]), x = "", y = "") +
  theme_minimal() + theme(plot.title = element_text(face = "bold", hjust = 0.5))

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

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

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

p_kappa <- ggplot(draws, aes(x = kappa)) +
  geom_density(fill = col7_t, alpha = 0.8, color = "white") +
  geom_vline(xintercept = kappa, 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 = col8_t, alpha = 0.8, color = "white") +
  geom_vline(xintercept = kappa, linetype = "dashed", color = "black", linewidth = 0.6) +
  labs(title = expression(kappa[eff]), x = "", y = "") +
  theme_minimal() + theme(plot.title = element_text(face = "bold", hjust = 0.5))

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

The dashed lines mark the generating values, so the panels show both recovery and posterior uncertainty.

The posterior of \(\kappa\) should not recover the generating von Mises concentration directly because the data were generated from a von Mises distribution but fitted with an LAvM distribution. At predictor value \(\eta\), the locally comparable LAvM quantity is the effective concentration \[ \kappa_{\text{eff}}(\eta)=\kappa\left(1 + \eta^{2}\right)^{2} + \frac{1}{2}\eta^{2}\left(1+\eta^{2}\right), \] which is approximately \(\kappa(1 + \eta^{2})^{2}\) for large \(\kappa\). The generated quantity averages this expression over the observed \(\eta_i^x\) values. Its posterior is therefore much closer to the generating concentration than the raw LAvM \(\kappa\), illustrating why the two scales should not be compared directly.

We compare posterior predictive replicates with both responses.

yrep <- fit$draws("y_rep", format = "matrix"); xrep <- fit$draws("x_rep", format = "matrix")

set.seed(520)
idx <- sample(nrow(yrep), 100)
yrep_sample <- yrep[idx, ]; xrep_sample <- xrep[idx, ]

# --- PPC for Linear Response (y) ---
p_ppc_y <- ppc_dens_overlay(y = as.numeric(y), yrep = yrep_sample) +
  labs(title = "Linear variable (y)", x = "", y = "Density") +
  theme_minimal() +
  theme(legend.position = "none", plot.title = element_text(face = "bold", size = 13, hjust = 0.5))
# --- PPC for Circular Response (x) ---
p_ppc_x <- ppc_dens_overlay(y = x, yrep = xrep_sample) +
  scale_x_continuous(limits = c(-pi, pi), breaks = c(-pi, 0, pi), 
                     labels = c(expression(-pi), "0", expression(pi))) +
  labs(title = "Circular variable (x)", x = "", y = "") +
  theme_minimal() +
  theme(legend.position = "none", plot.title = element_text(face = "bold", size = 13, hjust = 0.5))

(p_ppc_y | p_ppc_x) + 
  plot_annotation(
    title = "Posterior Predictive Checks",
    theme = theme(plot.title = element_text(face = "bold", size = 16, hjust = 0.5))
  )

The replicated and observed marginal densities agree closely in this example.

2.2 Simulation example 2 (AR2 random effects)

The second simulation replaces the observed driver \(z_i^x\) with a latent time series. This is the main advantage of the joint formulation: serial structure learned from the circular observations can enter the linear response without treating the raw angle as an ordinary covariate. The model is \[ \begin{align} y_{i} &\sim N(\beta_{0} + \beta_{1}z_{i} + b_{1}\left(\alpha_{0} + w_{i}\right) + s_{i}, \sigma_{y}^{2}) \\ x_{i} &\sim \text{LAvM}(\alpha_{0} + w_{i}, \kappa) \end{align} \] The process \(\boldsymbol w=(w_1,\ldots,w_n)\) is shared. It explains serial variation in \(x_i\) through \(\eta_i^x=\alpha_0+w_i\) and enters the mean of \(y_i\) after multiplication by \(b_1\). The process \(\boldsymbol s=(s_1,\ldots,s_n)\) is specific to \(y\) and absorbs serial variation that is not shared with \(x\). Keeping these components separate prevents every temporal pattern in \(y\) from being attributed to the circular response.

Both processes follow autoregressive models of order 2 (AR(2)), defined sequentially as \[ \begin{align} w_i &= \phi_{1,w} w_{i-1} + \phi_{2,w} w_{i-2} + \varepsilon_{w,i}, \quad \varepsilon_{w,i} \sim N(0, \sigma_w^2), \\ s_i &= \phi_{1,s} s_{i-1} + \phi_{2,s} s_{i-2} + \varepsilon_{s,i}, \quad \varepsilon_{s,i} \sim N(0, \sigma_s^2), \end{align} \] with both coefficient pairs restricted to the AR(2) stationarity triangle, \[ -1<\phi_2<1, \qquad \phi_2-1<\phi_1<1-\phi_2. \] Stationarity constrains the long-run dynamics, whereas centring each realised sequence at zero separates its level from the intercepts \(\beta_0\) and \(\alpha_0\). Conditional on \(z_i\) and the response-specific process \(s_i\), \(b_1\) remains the association between \(y_i\) and the shared circular predictor.

n <- 100
# Fixed Effects & Scales
beta0 <- 2; beta1 <- 1.2; b1 <- 0.7; alpha0 <- 0.1
sd_w <- 0.2; sd_s <- 0.3; kappa <- 200; sd_y <- 0.3
# AR(2) Parameters
phi1_w <- 1.2; phi2_w <- -0.3  # Shared process (w)
phi1_s <- 0.8; phi2_s <- -0.2  # Process (s) for y

set.seed(520)
# Generate AR2 processes w and s
z_w <- rnorm(n, 0, 1); z_s <- rnorm(n, 0, 1)
w_unc <- numeric(n); s_unc <- numeric(n)
w_unc[1:2] <- z_w[1:2] * 5.0; s_unc[1:2] <- z_s[1:2] * 5.0
for (i in 3:n) {
  w_unc[i] <- phi1_w * w_unc[i-1] + phi2_w * w_unc[i-2] + z_w[i] * sd_w
  s_unc[i] <- phi1_s * s_unc[i-1] + phi2_s * s_unc[i-2] + z_s[i] * sd_s
}
w <- w_unc - mean(w_unc); s <- s_unc - mean(s_unc)

# Generate Responses
eta_x <- alpha0 + w
z <- rnorm(n, mean = 0, sd = 1)
x <- rvm(n, mu = lfun(eta_x), kappa = kappa)
y <- rnorm(n, mean = beta0 + beta1*z + b1*eta_x + s, sd = sd_y)

par(mfrow = c(2, 3), pty = "s", mar = c(2, 2, 2, 1))
plot(1:n, s, type = "l", col = "black", lwd = 2, main = "Latent AR2 Path (s)", xlab = "", ylab = "")
hist(y, breaks = 15, col = col1_t, border = col1, main="y", xlab = "", ylab = "")
plot(1:n, y, type = "l", col = col1, lwd = 2, main="y", xlab = "", ylab = "")
points(1:n, y, pch = 16, cex = 0.6, col = "grey")
plot(1:n, w, type = "l", col = "black", lwd = 2, main = "Latent AR2 Path (w)", xlab = "", ylab = "")
rose.diag(x, bins = 20, col = col1_t, border = col1, prop = 1.5, main = "x")
plot(1:n, as.numeric(x), type = "l", col = col1, lwd = 2, main = "x", xlab = "", ylab = "")
points(1:n, as.numeric(x), pch = 16, cex = 0.6, col = "grey")

The histograms and rose diagram show the marginal distributions, while the ordered trajectories reveal serial persistence. The black curves are the latent AR(2) paths used to generate the data; the observed Gaussian and circular series contain only noisy, indirect information about them.

The Stan model constructs both paths from standard-normal innovations, scales the innovations by sd_w or sd_s, and then centres each realised path. This non-centred construction can reduce dependence between the latent states and their innovation scales. Parameter-dependent bounds restrict both AR(2) coefficient pairs to the stationary region; because the model assigns no additional density to these coefficients, their priors are uniform over that region.

stan_file <- file.path("stan", "joint_simulation_2.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;
  }
  
  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 = 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));
      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] x;
  vector[N] z;
  real<lower=0, upper=1> U;
  real<lower=0, upper=1> alpha;
}

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

parameters {
  real beta0;
  real beta1;
  real b1;
  real alpha0;
  
  // Sampled on unconstrained log scale for HMC geometry
  real log_kappa; 
  real<lower=0> sd_y;  
  
  // Dynamic Bounds for Shared Process (w)
  real<lower=-1, upper=1> phi2_w;
  real<lower=phi2_w-1, upper=1-phi2_w> phi1_w;
  real<lower=0> sd_w;
  vector[N] z_w; 
  
  // Dynamic Bounds for Specific Process (s)
  real<lower=-1, upper=1> phi2_s;
  real<lower=phi2_s-1, upper=1-phi2_s> phi1_s;
  real<lower=0> sd_s;
  vector[N] z_s; 
}

transformed parameters {
  // Exponentiate to get natural scale kappa
  real<lower=0> kappa = exp(log_kappa);
  
  vector[N] w_unc;
  vector[N] w;
  
  vector[N] s_unc;
  vector[N] s;
  
  w_unc[1] = z_w[1] * 5.0;
  w_unc[2] = z_w[2] * 5.0;
  for (i in 3:N) {
    w_unc[i] = phi1_w * w_unc[i-1] + phi2_w * w_unc[i-2] + z_w[i] * sd_w;
  }
  w = w_unc - mean(w_unc);
  
  s_unc[1] = z_s[1] * 5.0;
  s_unc[2] = z_s[2] * 5.0;
  for (i in 3:N) {
    s_unc[i] = phi1_s * s_unc[i-1] + phi2_s * s_unc[i-2] + z_s[i] * sd_s;
  }
  s = s_unc - mean(s_unc);
  
  vector[N] eta_x = alpha0 + w;
  vector[N] mu_y = beta0 + beta1 * z + b1 * eta_x + s;
}

model {
  // Jacobian adjustment for log_kappa
  target += log_kappa;

  // PC Prior Implementation
  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);
  
  // Priors
  beta0 ~ normal(0, 5);
  beta1 ~ normal(0, 5);
  b1 ~ normal(0, 5);
  alpha0 ~ normal(0, 5);
  sd_y ~ std_normal();
  
  sd_w ~ std_normal();
  sd_s ~ std_normal();
  z_w ~ std_normal();
  z_s ~ std_normal();
  
  // Joint model
  y ~ normal(mu_y, sd_y);
  x ~ lavm(eta_x, kappa);
}

generated quantities {
  vector[N] log_lik;     
  vector[N] log_lik_y;   
  vector[N] log_lik_x;   
  
  array[N] real y_rep;
  array[N] real x_rep;
  
  // Explicitly export effective concentration
  real kappa_eff; 
  
  { 
    vector[N] eta_sq = square(eta_x);
    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);
  
    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) {
      real ll_y = normal_lpdf(y[n] | mu_y[n], sd_y);
      log_lik_y[n] = ll_y;
      
      real tan_half_x = tan(x[n] / 2.0);
      real A = cos(2.0 * atan(tan_half_x - eta_x[n]));
      real D = 1.0 + square(eta_x[n]) - eta_x[n] * sin(x[n]) - square(eta_x[n]) * square(sin(x[n] / 2.0));
      real ll_x = (kappa * A) - log_bessel - log(D);
      log_lik_x[n] = ll_x;
      
      log_lik[n] = ll_y + ll_x;
      
      y_rep[n] = normal_rng(mu_y[n], sd_y);
      x_rep[n] = lavm_rng(eta_x[n], kappa);
    }
  }
}

mod <- cmdstan_model(stan_file)

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

init_fun <- function() {
  list(
    phi1_w = runif(1, -0.1, 0.1), phi2_w = runif(1, -0.1, 0.1),
    phi1_s = runif(1, -0.1, 0.1), phi2_s = runif(1, -0.1, 0.1),
    sd_w = runif(1, 0.1, 0.3), sd_s = runif(1, 0.1, 0.3),
    log_kappa = log(runif(1, 100, 300))
  )
}

set.seed(52003)
fit <- mod$sample(
  data = stan_data,
  init = init_fun,
  seed = 52003,
  refresh = 0,
  iter_warmup = 1500,
  iter_sampling = 5000
)

fit$summary(variables = c(
  "beta0", "beta1", "b1", "alpha0", 
  "sd_y", "kappa_eff", 
  "phi1_w", "phi2_w", "sd_w", 
  "phi1_s", "phi2_s", "sd_s"
))

## # A tibble: 12 × 10
##    variable     mean  median      sd     mad       q5      q95  rhat ess_bulk
##    <chr>       <dbl>   <dbl>   <dbl>   <dbl>    <dbl>    <dbl> <dbl>    <dbl>
##  1 beta0       1.97    1.97   0.0454  0.0441   1.89     2.04   1.00    13330.
##  2 beta1       1.25    1.25   0.0478  0.0477   1.17     1.33   1.00    19428.
##  3 b1          0.758   0.741  0.142   0.129    0.560    1.02   1.00     4318.
##  4 alpha0      0.131   0.132  0.0216  0.0214   0.0962   0.167  1.000   23840.
##  5 sd_y        0.369   0.368  0.0476  0.0451   0.295    0.449  1.00     3306.
##  6 kappa_eff 245.    239.    53.7    52.4    166.     340.     1.00     9338.
##  7 phi1_w      1.28    1.29   0.103   0.102    1.11     1.44   1.00     7551.
##  8 phi2_w     -0.369  -0.372  0.0962  0.0952  -0.520   -0.206  1.00     7693.
##  9 sd_w        0.159   0.157  0.0280  0.0274   0.116    0.208  1.00     6280.
## 10 phi1_s      0.862   0.867  0.0855  0.0825   0.715    0.994  1.00     5352.
## 11 phi2_s     -0.205  -0.205  0.0676  0.0668  -0.316   -0.0932 1.00    11097.
## 12 sd_s        0.243   0.239  0.0580  0.0538   0.157    0.343  1.00     2967.
## # ℹ 1 more variable: ess_tail <dbl>

The posterior intervals cover the generating values, with greater uncertainty than in the fixed-effect example because the model must separate a shared path, a response-specific path, and observation noise using only 100 time points.

draws <- fit$draws(
  variables = c("beta0", "beta1", "b1", "alpha0", "sd_y", "kappa_eff", 
                "phi1_w", "phi2_w", "sd_w", "phi1_s", "phi2_s", "sd_s"),
  format = "df")

plot_dens <- function(data, var, true_val, color_fill, label_y) {
  ggplot(data, aes(x = .data[[var]])) + 
    geom_density(fill = color_fill, color = "white", alpha = 0.8) +
    geom_vline(xintercept = true_val, linetype = "dashed", color = "black", linewidth = 0.6) +
    labs(title = label_y, x = "", y = "") + 
    theme_minimal() + 
    theme(plot.title = element_text(face = "bold", hjust = 0.5))
}

p_beta0 <- plot_dens(draws, "beta0", beta0, col1_t, expression(beta[0])) + labs(y = "Density")
p_beta1 <- plot_dens(draws, "beta1", beta1, col1_t, expression(beta[1]))
p_b1 <- plot_dens(draws, "b1", b1, col1_t, expression(b[1]))
p_alpha0 <- plot_dens(draws, "alpha0", alpha0, col1_t, expression(alpha[0])) + labs(y = "Density")
p_sd_y <- plot_dens(draws, "sd_y", sd_y, col2_t, expression(sigma[y]))
p_kappa_eff <- plot_dens(draws, "kappa_eff", kappa, col2_t, expression(kappa[eff]))
p_phi1_w <- plot_dens(draws, "phi1_w", phi1_w, col2_t, expression(phi[1]^w)) + labs(y = "Density")
p_phi2_w <- plot_dens(draws, "phi2_w", phi2_w, col2_t, expression(phi[2]^w))
p_sd_w <- plot_dens(draws, "sd_w", sd_w, col3_t, expression(sigma[w]))
p_phi1_s <- plot_dens(draws, "phi1_s", phi1_s, col3_t, expression(phi[1]^s)) + labs(y = "Density")
p_phi2_s <- plot_dens(draws, "phi2_s", phi2_s, col3_t, expression(phi[2]^s))
p_sd_s <- plot_dens(draws, "sd_s", sd_s, col3_t, expression(sigma[s]))

(p_beta0 | p_beta1 | p_b1 | p_alpha0) /
(p_sd_y | p_kappa_eff | p_phi1_w | p_phi2_w) /
(p_sd_w | p_phi1_s | p_phi2_s | p_sd_s) + 
  plot_annotation(
    title = "Posterior Distributions",
    theme = theme(plot.title = element_text(face = "bold", size = 15, hjust = 0.5))
  )

Dashed lines mark the generating values.

w_sum <- fit$summary("w")
s_sum <- fit$summary("s")

df_w <- data.frame(Index = 1:n, True = w, Est = w_sum$mean, Lower = w_sum$q5, Upper = w_sum$q95)
df_s <- data.frame(Index = 1:n, True = s, Est = s_sum$mean, Lower = s_sum$q5, Upper = s_sum$q95)

p_lat_w <- ggplot(df_w, aes(x = Index)) +
  geom_ribbon(aes(ymin = Lower, ymax = Upper), fill = col3_t, alpha = 0.4) +
  geom_line(aes(y = Est, color = "Posterior Mean"), linewidth = 1) +
  geom_line(aes(y = True, color = "True Latent Field"), linewidth = 0.8, linetype = "dashed") +
  scale_color_manual(values = c("Posterior Mean" = col3, "True Latent Field" = "black")) +
  labs(title = "w (Shared)", x = "index", y = "") +
  theme_minimal() + theme(legend.position = "none", plot.title = element_text(face = "bold", size = 13, hjust = 0.5))
p_lat_s <- ggplot(df_s, aes(x = Index)) +
  geom_ribbon(aes(ymin = Lower, ymax = Upper), fill = col4_t, alpha = 0.4) +
  geom_line(aes(y = Est, color = "Posterior Mean"), linewidth = 1) +
  geom_line(aes(y = True, color = "True Latent Field"), linewidth = 0.8, linetype = "dashed") +
  scale_color_manual(values = c("Posterior Mean" = col4, "True Latent Field" = "black")) +
  labs(title = "s (y Specific)", x = "index", y = "") +
  theme_minimal() + theme(legend.position = "none", plot.title = element_text(face = "bold", size = 13, hjust = 0.5))

(p_lat_w | p_lat_s) + 
  plot_annotation(
    title = "Latent processes",
    theme = theme(plot.title = element_text(face = "bold", size = 16, hjust = 0.5))
  )

The solid lines are posterior means, the ribbons are pointwise 90% intervals, and the dashed lines are the simulated paths. Recovery of \(w\) shows that the model can identify variation shared by \(x\) and \(y\); recovery of \(s\) shows that response-specific dynamics are not forced into that shared component.

yrep <- fit$draws("y_rep", format = "matrix")
xrep <- fit$draws("x_rep", format = "matrix")
set.seed(520)
idx <- sample(nrow(yrep), 100)
yrep_sample <- yrep[idx, ]
xrep_sample <- xrep[idx, ]

p_ppc_y <- ppc_dens_overlay(y = as.numeric(y), yrep = yrep_sample) +
  labs(title = "Linear response (y)", x = "", y = "Density") +
  theme_minimal() + theme(legend.position = "none", plot.title = element_text(face = "bold", size = 13, hjust = 0.5))
p_ppc_x <- ppc_dens_overlay(y = x, yrep = xrep_sample) +
  scale_x_continuous(limits = c(-pi, pi), breaks = c(-pi, 0, pi), labels = c(expression(-pi), "0", expression(pi))) +
  labs(title = "Circular response (x)", x = "", y = "") +
  theme_minimal() + theme(legend.position = "none", plot.title = element_text(face = "bold", size = 13, hjust = 0.5))

(p_ppc_y | p_ppc_x) + 
  plot_annotation(
    title = "Posterior predictive checks",
    theme = theme(plot.title = element_text(face = "bold", size = 16, hjust = 0.5))
  )

The posterior replicates reproduce the broad marginal shapes of both responses.

2.3 Wind application example

We now move from controlled simulations to hourly observations from the NOAA station at John F. Kennedy International Airport in New York. This tutorial uses the January subset: wind direction \(x\) is measured in radians, wind speed \(y\) in metres per second, and normal temperature \(z\) is included as a linear covariate. The aim is to model direction and speed jointly while allowing both long-run dependence and recurring hour-of-day patterns.

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")
# Subset the data to January
dat <- wind_data %>% filter(month %in% c(1))
n <- nrow(dat)

x <- as.numeric(dat$HLY.WIND.VCTDIR)
y <- dat$HLY.WIND.VCTSPD
z <- dat$HLY.TEMP.NORMAL

# Diurnal index (Stan uses 1-based indexing, so 0-23 becomes 1-24)
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

We rotate direction around its circular mean so that the dominant direction lies away from the branch cut at \(-\pi\) and \(\pi\). This changes only the coordinate origin, not the relative angles, and gives a more convenient representation for a concentrated LAvM response. Angular summaries for scientific interpretation should be rotated back to the original compass coordinates. The variable hour_idx maps repeated observations to 24 hour-of-day levels.

plot_dat <- data.frame(
  direction = x,
  direction_rotated = x_centered %% (2 * pi),
  speed = y,
  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("0", expression(pi / 2), expression(pi), expression(3 * pi / 2))) +
  labs(x = "direction relative to 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 = speed)) +
  geom_histogram(bins = 30, fill = "#E76F51", color = "white", alpha = 0.9) +
  scale_x_continuous(n.breaks = 4, guide = guide_axis(check.overlap = TRUE)) +
  labs(x = "wind speed (m/s)", y = NULL) +
  theme_minimal(base_size = 14)

p4 <- 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, speed = y)

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

p6 <- ggplot(trace_dat, aes(x = time, y = speed)) +
  geom_line(color = "#E76F51", alpha = 0.9, linewidth = 0.8) +
  scale_x_datetime(breaks = scales::breaks_pretty(n = 3), guide = guide_axis(check.overlap = TRUE)) +
  labs(x = NULL, y = "wind speed (m/s)") +
  theme_minimal(base_size = 14)

(p1 | p3 | p4) / (p2 | p5 | p6)

The polar plots verify the rotation, while the ordered series show that both direction and speed occur in temporal runs rather than as independent scatter. This motivates separate long-run and diurnal components in the joint model.

We model wind speed \(y\) and wind direction \(x\) jointly as \[ \begin{aligned} y_{i}\mid w_{i}, s_{i}, \boldsymbol{\phi},\boldsymbol{\psi} &\sim \operatorname{Gamma}\left( a = \rho, b = \rho \exp\left\{ -\eta^{y}_{i} \right\} \right) \\ x_{i}\mid w_{i}, \boldsymbol{\phi},\boldsymbol{\psi} &\sim \operatorname{LAvM}\left( a_{0} + a_{1}w_{i} + a_{2}w_{2i} + \alpha z_{i}, \kappa \right) \\ \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), \mathbf{s}\sim N\left( 0, \sigma_{s}\mathbf{Q}_{\text{AR2}}^{-1} \right), \mathbf{s}_{2}\sim N\left( 0, \sigma_{s_{2}}\mathbf{Q}_{\text{RW2}}^{-1} \right), \\ b_{1}, a_{0}, &b_{0}, \alpha, \beta \sim N\left(0, 1\right), \quad \kappa \sim \mathcal{PC}_{\kappa}\left(0.5,0.99\right), \\ a_{1}, a_{2}, \sigma_{s},& \sigma_{s_{2}} \sim \operatorname{Exp}\left(-\log(0.5)/0.5\right), \quad \rho \sim \operatorname{Gamma}\left(1,0.01\right), \\ (\phi_{1,w}, \phi_{2,w}), (\phi_{1,s}, \phi_{2,s}) &\sim \operatorname{Uniform}(\mathcal{S}_{\text{AR2}}), \\ \text{where} \quad \eta^{y}_{i} = &b_{0} + b_{1}\left(a_{0} + a_{1}w_{i} + a_{2}w_{2i} + \alpha z_{i}\right) + s_{i} + s_{2i} + \beta z_{i}. \end{aligned} \] The stationary triangle \(\mathcal{S}_{\text{AR2}}\) is \[\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\}.\] The direction predictor \[ \eta_i^x=a_0+a_1w_i+a_2w_{2i}+\alpha z_i \] is shared with wind speed through \(b_1\eta_i^x\). Thus \(b_1\) measures the conditional direction–speed association carried by the systematic direction component. Temperature affects direction through \(\alpha\), speed directly through \(\beta\), and speed indirectly through the product \(b_1\alpha\). The long-run process \(\mathbf w\) contributes to both responses through \(\eta_i^x\), whereas \(\mathbf s\) captures long-run variation specific to speed. Their AR(2) recursions are \[ \begin{align} w_i &= \phi_{1,w}w_{i-1}+\phi_{2,w}w_{i-2}+\varepsilon_{w,i}, && \varepsilon_{w,i}\sim N(0,1),\\ s_i &= \phi_{1,s}s_{i-1}+\phi_{2,s}s_{i-2}+\varepsilon_{s,i}, && \varepsilon_{s,i}\sim N(0,\sigma_s^2), \end{align} \] with both coefficient pairs restricted to \(\mathcal S_{\text{AR2}}\).

The 24-dimensional RW2 effects \(\mathbf w_2\) and \(\mathbf s_2\) represent smooth hour-of-day patterns for direction and speed. Sum and linear-trend constraints separate them from intercepts and trends. The recursion smooths the ordered levels 1 through 24 but does not join hour 24 back to hour 1, so the displayed effects are non-cyclic across midnight.

Under the Gamma shape–rate parameterisation, \[ E(y_i\mid\eta_i^y)=\exp(\eta_i^y), \qquad \operatorname{Var}(y_i\mid\eta_i^y)=\frac{\exp(2\eta_i^y)}{\rho}. \] Consequently, \(\rho\) controls conditional dispersion rather than the mean. This distinction matters when translating the code to software that uses a Gamma shape–scale convention.

stan_file <- file.path("stan", "joint_wind_application.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. 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 = 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));
      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] 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;
  
  real sum_log_y = sum(log(y));
}

parameters {
  // Coefficients
  real a0;
  real<lower=0> a1; 
  real<lower=0> a2; 
  real alpha;
  
  real b0;
  real b1;
  real beta;
  
  // Dispersion
  real log_kappa; 
  real<lower=1e-4> rho; // Bounded to prevent 0.0 underflow
  
  // AR(2) Shared Process (w)
  real<lower=-1, upper=1> phi2_w;
  real<lower=phi2_w-1, upper=1-phi2_w> phi1_w;
  vector[N] z_w; 
  
  // AR(2) Specific Process (s)
  real<lower=-1, upper=1> phi2_s;
  real<lower=phi2_s-1, upper=1-phi2_s> phi1_s;
  real<lower=0> sigma_s;
  vector[N] z_s; 
  
  // RW2 Diurnal Processes (w2 and s2)
  vector[24] w2;
  vector[24] s2;
  real<lower=0> sigma_s2;
}

transformed parameters {
  real<lower=0> kappa = exp(log_kappa);
  
  vector[N] w_unc;
  vector[N] s_unc;
  
  w_unc[1] = z_w[1] * 5.0;
  w_unc[2] = z_w[2] * 5.0;
  s_unc[1] = z_s[1] * 5.0;
  s_unc[2] = z_s[2] * 5.0;
  
  for (i in 3:N) {
    w_unc[i] = phi1_w * w_unc[i-1] + phi2_w * w_unc[i-2] + z_w[i];
    s_unc[i] = phi1_s * s_unc[i-1] + phi2_s * s_unc[i-2] + z_s[i] * sigma_s;
  }
  
  vector[N] w = w_unc - mean(w_unc);
  vector[N] s = s_unc - mean(s_unc);
  
  vector[N] eta_x = a0 + a1 * w + a2 * w2[hour_idx] + alpha * z;
  vector[N] eta_y = b0 + b1 * eta_x + s + s2[hour_idx] + beta * 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
  a0 ~ normal(0, 1);
  b0 ~ normal(0, 1);
  b1 ~ normal(0, 1);
  alpha ~ normal(0, 1);
  beta ~ normal(0, 1);
  
  a1 ~ exponential(lambda_tau);
  a2 ~ exponential(lambda_tau);
  sigma_s ~ exponential(lambda_tau);
  sigma_s2 ~ exponential(lambda_tau);
  rho ~ gamma(1, 0.01);
  
  z_w ~ std_normal();
  z_s ~ std_normal();
  
  // RW2 Generative Models
  w2[1] ~ normal(0, 10);
  w2[2] ~ normal(0, 10);
  for (h in 3:24) w2[h] ~ normal(2.0 * w2[h-1] - w2[h-2], 1.0); 
  sum(w2) ~ normal(0, 0.001 * 24); 
  dot_product(w2, seq_idx_24) ~ normal(0, 0.001 * 24);
  
  s2[1] ~ normal(0, 10);
  s2[2] ~ normal(0, 10);
  for (h in 3:24) s2[h] ~ normal(2.0 * s2[h-1] - s2[h-2], sigma_s2); 
  sum(s2) ~ normal(0, 0.001 * 24);
  dot_product(s2, seq_idx_24) ~ normal(0, 0.001 * 24);
  
  // --- Model for x and likelihood given eta_y ---
  x ~ lavm(eta_x, kappa);
  target += N * (rho * log(rho) - lgamma(rho)) 
          - rho * sum(eta_y) 
          + (rho - 1.0) * sum_log_y 
          - rho * sum(y .* exp(-eta_y));
}

generated quantities {
  array[N] real x_rep;
  array[N] real y_rep;
  
  for (i in 1:N) {
    x_rep[i] = lavm_rng(eta_x[i], kappa);
    y_rep[i] = gamma_rng(rho, rho * exp(-eta_y[i])); 
  }
}

mod <- cmdstan_model(
  stan_file = stan_file,
  cpp_options = list(
    ## Speedup for vectorized math
    STAN_CPP_OPTIMS = TRUE,
    ## disable indexing checks
    STAN_NO_RANGE_CHECKS = TRUE
  )
)

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

init_fun <- function() {
  list(
    # --- Fixed Effects ---
    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),
    b0 = rnorm(1, mean = 0, sd = 0.1),
    b1 = rnorm(1, mean = 0, sd = 0.1),
    beta = rnorm(1, mean = 0, sd = 0.1),
    
    # --- Scaling parameters ---
    log_kappa = log(runif(1, min = 3.0, max = 7.0)),
    rho = runif(1, min = 1.5, max = 2.5),
    
    # --- AR(2) Processes ---
    phi2_w = runif(1, min = -0.1, max = 0.1),
    phi1_w = runif(1, min = 0.3, max = 0.7),
    z_w = rnorm(n, mean = 0, sd = 0.1),

    phi2_s = runif(1, min = -0.1, max = 0.1),
    phi1_s = runif(1, min = 0.3, max = 0.7),
    sigma_s = runif(1, min = 0.05, max = 0.2),
    z_s = rnorm(n, mean = 0, sd = 0.1),
    
    # --- RW2 Processes ---
    w2 = rnorm(24, mean = 0, sd = 0.1),
    s2 = rnorm(24, mean = 0, sd = 0.1),
    sigma_s2 = runif(1, min = 0.05, max = 0.2)
  )
}

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

structural_vars <- c(
  "a0", "a1", "a2", "alpha",             # Direction fixed effects
  "b0", "b1", "beta",                    # Speed fixed effects & coupling
  "rho", "log_kappa",                    # Dispersion and concentration
  "phi1_w", "phi2_w",                    # Shared AR(2)
  "phi1_s", "phi2_s",                    # AR(2)
  "sigma_s", "sigma_s2"                  # Latent-process scales
)

fit$summary(variables = structural_vars)

## # A tibble: 15 × 10
##    variable       mean   median      sd     mad       q5      q95  rhat ess_bulk
##    <chr>         <dbl>    <dbl>   <dbl>   <dbl>    <dbl>    <dbl> <dbl>    <dbl>
##  1 a0          6.53e-1  6.54e-1 2.60e-2 2.67e-2  6.10e-1  6.94e-1  1.12    23.6 
##  2 a1          5.81e-3  5.70e-3 7.84e-4 8.26e-4  4.68e-3  7.23e-3  1.38     8.63
##  3 a2          2.27e-2  2.21e-2 3.59e-3 3.06e-3  1.82e-2  2.84e-2  1.02    70.8 
##  4 alpha      -1.98e-2 -1.98e-2 7.89e-4 8.09e-4 -2.11e-2 -1.85e-2  1.12    23.6 
##  5 b0          1.02e+0  1.02e+0 1.05e-1 1.07e-1  8.52e-1  1.20e+0  1.04   102.  
##  6 b1          4.35e-1  4.36e-1 1.32e-1 1.35e-1  2.17e-1  6.48e-1  1.07    42.8 
##  7 beta        2.58e-2  2.58e-2 3.17e-3 3.25e-3  2.04e-2  3.09e-2  1.04   101.  
##  8 rho         2.72e+3  2.71e+3 2.45e+2 2.38e+2  2.34e+3  3.15e+3  1.06    47.6 
##  9 log_kappa   8.66e+0  8.64e+0 1.28e-1 1.24e-1  8.47e+0  8.89e+0  1.28    10.6 
## 10 phi1_w      1.31e+0  1.32e+0 1.17e-1 1.29e-1  1.10e+0  1.48e+0  1.36     8.98
## 11 phi2_w     -3.25e-1 -3.41e-1 1.15e-1 1.24e-1 -4.89e-1 -1.23e-1  1.33     9.60
## 12 phi1_s      1.18e+0  1.16e+0 1.47e-1 1.42e-1  9.59e-1  1.44e+0  1.09    33.5 
## 13 phi2_s     -1.83e-1 -1.70e-1 1.46e-1 1.41e-1 -4.50e-1  3.33e-2  1.09    32.5 
## 14 sigma_s     1.03e-2  1.04e-2 1.88e-3 1.90e-3  7.02e-3  1.32e-2  1.10    27.0 
## 15 sigma_s2    2.63e-2  2.56e-2 4.85e-3 4.49e-3  1.94e-2  3.50e-2  1.00  2016.  
## # ℹ 1 more variable: ess_tail <dbl>

Diagnostic warning. This wind fit has not converged adequately: several split-\(\widehat R\) values are well above 1.01 and some effective sample sizes are very small. Its summaries and figures illustrate the workflow only and should not be used for scientific conclusions.

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_b0 <- plot_density(draws, "b0", expression(b[0])) + labs(y = "Density")
p_b1 <- plot_density(draws, "b1", expression(b[1]))
p_beta <- plot_density(draws, "beta", expression(beta))
p_rho <- plot_density(draws, "rho", expression(rho))

p_log_kappa <- plot_density(draws, "log_kappa", expression(log(kappa))) + labs(y = "Density")
p_phi1_w <- plot_density(draws, "phi1_w", expression(phi[1]^w))
p_phi2_w <- plot_density(draws, "phi2_w", expression(phi[2]^w))
p_phi1_s <- plot_density(draws, "phi1_s", expression(phi[1]^s))

p_phi2_s <- plot_density(draws, "phi2_s", expression(phi[2]^s)) + labs(y = "Density")
p_sigma_s <- plot_density(draws, "sigma_s", expression(sigma[s]))
p_sigma_s2 <- plot_density(draws, "sigma_s2", expression(sigma[s2]))

(p_a0 | p_a1 | p_a2 | p_alpha) /
(p_b0 | p_b1 | p_beta | p_rho) /
(p_log_kappa | p_phi1_w | p_phi2_w | p_phi1_s) /
(p_phi2_s | p_sigma_s | p_sigma_s2 | plot_spacer()) + 
  plot_annotation(
    title = "Posterior Marginal Distributions",
    theme = theme(plot.title = element_text(face = "bold", size = 18, hjust = 0.5))
  )

Once a refitted model has converged, \(\alpha\) and \(\beta\) describe the temperature associations with the direction predictor and log mean speed, while \(b_1\) links the latent direction predictor to log mean speed. The left-hand latent plots show long-run AR(2) paths over observation time; the right-hand plots show the 24 hour-of-day effects reused across days.

get_latent_summary <- function(fit, param_name, n_state) {
  draws_mat <- fit$draws(param_name, format = "matrix")
  data.frame(
    idx = seq_len(n_state),
    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)
s_summary <- get_latent_summary(fit, "s", n)
w2_summary <- get_latent_summary(fit, "w2", 24)
s2_summary <- get_latent_summary(fit, "s2", 24)

w_summary$time <- dat$datetime
s_summary$time <- dat$datetime

# --- Shared AR(2) Process (w) ---
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("Shared AR(2) Process (" * w * ")"), x = "", y = "Value") +
  theme_minimal(base_size = 14) + theme(plot.title = element_text(face = "bold"))

# --- AR(2) Process (s) ---
p_s <- ggplot(s_summary, aes(x = time, y = mean)) +
  geom_ribbon(aes(ymin = lower, ymax = upper), fill = col2_t, alpha = 0.3) +
  geom_line(color = col2, linewidth = 0.8) +
  labs(title = expression("AR(2) Process (" * s * ")"), x = "Time", y = "Value") +
  theme_minimal(base_size = 14) + theme(plot.title = element_text(face = "bold"))

# --- Shared Diurnal RW2 (w2) ---
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("Shared Diurnal RW2 (" * w[2] * ")"), x = "", y = "") +
  theme_minimal(base_size = 14) + theme(plot.title = element_text(face = "bold"))

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

(p_w | p_w2) / (p_s | p_s2) +
  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))
  )

Finally, we compare the observed marginals with posterior predictive replicates.

y_rep <- fit$draws("y_rep", format = "matrix")
x_rep <- fit$draws("x_rep", format = "matrix")

set.seed(123)
sample_idx <- sample(nrow(y_rep), 100)
y_rep_sample <- y_rep[sample_idx, ]
x_rep_sample <- x_rep[sample_idx, ]

# --- Wind Speed ---
p_ppc_y <- ppc_dens_overlay(y, y_rep_sample) +
  labs(
    title = "Posterior Predictive Check: Wind Speed",
    subtitle = "Gamma Distribution",
    x = "Wind Speed (m/s)",
    y = "Density"
  ) +
  theme_minimal(base_size = 14) +
  theme(plot.title = element_text(face = "bold"), legend.position = "none")

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

p_ppc_y | p_ppc_x

Across the examples, the circular observation is modeled in an LAvM submodel and its real-valued predictor—rather than the raw angle—is shared with the second response. This preserves the circular geometry while allowing \(b_1\) to describe cross-response association, and it extends naturally to shared and response-specific temporal effects.

Tutorial: Joint regression with circular covariates

Xiang Ye

2026-09-07