Tutorial: Penalized Complexity Priors for Circular Distributions

In the circular models, sparsity is controlled by concentration parameters. For example, for the von Mises distribution, which also known as the circular normal distribution, the parameter \(\kappa\) controls the concentration of the distribution around its dominant direction \(\mu\). However, 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 address this by starting from a scientifically meaningful simpler model and controlling 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.

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.

1 Preliminary

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 flexible model. The right panel shows the pointwise integrand whose signed area is the KLD.

1.3 Penalizing Complexity (PC) Prior

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 direction of the KLD is fixed by the selected base model. The PC construction assigns an exponential prior to the resulting distance,

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

Thus the prior density is largest at the base model, \(d=0\), and decreases at a constant rate on the distance scale. Larger \(\lambda\) gives stronger shrinkage; smaller \(\lambda\) allows larger departures.

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 Model 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 choice is scientific rather than automatic. 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

Rather than selecting \(\lambda\) on an unfamiliar parameter scale, we calibrate it 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 represent angular spread. The calibration \(P_{\mathcal{PC}}\{Q(\xi)>U\}=\alpha\) assigns probability \(\alpha\) to spread beyond \(U\) around the mean direction. 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.

2 PC priors for 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)

# Reusable plotting palette
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,

\[ 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 two densities 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 PC prior for \(\kappa\) with point mass base model (\(\kappa_{0} \rightarrow \infty\)) has density

\[ 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, the continuous distance decreases from 1 at \(\kappa=0\) to 0 as \(\kappa\rightarrow\infty\). The full prior also has boundary mass at \(\kappa=0\); the Stan examples below use its continuous component.

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 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 for the PC prior with circular uniform base model through mean resultant length
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 for the PC prior with point mass base model through mean resultant length
lambda.vminf.mrl <- function(u, alpha) {-log(1-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 is an interpretable calibration of angular spread, rather than the literal support width of a von Mises density. 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 for the PC prior with circular uniform base model through circle partitioning
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
  lambda <- -log(1 - alpha) / sqrt(d_squared)
  return(lambda)
}
# lambda for the PC prior with point mass base model through circle partitioning
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
  lambda <- -log(alpha) / sqrt(1 - ratio)
  return(lambda)
}

An INLA implementation is available in INLAcircular.

2.4 Simulated example (weakly informative scaling)

We now fit the model to the three datasets generated with \(\kappa=0.5\), \(3\), and \(10\), representing low, moderate, and high concentration. Because the generating values are known, the example checks the implementation and shows how the same prior behaves at different concentration levels.

We use the point-mass base to demonstrate its Stan implementation. The code evaluates Bessel-function ratios in log space and uses asymptotic expressions for very large \(\kappa\).

library(cmdstanr)

stan_file <- file.path(pc_prior_dir, "stan", "pc_vminf_simulation.stan")

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

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

model {
  // Uniform prior for mu
  target += -log(2 * pi());

  // PC prior for kappa (Point Mass Base Model limit)
  real Q;
  real dQ;
  real log_I0_val;

  if (kappa < 10000) {
    // Calculate strictly in log-space to prevent overflow
    real log_I0 = log_modified_bessel_first_kind(0, kappa);
    real log_I1 = log_modified_bessel_first_kind(1, kappa);
    real log_I2 = log_modified_bessel_first_kind(2, kappa);

    // Safe ratios using exponentiated differences
    real I1_over_I0 = exp(log_I1 - log_I0);
    real I2_over_I0 = exp(log_I2 - log_I0);

    Q = I1_over_I0;
    dQ = 0.5 * (1.0 + I2_over_I0) - (I1_over_I0^2);
    log_I0_val = log_I0;
  } else {
    // Asymptotic expansion limits for large kappa values
    Q = 1.0 - 0.5 / kappa - 0.125 / (kappa^2);
    dQ = 0.5 / (kappa^2) + 0.25 / (kappa^3);
    log_I0_val = kappa - 0.5 * log(2 * pi() * kappa) + 1.0 / (8.0 * kappa);
  }

  real sqrt_term = 1.0 - Q;

  target += log(lambda) + log(dQ) - log(2.0) - 0.5 * log(sqrt_term) - lambda * sqrt(sqrt_term);

  // Manual Likelihood
  for (n in 1:N) {
    target += kappa * cos(x[n] - mu) - log_I0_val;
  }
}

generated quantities {
  vector[N] log_lik;
  array[N] real x_rep;

  real log_I0_gen = kappa < 10000 ? log_modified_bessel_first_kind(0, kappa)
                                   : (kappa - 0.5 * log(2 * pi() * kappa)
                                      + 1.0 / (8.0 * kappa));

  for (n in 1:N) {
    log_lik[n] = kappa * cos(x[n] - mu) - log(2 * pi()) - log_I0_gen;
    x_rep[n] = von_mises_rng(mu, kappa);
  }
}

mod <- cmdstan_model(stan_file)

We begin with the midpoint calibration \(P\{\psi(\kappa)>0.5\}=0.5\). It assigns equal prior probability to mean resultant lengths below and above 0.5 and serves as a weak-information reference for these examples.

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

## [1] 0.9802581

fit1 <- mod$sample(data = list(N = n, x = x1, lambda = lambda_val), refresh = 0)
fit2 <- mod$sample(data = list(N = n, x = x2, lambda = lambda_val), refresh = 0)
fit3 <- mod$sample(data = list(N = n, x = x3, lambda = lambda_val), refresh = 0)

# Summary for kappa = 0.5 model
fit1$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.819  0.932 0.970 0.701 -1.12   2.15   1.00    1706.     773.
## 2 kappa    0.430  0.407 0.270 0.292  0.0436 0.919  1.00     701.     619.

# Summary for kappa = 3 model
fit2$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.827  0.827 0.114 0.115 0.644  1.01  1.00    3050.    2378.
## 2 kappa    3.26   3.19  0.761 0.730 2.13   4.62  1.00    2608.    2231.

# Summary for kappa = 10 model
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.763  0.763 0.0670 0.0647 0.653  0.873  1.00    3816.    2786.
## 2 kappa    8.26   8.09  2.12   2.07   5.13  12.1    1.00    2693.    2421.

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 prior–posterior overlays show how much each dataset updates the common prior. When \(\kappa=0.5\), weak concentration leaves substantial uncertainty about \(\mu\). As concentration increases, the posterior for \(\mu\) narrows, while the posterior for \(\kappa\) remains centred near the generating values.

2.5 Simulated example (informative choices for \(\lambda\))

We next use three calibrations representing low, moderate, and high expected concentration. Each statement is matched to one simulated setting to make its effect visible; in an application, \((U,\alpha)\) should instead 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 = x1, lambda = lambda4), refresh = 0)
fit5 <- mod$sample(data = list(N = n, x = x2, lambda = lambda5), refresh = 0)
fit6 <- mod$sample(data = list(N = n, x = x3, lambda = lambda6), refresh = 0)

