For the von Mises distribution, also known as the circular normal distribution, \(\kappa\) controls concentration around the mean direction \(\mu\). A prior stated directly on \(\kappa\) can be difficult to interpret because changes in \(\kappa\) do not correspond uniformly to changes in angular spread. Penalized-complexity (PC) priors start from a scientifically meaningful simpler model and control departure from it through an interpretable probability statement.

This tutorial develops that construction using circular-uniform and point-mass base models, calibrates the prior through mean resultant length or angular spread, and implements both choices in Stan. A basic Stan workflow is assumed; no prior familiarity with Kullback–Leibler divergence or PC priors is required.

The examples illustrate the workflow; their models and calibrations are teaching choices, not optimized analyses.

After completing the tutorial, you will be able to:

explain how KLD and a base model define a PC prior;
choose a base model and calibrate the prior on an interpretable circular scale; and
implement the prior in Stan and examine sensitivity to its calibration.

To reproduce this tutorial, run Rscript pc_prior/R/run_all.R from the repository root, then knit the document. The script saves twelve fits and their input data; knitting loads these results. Stan programs are displayed from the files in stan/. Rerun the relevant fitting script after changing a model, dataset, or calibration.

1 Background

Companion article. The derivations follow Ye, Van Niekerk, and Rue (2026), “Penalizing complexity priors for Bayesian inference of circular models”, which also treats the cardioid and wrapped Cauchy distributions.

1.1 Circular data

Circular data arise whenever the measurement scale wraps around, as with wind direction, animal headings, or time of day. The values \(0\) and \(2\pi\) represent the same direction, so ordinary linear summaries can be misleading: for example, the arithmetic mean of \(1^\circ\) and \(359^\circ\) is \(180^\circ\), whereas their circular mean is \(0^\circ\).

Circular models respect this geometry through distributions such as the von Mises distribution. Their dispersion is commonly described by a concentration parameter or by the mean resultant length, and prior information should be expressed on similarly meaningful circular scales.

1.2 Kullback–Leibler Divergence

The Kullback–Leibler divergence (KLD) measures the expected log-density contrast between two distributions. Here it quantifies how far a model \(p(x \mid \xi)\) departs from a deliberately simpler base model \(p(x \mid \xi_0)\):

\[\text{KLD}\left(p(x \mid \xi) \parallel p(x \mid \xi_0)\right) = \int p(x \mid \xi) \log\left( \frac{p(x \mid \xi)}{p(x \mid \xi_0)} \right) dx.\]

The expectation is taken under \(p(x\mid\xi)\), so the order matters: in general, \(\operatorname{KLD}(p\parallel q)\neq\operatorname{KLD}(q\parallel p)\). Although KLD is not a metric, it is non-negative and equals zero when the two distributions agree. A PC prior uses it as a one-dimensional measure of departure from the base model.

library(ggplot2)
library(patchwork)

The left panel below compares a base density with a more dispersed model. The right panel shows the KLD integrand, which can be negative locally even though its integral over the full support is non-negative.

1.3 Penalized complexity priors

A flexible parameter—for example, a scale, smoothness, or circular concentration—allows a model to move away from a simpler limiting case. A PC prior makes that base model explicit and penalizes departures from it, while still allowing the data to support additional complexity.

For circular distributions, the KLD is mapped to the non-negative distance

\[ \begin{aligned} d\left(\xi\right) &= \sqrt{ \operatorname{KLD}\left( p\left(x\mid\xi\right) \,\Vert\, p\left(x\mid\xi_0\right) \right) } \\ &= \sqrt{ \int_0^{2\pi} p\left(x\mid\xi\right) \log\left\{ \frac{ p\left(x\mid\xi\right) }{ p\left(x\mid\xi_0\right) } \right\} \,dx }. \end{aligned} \]

The base model occupies the second argument of the KLD. For a singular base, such as a point mass, the distance is obtained through a rescaled limit, as in the construction below. The PC construction assigns an exponential prior to distance,

\[ \pi(d)=\lambda\exp(-\lambda d), \qquad d\geq 0,\quad \lambda>0. \]

Thus the density is largest at the base model, \(d=0\), and its logarithm decreases linearly with distance. Larger \(\lambda\) gives stronger shrinkage; smaller \(\lambda\) allows larger departures. When the available distance is bounded, its boundary probability also needs to be accounted for.

The exponential shape applies to \(d\), not necessarily to the original parameter. Transforming back to \(\xi\) introduces the Jacobian,

\[ p\left(d\right) = \lambda \exp\left\{ -\lambda d \right\} \Longleftrightarrow p\left(\xi\right) = \lambda \exp\left\{ -\lambda d\left(\xi\right) \right\} \left\lvert \frac{\partial d\left(\xi\right)}{\partial \xi}\right\rvert \]

so the induced density on \(\xi\) can have a different visual shape.

1.4 Base models for circular distributions

Two natural base models occur at opposite boundaries of many circular families:

The circular-uniform distribution represents no preferred direction and shrinks toward complete dispersion.
The point-mass limit represents complete concentration at the mean direction and shrinks toward vanishing angular dispersion.

The uniform base is natural when “no directional signal” is the simpler state; the point-mass base is useful when an almost deterministic direction is the baseline and dispersion represents additional complexity.

1.5 User-defined scaling

We calibrate \(\lambda\) through an interpretable transformation \(Q(\xi)\):

\[P_{\mathcal{PC}}\left( Q\left(\xi\right) > U \right) = \alpha\]

Here \(U\) is a meaningful threshold and \(\alpha\) is the prior probability of exceeding it. The event direction must match the monotonicity of \(Q\): a quantity that increases with dispersion corresponds to a different tail of a concentration parameter than one that increases with concentration.

For comparison, a PC prior for Gaussian precision \(\tau\) can be calibrated through the standard deviation \(Q(\tau)=1/\sqrt{\tau}\). The statement \(P(\sigma>3)=0.05\) directly expresses that only 5% prior probability is assigned to standard deviations above 3. The same principle is used below on circular scales.

1.5.1 Mean resultant length

For observations \(x_1,\ldots,x_n\), the mean resultant length is \[ \widehat\psi = \frac{1}{n}\sqrt{\left( \sum_{i=1}^{n}\cos\left(x_{i}\right) \right)^{2} + \left( \sum_{i=1}^{n}\sin\left(x_{i}\right) \right)^{2}} \in \left[0,1\right] \] Values near zero indicate little first-order directional concentration, whereas values near one indicate observations clustered around a common direction. For prior calibration we use the population counterpart \(\psi(\xi)\) and set \(P_{\mathcal{PC}}\{\psi(\xi)>U\}=\alpha\).

