# Sampling from the von-Mises Fisher distribution?

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!

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

• Thanks. I was aware of this paper. But, to my understanding, it only concerns sampling over the sphere $S^2$ while I am interested in sampling over $S^d$ for arbitrary d > 2. – mic Jun 16 '15 at 7:44
• @mic Oh, I see. Nice that you put your implementation as a reference! – MBW Jun 16 '15 at 20:49

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)