Échantillonnage à partir de la distribution de von Mises-Fisher en Python?

Je recherche un moyen simple d'échantillonner à partir d'une distribution multivariée de von Mises-Fisher en Python. J'ai regardé dans le module stats dans scipy et le module numpy mais j'ai seulement trouvé la distribution univariée de von Mises. Y a-t-il un code disponible? Je n'ai pas encore trouvé.

Apparemment, Wood (1994) a conçu un algorithme d'échantillonnage à partir de la distribution vMF selon ce lien , mais je ne trouve pas l'article.

- edit Pour plus de précision, je suis intéressé par l'algorithme qui est difficile à trouver dans la littérature (la plupart des articles se concentrent sur ). L'article fondateur (Wood, 1994) est introuvable à ma connaissance.$S^2$

Je l'ai finalement eu. Voici ma réponse.

J'ai enfin mis la main sur les statistiques directionnelles (Mardia et Jupp, 1999) et sur l'algorithme d'échantillonnage d'Ulrich-Wood. Je poste ici ce que j'en ai compris, c'est-à-dire mon code (en Python).

Le schéma d'échantillonnage de rejet:

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

$v \sqrt{1-w^2} + w \mu$$w$$v$

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

    result = []
    for sample in range(0,n):
        w = rW(n, 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

Et, pour un échantillonnage efficace avec ce code, voici un exemple:

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)

(Je m'excuse pour la mise en forme ici, j'ai créé un compte juste pour répondre à cette question, car j'essayais également de comprendre cela récemment).

La réponse du micro n'est pas tout à fait juste, le vecteur $v$ doit venir de $S^{p-2}$ dans l'espace tangent à $\mu$, C'est, $v$ devrait être un vecteur unitaire orthogonal à $\mu$. Sinon, le vecteur$v\sqrt{1-w^2}+w\mu$n'aura pas de norme un. Vous pouvez le voir dans l'exemple fourni par mic. Pour résoudre ce problème, utilisez quelque chose comme:

import scipy.linalg as la
def sample_tangent_unit(mu):
    mat = np.matrix(mu)

    if mat.shape[1]>mat.shape[0]:
        mat = mat.T

    U,_,_ = la.svd(mat)
    nu = np.matrix(np.random.randn(mat.shape[0])).T
    x = np.dot(U[:,1:],nu[1:,:])
    return x/la.norm(x)

et remplacer

v = np.random.randn(dim)
v = v / np.linalg.norm(v)

dans l'exemple du micro avec un appel à

v = sample_tangent_unit(mu)