1.5.2 Circle partitioning

Alternatively, a decreasing \(Q(\xi)\in(0,2\pi]\) can provide an angular scale for calibration. For the von Mises model we use \(Q(\kappa)=2\pi/(1+\kappa)\), which decreases from \(2\pi\) at the uniform model to zero at the point-mass limit. This is a chosen index of spread, rather than an interval containing a specified probability of observations.

2 PC priors for the von Mises distribution

2.1 Von Mises distribution

The von Mises distribution is a standard unimodal model on the circle. Its density is

\[ p_{\mathcal{VM}}\left(x \mid \mu, \kappa\right) = \frac{1}{2\pi \mathcal{I}_{0} \left(\kappa\right)} \exp \left\{ \kappa \cos \left( x - \mu \right) \right\}, \quad x,\mu \in \left[-\pi,\pi\right), \quad \kappa \in \left[0,\infty\right), \] where \(\mu\) is the mean direction, \(\mathcal{I}_{a}(\cdot)\) is the modified Bessel function of the first kind, and \(\kappa\) controls concentration. At \(\kappa=0\) the density is circular uniform and \(\mu\) is not identified. As \(\kappa\) increases, observations cluster more tightly around \(\mu\); as \(\kappa\rightarrow\infty\), the distribution approaches a point mass at \(\mu\).

library(circular)

n <- 30
set.seed(123)
x1 <- rvonmises(n = n, mu = pi/4, kappa = 0.5)
x2 <- rvonmises(n = n, mu = pi/4, kappa = 3)
x3 <- rvonmises(n = n, mu = pi/4, kappa = 10)

col1 <- "#77AFA9"
col2 <- "#DD9871"
col3 <- "#839BB2" 
col1_t <- adjustcolor(col1, alpha.f = 0.6)
col2_t <- adjustcolor(col2, alpha.f = 0.6)
col3_t <- adjustcolor(col3, alpha.f = 0.6)

par(mfrow = c(1, 4), pty = "s", mar = c(0.5, 0.5, 2, 0.5))
plot.circular(x1, stack = TRUE, shrink = 1, col = col1, main = expression(kappa == 0.5))
plot.circular(x2, stack = TRUE, shrink = 1, col = col2, main = expression(kappa == 3))
plot.circular(x3, stack = TRUE, shrink = 1, col = col3, main = expression(kappa == 10))
rose.diag(x1, bins = 16, main = "Rose Diagram", col = col1_t, border = col1, prop = 1.5)
rose.diag(x2, bins = 16, add = TRUE, col = col2_t, border = col2, prop = 1.5)
rose.diag(x3, bins = 16, add = TRUE, col = col3_t, border = col3, prop = 1.5)

The first three panels show the individual samples and the final panel overlays their rose diagrams. Higher values of \(\kappa\) produce stronger concentration around the common mean direction \(\mu=\pi/4\).

2.2 PC priors

Write

\[ A(\kappa)=\frac{\mathcal I_1(\kappa)}{\mathcal I_0(\kappa)}. \]

For the circular-uniform base, the distance is

\[ d_0(\kappa)=\sqrt{\kappa A(\kappa)-\log\{\mathcal I_0(\kappa)\}}, \]

which is zero at \(\kappa=0\) and increases with concentration. For the point-mass base, a rescaled limit of the KLD gives

\[ d_\infty(\kappa)=\sqrt{1-A(\kappa)}, \]

which approaches zero as \(\kappa\rightarrow\infty\). Substituting either distance into \(\lambda\exp\{-\lambda d(\kappa)\}|d'(\kappa)|\) gives the expressions below.

Proposition 1

The PC prior for \(\kappa\) with circular uniform base model (\(\kappa_{0}=0\)) has density

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

Proposition 2

The continuous component of the PC prior for \(\kappa\) with point mass base model (\(\kappa_{0} \rightarrow \infty\)) has density, for \(\kappa>0\),

\[ 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)}}}, \quad \lambda>0. \]

The two expressions use the same construction but shrink in opposite directions. Under the circular-uniform base, increasing \(\kappa\) increases the distance from the base. Under the point-mass base, increasing \(\kappa\) decreases that distance because the base is reached as \(\kappa\rightarrow\infty\). This reversal is important when interpreting both \(\lambda\) and the calibration probability.

The figure shows how the base model changes the meaning of \(\lambda\). For the uniform base, larger \(\lambda\) places more mass near \(\kappa=0\). For the point-mass base, larger \(\lambda\) favors higher concentration, so probability can move beyond the displayed range of \(\kappa\).

For the point-mass construction, distance decreases from 1 at \(\kappa=0\) to 0 as \(\kappa\rightarrow\infty\). The full prior has CDF \(F_\infty(\kappa)=\exp\{-\lambda d_\infty(\kappa)\}\) and probability \(\exp(-\lambda)\) at \(\kappa=0\). Its continuous component, shown in the curves, integrates to \(1-\exp(-\lambda)\).

The Stan examples use this prior conditional on \(\kappa>0\). Its normalizing factor depends only on the fixed calibration rate \(\lambda\), so it can be omitted from the Stan target. The calibration equations below describe the full prior; conditioning on positive concentration changes those tail probabilities.

2.3 User-defined scaling \(\lambda\)

2.3.1 Through mean resultant length

For a von Mises distribution, the population mean resultant length is

\[ \psi(\kappa)=\frac{\mathcal I_1(\kappa)}{\mathcal I_0(\kappa)}. \]

Because \(\psi(\kappa)\) increases monotonically from 0 to 1, we calibrate \(\lambda\) through

\[P\left( \psi(\kappa) > U \right) = \alpha,\]

This states that the full prior probability of a population mean resultant length above \(U\) is \(\alpha\). For \(U,\alpha\in(0,1)\), the two base models give \[ \lambda = \begin{cases} -\frac{\log(\alpha)}{\sqrt{\kappa_U U - \log \left( \mathcal{I}_{0} (\kappa_U) \right)}} \quad &\text{for the base model at } \kappa_{0} = 0; \\[8pt] -\frac{\log(1 - \alpha)}{\sqrt{1 - U}} \quad &\text{for the base model at } \kappa_{0} \to \infty, \end{cases} \] where \(\kappa_U\) solves \(\psi(\kappa_U)=U\). The logarithm changes from \(\log(\alpha)\) to \(\log(1-\alpha)\) because the distance increases with \(\kappa\) under the uniform base but decreases under the point-mass base. The helper functions solve for \(\kappa_U\) and use scaled Bessel functions for numerical stability.

