Numerically stable sampling of the von Mises
Fisher distribution on S 2 (and other tricks)
Wenzel Jakob
July 21, 2012 (revised Jun 28, 2015)
1 Introduction
The von Mises-Fisher (vMF) distribution is a popular all-purpose distribution [1, 4] for statisticalinference involving directional data. On the 2-sphere, it is defined asκ
fvMF (ω) =exp(κµT ω)
4π sinh κ
where µ ∈ S 2 is the mean direction and κ denotes the concentration parameter (κ → 0 approaching
the uniform distribution). A recent application [2] of this distribution in computer graphics entailed
fitting mixture models composed of vMF functions to arbitrary spherical data using the expectation
maximization procedure.
Unfortunately, many basic operations involving this distribution are prone to severe numericalissues when implemented in finite precision computer arithmetic. There is a surprising lack on
information on how these can be circumvented, and hence the purpose of this document is to serve
as a collection of numerically-well behaved recipes for common operations.

2 Evaluation
Evaluation of the vMF distribution easily overflows single precision arithmetic even for moderateconcentration values (for instance, sinh 100 = 1.34406 · 1043 ), and double precision fails shortly
thereafter. The following equivalent expression follows from exponential function identities andworks reliably over a much larger range of concentrations: 1 ,κ=0fvMF (ω) = 4πκTeκ(µ ω−1) , κ > 0
2π(1 − exp(−2κ))

3 Sample generation
Several prior works have investigated how independent samples can be drawn so that they aredistributed according to the vMF distribution [5, 6, 3]. The following is a brief summary of [3],which leads to a simple but numerically ill-behaved method:1

Observe that the following random vector with mean direction µ = (0, 0, 1) is distributed according to fvMF [5]:pωκ = ( 1 − W 2 V, W )T
where V and W are independent random variables, V ∈ R2 is a uniformly distributed vector on
the unit circle, and W ∈ [−1, 1] follows the densityκfW (w) =exp(κw).
2 sinh κ
All that is needed for a computer implementation is a way to generate realizations of W . Applying
the inversion method results in−1FW(ξ) = κ−1 log exp−κ +2 ξ sinh κTo handle other values of µ, one can simply apply a rotation to directions obtained in this manner.

3.1 Numerically stable variant
Again, we can apply exponential function identities to arrive at an equivalent and numericallywell-behaved expression, which avoids overflow for large values of κ:−1(ξ) = 1 + κ−1 log ξ + (1 − ξ)e−2κ
FW4 Finding κ such that fvMF (µ) = c
One very useful tool is the ability to create distributions that have a specified solid angle density
into a certain direction. In the case of the van Mises-Fisher distribution, we can see that fvMF
takes on its maximum into direction µ, whereκ(1 + coth κ).
g(κ) :=4π
gives the maximum as a function of the concentration. Unfortunately, it is not convenient toanalytically invert this expression for κ. However, note thatcoth κ =e2κ + 1e2κ − 1

rapidly approaches 1. For instance, coth 5 is already approximately equal to 1.0009. Assuming that
there are no particularly stringent accuracy requirements on the inversion, we can use the following
approximate scheme:(2πx,x > g(5) ≈ 0.795−1g (x) ≈−1
grat (x), otherwise
where we have approximated coth κ ≈ 1 for κ > 5 and make use the following rational interpolant
elsewhere:
2−1−5 168.479x + 16.4585x − 2.39942
grat (x) := max 10 ,.
−1.12718x2 + 29.1433x + 1−1
On the interval [1/4π, g(5)], the function grat
has an absolute error of < 0.007 The relative error is−1
1infinite, as g (x) → 0 (x → /4π). Figure 1 shows an illustration of the fit.

2

540.2

0.3

0.4

0.5

0.6

0.7

0.80.1320.001

110 - 50.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

(a) Linear plot(b) Logarithmic plot

−1
Figure 1: Fit of the rational function grat(blue) to g −1 (red).

5 Convolution
The convolution of two vMF distributions does not generally produce another vMF distribution.
However, the result of this operation can be well-approximated by a vMF distribution with a suitably chosen value of κ. Mardia and Jupp [4] describe one approach to obtain this parameter, which
entails approximating the distributions to be convolved by wrapped normal distributions, convolving them instead, and transforming the result back into a vMF distribution. A C implementationof this is given below:
float A3(float kappa) {return 1 / std::tanh(kappa) - 1 / kappa;}
float dA3(float kappa) {float csch = 2.0f / (std::exp(kappa) - std::exp(-kappa));
return 1 / (kappa*kappa) - csch*csch;}float A3inv(float y, float guess) {/* Initial guess */float x = guess, residual = 0;/* Invert using Newton’s method */do {residual = A3(x)-y, deriv = dA3(x);x -= residual/deriv;} while (std::abs(residual) > 1e-5f);
return x;}
float convolve(float kappa1, float kappa2) {return A3inv(A3(kappa1) * A3(kappa2), std::min(kappa1, kappa2));}

3

References[1] Fisher, N., Lewis, T., and Embleton, B. Statistical analysis of spherical data. Cambridge
University Press, 1993.
[2] Han, C., Sun, B., Ramamoorthi, R., and Grinspun, E. Frequency domain normal mapfiltering. ACM Transactions on Graphics (Proceedings of SIGGRAPH 2007) 26, 3 (2007), 28:1–
28:12.
[3] Jung, S. Generating von mises fisher distribution on the unit sphere (s2 ). Tech. rep., Universityof Pittsburgh, 10 2009.
[4] Mardia, K., and Jupp, P. Directional statistics. Wiley, 2000.
[5] Ulrich, G. Computer generation of distributions on the m-sphere. Applied Statistics (1984),158–163.
[6] Wood, A. Simulation of the von mises fisher distribution. Communications in statisticssimulation and computation 23, 1 (1994), 157–164.

4