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

Estoy buscando una forma simple de muestrear de una distribución multivariada de von Mises-Fisher en Python. He buscado en el módulo de estadísticas en scipy y en el módulo numpy, pero solo encontré la distribución univariante de von Mises. ¿Hay algún código disponible? No he encontrado aún.

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

- editar Para mayor precisión, me interesa el algoritmo que es difícil de encontrar en la literatura (la mayoría de los artículos se centran en ). El artículo seminal (Wood, 1994) no se puede encontrar de forma gratuita, que yo sepa. $S^2$

distributions sampling python mic
fuente

La entrada a scipy.stats.vonmisespuede ser de tipo matriz, por lo que puede especificar la distribución como un array. Vea este ejemplo

rightskewed

Gracias por tu respuesta. Sin embargo, parece que es más un producto de 1-D von Mises que una verdadera nD von Mises-Fisher: K = vonmises.pdf([x,x], kappa=[[1],[10]]). Un vMF 2-D debe tener solo un

real como parámetro. ¿Estás de acuerdo?

κ

$\kappa$

micrófono

Estoy buscando el algoritmo VM * originalmente en Simulación de la distribución de von Mises Fisher (Wood, 1994). ¿Nadie?

mic

Encontré las respuestas en este hilo aquí realmente útiles. He proporcionado una función de utilidad ligeramente limpia para hacer esto como parte de este paquete: https://github.com/clara-labs/spherecluster/blob/develop/spherecluster/util.py , para aquellos que todavía buscan generar esto datos.

Jaska

Respuestas:

Finalmente lo tengo. Aquí está mi respuesta.

Finalmente puse mis manos en Estadísticas Direccionales (Mardia y Jupp, 1999) y en el algoritmo de Ulrich-Wood para el muestreo. Publico aquí lo que entendí de él, es decir, mi código (en Python).

El esquema de muestreo de rechazo:

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

Entonces, el muestreo deseado es , dondees el resultado del esquema de muestreo de rechazo, yse muestrea de manera uniforme sobre la hiperesfera. $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

Y, para un muestreo efectivo con este código, aquí hay 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)

mic
fuente

(+1) ¡Gracias por compartir su respuesta (especialmente a pesar del desaliento potencial de tener su pregunta inicialmente cerrada)!

whuber

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

$v$ $S^{p-2}$ $\mu$ $v$ $\mu$ $v\sqrt{1-w^2}+w\mu$ No tendrá la norma uno. Puede ver esto en el ejemplo proporcionado por mic. Para solucionar esto, use 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 el ejemplo del micrófono con una llamada a

v = sample_tangent_unit(mu)

Kevin
fuente