lambda.vm0.mrl <- function(u, alpha) {
  obj_fn <- function(kappa) {
    mrl <- besselI(kappa, 1, expon.scaled = TRUE) / besselI(kappa, 0, expon.scaled = TRUE)
    return((mrl - u)^2)
  }
  opt <- optim(par = 2 * u, fn = obj_fn, method = "L-BFGS-B", lower = 1e-6)
  kappa_U <- opt$par
  log_I0 <- log(besselI(kappa_U, 0, expon.scaled = TRUE)) + kappa_U
  lambda <- -log(alpha) / sqrt((kappa_U * u) - log_I0)
  return(lambda)
}
lambda.vminf.mrl <- function(u, alpha) {
  -log1p(-alpha) / sqrt(1 - u)
}

2.3.2 Through circle partitioning

The transformation

\[ Q(\kappa)=\frac{2\pi}{1+\kappa}\in(0,2\pi] \]

decreases from the circular-uniform model to the point-mass limit. It summarizes angular spread; the von Mises density still has support on the full circle. For \(0<U<2\pi\),

\[ P\left( Q\left(\kappa\right) > U \right) = P\left( \frac{2\pi}{1 + \kappa} > U \right) = P\left( \kappa < \frac{2\pi}{U} - 1 \right) = F\left( \frac{2\pi}{U} - 1 \right) = \alpha. \]

Let \(\kappa_U=2\pi/U-1\). Solving this CDF equation gives \[ \begin{aligned} \lambda = \begin{cases} -\frac{\log\left(1-\alpha\right)}{\sqrt{\frac{ \left(\frac{2\pi}{U} - 1\right) \mathcal{I}_{1}( \frac{2\pi}{U} - 1)}{\mathcal{I}_{0}( \frac{2\pi}{U} - 1)} - \log \left( \mathcal{I}_{0} \left( \frac{2\pi}{U} - 1\right) \right)}} \quad &\text{for base model at} \quad \kappa_{0} = 0; \\ -\frac{\log\left(\alpha\right)}{\sqrt{1 - \frac{\mathcal{I}_{1}(\frac{2\pi}{U} - 1)}{\mathcal{I}_{0}(\frac{2\pi}{U} - 1)}}} \quad &\text{for base model at} \quad \kappa_{0} \rightarrow \infty. \end{cases} \end{aligned} \]

For the point-mass base, the \(\log(\alpha)\) term follows from \(F_\infty(\kappa)=\exp\{-\lambda d_\infty(\kappa)\}\). Reporting \((U,\alpha)\) alongside the fitted model is more informative than reporting the resulting numerical \(\lambda\) alone.

lambda.vm0.cp <- function(u, alpha) {
  k_hat <- (2 * pi / u) - 1
  I0_scaled <- besselI(k_hat, 0, expon.scaled = TRUE)
  I1_scaled <- besselI(k_hat, 1, expon.scaled = TRUE)
  ratio <- I1_scaled / I0_scaled
  log_I0 <- log(I0_scaled) + k_hat
  d_squared <- (k_hat * ratio) - log_I0
  -log1p(-alpha) / sqrt(d_squared)
}
lambda.vminf.cp <- function(u, alpha) {
  k_hat <- (2 * pi / u) - 1
  I0_scaled <- besselI(k_hat, 0, expon.scaled = TRUE)
  I1_scaled <- besselI(k_hat, 1, expon.scaled = TRUE)
  ratio <- I1_scaled / I0_scaled
  -log(alpha) / sqrt(1 - ratio)
}

An INLA implementation is available in INLAcircular.

2.4 Simulated example (a common calibration)

We now fit the model to the three datasets generated with \(\kappa=0.5\), \(3\), and \(10\), representing low, moderate, and high concentration. The known generating values provide a reference for comparing posterior uncertainty across these three samples.

PC_vm_inf(kappa, lambda) gives the log PC-prior contribution for the point-mass base. The circular-uniform counterpart, PC_vm_0(kappa, lambda), is used in the wind example.

The functions evaluate Bessel-function ratios in log space and use asymptotic approximations for \(\kappa\geq10{,}000\).

library(cmdstanr)

stan_file <- file.path(pc_prior_dir, "stan", "pc_vminf_simulation.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 vm_log_densities(vector x, real mu, real kappa) {
    if (kappa < 10000) {
      return kappa * cos(x - mu) - log(2 * pi())
             - log_modified_bessel_first_kind(0, kappa);
    } else {
      return -2 * (kappa * square(sin(0.5 * (x - mu))))
             + 0.5 * (log(kappa) - log(2 * pi())) - 0.125 / kappa;
    }
  }

  real vm_lpdf(vector x, real mu, real kappa) {
    return sum(vm_log_densities(x, mu, kappa));
  }
}

data {
  int<lower=0> N;
  vector[N] x;
  real<lower=0> lambda;
}

parameters {
  real<lower=-pi(), upper=pi()> mu;
  real<lower=0> kappa;
}

model {
  target += -log(2 * pi());
  target += PC_vm_inf(kappa, lambda);
  x ~ vm(mu, kappa);
}

generated quantities {
  vector[N] log_lik = vm_log_densities(x, mu, kappa);
  array[N] real x_rep;

  for (n in 1:N) {
    x_rep[n] = von_mises_rng(mu, kappa);
  }
}

mod <- cmdstan_model(stan_file)

We begin with the full-prior calibration \(P\{\psi(\kappa)>0.5\}=0.5\), giving \(\lambda\approx0.98\). This provides a common reference across the three datasets. In the positive-concentration prior used by Stan, the corresponding probability is \(0.5/\{1-\exp(-\lambda)\}\approx0.80\).

lambda_val <- lambda.vminf.mrl(u = 0.5, alpha = 0.5)
lambda_val

## [1] 0.9802581

fit1 <- mod$sample(data = list(N = n, x = as.numeric(x1), lambda = lambda_val), seed = 12301, chains = 4, parallel_chains = 4, refresh = 0)
fit2 <- mod$sample(data = list(N = n, x = as.numeric(x2), lambda = lambda_val), seed = 12302, chains = 4, parallel_chains = 4, refresh = 0)
fit3 <- mod$sample(data = list(N = n, x = as.numeric(x3), lambda = lambda_val), seed = 12303, chains = 4, parallel_chains = 4, refresh = 0)

