Próbkowanie z dystrybucji von Misesa-Fishera w Pythonie?

Szukam prostego sposobu na pobranie próbki z wielowymiarowej dystrybucji von Misesa-Fishera w Pythonie. Mam spojrzał w statystyki modułu w scipy i modułem numpy ale tylko znaleźć jednoczynnikowej dystrybucję von Misesa. Czy jest dostępny kod? Jeszcze nie znalazłem

Najwyraźniej Wood (1994) opracował algorytm próbkowania z rozkładu vMF zgodnie z tym linkiem , ale nie mogę znaleźć papieru.

- edycja Dla precyzji interesuje mnie algorytm, który jest trudny do znalezienia w literaturze (większość artykułów dotyczy ). Według mojej wiedzy nie można znaleźć tego przełomowego artykułu (Wood, 1994). $S^2$

distributions sampling python

— mikrofon
źródło

Dane wejściowe scipy.stats.vonmisesmogą być podobne do tablic, dzięki czemu można określić rozkład jako array. Zobacz ten przykład

— przesunięty w prawo

Dziękuję za odpowiedź. Ale wydaje się, że jest to bardziej produktem 1-D von Mises niż prawdziwy nD von Misesa-Fishera: K = vonmises.pdf([x,x], kappa=[[1],[10]]). 2-D vMF powinien mieć tylko jeden rzeczywisty

jako parametr. Czy sie zgadzasz?

κ

$\kappa$

— mic

Szukam algorytmu VM * pierwotnie w Symulacji rozkładu von Misesa Fishera (Wood, 1994). Ktoś?

— mic

Znalazłem odpowiedzi w tym wątku tutaj bardzo przydatne. Jako część tego pakietu podałem nieco oczyszczoną funkcję narzędzia: https://github.com/clara-labs/spherecluster/blob/develop/spherecluster/util.py , dla tych, którzy nadal chcą to wygenerować dane.

— Jaska

Odpowiedzi:

W końcu to dostałam. Oto moja odpowiedź.

W końcu zdecydowałem się na statystyki kierunkowe (Mardia i Jupp, 1999) oraz algorytm Ulricha-Wooda do próbkowania. Publikuję tutaj to, co z niego zrozumiałem, tj. Mój kod (w języku Python).

Odrzucony schemat próbkowania:

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

Następnie pożądanym próbkowaniem jest , gdziejest wynikiem ze schematu prób odrzucenia, ajest równomiernie próbkowany na hipersferze. $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

I, dla skutecznego próbkowania z tym kodem, oto przykład:

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)

— mikrofon
źródło

(+1) Dziękujemy za podzielenie się odpowiedzią (szczególnie pomimo potencjalnego zniechęcenia do początkowego zamknięcia pytania)!

— whuber

(Przepraszam za formatowanie tutaj, utworzyłem konto tylko po to, aby odpowiedzieć na to pytanie, ponieważ ostatnio próbowałem to rozgryźć).

$v$ $S^{p-2}$ w przestrzeni stycznej do $\mu$ , to jest, $v$ powinien być wektorem jednostkowym prostopadłym do $\mu$ . W przeciwnym razie wektor $v\sqrt{1-w^2}+w\mu$ nie będzie miał normy jeden. Możesz to zobaczyć w przykładzie dostarczonym przez mikrofon. Aby to naprawić, użyj czegoś takiego:

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)

i zamień

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

w przykładzie mikrofonu z wezwaniem do

v = sample_tangent_unit(mu)

— Kevin
źródło