Pythonの多変量Von Mises-Fisher分布からのサンプリング

Question

私はからサンプリングする簡単な方法を探しています 多変量フォンミーゼスフィッシャー Pythonの分布。私は調べました Scipyの統計モジュール そしてその Numpyモジュール しかし、単変量のフォンミーゼス分布のみが見つかりました。利用可能なコードはありますか？まだ見つかりません。

- 編集。どうやら、Wood（1994）は、VMF分布からサンプリングするためのアルゴリズムを設計しています。 このリンク, 、しかし、私は紙が見つかりません。

Answer 1

あなたの助けに感謝します。ついにコードが機能し、さらに書誌が機能しました。

手に入れました 方向統計 （Mardia and Jupp、1999）およびサンプリングのためのUlrich-Woodのアルゴリズムについて。ここに私が理解したこと、つまり私のコード（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)

Answer 2

そのRパッケージでVon Mises-Fisher Distributionをサンプリングできるようです。 Python内からRを呼び出すことを考えましたか 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)))