fit1$summary(variables = c("mu", "kappa"))
fit2$summary(variables = c("mu", "kappa"))
fit3$summary(variables = c("mu", "kappa"))

## # A tibble: 2 × 10
##   variable  mean median    sd   mad      q5   q95  rhat ess_bulk ess_tail
##   <chr>    <dbl>  <dbl> <dbl> <dbl>   <dbl> <dbl> <dbl>    <dbl>    <dbl>
## 1 mu       0.839  0.915 0.780 0.601 -0.512  1.92   1.00    1342.     905.
## 2 kappa    0.470  0.458 0.260 0.274  0.0619 0.921  1.00     712.     750.

## # A tibble: 2 × 10
##   variable  mean median    sd   mad    q5   q95  rhat ess_bulk ess_tail
##   <chr>    <dbl>  <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>    <dbl>    <dbl>
## 1 mu       0.828  0.828 0.115 0.112 0.641  1.02  1.00    2881.    2595.
## 2 kappa    3.27   3.21  0.754 0.723 2.16   4.62  1.00    2902.    2266.

## # A tibble: 2 × 10
##   variable  mean median     sd    mad    q5    q95  rhat ess_bulk ess_tail
##   <chr>    <dbl>  <dbl>  <dbl>  <dbl> <dbl>  <dbl> <dbl>    <dbl>    <dbl>
## 1 mu       0.763  0.763 0.0692 0.0678 0.647  0.878  1.00    3083.    2536.
## 2 kappa    8.24   8.02  2.13   2.12   5.15  12.0    1.00    2912.    2670.

library(bayesplot)

draws1 <- fit1$draws(variables = c("kappa", "mu"), format = "df")
draws2 <- fit2$draws(variables = c("kappa", "mu"), format = "df")
draws3 <- fit3$draws(variables = c("kappa", "mu"), format = "df")

plot_k <- function(draws, true_k, max_x, current_lambda) {
  prior_x <- seq(0.001, max_x, length.out = 500)
  prior_y <- sapply(prior_x, function(x) dpc.vminf(x, lambda = current_lambda))
  prior_df <- data.frame(kappa = prior_x, density = prior_y)
  
  ggplot() +
    geom_area(data = prior_df, aes(x = kappa, y = density, fill = "Prior")) +
    geom_density(data = draws, aes(x = kappa, fill = "Posterior"), color = "white") +
    geom_vline(xintercept = true_k, linetype = "dashed", color = "black", linewidth = 0.6) +
    scale_fill_manual(values = c("Prior" = col1_t, "Posterior" = col2_t)) +
    coord_cartesian(xlim = c(0, max_x)) +
    labs(x = expression(kappa), y = "Density") +
    theme_minimal() +
    theme(plot.title = element_text(face = "bold", hjust = 0.5, size = 11))
}

plot_mu <- function(draws, true_mu) {
  ggplot() +
    geom_rect(data = data.frame(1), aes(xmin = -pi, xmax = pi, ymin = 0, ymax = 1/(2*pi), fill = "Prior")) +
    geom_density(data = draws, aes(x = mu, fill = "Posterior"), color = "white") +
    geom_vline(xintercept = true_mu, linetype = "dashed", color = "black", linewidth = 0.6) +
    scale_fill_manual(values = c("Prior" = col1_t, "Posterior" = col2_t)) +
    coord_cartesian(xlim = c(0.2, 1.4)) +
    labs(x = expression(mu), y = "") +
    theme_minimal() +
    theme(plot.title = element_text(face = "bold", hjust = 0.5, size = 11))
}

pk1 <- plot_k(draws1, 0.5, 3, lambda_val) + labs(title = expression("Study 1: True "~kappa == 0.5))
pmu1 <- plot_mu(draws1, pi/4) + labs(title = expression("Study 1: True "~mu == pi/4)) + guides(fill = "none")

pk2 <- plot_k(draws2, 3, 6, lambda_val) + labs(title = expression("Study 2: True "~kappa == 3)) + guides(fill = "none")
pmu2 <- plot_mu(draws2, pi/4) + labs(title = expression("Study 2: True "~mu == pi/4)) + guides(fill = "none")

pk3 <- plot_k(draws3, 10, 15, lambda_val) + labs(title = expression("Study 3: True "~kappa == 10)) + guides(fill = "none")
pmu3 <- plot_mu(draws3, pi/4) + labs(title = expression("Study 3: True "~mu == pi/4)) + guides(fill = "none")

(pk1 | pmu1) /
(pk2 | pmu2) /
(pk3 | pmu3) +
  plot_annotation(
    title = "Prior vs. Posterior Distributions",
    theme = theme(plot.title = element_text(face = "bold", size = 15, hjust = 0.5))
  ) +
  plot_layout(guides = "collect") & 
  theme(legend.position = "top",
        legend.direction = "horizontal",
        legend.text = element_text(size = 12),
        legend.title = element_blank())

The overlays compare the continuous prior component with the posterior. When \(\kappa=0.5\), weak concentration leaves substantial uncertainty about \(\mu\); the summary intervals above show its full extent beyond the common plotting window. The higher-concentration samples give narrower intervals for \(\mu\), and all three generating values of \(\kappa\) lie within their respective 90% posterior intervals.

2.5 Simulated example (alternative calibrations)

We next compare three calibration choices, one for each simulated dataset. They illustrate how changing \(\lambda\) affects the fit; in an application, \((U,\alpha)\) should come from substantive prior knowledge.

# P(mean resultant length > 0.1) = 0.1
lambda4 <- lambda.vminf.mrl(u = 0.1, alpha = 0.1)
# P(mean resultant length > 0.5) = 0.5
lambda5 <- lambda.vminf.mrl(u = 0.5, alpha = 0.5)
# P(mean resultant length > 0.9) = 0.9
lambda6 <- lambda.vminf.mrl(u = 0.9, alpha = 0.9)
lambda4; lambda5; lambda6

## [1] 0.1110597

## [1] 0.9802581

## [1] 7.281413

