Sampling aus der von Mises-Fisher-Distribution in Python?

Ich suche nach einer einfachen Möglichkeit, eine multivariate von Mises-Fisher-Distribution in Python abzutasten. Ich habe im Statistikmodul in Scipy und im Numpy-Modul nachgesehen, aber nur die univariate von Mises-Verteilung gefunden. Gibt es einen Code zur Verfügung? Habe ich noch nicht gefunden.

Anscheinend hat Wood (1994) einen Algorithmus zum Abtasten von der vMF-Verteilung gemäß diesem Link entworfen , aber ich kann das Papier nicht finden.

- edit Aus Gründen der Genauigkeit interessiert mich der Algorithmus, der in der Literatur schwer zu finden ist (die meisten Arbeiten konzentrieren sich auf ). Meines Wissens kann der wegweisende Artikel (Wood, 1994) nicht kostenlos gefunden werden.$S^2$

Endlich hab ich es. Hier ist meine Antwort.

Ich habe schließlich die Richtungsstatistik (Mardia und Jupp, 1999) und den Ulrich-Wood-Algorithmus für die Stichprobenentnahme in die Hand genommen . Ich poste hier, was ich davon verstanden habe, dh meinen Code (in Python).

Stichprobenplan für die Ablehnung:

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

Und um effektiv mit diesem Code zu arbeiten, hier ein Beispiel:

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)

(Ich entschuldige mich für die Formatierung hier. Ich habe ein Konto erstellt, um auf diese Frage zu antworten, da ich kürzlich auch versucht habe, dies herauszufinden.)

Die Antwort von mic ist nicht ganz richtig, der Vektor $v$ muss kommen aus $S^{p-2}$ im tangentialen Raum zu $\mu$, das ist, $v$ sollte ein Einheitsvektor orthogonal zu sein $\mu$. Ansonsten der Vektor$v\sqrt{1-w^2}+w\mu$wird nicht norm eins haben. Dies können Sie dem Beispiel von mic entnehmen. Verwenden Sie zum Beheben dieses Problems Folgendes:

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)

und ersetzen

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

in mic's beispiel mit einem anruf an

v = sample_tangent_unit(mu)