¿Muestreo de la distribución de von Mises-Fisher en Python?

Question

Estoy buscando una manera simple muestra de un multivariante de von Mises-Fisher distribución en Python. Me he mirado en el módulo de estadísticas en scipy y la numpy módulo , pero sólo encontró la univariante de von Mises de distribución. ¿Hay algún código? No he encontrado todavía.

Al parecer, Wood (1994) ha diseñado un algoritmo para el muestreo de la vMF distribución de acuerdo a este enlace, pero no puedo encontrar el papel.

-- edit De precisión, estoy interesado por el algoritmo de la que es difícil encontrar en la literatura (la mayoría de los documentos se centran en $S^2$). El artículo seminal (Wood, 1994) no se puede encontrar de forma gratuita, para mi conocimiento.

Answer 1

Finalmente, lo conseguí. Aquí está mi respuesta.

Finalmente me puse las manos en las Estadísticas de dirección (Mardia y Jupp, 1999) y en el Ulrich-Madera del algoritmo de muestreo. He puesto aquí lo que he entendido, es decir, mi código (en Python).

El rechazo de un esquema de muestreo:

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

A continuación, el muestreo deseado es $v \sqrt{1-w^2} + w \mu$ donde $w$ es el resultado del rechazo de un esquema de muestreo, y $v$ es con muestreo uniforme sobre la 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

Y, para la eficacia de muestreo con este código, he aquí un ejemplo:

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

(Pido disculpas por el formato aquí, he creado una cuenta solo para responder a esta pregunta, ya que yo también estaba tratando de averiguar esto recientemente).

La respuesta de mic no está bien, el vector $v$ debe provenir de $S^{p-2}$ en el espacio de la tangente a$\mu$, $v$ debe ser un vector unitario ortogonal a $\mu$. De lo contrario, el vector $v\sqrt{1-w^2}+w\mu$ no tiene norma. Usted puede ver esto en el ejemplo proporcionado por el micrófono. Para solucionar este problema, utilice algo como:

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)

y reemplazar

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

en mic ejemplo con una llamada a

v = sample_tangent_unit(mu)