fit4 <- mod$sample(data = list(N = n, x = as.numeric(x1), lambda = lambda4), seed = 12304, chains = 4, parallel_chains = 4, refresh = 0)
fit5 <- mod$sample(data = list(N = n, x = as.numeric(x2), lambda = lambda5), seed = 12305, chains = 4, parallel_chains = 4, refresh = 0)
fit6 <- mod$sample(data = list(N = n, x = as.numeric(x3), lambda = lambda6), seed = 12306, chains = 4, parallel_chains = 4, refresh = 0)

fit4$summary(variables = c("mu", "kappa"))
fit5$summary(variables = c("mu", "kappa"))
fit6$summary(variables = c("mu", "kappa"))

## # A tibble: 2 × 10
##   variable  mean median    sd   mad      q5   q95  rhat ess_bulk ess_tail
##   <chr>    <dbl>  <dbl> <dbl> <dbl>   <dbl> <dbl> <dbl>    <dbl>    <dbl>
## 1 mu       0.847  0.922 0.891 0.638 -0.667  2.06   1.00    1700.     957.
## 2 kappa    0.448  0.428 0.264 0.281  0.0539 0.914  1.01     810.     707.

## # A tibble: 2 × 10
##   variable  mean median    sd   mad    q5   q95  rhat ess_bulk ess_tail
##   <chr>    <dbl>  <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>    <dbl>    <dbl>
## 1 mu       0.829  0.828 0.114 0.109 0.642  1.02  1.00    2921.    2432.
## 2 kappa    3.27   3.20  0.771 0.767 2.13   4.66  1.00    2377.    2171.

## # A tibble: 2 × 10
##   variable  mean median     sd    mad    q5    q95  rhat ess_bulk ess_tail
##   <chr>    <dbl>  <dbl>  <dbl>  <dbl> <dbl>  <dbl> <dbl>    <dbl>    <dbl>
## 1 mu       0.761  0.761 0.0656 0.0657 0.652  0.866  1.00    3521.    2269.
## 2 kappa    8.69   8.50  2.11   2.09   5.59  12.5    1.00    3609.    2748.

draws4 <- fit4$draws(variables = c("kappa", "mu"), format = "df")
draws5 <- fit5$draws(variables = c("kappa", "mu"), format = "df")
draws6 <- fit6$draws(variables = c("kappa", "mu"), format = "df")

pk4 <- plot_k(draws4, 0.5, 3, lambda4) + labs(title = expression("Study 1: True "~kappa == 0.5))
pmu4 <- plot_mu(draws4, pi/4) + labs(title = expression("Study 1: True "~mu == pi/4)) + guides(fill = "none")

pk5 <- plot_k(draws5, 3, 6, lambda5) + labs(title = expression("Study 2: True "~kappa == 3)) + guides(fill = "none")
pmu5 <- plot_mu(draws5, pi/4) + labs(title = expression("Study 2: True "~mu == pi/4)) + guides(fill = "none")

pk6 <- plot_k(draws6, 10, 15, lambda6) + labs(title = expression("Study 3: True "~kappa == 10)) + guides(fill = "none")
pmu6 <- plot_mu(draws6, pi/4) + labs(title = expression("Study 3: True "~mu == pi/4)) + guides(fill = "none")

wrap_plots(pk4, pmu4, pk5, pmu5, pk6, pmu6, ncol = 2) +
  plot_annotation(
    title = "Prior vs. Posterior Distributions",
    theme = theme(plot.title = element_text(face = "bold", size = 15, hjust = 0.5))
  ) +
  plot_layout(guides = "collect") & 
  theme(legend.position = "top",
        legend.direction = "horizontal",
        legend.text = element_text(size = 12),
        legend.title = element_blank())

For the full point-mass-base prior, the three calibrations imply increasing preference for concentration:

\(P\{\psi(\kappa)>0.1\}=0.1\) permits substantial dispersion;
\(P\{\psi(\kappa)>0.5\}=0.5\) is the midpoint calibration; and
\(P\{\psi(\kappa)>0.9\}=0.9\) strongly favors concentrated distributions.

The Stan fits use the corresponding positive-concentration components. Their ordering in \(\lambda\) is preserved, while the probabilities listed above include the boundary mass at \(\kappa=0\).

2.6 Application: wind data

We finish with the 310 wind directions in circular::wind, recorded at Col de la Roa in the Italian Alps. The circular plot and rose diagram show a preferred direction together with appreciable spread.

data("wind", package = "circular")
ys <- wind
N <- length(ys)

To examine prior sensitivity, we fit both base-model choices with angular thresholds \(U\in\{\pi,\pi/2,\pi/4\}\) and \(\alpha=0.5\). Smaller \(U\) expresses a tighter prior scale for angular spread.

u_vals <- c(pi, pi/2, pi/4)
alpha_val <- 0.5

lambdas_vm0 <- sapply(u_vals, lambda.vm0.cp, alpha = alpha_val)
lambdas_vminf <- sapply(u_vals, lambda.vminf.cp, alpha = alpha_val)

data.frame(
  U_Value = c("pi", "pi/2", "pi/4"),
  Lambda_vm0 = lambdas_vm0,
  Lambda_vminf = lambdas_vminf
)

##   U_Value Lambda_vm0 Lambda_vminf
## 1      pi  1.5108612    0.9315871
## 2    pi/2  0.7542017    1.5901272
## 3    pi/4  0.5962940    2.5400439

The resulting rates are not directly comparable across base models because their distance functions differ. The pair \((U,\alpha)\) describes the calibration before conditioning the point-mass-base prior on \(\kappa>0\). The two Stan programs differ only in their prior contributions: PC_vm_0(kappa, lambda) uses \(d_0(\kappa)\) and PC_vm_inf(kappa, lambda) uses \(d_\infty(\kappa)\).

The circular-uniform-base program is:

file_vm0 <- file.path(pc_prior_dir, "stan", "pc_vm0_wind_application.stan")

functions {
  real PC_vm_0(real kappa, real lambda) {
    real 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;
      d_sq = kappa * exp(log_Q) - log_I0;
      log_numerator = log(kappa) + log(numerator);
    } else {
      real r = inv(kappa);
      d_sq = 0.5 * (log(2 * pi()) + log(kappa) - 1) - 0.25 * r;
      log_numerator = -log(2.0) - log(kappa) + log1p(0.5 * r);
    }

    return log(lambda) + log_numerator - log(2.0) - 0.5 * log(d_sq)
           - lambda * sqrt(d_sq);
  }

  vector vm_log_densities(vector x, real mu, real kappa) {
    if (kappa < 10000) {
      return kappa * cos(x - mu) - log(2 * pi())
             - log_modified_bessel_first_kind(0, kappa);
    } else {
      return -2 * (kappa * square(sin(0.5 * (x - mu))))
             + 0.5 * (log(kappa) - log(2 * pi())) - 0.125 / kappa;
    }
  }

  real vm_lpdf(vector x, real mu, real kappa) {
    return sum(vm_log_densities(x, mu, kappa));
  }
}