# Summary for kappa = 0.5 model
fit4$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.822  0.915 0.917 0.641 -0.864  2.10   1.00    1072.     686.
## 2 kappa    0.439  0.421 0.260 0.282  0.0470 0.895  1.00     768.     904.

# Summary for kappa = 3 model
fit5$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.828  0.826 0.113 0.115 0.643  1.01 1.000    3004.    2646.
## 2 kappa    3.27   3.21  0.725 0.697 2.21   4.57 1.00     2622.    2714.

# Summary for kappa = 10 model
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.763  0.764 0.0646 0.0623 0.657  0.869  1.00    2809.    2455.
## 2 kappa    8.69   8.51  2.18   2.16   5.48  12.6    1.00    2508.    2418.

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

# Prior and posterior densities for each calibration
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")

(pk4 | pmu4) /
(pk5 | pmu5) /
(pk6 | pmu6) +
  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())

Under the point-mass base, the three statements imply progressively stronger contraction toward high 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 overlays show how the likelihood updates each of these prior specifications.

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.

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.

# Prior probability of angular spread beyond U is 0.5
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 have different directions. The probability statement \((U,\alpha)\) is the meaningful comparison. The two Stan programs differ only in their prior contributions: pc_vm0_wind_application.stan uses \(d_0(\kappa)\) and pc_vminf_wind_application.stan uses \(d_\infty(\kappa)\).

The circular-uniform-base program is:

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

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

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

model {
  // Uniform prior for mu
  target += -log(2 * pi());

  real numerator;
  real d_sq;
  real log_I0_val;

  if (kappa < 10000) {
    // Calculate strictly in log-space to prevent overflow
    real log_I0 = log_modified_bessel_first_kind(0, kappa);
    real log_I1 = log_modified_bessel_first_kind(1, kappa);
    real log_I2 = log_modified_bessel_first_kind(2, kappa);

    // Safe ratios using exponentiated differences
    real I1_over_I0 = exp(log_I1 - log_I0);
    real I2_over_I0 = exp(log_I2 - log_I0);

    d_sq = kappa * I1_over_I0 - log_I0;
    numerator = kappa * (0.5 * (1.0 + I2_over_I0) - (I1_over_I0^2));
    log_I0_val = log_I0;
  } else {
    // Asymptotic expansion limits for large kappa values
    real Q = 1.0 - 0.5 / kappa - 0.125 / (kappa^2);
    real dQ = 0.5 / (kappa^2) + 0.25 / (kappa^3);
    log_I0_val = kappa - 0.5 * log(2 * pi() * kappa) + 1.0 / (8.0 * kappa);

    d_sq = kappa * Q - log_I0_val;
    numerator = kappa * dQ;
  }

  // Prevent log(0) during warmup initialization
  if (d_sq < 1e-12) d_sq = 1e-12;
  if (numerator < 1e-12) numerator = 1e-12;

  // PC prior with circular uniform base model
  target += log(lambda) + log(numerator) - log(2.0) - 0.5 * log(d_sq) - lambda * sqrt(d_sq);

  // Manual safe von Mises likelihood
  for (n in 1:N) {
    target += kappa * cos(x[n] - mu) - log_I0_val;
  }
}

