Circular data and the von Mises distribution

R
circular statistics
directional data
ecology tutorial
ggplot2
Why arithmetic means fail on the circle: the mean direction and resultant length for animal orientations, and fitting a von Mises distribution in base R.
Author

Tidy Ecology

Published

2026-08-18

Ecologists collect angles more often than they realise. Compass bearings of departing birds, the aspect of a nest slope, the time of day an animal is active, the flowering date on the annual cycle: all of these live on a circle, not on a line. The awkward part is that the number 359 and the number 1 are two degrees apart, not 358, and every ordinary summary you reach for assumes the opposite. This post sets out the base-R toolkit for angular data: the mean direction, the mean resultant length as a measure of concentration, and the von Mises distribution, the circular counterpart of the normal.

Everything here is hand-coded with base stats plus ggplot2. There is no circular-statistics package to install, and building the pieces yourself makes the assumptions visible.

The wrap-around trap

Take five bearings clustered around north: 350, 10, 20, 340 and 0 degrees. By eye they point north. The arithmetic mean does not.

# work in radians throughout; degrees only for display
deg2rad <- function(d) d * pi / 180
rad2deg <- function(r) r * 180 / pi

circ_mean <- function(theta) atan2(mean(sin(theta)), mean(cos(theta)))   # in (-pi, pi]

bearings  <- c(350, 10, 20, 340, 0)
arith     <- mean(bearings)                                              # naive linear mean
circular  <- rad2deg(circ_mean(deg2rad(bearings)))                       # correct mean direction, (-180, 180]
circular  <- ifelse(abs(circular) < 0.5, 0, circular %% 360)            # tidy the exact-north case
c(arithmetic = arith, circular = circular)
arithmetic   circular 
       144          0 

The arithmetic mean returns 144 degrees, pointing south-east, because it cannot see that 350 and 10 are neighbours. The circular mean returns 0 degrees, which is north. The fix is to average the unit vectors, not the raw numbers: replace each angle by its point on the unit circle, sum the horizontal and vertical components, and take the direction of the total. That single idea drives everything below.

Mean direction and resultant length

For a sample of angles the two building blocks are the mean cosine and mean sine. Their vector length is the mean resultant length, written R, and its direction is the mean direction.

circ_R <- function(theta) sqrt(mean(cos(theta))^2 + mean(sin(theta))^2)  # resultant length in [0, 1]

# synthetic orientations: 60 departure bearings, illustrative only
set.seed(4250)
theta <- rvm(60, mu = deg2rad(40), kappa = 3.0)

Rbar     <- circ_R(theta)
mean_dir <- rad2deg(circ_mean(theta)) %% 360
circ_var <- 1 - Rbar                       # circular variance in [0, 1]
circ_sd  <- rad2deg(sqrt(-2 * log(Rbar)))  # circular standard deviation, degrees
round(c(Rbar = Rbar, mean_dir = mean_dir, circ_var = circ_var, circ_sd = circ_sd), 3)
    Rbar mean_dir circ_var  circ_sd 
   0.830   45.127    0.170   34.933 

The resultant length is 0.830, a number between zero and one. At one, every angle is identical; at zero, the angles cancel out and there is no preferred direction. Here it is 0.83, a moderately tight cluster, with a mean direction of 45.1 degrees and a circular standard deviation of 34.9 degrees.

The trap from the first section is not a toy. The arithmetic mean of these same 60 bearings is 81.1 degrees, off by more than thirty degrees, because a handful of angles fall just below 360 and drag the linear average upward. The circular mean is immune to where you place the zero.

Two compass plots. The left shows five points near the top of the circle with a green arrow pointing up and a rust arrow pointing down-left. The right shows sixty points clustered top-right with a single arrow of length 0.83 pointing towards forty-five degrees.
Figure 1: Left: five bearings near north, with the arithmetic mean (rust) pointing the wrong way and the circular mean (green) pointing north. Right: 60 synthetic orientations with the mean resultant vector; its length R measures concentration.

The von Mises distribution

The normal distribution does not wrap, so it cannot describe angles. Its circular replacement is the von Mises distribution, with a mean direction mu and a concentration kappa. The density is

f(theta) = exp(kappa * cos(theta - mu)) / (2 * pi * I0(kappa))

where I0 is the modified Bessel function of the first kind, order zero, available in base R as besselI(kappa, 0). When kappa is zero the distribution is uniform on the circle; as kappa grows it concentrates around mu, and for large kappa it looks almost normal.

Maximum likelihood is unusually clean. The estimate of mu is simply the mean direction, and the estimate of kappa solves A(kappa) = R, where A is the ratio of two Bessel functions. There is no closed form, so solve it numerically.