data {
  int<lower=0> N;
  vector[N] x;
  real<lower=0> lambda;
}

parameters {
  real<lower=-pi(), upper=pi()> mu;
  real<lower=0> kappa;
}

model {
  target += -log(2 * pi());
  target += PC_vm_0(kappa, lambda);
  x ~ vm(mu, kappa);
}

generated quantities {
  vector[N] log_lik = vm_log_densities(x, mu, kappa);
  array[N] real x_rep;

  for (n in 1:N) {
    x_rep[n] = von_mises_rng(mu, kappa);
  }
}

The point-mass-base program is:

file_vminf <- file.path(pc_prior_dir, "stan", "pc_vminf_wind_application.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 vm_log_densities(vector x, real mu, real kappa) {
    if (kappa < 10000) {
      return kappa * cos(x - mu) - log(2 * pi())
             - log_modified_bessel_first_kind(0, kappa);
    } else {
      return -2 * (kappa * square(sin(0.5 * (x - mu))))
             + 0.5 * (log(kappa) - log(2 * pi())) - 0.125 / kappa;
    }
  }

  real vm_lpdf(vector x, real mu, real kappa) {
    return sum(vm_log_densities(x, mu, kappa));
  }
}

data {
  int<lower=0> N;
  vector[N] x;
  real<lower=0> lambda;
}

parameters {
  real<lower=-pi(), upper=pi()> mu;
  real<lower=0> kappa;
}

model {
  target += -log(2 * pi());
  target += PC_vm_inf(kappa, lambda);
  x ~ vm(mu, kappa);
}

generated quantities {
  vector[N] log_lik = vm_log_densities(x, mu, kappa);
  array[N] real x_rep;

  for (n in 1:N) {
    x_rep[n] = von_mises_rng(mu, kappa);
  }
}

mod_vm0 <- cmdstan_model(file_vm0)
mod_vminf <- cmdstan_model(file_vminf)

Each model is fitted at all three calibrations. This six-model grid separates sensitivity to the base model from sensitivity to \(U\) while holding the likelihood and data fixed.

fit_vm0_pi  <- mod_vm0$sample(data = list(N = N, x = as.numeric(ys), lambda = lambdas_vm0[1]), seed = 12401, chains = 4, parallel_chains = 4, refresh = 0)
fit_vm0_pi2 <- mod_vm0$sample(data = list(N = N, x = as.numeric(ys), lambda = lambdas_vm0[2]), seed = 12402, chains = 4, parallel_chains = 4, refresh = 0)
fit_vm0_pi4 <- mod_vm0$sample(data = list(N = N, x = as.numeric(ys), lambda = lambdas_vm0[3]), seed = 12403, chains = 4, parallel_chains = 4, refresh = 0)

fit_vminf_pi  <- mod_vminf$sample(data = list(N = N, x = as.numeric(ys), lambda = lambdas_vminf[1]), seed = 12404, chains = 4, parallel_chains = 4, refresh = 0)
fit_vminf_pi2 <- mod_vminf$sample(data = list(N = N, x = as.numeric(ys), lambda = lambdas_vminf[2]), seed = 12405, chains = 4, parallel_chains = 4, refresh = 0)
fit_vminf_pi4 <- mod_vminf$sample(data = list(N = N, x = as.numeric(ys), lambda = lambdas_vminf[3]), seed = 12406, chains = 4, parallel_chains = 4, refresh = 0)

fit_vm0_pi$summary(variables = c("mu", "kappa"))
fit_vm0_pi2$summary(variables = c("mu", "kappa"))
fit_vm0_pi4$summary(variables = c("mu", "kappa"))

fit_vminf_pi$summary(variables = c("mu", "kappa"))
fit_vminf_pi2$summary(variables = c("mu", "kappa"))
fit_vminf_pi4$summary(variables = c("mu", "kappa"))

## # A tibble: 2 × 10
##   variable  mean median     sd    mad    q5   q95  rhat ess_bulk ess_tail
##   <chr>    <dbl>  <dbl>  <dbl>  <dbl> <dbl> <dbl> <dbl>    <dbl>    <dbl>
## 1 mu       0.291  0.293 0.0537 0.0565 0.203 0.377  1.00    3438.    2550.
## 2 kappa    1.75   1.75  0.126  0.127  1.55  1.96   1.00    2969.    2421.

## # A tibble: 2 × 10
##   variable  mean median     sd    mad    q5   q95  rhat ess_bulk ess_tail
##   <chr>    <dbl>  <dbl>  <dbl>  <dbl> <dbl> <dbl> <dbl>    <dbl>    <dbl>
## 1 mu       0.293  0.292 0.0533 0.0530 0.207 0.381  1.00    3735.    2581.
## 2 kappa    1.76   1.75  0.129  0.124  1.55  1.97   1.00    4111.    2497.

## # A tibble: 2 × 10
##   variable  mean median     sd    mad    q5   q95  rhat ess_bulk ess_tail
##   <chr>    <dbl>  <dbl>  <dbl>  <dbl> <dbl> <dbl> <dbl>    <dbl>    <dbl>
## 1 mu       0.292  0.292 0.0529 0.0528 0.206 0.378  1.00    3402.    2501.
## 2 kappa    1.76   1.75  0.124  0.129  1.56  1.96   1.00    3289.    2617.

## # A tibble: 2 × 10
##   variable  mean median     sd    mad    q5   q95  rhat ess_bulk ess_tail
##   <chr>    <dbl>  <dbl>  <dbl>  <dbl> <dbl> <dbl> <dbl>    <dbl>    <dbl>
## 1 mu       0.290  0.290 0.0533 0.0517 0.203 0.379  1.00    3102.    2482.
## 2 kappa    1.77   1.76  0.128  0.125  1.56  1.98   1.00    2643.    2225.

## # A tibble: 2 × 10
##   variable  mean median     sd    mad    q5   q95  rhat ess_bulk ess_tail
##   <chr>    <dbl>  <dbl>  <dbl>  <dbl> <dbl> <dbl> <dbl>    <dbl>    <dbl>
## 1 mu       0.292  0.292 0.0528 0.0532 0.206 0.381 1.00     2771.    2506.
## 2 kappa    1.77   1.77  0.127  0.126  1.56  1.98  1.000    3530.    2727.