generated quantities {
  vector[N] log_lik;
  array[N] real x_rep;

  real log_I0_gen = kappa < 10000 ? log_modified_bessel_first_kind(0, kappa)
                                  : (kappa - 0.5 * log(2 * pi() * kappa)
                                     + 1.0 / (8.0 * kappa));

  for (n in 1:N) {
    // Pointwise log-likelihood for LOO (stable manual formulation)
    log_lik[n] = kappa * cos(x[n] - mu) - log(2 * pi()) - log_I0_gen;
    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")

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

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

model {
  // Uniform prior for mu
  target += -log(2 * pi());

  real numerator;
  real d_sq;
  real log_I0_val;

  if (kappa < 10000) {
    // Calculate strictly in log-space to prevent overflow
    real log_I0 = log_modified_bessel_first_kind(0, kappa);
    real log_I1 = log_modified_bessel_first_kind(1, kappa);
    real log_I2 = log_modified_bessel_first_kind(2, kappa);

    // Safe ratios using exponentiated differences
    real I1_over_I0 = exp(log_I1 - log_I0);
    real I2_over_I0 = exp(log_I2 - log_I0);

    d_sq = 1.0 - I1_over_I0;
    numerator = 0.5 * (1.0 + I2_over_I0) - (I1_over_I0^2);
    log_I0_val = log_I0;
  } else {
    // Asymptotic expansion limits for large kappa values
    d_sq = 0.5 / kappa + 0.125 / (kappa^2);
    numerator = 0.5 / (kappa^2) + 0.25 / (kappa^3);
    log_I0_val = kappa - 0.5 * log(2 * pi() * kappa) + 1.0 / (8.0 * kappa);
  }

  // Prevent log(0) during warmup initialization
  if (d_sq < 1e-12) d_sq = 1e-12;
  if (numerator < 1e-12) numerator = 1e-12;

  // PC prior with point mass base model
  target += log(lambda) + log(numerator) - log(2.0) - 0.5 * log(d_sq) - lambda * sqrt(d_sq);

  // Manual safe von Mises likelihood
  for (n in 1:N) {
    target += kappa * cos(x[n] - mu) - log_I0_val;
  }
}

generated quantities {
  vector[N] log_lik;
  array[N] real x_rep;

  real log_I0_gen = kappa < 10000 ? log_modified_bessel_first_kind(0, kappa)
                                  : (kappa - 0.5 * log(2 * pi() * kappa)
                                     + 1.0 / (8.0 * kappa));

  for (n in 1:N) {
    // Pointwise log-likelihood for LOO (stable manual formulation)
    log_lik[n] = kappa * cos(x[n] - mu) - log(2 * pi()) - log_I0_gen;
    x_rep[n] = von_mises_rng(mu, kappa);
  }
}

# Compile the models using cmdstanr
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 = ys, lambda = lambdas_vm0[1]), refresh = 0)
fit_vm0_pi2 <- mod_vm0$sample(data = list(N = N, x = ys, lambda = lambdas_vm0[2]), refresh = 0)
fit_vm0_pi4 <- mod_vm0$sample(data = list(N = N, x = ys, lambda = lambdas_vm0[3]), refresh = 0)

fit_vminf_pi  <- mod_vminf$sample(data = list(N = N, x = ys, lambda = lambdas_vminf[1]), refresh = 0)
fit_vminf_pi2 <- mod_vminf$sample(data = list(N = N, x = ys, lambda = lambdas_vminf[2]), refresh = 0)
fit_vminf_pi4 <- mod_vminf$sample(data = list(N = N, x = ys, lambda = lambdas_vminf[3]), refresh = 0)

fit_vm0_pi$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.289  0.289 0.0516 0.0514 0.204 0.373  1.00    3580.    2474.
## 2 kappa    1.76   1.76  0.124  0.127  1.56  1.96   1.00    3556.    2968.

fit_vm0_pi2$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.293  0.293 0.0528 0.0540 0.208 0.380  1.00    3411.    2772.
## 2 kappa    1.76   1.76  0.126  0.124  1.55  1.96   1.00    3148.    2841.

fit_vm0_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.292  0.292 0.0541 0.0552 0.205 0.381  1.00    3316.    2715.
## 2 kappa    1.76   1.75  0.124  0.124  1.55  1.96   1.00    3051.    2635.