dvm     <- function(x, mu, kappa) exp(kappa * cos(x - mu)) / (2 * pi * besselI(kappa, 0))
A_kappa <- function(kappa) besselI(kappa, 1) / besselI(kappa, 0)

kappa_mle <- function(Rbar) {
  if (Rbar < 1e-8)      return(0)
  if (Rbar > 0.999999)  return(Inf)
  uniroot(function(k) A_kappa(k) - Rbar, interval = c(1e-6, 700))$root
}

kappa_hat <- kappa_mle(Rbar)
kappa_hat
[1] 3.307571

The fitted concentration is 3.31. That is the value that makes the theoretical resultant length A(kappa) match the observed 0.83.

One caution worth carrying forward: the maximum-likelihood kappa is biased upward in small samples, sometimes badly. The standard correction (Best and Fisher 1981) nudges it back:

kappa_bc <- function(k, n) {
  if (!is.finite(k)) return(k)
  if (k < 2) return(max(k - 2 / (n * k), 0))
  (n - 1)^3 * k / (n^3 + n)
}
kappa_star <- kappa_bc(kappa_hat, length(theta))
kappa_star
[1] 3.14406

The corrected value is 3.14, a little lower. With sixty angles the difference is small; with fifteen it would not be, and the checking post in this series returns to how far the reported concentration can drift from the truth.

Warning: Removed 1 row containing missing values or values outside the scale range
(`geom_col()`).
Left panel is a rose diagram with the longest wedges near north-east. Right panel is a histogram with three density curves: the fitted kappa of 3.3 in green tracks the histogram, a diffuse kappa of 1 in amber is nearly flat, and a tight kappa of 8 in rust is a narrow spike.
Figure 2: Left: a rose diagram of the 60 orientations in 30-degree sectors. Right: the same data unrolled and centred on the mean direction, with the fitted von Mises density and two reference concentrations. Larger kappa means a tighter peak.

Grouped data

Field bearings often arrive already binned into sectors, and the same summaries work from counts. Weight each sector midpoint by its count, average the components, and read off the resultant.

Cg <- sum(cnt * cos(mids)) / sum(cnt)
Sg <- sum(cnt * sin(mids)) / sum(cnt)
Rbar_g <- sqrt(Cg^2 + Sg^2)
mean_dir_g <- rad2deg(atan2(Sg, Cg)) %% 360
round(c(Rbar_grouped = Rbar_g, mean_dir_grouped = mean_dir_g), 3)
    Rbar_grouped mean_dir_grouped 
           0.822           46.007 

The grouped resultant length is 0.82 against 0.83 for the raw angles, and the grouped mean direction is 46.0 degrees against 45.1. Binning into twelve sectors costs almost nothing here, though very coarse sectors would smear a tight cluster.

What R does not tell you

The mean direction and resultant length describe the data as though a single preferred direction exists. That assumption is doing quiet work. If animals move along a two-way axis, say up and down a valley, the angles split into two opposite clusters, the components cancel, and R collapses towards zero even though the pattern is strongly structured. A near-zero R can mean “no preferred direction” or “two opposite directions”, and the mean direction becomes meaningless. Reporting R alone would hide the difference.

That is the reason the next post in this series asks a sharper question before any of these summaries: is there a preferred direction at all, and which test can actually see the structure that is present. The checking post at the end returns to how modality, sample size and the choice of test decide what a circular summary is allowed to claim.

References

Batschelet E 1981 Circular Statistics in Biology. Academic Press. ISBN 978-0-12-081050-9.

Fisher NI 1993 Statistical Analysis of Circular Data. Cambridge University Press. ISBN 978-0-521-56890-6.

Mardia KV, Jupp PE 2000 Directional Statistics. Wiley. ISBN 978-0-471-95333-3.

Jammalamadaka SR, SenGupta A 2001 Topics in Circular Statistics. World Scientific. ISBN 978-981-02-3778-3.

Pewsey A, Neuhauser M, Ruxton GD 2013 Circular Statistics in R. Oxford University Press. ISBN 978-0-19-967113-7.

Best DJ, Fisher NI 1979 Applied Statistics 28(2):152-157 (10.2307/2346732).

Best DJ, Fisher NI 1981 Communications in Statistics - Simulation and Computation 10(5):493-502 (10.1080/03610918108812225).

Newsletter

Get new tutorials by email

New R and QGIS tutorials for ecologists, straight to your inbox. No spam; unsubscribe anytime.

By subscribing you agree to receive these emails and confirm your address once. See the privacy policy.