Выборка из многомерного распределения фон Мизес-Фишера в Python
-
16-10-2019 - |
Вопрос
Я ищу простой способ попробовать из Многомерный фон Мизес-Фишер Распределение в Python. Я посмотрел модуль статистики в Scipy и Numpy Module но нашел только одномерное распределение фон Мизеса. Есть ли код доступен? Я еще не нашел.
-- редактировать. По -видимому, Вуд (1994) разработал алгоритм для отбора проб из распределения VMF в соответствии с эта ссылка, но я не могу найти бумагу.
Решение 2
Благодаря вашей помощи. Я наконец -то получил свой код, плюс библиография.
Я положил руки на Направленная статистика (Mardia and Jupp, 1999) и об алгоритме Ульрих-Вуда для отбора проб. Я публикую здесь то, что я понял из этого, то есть мой код (в Python), со вкусом «Movmf».
Схема отбора проб отказа:
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(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
И для эффективной выборки с этим кодом, вот пример:
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)
Другие советы
Похоже, вы можете попробовать распределение фон Мизес-Фишер с этим пакетом R. Вы думали о том, чтобы позвонить R из Python, используя RPY2 упаковка? Я не пробовал это для себя, но не могли бы вы попробовать следующее?
from numpy import *
import scipy as sp
from pandas import *
from rpy2.robjects.packages import importr
import rpy2.robjects as ro
import pandas.rpy.common as com
from rpy2.robjects.packages import importr
# import the movMF R package
movMF = importr('movMF')
# call the rmovMF sampling function from the R package
print(movMF.rmovMF(10, 3 * c(1, -1) / sqrt(2)))