fit_vminf_pi$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.292  0.291 0.0529 0.0532 0.207 0.381  1.00    3537.    2771.
## 2 kappa    1.76   1.76  0.127  0.127  1.56  1.98   1.00    3472.    2350.

fit_vminf_pi2$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.292  0.292 0.0523 0.0521 0.208 0.378  1.00    3956.    2692.
## 2 kappa    1.76   1.76  0.127  0.127  1.56  1.98   1.00    3433.    2325.

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.291 0.0518 0.0514 0.205 0.375  1.00    2830.    2476.
## 2 kappa    1.77   1.77  0.128  0.125  1.56  1.98   1.00    3554.    2877.

## Plotting posteriors
# Extract posterior draws
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))

# Uniform Base Plots
p_kappa_vm0 <- ggplot(posteriors_vm0, aes(x = Kappa_Estimate, fill = Prior)) +
  geom_density(color = "white") +
  scale_fill_manual(values = my_colors_t, labels = sim_labels) +
  labs(title = expression("Uniform Base"), x ="", y = "Posterior Density") +
  theme(plot.title = element_text(face = "bold", hjust = 0.5))

p_mu_vm0 <- ggplot(posteriors_vm0, aes(x = Mu_Estimate, fill = Prior)) +
  geom_density(color = "white") +
  scale_fill_manual(values = my_colors_t, labels = sim_labels) +
  labs(title = expression("Uniform Base"), x = "", y = "") +
  theme(plot.title = element_text(face = "bold", hjust = 0.5))

# Point Mass Base Plots
p_kappa_vminf <- ggplot(posteriors_vminf, aes(x = Kappa_Estimate, fill = Prior)) +
  geom_density(color = "white") +
  scale_fill_manual(values = my_colors_t, labels = sim_labels) +
  labs(title = expression("Point Mass"), x = expression(kappa), y = "Posterior Density") +
  theme(plot.title = element_text(face = "bold", hjust = 0.5))

p_mu_vminf <- ggplot(posteriors_vminf, aes(x = Mu_Estimate, fill = Prior)) +
  geom_density(color = "white") +
  scale_fill_manual(values = my_colors_t, labels = sim_labels) +
  labs(title = expression("Point Mass"), x = expression(mu), y = "") +
  theme(plot.title = element_text(face = "bold", hjust = 0.5))

(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. This is useful sensitivity evidence for this example, not a general claim that the base-model choice is unimportant.

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

library(loo)

# Extract Log-Likelihood Matrices
# With circular uniform base model PC prior
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")
# With point mass base model PC prior
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")

# Calculate LOO
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
)
# Calculate WAIC
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
df$SE_dWAIC <- df$SE_dELPD * 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)",
             "k > 0.7", "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.08    -419.17    -419.15 -419.17   -419.17   -419.15
## 2    SE(ELPD)    16.63      16.61      16.59   16.68     16.69     16.72
## 3       dELPD     0.00      -0.09      -0.07   -0.09     -0.10     -0.08
## 4   SE(dELPD)     0.00       0.06       0.06    0.07      0.07      0.09
## 5       LOOIC   838.15     838.34     838.30  838.34    838.34    838.30
## 6   SE(LOOIC)    33.26      33.23      33.18   33.36     33.37     33.44
## 7      dLOOIC     0.00       0.19       0.15    0.19      0.19      0.15
## 8  SE(dLOOIC)     0.00       0.11       0.13    0.15      0.14      0.18
## 9        WAIC   838.15     838.34     838.30  838.34    838.34    838.30
## 10   SE(WAIC)    33.26      33.23      33.18   33.36     33.37     33.44
## 11      dWAIC     0.00       0.19       0.15    0.19      0.19      0.15
## 12  SE(dWAIC)     0.00       0.11       0.13    0.15      0.14      0.18
## 13    k > 0.7        0          0          0       0         0         0
## 14    k > 1.0        0          0          0       0         0         0

Differences should be interpreted relative to their standard errors rather than by rank alone. The predictive scores are essentially the same across the six fits, so none of these prior specifications has a meaningful predictive advantage for the wind data.

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.

# Calculate the Mean Resultant Length function
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, indicating that the fitted von Mises model reproduces the overall concentration reasonably well.

The main workflow is therefore:

1. choose a scientifically meaningful base model;
2. express prior knowledge through \((U,\alpha)\) on an interpretable circular scale;
3. inspect the induced prior on \(\kappa\) before fitting;
4. verify computation and assess sensitivity to reasonable calibrations; and
5. evaluate the likelihood with circular posterior predictive checks.

The base model determines the direction of regularization, while \((U,\alpha)\) determines its strength on an interpretable circular scale. PC priors do not remove prior judgment; they express it through these two explicit choices.