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

  • $\begingroup$ 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. $\endgroup$ – mic Jun 16 '15 at 7:44
  • $\begingroup$ @mic Oh, I see. Nice that you put your implementation as a reference! $\endgroup$ – 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
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)

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.