Sampling from the von-Mises Fisher distribution?

Question

This topic has already been tackled on this website (here). But, unfortunately, no clear cut answers were given. In (Wood,1994), there is apparently a rejection algorithm for sampling from this distribution, but I can't find it on the web. I am thus looking for a pseudo-code (I have not much experience with sampling random variables). Many thanks!

Answer 1

I finally put my hands on Directional Statistics (Mardia and Jupp, 1999) and on the Ulrich-Wood's algorithm for sampling. I post here what I understood from it, i.e. my code (in Python).

The rejection sampling scheme:

def rW(n,kappa,m):
dim = m-1
b = dim / (np.sqrt(4*kappa*kappa + dim*dim) + 2*kappa)
x = (1-b) / (1+b)
c = kappa*x + dim*np.log(1-x*x)

y = []
for i in range(0,n):
    done = False
    while not done:
        z = sc.stats.beta.rvs(dim/2,dim/2)
        w = (1 - (1+b)*z) / (1 - (1-b)*z)
        u = sc.stats.uniform.rvs()
        if kappa*w + dim*np.log(1-x*w) - c >= np.log(u):
            done = True
    y.append(w)
return y

Then, the desired sampling is $v \sqrt{1-w^2} + w \mu$, where $w$ is the result from the rejection sampling scheme, and $v$ is uniformly sampled over the hypersphere.

def rvMF(n,theta):
dim = len(theta)
kappa = np.linalg.norm(theta)
mu = theta / kappa

result = []
for sample in range(0,n):
    w = rW(kappa,dim)
    v = np.random.randn(dim)
    v = v / np.linalg.norm(v)

    result.append(np.sqrt(1-w**2)*v + w*mu)

return result

And, for effectively sampling with this code, here is an example:

import numpy as np
import scipy as sc
import scipy.stats

n = 10
kappa = 100000
direction = np.array([1,-1,1])
direction = direction / np.linalg.norm(direction)

res_sampling = rvMF(n, kappa * direction)

Answer 2

Here you can find a relatively simple way to sample the von Mises-Fisher distribution. More specifically, on section 3:

Observe that the following random vector with mean direction $\mu = (0, 0, 1)$ is distributed according to $f_{\text{vMF}}$ [5]: $$\omega_{\kappa} = (\sqrt{1 − W^2} V, W)^T$$ where $V$ and $W$ are independent random variables, $V \in \mathbb{R}^2$ is a uniformly distributed vector on the unit circle, and $W \in [−1, 1]$ follows the density $$f_W (w) = \frac{\kappa}{ 2 \sinh \kappa} \exp(κw).$$ All that is needed for a computer implementation is a way to generate realizations of $W$. Applying the inversion method results in $$F^{−1}_W (\xi) = \kappa^{-1}\log\left(\exp(−\kappa) +2 \xi \sinh \kappa\right)$$ To handle other values of $\mu$, one can simply apply a rotation to directions obtained in this manner

Answer 3

The above answer by @mic does not work if $\mu$ and $v$ are not orthogonal, because the resulting $v\sqrt{1 - w^2} + w \mu$ does not lie on the sphere in that case. His procedure can be used with $\mu=(1,0,\ldots,0)$ and $v=(0, \tilde v)$ to enforce orthogonality, then the samples can be mapped around an arbitrary mean $\mu$ by applying a rotation matrix (this does not change the distribution).

For info, the sampling method is available in the open source Python package geomstats, and can be used with the following:

from geomstats.geometry.hypersphere import Hypersphere

sphere = Hypersphere(dim=4)

kappa = 10.
mean = sphere.random_uniform()

samples = sphere.random_von_mises_fisher(
    mu=mean, kappa=kappa, n_samples=1000)

The python code is available here.