## # A tibble: 2 × 10
##   variable  mean median     sd    mad    q5   q95  rhat ess_bulk ess_tail
##   <chr>    <dbl>  <dbl>  <dbl>  <dbl> <dbl> <dbl> <dbl>    <dbl>    <dbl>
## 1 mu       0.292  0.291 0.0511 0.0510 0.206 0.375 1.00     3491.    2765.
## 2 kappa    1.77   1.76  0.128  0.128  1.57  1.98  1.000    3612.    2632.

draws1 <- fit_vm0_pi$draws(variables = c("kappa", "mu"), format = "df")
draws2 <- fit_vm0_pi2$draws(variables = c("kappa", "mu"), format = "df")
draws3 <- fit_vm0_pi4$draws(variables = c("kappa", "mu"), format = "df")

posteriors_vm0 <- data.frame(
  Kappa_Estimate = c(draws1$kappa, draws2$kappa, draws3$kappa),
  Mu_Estimate = c(draws1$mu, draws2$mu, draws3$mu),
  Prior = factor(rep(c("U_pi", "U_pi2", "U_pi4"), each = nrow(draws1)),
                 levels = c("U_pi", "U_pi2", "U_pi4"))
)

draws4 <- fit_vminf_pi$draws(variables = c("kappa", "mu"), format = "df")
draws5 <- fit_vminf_pi2$draws(variables = c("kappa", "mu"), format = "df")
draws6 <- fit_vminf_pi4$draws(variables = c("kappa", "mu"), format = "df")

posteriors_vminf <- data.frame(
  Kappa_Estimate = c(draws4$kappa, draws5$kappa, draws6$kappa),
  Mu_Estimate = c(draws4$mu, draws5$mu, draws6$mu),
  Prior = factor(rep(c("U_pi", "U_pi2", "U_pi4"), each = nrow(draws4)),
                 levels = c("U_pi", "U_pi2", "U_pi4"))
)

my_colors_t <- c(col1_t, col2_t, col3_t)
sim_labels <- c(expression(U == pi), expression(U == pi/2), expression(U == pi/4))

plot_prior_density <- function(draws, variable, title, x_label = "", y_label = "") {
  ggplot(draws, aes(x = .data[[variable]], fill = Prior)) +
    geom_density(color = "white") +
    scale_fill_manual(values = my_colors_t, labels = sim_labels) +
    labs(title = title, x = x_label, y = y_label) +
    theme(plot.title = element_text(face = "bold", hjust = 0.5))
}

p_kappa_vm0 <- plot_prior_density(posteriors_vm0, "Kappa_Estimate", "Uniform Base",
                                  y_label = "Posterior Density")
p_mu_vm0 <- plot_prior_density(posteriors_vm0, "Mu_Estimate", "Uniform Base")
p_kappa_vminf <- plot_prior_density(posteriors_vminf, "Kappa_Estimate", "Point Mass",
                                    x_label = expression(kappa), y_label = "Posterior Density")
p_mu_vminf <- plot_prior_density(posteriors_vminf, "Mu_Estimate", "Point Mass",
                                 x_label = expression(mu))

(p_kappa_vm0 | p_mu_vm0) / 
(p_kappa_vminf | p_mu_vminf) + 
  plot_annotation(
    title = "Posterior Distributions",
    theme = theme(plot.title = element_text(face = "bold", size = 16, hjust = 0.5))
  ) +
  plot_layout(guides = "collect") & 
  theme(legend.position = "bottom", legend.direction = "horizontal", 
        legend.text = element_text(size = 12), legend.title = element_blank())

The six posterior distributions are nearly indistinguishable, with \(\mu\) around \(0.29\) radians and \(\kappa\) around \(1.76\). For these 310 observations, the likelihood is sufficiently informative that the tested base models and thresholds have little effect.

We also compare predictive criteria across the six prior specifications using the pointwise log-likelihood draws.

library(loo)

ll_vm0_pi  <- fit_vm0_pi$draws("log_lik", format = "matrix")
ll_vm0_pi2 <- fit_vm0_pi2$draws("log_lik", format = "matrix")
ll_vm0_pi4 <- fit_vm0_pi4$draws("log_lik", format = "matrix")
ll_vminf_pi  <- fit_vminf_pi$draws("log_lik", format = "matrix")
ll_vminf_pi2 <- fit_vminf_pi2$draws("log_lik", format = "matrix")
ll_vminf_pi4 <- fit_vminf_pi4$draws("log_lik", format = "matrix")

loo_vm0_pi  <- loo(ll_vm0_pi)
loo_vm0_pi2 <- loo(ll_vm0_pi2)
loo_vm0_pi4 <- loo(ll_vm0_pi4)
loo_vminf_pi  <- loo(ll_vminf_pi)
loo_vminf_pi2 <- loo(ll_vminf_pi2)
loo_vminf_pi4 <- loo(ll_vminf_pi4)
loo_list <- list(
  "Uni_U=pi" = loo_vm0_pi,  "Uni_U=pi/2" = loo_vm0_pi2,  "Uni_U=pi/4" = loo_vm0_pi4,
  "Pm_U=pi"  = loo_vminf_pi, "Pm_U=pi/2" = loo_vminf_pi2, "Pm_U=pi/4" = loo_vminf_pi4
)
waic_vm0_pi    <- waic(ll_vm0_pi);    waic_vm0_pi2   <- waic(ll_vm0_pi2)
waic_vm0_pi4   <- waic(ll_vm0_pi4);   waic_vminf_pi  <- waic(ll_vminf_pi)
waic_vminf_pi2 <- waic(ll_vminf_pi2); waic_vminf_pi4 <- waic(ll_vminf_pi4)
waic_list <- list(
  "Uni_U=pi" = waic_vm0_pi,  "Uni_U=pi/2" = waic_vm0_pi2,   "Uni_U=pi/4" = waic_vm0_pi4,
  "Pm_U=pi"  = waic_vminf_pi, "Pm_U=pi/2" = waic_vminf_pi2, "Pm_U=pi/4"  = waic_vminf_pi4
)

