Tutorial: Joint regression with circular covariates
Xiang Ye
2026-09-11
This tutorial uses the circular model’s linear predictor as a local linear representation of circular data in a joint response model. We begin with fixed effects, add shared temporal effects, and finish with a joint analysis of 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.
The rendered tutorial uses saved fits and their matching datasets. To
recompute them, run Rscript R/run_all.R from this
tutorial’s directory; rendering then regenerates summaries and figures
without sampling.
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
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 transformed location is \(g(\eta)\), and the density concentrates near it as \(\kappa\) increases. Its mode can differ from \(g(\eta)\) because of the Jacobian in the denominator.
The link adjustment compares \(g^{-1}(x)\) with \(\eta\) on the real line and maps the difference back to the circle. Regression coefficients act linearly on the \(\eta\) scale; applying \(g\) gives their nonlinear effect on the angular location.
When \(\eta=0\), the LAvM density is the corresponding von Mises density centred at zero. Changing \(\eta\) changes both angular location and local concentration, while increasing \(\kappa\) sharpens the density. Thus the fitted LAvM \(\kappa\) alone does not describe concentration on the angular scale.
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}) \] A first-harmonic model avoids this discontinuity: \[ 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. To share temporal, spatial, or other latent structure underlying \(x\) with \(y\), we can instead 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\) provides a local linear representation of the systematic circular variation on the chosen angular interval. 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\).
A positive \(b_1\) means that larger values of the linear predictor for \(x_i\) are associated with larger expected \(y_i\), conditional on the other terms in the model. The coefficient acts on the predictor scale, rather than per radian of observed angle. The fitted angular location is \(g(\eta)=2\arctan(\eta)\); predictive angles also include LAvM variation around that location.
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 contribution \(b_1\alpha_1z_i\) through the circular predictor without further constraints or prior information. With distinct, noncollinear covariates and \(\alpha_1\ne0\), the circular data identify \(\alpha_1\), while the Gaussian data distinguish \(\beta_1\) from \(b_1\).
library(circular)We simulate \(n=100\) observations using a von Mises distribution with high concentration and mean direction \(g(\eta_i^x)\). Fitting an LAvM model illustrates its local approximation to this concentrated distribution. We compare the regression estimates with the generating values and use an effective concentration to compare the two response families.
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 <- 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 \(\kappa\), we use the
penalised-complexity (PC) prior derived in the PC-prior tutorial, with the
point-mass base model \(\kappa\rightarrow\infty\). We set \(\lambda=-\log(1-\alpha)/\sqrt{1-U}\) using
its mean-resultant-length calibration, with \(U=\alpha=0.5\). The function
PC_vm_inf(kappa, lambda) supplies the log-density kernel
for the positive-concentration component of this prior.
library(cmdstanr)The model estimates concentration on the log scale. Adding
log_kappa to the target accounts for the PC prior being
defined on \(\kappa\). The function
lavm_rng samples predictive angles on a 400-point grid;
check the grid resolution when fine angular detail matters.
stan_file <- file.path("stan", "joint_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));
}
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=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 = -log1m(alpha) / sqrt(1.0 - U);
}
parameters {
real beta0;
real beta1;
real b1;
real alpha0;
real alpha1;
real log_kappa;
real<lower=0> sd_y;
}
transformed parameters {
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;
target += PC_vm_inf(kappa, lambda);
beta0 ~ normal(0, 5);
beta1 ~ normal(0, 5);
b1 ~ normal(0, 5);
alpha0 ~ normal(0, 5);
alpha1 ~ normal(0, 5);
sd_y ~ std_normal();
y ~ normal(mu_y, sd_y);
x ~ lavm(eta_x, kappa);
}
generated quantities {
vector[N] log_lik = lavm_log_densities(x, eta_x, kappa);
array[N] real y_rep = normal_rng(mu_y, sd_y);
array[N] real x_rep = lavm_rng(eta_x, kappa);
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);
for (n in 1:N) {
log_lik[n] += normal_lpdf(y[n] | mu_y[n], sd_y);
}
}
}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.0566 0.0567 1.89 2.08 1.000 4934.
## 2 beta1 1.19 1.19 0.0576 0.0551 1.09 1.28 1.00 5702.
## 3 b1 0.692 0.692 0.102 0.102 0.525 0.861 1.00 5328.
## 4 alpha0 0.102 0.102 0.00451 0.00434 0.0950 0.110 1.00 8417.
## 5 alpha1 0.500 0.500 0.00422 0.00410 0.493 0.507 1.000 6449.
## 6 sd_y 0.541 0.539 0.0391 0.0390 0.479 0.609 1.00 4876.
## 7 kappa 133. 131. 18.7 18.0 104. 166. 1.00 4793.
## 8 kappa_eff 252. 250. 35.6 34.5 197. 315. 1.00 4772.
## # ℹ 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)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", "b1", "alpha0", "alpha1", "sd_y", "kappa", "kappa_eff"), 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_b1 <- plot_density(draws, "b1", expression(b[1]), b1, col3_t)
p_alpha0 <- plot_density(draws, "alpha0", expression(alpha[0]), alpha0, col4_t)
p_alpha1 <- plot_density(draws, "alpha1", expression(alpha[1]), alpha1, col5_t,
y_label = "Density")
p_sd_y <- plot_density(draws, "sd_y", expression(sigma[y]), sd_y, col6_t)
p_kappa <- plot_density(draws, "kappa", expression(kappa), kappa, col7_t)
p_kappa_eff <- plot_density(draws, "kappa_eff", expression(kappa[eff]), kappa, col8_t)
(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))
)
To compare the fitted LAvM concentration with the generating von
Mises concentration, we use the negative log-density curvature at \(x=g(\eta)\): \[
\kappa_{\text{eff}}(\eta)=\kappa\left(1 + \eta^{2}\right)^{2} +
\frac{1}{2}\eta^{2}\left(1+\eta^{2}\right),
\] For large \(\kappa\), its
leading term is \(\kappa(1 +
\eta^{2})^{2}\). This is a local concentration comparison near
the transformed location, rather than an exact equivalence between the
two distributions. The reported kappa_eff averages the full
expression over observation-specific predictors within each posterior
draw; in this example, its posterior centre is closer to the generating
concentration than that of the baseline LAvM \(\kappa\).
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, ]
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))
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 uses a latent time series for the circular predictor, allowing serial structure learned from circular observations to inform the linear response. As in the first simulation, the circular observations are generated from a von Mises distribution and fitted with LAvM. The fitted 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)\) allows additional serial variation specific to \(y\).
Both paths are mean-centered AR(2) processes, with independent first states \(\widetilde w_{1:2},\widetilde s_{1:2}\sim N(0,5^2)\) and subsequent innovation standard deviations \(\sigma_w\) and \(\sigma_s\). Each coefficient pair is uniform on the stationary triangle \[ \mathcal S_{\mathrm{AR2}}=\{(\phi_1,\phi_2):-1<\phi_2<1,\ \phi_2-1<\phi_1<1-\phi_2\}. \] The displayed Stan code specifies the recursions and priors.
n <- 100
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
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)
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)
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 implementation constructs each path from standard-normal innovations.
stan_file <- file.path("stan", "joint_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] x;
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 b1;
real alpha0;
real log_kappa;
real<lower=0> sd_y;
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;
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 {
real<lower=0> kappa = exp(log_kappa);
vector[N] w_unc;
vector[N] w;
vector[N] s_unc;
vector[N] s;
w_unc[1:2] = 5 * z_w[1:2];
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:2] = 5 * z_s[1:2];
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;
target += PC_vm_inf(kappa, lambda);
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();
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 = lavm_log_densities(x, eta_x, kappa);
array[N] real y_rep = normal_rng(mu_y, sd_y);
array[N] real x_rep = lavm_rng(eta_x, kappa);
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);
for (n in 1:N) {
log_lik_y[n] = normal_lpdf(y[n] | mu_y[n], sd_y);
}
log_lik = log_lik_x + log_lik_y;
}
}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.0449 0.0441 1.89 2.04 1.00 10214.
## 2 beta1 1.25 1.25 0.0476 0.0473 1.17 1.33 1.000 19040.
## 3 b1 0.761 0.745 0.142 0.129 0.559 1.02 1.00 3132.
## 4 alpha0 0.131 0.131 0.0218 0.0215 0.0955 0.167 1.00 21730.
## 5 sd_y 0.367 0.367 0.0491 0.0441 0.294 0.446 1.00 2176.
## 6 kappa_eff 244. 239. 54.1 52.7 165. 342. 1.00 8475.
## 7 phi1_w 1.28 1.29 0.104 0.103 1.11 1.45 1.00 7646.
## 8 phi2_w -0.369 -0.373 0.0969 0.0956 -0.522 -0.205 1.00 7883.
## 9 sd_w 0.159 0.157 0.0278 0.0272 0.116 0.207 1.00 5734.
## 10 phi1_s 0.862 0.868 0.0890 0.0843 0.707 0.996 1.00 2922.
## 11 phi2_s -0.204 -0.205 0.0676 0.0665 -0.313 -0.0926 1.00 9234.
## 12 sd_s 0.245 0.239 0.0595 0.0550 0.158 0.345 1.00 1992.
## # ℹ 1 more variable: ess_tail <dbl>The displayed 90% posterior intervals contain the generating values. Their widths reflect uncertainty in separating the shared path, the response-specific path, and observation noise using 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")
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, col1_t)
p_b1 <- plot_density(draws, "b1", expression(b[1]), b1, col1_t)
p_alpha0 <- plot_density(draws, "alpha0", expression(alpha[0]), alpha0, col1_t,
y_label = "Density")
p_sd_y <- plot_density(draws, "sd_y", expression(sigma[y]), sd_y, col2_t)
p_kappa_eff <- plot_density(draws, "kappa_eff", expression(kappa[eff]), kappa, col2_t)
p_phi1_w <- plot_density(draws, "phi1_w", expression(phi[1]^w), phi1_w, col2_t,
y_label = "Density")
p_phi2_w <- plot_density(draws, "phi2_w", expression(phi[2]^w), phi2_w, col2_t)
p_sd_w <- plot_density(draws, "sd_w", expression(sigma[w]), sd_w, col3_t)
p_phi1_s <- plot_density(draws, "phi1_s", expression(phi[1]^s), phi1_s, col3_t,
y_label = "Density")
p_phi2_s <- plot_density(draws, "phi2_s", expression(phi[2]^s), phi2_s, col3_t)
p_sd_s <- plot_density(draws, "sd_s", expression(sigma[s]), sd_s, col3_t)
(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. These panels compare inference on shared and response-specific variation with the paths used to generate the data.
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 analyse the January hourly climate normals from the NOAA station at John F. Kennedy International Airport in New York. The variables are mean wind-vector direction \(x\), converted to radians, mean wind-vector magnitude \(y\), and normal temperature \(z\). These are climatological summaries for each calendar day and hour. The joint model describes their smooth calendar variation and recurring hour-of-day patterns. See the NOAA hourly normals documentation for the variable definitions.
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)
y <- dat$HLY.WIND.VCTSPD
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) - piWe rotate direction around its circular mean so that the dominant
direction lies away from the branch cut at \(-\pi\) and \(\pi\). The rotation preserves relative
angles and gives a 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 = 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 = 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 vector magnitude", 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 = "normal 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 vector magnitude") +
theme_minimal(base_size = 14)
(p1 | p3 | p4) / (p2 | p5 | p6)
The polar plots show the effect of rotation, while the ordered series reveal smooth calendar and daily patterns. We use an AR(1) process for direction, an AR(2) process for additional speed variation, and separate diurnal components.
We model wind speed \(y\) and wind direction \(x\) jointly as \[ \begin{aligned} y_{i}\mid\eta_i^y,\rho &\sim \operatorname{Gamma}\left( a = \rho, b = \rho \exp\left\{ -\eta^{y}_{i} \right\} \right), \\ x_{i}\mid\eta_i^x,\kappa &\sim \operatorname{LAvM}\left(\eta_i^x,\kappa\right), \\ \mathbf{w}_{2} \sim N\left( 0, \mathbf{Q}_{\text{cyclic RW2}}^{-1} \right), \quad &\mathbf{s}_{2}\sim N\left( 0, \sigma_{s_{2}}^{2}\mathbf{Q}_{\text{cyclic 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} \sim \operatorname{Uniform}(-1,1), \quad (\phi_{1,s}, \phi_{2,s}) &\sim \operatorname{Uniform}(\mathcal{S}_{\text{AR2}}), \\ \eta_i^x &= a_0+a_1w_i+a_2w_{2,h(i)}+\alpha z_i, \\ \eta^{y}_{i} &= b_{0} + b_{1}\eta_i^x + s_{i} + s_{2,h(i)} + \beta z_{i}. \end{aligned} \] Here \(h(i)\) is the hour-of-day index, and \(\mathcal S_{\mathrm{AR2}}\) is the stationary triangle defined above. The direction predictor enters the speed model through \(b_1\eta_i^x\), sharing its local linear representation. Temperature therefore contributes through \(\alpha\) to direction and through \(\beta+b_1\alpha\) to log mean speed.
The calendar-time paths are mean-centered stationary processes: AR(1) for direction and AR(2) for additional speed variation, initialized from their stationary distributions. Before centering, the direction path has variance 1 and the speed path has variance \(\sigma_s^2\), so \(a_1\) and \(\sigma_s\) describe their marginal scales. The hourly effects follow cyclic RW2 priors on the zero-sum subspace, smoothing second differences across midnight. Their precision inverses above are restricted to that subspace; the Stan code implements these priors.
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.
The implementation integrates Gaussian coefficients and hourly effects analytically, then recovers their conditional posterior draws for summaries. Scaled predictor coordinates keep the circular likelihood evaluation stable at large concentration.
stan_file <- file.path("stan", "joint_wind_application.stan")functions {
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));
}
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[8] 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]);
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';
grams[8] = Xc[1]' * Xc[2] + Xc[2]' * Xc[1];
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];
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);
real r1 = phi1 / (1 - phi2);
real initial_variance_factor = 1 - square(r1);
real innovation_variance_factor = initial_variance_factor * (1 - square(phi2));
real initial_sd = sigma * sqrt(initial_variance_factor);
real innovation_sd = sigma * sqrt(innovation_variance_factor);
vector[N] innovations;
// sigma is the stationary marginal SD before mean-centering.
innovations[1] = path[1] / sigma;
innovations[2] = (path[2] - r1 * path[1]) / initial_sd;
innovations[3:N] = (path[3:N] - phi1 * path[2:(N - 1)]
- phi2 * path[1:(N - 2)]) / innovation_sd;
return innovations;
}
vector ar2_static_cross(data matrix Xc, vector innovations,
real phi1, real phi2, real sigma) {
int N = rows(Xc);
real r1 = phi1 / (1 - phi2);
real initial_variance_factor = 1 - square(r1);
real innovation_variance_factor = initial_variance_factor * (1 - square(phi2));
real initial_sd = sigma * sqrt(initial_variance_factor);
real innovation_sd = sigma * sqrt(innovation_variance_factor);
return Xc[1]' * (innovations[1] / sigma - r1 * innovations[2] / initial_sd)
+ Xc[2]' * (innovations[2] / initial_sd)
+ (Xc[3:N]' * innovations[3:N]
- phi1 * (Xc[2:(N - 1)]' * innovations[3:N])
- phi2 * (Xc[1:(N - 2)]' * innovations[3:N])) / innovation_sd;
}
vector ar2_static_transpose(vector innovations, real phi1, real phi2, real sigma) {
int N = num_elements(innovations);
real r1 = phi1 / (1 - phi2);
real initial_variance_factor = 1 - square(r1);
real innovation_variance_factor = initial_variance_factor * (1 - square(phi2));
real initial_sd = sigma * sqrt(initial_variance_factor);
real innovation_sd = sigma * sqrt(innovation_variance_factor);
vector[N] result = rep_vector(0, N);
result[1] = innovations[1] / sigma - r1 * innovations[2] / initial_sd;
result[2] = innovations[2] / initial_sd;
result[3:N] += innovations[3:N] / innovation_sd;
result[2:(N - 1)] -= phi1 * innovations[3:N] / innovation_sd;
result[1:(N - 2)] -= phi2 * innovations[3:N] / innovation_sd;
return result;
}
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 r1 = phi1 / (1 - phi2);
real initial_variance_factor = 1 - square(r1);
real innovation_variance_factor = initial_variance_factor * (1 - square(phi2));
real initial_precision = inv(square(sigma) * initial_variance_factor);
real innovation_precision = inv(square(sigma) * innovation_variance_factor);
vector[N] v = ar2_static_standardize(rep_vector(1, N), phi1, phi2, sigma);
real precision_mean = dot_self(v);
vector[P] mean_cross = ((1 - r1) * initial_precision * sums[1]
+ ar_constant * innovation_precision
* (sums[2] - phi1 * sums[3] - phi2 * sums[4]))';
matrix[P, P] gram_raw = initial_precision * (grams[1] - r1 * grams[8])
+ innovation_precision
* (grams[2] + square(phi1) * grams[3] + square(phi2) * grams[4]
- phi1 * grams[5] - phi2 * grams[6] + phi1 * phi2 * grams[7]);
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 = dot_product(v, r);
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 = dot_product(v, residual) / precision_mean;
residual -= v * residual_mean;
real residual_ss = dot_self(residual)
+ square(mean_predictor - mean_design * conditional_mean)
+ dot_self(conditional_mean);
real log_det_R = -N * log(sigma) - 0.5 * log(initial_variance_factor)
- 0.5 * (N - 2) * log(innovation_variance_factor);
real log_density = -0.5 * N * log(2 * pi()) + log_det_R
- (N - 1) * log(loading)
- 0.5 * (log(precision_mean) + residual_ss)
- sum(log(diagonal(L)));
vector[N] score = ar2_static_transpose(residual, phi1, phi2, 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) ar2_observation_coordinates(
vector coordinates, vector score, real phi1, real phi2,
real sigma, real loading, real log_scale,
real curvature, real curvature_offset) {
int N = num_elements(coordinates);
real inverse_scale = exp(-log_scale);
real inverse_sqrt_scale = exp(-0.5 * log_scale);
real r1 = phi1 / (1 - phi2);
real initial_variance_factor = 1 - square(r1);
real innovation_variance_factor = initial_variance_factor * (1 - square(phi2));
real initial_precision = inverse_scale / (square(loading * sigma) * initial_variance_factor);
real transition_precision = inverse_scale / (square(loading * sigma) * innovation_variance_factor);
vector[N] diagonal = rep_vector(curvature + curvature_offset * inverse_scale, N);
vector[N - 1] subdiagonal = rep_vector(0, N - 1);
vector[N - 2] second_subdiagonal = rep_vector(-phi2 * transition_precision, N - 2);
matrix[N, 3] L = rep_matrix(0, N, 3);
vector[N] filtered_score;
vector[N] residual;
diagonal[1:2] += rep_vector(initial_precision, 2);
subdiagonal[1] -= r1 * initial_precision;
diagonal[3:N] += rep_vector(transition_precision, N - 2);
diagonal[2:(N - 1)] += rep_vector(square(phi1) * transition_precision, N - 2);
diagonal[1:(N - 2)] += rep_vector(square(phi2) * transition_precision, N - 2);
subdiagonal[1:(N - 2)] += rep_vector(phi1 * phi2 * transition_precision, N - 2);
subdiagonal[2:(N - 1)] -= rep_vector(phi1 * transition_precision, N - 2);
// AR(2) precision has two nonzero bands on either side of its diagonal.
for (i in 1:N) {
if (i > 2) L[i, 3] = second_subdiagonal[i - 2] / L[i - 2, 1];
if (i > 1) L[i, 2] = (subdiagonal[i - 1] - L[i, 3] * L[i - 1, 2]) / L[i - 1, 1];
L[i, 1] = sqrt(diagonal[i] - square(L[i, 2]) - square(L[i, 3]));
filtered_score[i] = score[i];
if (i > 1) filtered_score[i] -= L[i, 2] * filtered_score[i - 1];
if (i > 2) filtered_score[i] -= L[i, 3] * filtered_score[i - 2];
filtered_score[i] /= L[i, 1];
}
for (j in 1:N) {
int i = N + 1 - j;
residual[i] = coordinates[i] + inverse_sqrt_scale * filtered_score[i];
if (i < N) residual[i] -= L[i + 1, 2] * residual[i + 1];
if (i < N - 1) residual[i] -= L[i + 2, 3] * residual[i + 2];
residual[i] /= L[i, 1];
}
return (residual, -sum(log(L[, 1])));
}
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)';
}
tuple(real, vector, matrix, vector) gaussian_shared_regression_distribution(
vector predictor, vector shared_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);
tuple(real, vector[P], matrix[P, P], vector[N]) base = gaussian_regression_distribution(
predictor, Xc, mean_X, grams, sums, scales, loading, phi1, phi2, sigma);
real mean_predictor = mean(predictor);
real mean_shared = mean(shared_predictor);
vector[N] shared_centered = shared_predictor - mean_shared;
row_vector[P] mean_design = mean_X .* scales';
real r1 = phi1 / (1 - phi2);
real initial_variance_factor = 1 - square(r1);
real innovation_variance_factor = initial_variance_factor * (1 - square(phi2));
vector[N] v = ar2_static_standardize(rep_vector(1, N), phi1, phi2, sigma);
real precision_mean = dot_self(v);
vector[N] h = ar2_static_standardize(shared_centered / loading, phi1, phi2, sigma);
h -= v * (dot_product(v, h) / precision_mean);
vector[P] cross = mean_design' * mean_shared
+ (scales / loading) .* ar2_static_cross(Xc, h, phi1, phi2, sigma);
vector[P] L_cross = mdivide_left_tri_low(base.3, cross);
vector[P] shared_coef = mdivide_right_tri_low(L_cross', base.3)';
vector[N] shared_residual = ar2_static_standardize(
(shared_centered - Xc * (scales .* shared_coef)) / loading, phi1, phi2, sigma);
shared_residual -= v * (dot_product(v, shared_residual) / precision_mean);
vector[N] residual = ar2_static_standardize(
(predictor - mean_predictor - Xc * (scales .* base.2)) / loading, phi1, phi2, sigma);
residual -= v * (dot_product(v, residual) / precision_mean);
real mean_shared_residual = mean_shared - mean_design * shared_coef;
real mean_residual = mean_predictor - mean_design * base.2;
real shared_precision = 1 + dot_self(shared_residual) + square(mean_shared_residual)
+ dot_self(shared_coef);
real shared_mean = (dot_product(residual, shared_residual)
+ mean_residual * mean_shared_residual
+ dot_product(base.2, shared_coef)) / shared_precision;
vector[P + 1] conditional_mean = append_row(base.2 - shared_coef * shared_mean, shared_mean);
matrix[P + 1, P + 1] L = rep_matrix(0, P + 1, P + 1);
L[1:P, 1:P] = base.3;
L[P + 1, 1:P] = L_cross';
L[P + 1, P + 1] = sqrt(shared_precision);
real residual_ss = dot_self(residual - shared_residual * shared_mean)
+ square(mean_residual - mean_shared_residual * shared_mean)
+ dot_self(conditional_mean);
real log_det_R = -N * log(sigma) - 0.5 * log(initial_variance_factor)
- 0.5 * (N - 2) * log(innovation_variance_factor);
real log_density = -0.5 * N * log(2 * pi()) + log_det_R
- (N - 1) * log(loading)
- 0.5 * (log(precision_mean) + residual_ss)
- sum(log(diagonal(L)));
vector[N] final_residual = residual - shared_residual * shared_mean;
vector[N] score = ar2_static_transpose(final_residual, phi1, phi2, sigma);
score = -(score - mean(score)) / loading
- rep_vector((mean_residual - mean_shared_residual * shared_mean) / N, N);
return (log_density, conditional_mean, L, score);
}
real gaussian_shared_regression_lpdf(vector predictor, vector shared_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_shared_regression_distribution(predictor, shared_predictor, Xc, mean_X, grams, sums,
scales, loading, phi1, phi2, sigma).1;
}
vector gaussian_shared_regression_coefficients_rng(vector predictor, vector shared_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) + 1;
tuple(real, vector[P], matrix[P, P], vector[num_elements(predictor)]) conditional = gaussian_shared_regression_distribution(
predictor, shared_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)';
}
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;
}
}
data {
int<lower=3> 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 = -log1m(alpha_pc) / sqrt(1 - U);
real lambda_tau = -log(0.5) / 0.5;
real sum_log_y = sum(log(y));
vector[N] x_location = tan(0.5 * x);
vector[N] log_y = log(y);
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_design = rep_row_vector(1.0 / N, N) * fixed_design;
matrix[N, 24] centered_design = fixed_design - rep_matrix(mean_design, N);
array[8] matrix[24, 24] design_grams = ar2_static_grams(centered_design);
matrix[4, 24] design_sums = ar2_static_sums(centered_design);
}
parameters {
vector[N] x_residual;
vector[N] y_residual;
real<lower=0> a1;
real<lower=0> a2;
real log_kappa;
real<lower=1e-4> rho;
real<lower=-1, upper=1> phi1_w;
real<lower=-1, upper=1> phi2_s;
real<lower=phi2_s-1, upper=1-phi2_s> phi1_s;
real<lower=0> sigma_s;
real<lower=0> sigma_s2;
}
transformed parameters {
real phi2_w = 0;
real<lower=0> kappa = exp(log_kappa);
vector[N] eta_x;
vector[N] effective_residual;
vector[N] eta_y;
real residual_log_jacobian;
real speed_log_jacobian;
{
vector[24] scales = append_row(1, rep_vector(a2, 23));
vector[N] score = gaussian_regression_distribution(
x_location, centered_design, mean_design, design_grams, design_sums,
scales, a1, phi1_w, phi2_w, 1).4;
tuple(vector[N], real) coordinates = ar2_observation_coordinates(
x_residual, score, phi1_w, phi2_w, 1, a1, log_kappa, 4, 2);
effective_residual = coordinates.1;
residual_log_jacobian = coordinates.2;
eta_x = x_location + exp(-0.5 * log_kappa) * effective_residual;
}
{
vector[24] scales = append_row(1, rep_vector(sigma_s2, 23));
vector[N] score = gaussian_shared_regression_distribution(
log_y, eta_x, centered_design, mean_design, design_grams, design_sums,
scales, 1, phi1_s, phi2_s, sigma_s).4;
tuple(vector[N], real) coordinates = ar2_observation_coordinates(
y_residual, score, phi1_s, phi2_s, sigma_s, 1, log(rho), 1, 0);
eta_y = log_y + coordinates.1 / sqrt(rho);
speed_log_jacobian = coordinates.2;
}
}
model {
target += log_kappa;
target += PC_vm_inf(kappa, lambda);
target += -0.5 * N * log(rho) + speed_log_jacobian;
a1 ~ exponential(lambda_tau);
a2 ~ exponential(lambda_tau);
sigma_s ~ exponential(lambda_tau);
sigma_s2 ~ exponential(lambda_tau);
rho ~ gamma(1, 0.01);
vector[24] direction_scales = append_row(1, rep_vector(a2, 23));
vector[24] speed_scales = append_row(1, rep_vector(sigma_s2, 23));
target += gaussian_regression_lpdf(eta_x | centered_design, mean_design, design_grams, design_sums,
direction_scales, a1, phi1_w, phi2_w, 1);
target += gaussian_shared_regression_lpdf(eta_y | eta_x, centered_design, mean_design, design_grams, design_sums,
speed_scales, 1, phi1_s, phi2_s, sigma_s);
target += lavm_scaled_residual_lpdf(effective_residual | x, log_kappa);
target += residual_log_jacobian;
target += N * (rho * log(rho) - lgamma(rho))
- rho * sum(eta_y)
+ (rho - 1) * sum_log_y
- rho * sum(y .* exp(-eta_y));
}
generated quantities {
real alpha;
real a0;
vector[N] w;
vector[24] w2;
real beta;
real b0;
real b1;
vector[N] s;
vector[24] s2;
{
vector[24] direction_scales = append_row(1, rep_vector(a2, 23));
vector[24] speed_scales = append_row(1, rep_vector(sigma_s2, 23));
vector[24] direction_coef = gaussian_regression_coefficients_rng(eta_x, centered_design,
mean_design, design_grams, design_sums, direction_scales, a1, phi1_w, phi2_w, 1);
vector[25] speed_coef = gaussian_shared_regression_coefficients_rng(eta_y, eta_x, centered_design,
mean_design, design_grams, design_sums, speed_scales, 1, phi1_s, phi2_s, sigma_s);
vector[N] direction_fit = fixed_design * (direction_scales .* direction_coef);
vector[N] speed_fit = fixed_design * (speed_scales .* speed_coef[1:24]) + eta_x * speed_coef[25];
alpha = direction_coef[1];
w2 = hour_basis * direction_coef[2:24];
a0 = mean(eta_x - direction_fit);
w = (eta_x - a0 - direction_fit) / a1;
beta = speed_coef[1];
b1 = speed_coef[25];
s2 = sigma_s2 * hour_basis * speed_coef[2:24];
b0 = mean(eta_y - speed_fit);
s = eta_y - b0 - speed_fit;
}
array[N] real x_rep = lavm_rng(eta_x, kappa);
array[N] real y_rep = gamma_rng(rho, rho * exp(-eta_y));
}mod <- cmdstan_model(
stan_file = stan_file,
stanc_options = list(O1 = TRUE),
cpp_options = list(
STAN_CPP_OPTIMS = TRUE,
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(
x_residual = rnorm(n, sd = 0.1),
y_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)),
rho = runif(1, min = 1.5, max = 2.5),
phi1_w = runif(1, min = 0.3, max = 0.7),
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),
sigma_s2 = runif(1, min = 0.05, max = 0.2)
)
}
output_dir <- file.path("results", "wind_csv")
dir.create(output_dir, recursive = TRUE, showWarnings = FALSE)
set.seed(520)
inits <- lapply(seq_len(4), function(i) init_fun())
fit <- mod$sample(
data = stan_data,
output_dir = output_dir,
init = inits,
seed = 520,
refresh = 100,
iter_warmup = 500,
iter_sampling = 3000,
adapt_delta = 0.95,
chains = 4,
parallel_chains = 4
)The saved fit used below contains 3,000 post-warmup draws from each of 4 chains.
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", # Shared AR(1)
"phi1_s", "phi2_s", # AR(2)
"sigma_s", "sigma_s2" # Latent-process scales
)
fit$summary(variables = structural_vars)## # A tibble: 14 × 10
## variable mean median sd mad q5 q95 rhat ess_bulk
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 a0 0.431 4.33e-1 1.05e-1 1.06e-1 2.55e-1 6.00e-1 1.00 8116.
## 2 a1 0.0322 2.95e-2 1.96e-2 5.40e-3 2.28e-2 4.62e-2 1.00 4908.
## 3 a2 0.0222 2.18e-2 3.68e-3 3.51e-3 1.70e-2 2.88e-2 1.00 19715.
## 4 alpha -0.0131 -1.31e-2 3.19e-3 3.21e-3 -1.82e-2 -7.75e-3 1.00 8117.
## 5 b0 2.84 2.83e+0 1.84e-1 1.82e-1 2.54e+0 3.14e+0 1.00 4088.
## 6 b1 0.143 1.42e-1 1.12e-1 1.12e-1 -4.08e-2 3.29e-1 1.00 5397.
## 7 beta -0.0294 -2.92e-2 5.57e-3 5.52e-3 -3.87e-2 -2.04e-2 1.00 4086.
## 8 rho 3626. 3.61e+3 3.53e+2 3.48e+2 3.08e+3 4.23e+3 1.00 1908.
## 9 log_kappa 9.21 9.19e+0 1.93e-1 1.84e-1 8.92e+0 9.55e+0 1.00 10998.
## 10 phi1_w 0.957 9.59e-1 1.79e-2 1.68e-2 9.25e-1 9.84e-1 1.00 4489.
## 11 phi1_s 1.17 1.16e+0 9.05e-2 8.86e-2 1.02e+0 1.32e+0 1.01 1190.
## 12 phi2_s -0.198 -1.97e-1 8.82e-2 8.65e-2 -3.45e-1 -5.36e-2 1.00 1270.
## 13 sigma_s 0.0541 5.18e-2 1.15e-2 7.46e-3 4.25e-2 7.28e-2 1.00 5931.
## 14 sigma_s2 0.0371 3.62e-2 6.34e-3 5.87e-3 2.82e-2 4.86e-2 1.00 11175.
## # ℹ 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")
draws$log_a1 <- log(draws$a1)
p_a0 <- plot_density(draws, "a0", expression(a[0]), alpha = 0.6) + labs(y = "Density")
p_a1 <- plot_density(draws, "log_a1", expression(log(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)
p_b0 <- plot_density(draws, "b0", expression(b[0]), alpha = 0.6) + labs(y = "Density")
p_b1 <- plot_density(draws, "b1", expression(b[1]), alpha = 0.6)
p_beta <- plot_density(draws, "beta", expression(beta), alpha = 0.6)
p_rho <- plot_density(draws, "rho", expression(rho), alpha = 0.6)
p_log_kappa <- plot_density(draws, "log_kappa", expression(log(kappa)), alpha = 0.6) + labs(y = "Density")
p_phi1_w <- plot_density(draws, "phi1_w", expression(phi[1]^w), alpha = 0.6)
p_phi1_s <- plot_density(draws, "phi1_s", expression(phi[1]^s), alpha = 0.6)
p_phi2_s <- plot_density(draws, "phi2_s", expression(phi[2]^s), alpha = 0.6)
p_sigma_s <- plot_density(draws, "sigma_s", expression(sigma[s]), alpha = 0.6) + labs(y = "Density")
p_sigma_s2 <- plot_density(draws, "sigma_s2", expression(sigma[s2]), alpha = 0.6)
(p_a0 | p_a1 | p_a2 | p_alpha) /
(p_b0 | p_b1 | p_beta | p_rho) /
(p_log_kappa | p_phi1_w | p_phi1_s | p_phi2_s) /
(p_sigma_s | p_sigma_s2 | plot_spacer() | plot_spacer()) +
plot_annotation(
title = "Posterior Marginal Distributions",
theme = theme(plot.title = element_text(face = "bold", size = 18, hjust = 0.5))
)
The marginal plots show \(a_1\) and \(\kappa\) on log scales.
The left-hand latent plots show the centred shared AR(1) direction path and additional AR(2) speed path over the January calendar; the right-hand plots show the 24 hourly effects repeated across days. Their intervals describe posterior uncertainty under the model for these climate normals.
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
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(1) Process (" * w * ")"), x = "", y = "Value") +
theme_minimal(base_size = 14) + theme(plot.title = element_text(face = "bold"))
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"))
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"))
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, ]
p_ppc_y <- ppc_dens_overlay(y, y_rep_sample) +
labs(
title = "Posterior Predictive Check: Wind Speed",
subtitle = "Gamma Distribution",
x = "Wind Vector Magnitude",
y = "Density"
) +
theme_minimal(base_size = 14) +
theme(plot.title = element_text(face = "bold"), legend.position = "none")
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
3 Conclusion
The circular model’s linear predictor provides a local linear representation of circular data. Sharing it as a covariate in another response model estimates their association while propagating uncertainty from the circular observations.