random - SciPy:半圆上的 von Mises 分布？

Question

我正在尝试找出定义包裹在半圆上的 von-Mises 分布的最佳方法(我用它来绘制不同浓度的无方向线)。我目前正在使用 SciPy 的 vonmises.rvs()。本质上，我希望能够输入 pi/2 的平均方向，并将分布截断到任一侧不超过 pi/2。

我可以使用截断的正态分布，但我会失去 von-mises 的环绕(如果我想要平均方向为 0)

我在研究映射纤维方向的研究论文中看到过这样做，但我不知道如何实现它(在 python 中)。我有点不知道从哪里开始。

如果我的 von Mesis 定义为(来自 numpy.vonmises):

np.exp(kappa*np.cos(x-mu))/(2*np.pi*i0(kappa))

与:

mu, kappa = 0, 4.0

x = np.linspace(-np.pi, np.pi, num=51)

我如何更改它以使用环绕半圆来代替？

有这方面经验的人可以提供一些指导吗？

Answer 1

直接数值逆 CDF 采样很有用，它应该非常适合有界域的分布。这是代码示例，构建 PDF 和 CDF 表并使用逆 CDF 方法进行采样。当然可以优化和矢量化

代码，Python 3.8，x64 Windows 10

import numpy as np
import matplotlib.pyplot as plt
import scipy.integrate as integrate

def PDF(x, μ, κ):
    return np.exp(κ*np.cos(x - μ))

N = 201

μ = np.pi/2.0
κ = 4.0

xlo = μ - np.pi/2.0
xhi = μ + np.pi/2.0

# PDF normaliztion

I = integrate.quad(lambda x: PDF(x, μ, κ), xlo, xhi)
print(I)
I = I[0]

x = np.linspace(xlo, xhi, N, dtype=np.float64)
step = (xhi-xlo)/(N-1)

p = PDF(x, μ, κ)/I # PDF table

# making CDF table
c = np.zeros(N, dtype=np.float64)

for k in range(1, N):
    c[k] = integrate.quad(lambda x: PDF(x, μ, κ), xlo, x[k])[0] / I

c[N-1] = 1.0 # so random() in [0...1) range would work right

#%%
# sampling from tabular CDF via insverse CDF method

def InvCDFsample(c, x, gen):
    r = gen.random()
    i = np.searchsorted(c, r, side='right')
    q = (r - c[i-1]) / (c[i] - c[i-1])
    return (1.0 - q) * x[i-1] + q * x[i]

# sampling test
RNG = np.random.default_rng()

s = np.empty(20000)

for k in range(0, len(s)):
    s[k] = InvCDFsample(c, x, RNG)

# plotting PDF, CDF and sampling density
plt.plot(x, p, 'b^') # PDF
plt.plot(x, c, 'r.') # CDF
n, bins, patches = plt.hist(s, x, density = True, color ='green', alpha = 0.7)
plt.show()

以及带有 PDF、CDF 和采样直方图的图形