model_names <- names(loo_list)
metrics <- lapply(model_names, function(m) {
  l <- loo_list[[m]]
  w <- waic_list[[m]]
  elpd <- l$estimates["elpd_loo", "Estimate"]
  se_elpd <- l$estimates["elpd_loo", "SE"]
  looic <- l$estimates["looic", "Estimate"]
  se_looic <- l$estimates["looic", "SE"]
  waic_est <- w$estimates["waic", "Estimate"]
  se_waic <- w$estimates["waic", "SE"]
  pareto_k <- l$diagnostics$pareto_k
  k_07 <- sum(pareto_k > 0.7 & pareto_k <= 1.0)
  k_10 <- sum(pareto_k > 1.0)
  data.frame(
    Model = m, ELPD = elpd, SE_ELPD = se_elpd, LOOIC = looic, SE_LOOIC = se_looic,
    WAIC = waic_est, SE_WAIC = se_waic, k_07 = k_07, k_10 = k_10
  )
})
df <- do.call(rbind, metrics)
best_elpd <- max(df$ELPD); best_looic <- min(df$LOOIC); best_waic <- min(df$WAIC)
df$dELPD <- df$ELPD - best_elpd
df$dLOOIC <- df$LOOIC - best_looic
df$dWAIC <- df$WAIC - best_waic
comp <- as.data.frame(loo_compare(loo_list))
df$SE_dELPD <- comp[match(df$Model, rownames(comp)), "se_diff"]
df$SE_dLOOIC <- df$SE_dELPD * 2
comp_waic <- as.data.frame(loo_compare(waic_list))
df$SE_dWAIC <- comp_waic[match(df$Model, rownames(comp_waic)), "se_diff"] * 2
formatted_df <- data.frame(
  Metric = c("ELPD", "SE(ELPD)", "dELPD", "SE(dELPD)",
             "LOOIC", "SE(LOOIC)", "dLOOIC", "SE(dLOOIC)",
             "WAIC", "SE(WAIC)", "dWAIC", "SE(dWAIC)",
             "0.7 < k <= 1.0", "k > 1.0")
)
for (i in 1:nrow(df)) {
  formatted_df[[df$Model[i]]] <- c(
    sprintf("%.2f", df$ELPD[i]), sprintf("%.2f", df$SE_ELPD[i]),
    sprintf("%.2f", df$dELPD[i]), sprintf("%.2f", df$SE_dELPD[i]),
    sprintf("%.2f", df$LOOIC[i]), sprintf("%.2f", df$SE_LOOIC[i]),
    sprintf("%.2f", df$dLOOIC[i]), sprintf("%.2f", df$SE_dLOOIC[i]),
    sprintf("%.2f", df$WAIC[i]), sprintf("%.2f", df$SE_WAIC[i]),
    sprintf("%.2f", df$dWAIC[i]), sprintf("%.2f", df$SE_dWAIC[i]),
    df$k_07[i], df$k_10[i]
  )
}
print(formatted_df)

##            Metric Uni_U=pi Uni_U=pi/2 Uni_U=pi/4 Pm_U=pi Pm_U=pi/2 Pm_U=pi/4
## 1            ELPD  -419.19    -419.22    -419.12 -419.20   -419.17   -419.12
## 2        SE(ELPD)    16.56      16.59      16.62   16.70     16.72     16.72
## 3           dELPD    -0.07      -0.10       0.00   -0.08     -0.05     -0.00
## 4       SE(dELPD)     0.05       0.02       0.00    0.09      0.10      0.10
## 5           LOOIC   838.38     838.44     838.24  838.41    838.35    838.25
## 6       SE(LOOIC)    33.13      33.19      33.23   33.40     33.44     33.43
## 7          dLOOIC     0.14       0.20       0.00    0.16      0.11      0.00
## 8      SE(dLOOIC)     0.11       0.05       0.00    0.18      0.21      0.20
## 9            WAIC   838.37     838.44     838.24  838.40    838.34    838.24
## 10       SE(WAIC)    33.13      33.19      33.23   33.40     33.44     33.43
## 11          dWAIC     0.13       0.20       0.00    0.16      0.11      0.00
## 12      SE(dWAIC)     0.11       0.05       0.00    0.18      0.21      0.20
## 13 0.7 < k <= 1.0        0          0          0       0         0         0
## 14        k > 1.0        0          0          0       0         0         0

Compare predictive-score differences with their paired standard errors. Here the largest ELPD difference is about 0.10 across 310 observations, so the practical predictive differences among the tested specifications are small. The Pareto-\(k\) counts summarize the importance-sampling diagnostics for the LOO approximation.

As a direct circular posterior predictive check, we compare the observed mean resultant length with values from replicated datasets. This statistic targets overall first-order concentration.

mrl <- function(x) {
  cos_sum <- sum(cos(x))
  sin_sum <- sum(sin(x))
  sqrt(cos_sum^2 + sin_sum^2) / length(x)
}

yrep_matrix <- fit_vm0_pi$draws("x_rep", format = "matrix")

p_mrl <- ppc_stat(y = as.numeric(ys), 
                  yrep = yrep_matrix, 
                  stat = "mrl") +
  labs(title = "Mean resultant length") +
  scale_color_manual(name = "", 
                     values = c("y" = "darkred", "yrep" = "#839BB2"), 
                     labels = c("y" = "Observed", "yrep" = "Posterior predictive")) +
  scale_fill_manual(name = "", 
                    values = c("y" = "darkred", "yrep" = "#B3CDE3"), 
                    labels = c("y" = "Observed", "yrep" = "Posterior predictive")) +
  guides(color = guide_legend(title = NULL), fill = guide_legend(title = NULL)) +
  theme(legend.position = "bottom", legend.text = element_text(size = 11))

p_mrl

The observed mean resultant length lies within the replicated distribution, so the observed first-order concentration is consistent with this fitted model. Further checks could examine other features of the angular distribution.

3 Conclusion

A PC prior penalizes distance from a chosen base model, discouraging unsupported complexity and helping to limit overfitting. The base determines the direction of shrinkage, while \(\lambda\) controls its strength.

This tutorial calibrates \(\lambda\) through prior tail probabilities for either mean resultant length \(\psi(\kappa)\) or the angular-spread proxy \(Q(\kappa)=2\pi/(1+\kappa)\). A meaningful threshold \(U\) and prior probability \(\alpha\) express uncertainty about the parameter on an interpretable circular scale and provide a basis for sensitivity checks.

Tutorial: Penalized Complexity Priors for Circular Distributions

Xiang Ye

2026-09-11