SciPy：半圆上的 von Mises 分布？

Posted 2023-03-12

【中文标题】SciPy：半圆上的 von Mises 分布？

【英文标题】：SciPy: von Mises distribution on a half circle?

【发布时间】：2020-12-06 11:49:44

【问题描述】：

我试图找出定义包裹在半圆上的 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)

我将如何更改它以使用半圆环绕？

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

【参考方案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 和采样直方图的图形

【讨论】：

感谢您如此详细的回复。这对我的情况非常有用，我将努力优化它（尽管我只需要每 3 秒执行一次操作！）。我在论坛 (link) 上看到，并通过询问其他人，我可以将我的角度乘以 2 以使其包裹在一个半圆上。我不确定我是否理解这一点，我肯定会将角度除以二吗？我想知道您对此是否有任何意见，但我承认您已经给出了很好的答案！

@mscone Thanks for such a detailed response 好吧，您可以接受答案，也可以颠倒它或两者兼而有之。我阅读了这些类比的讨论，并且总是如此。你没有冯米塞斯分布，时期。任何你认为可能来自 von Mises 的东西，或者任何其他人声称它会因为它为 von Mises 工作而声称它会起作用，都需要进行测试。这实际上是我编写的代码的美妙之处。 F.e. JohanC 建议镜像角度 - 好的，镜像它并针对数字 PDF 绘图，做 QQ 绘图，KS 测试等。如果有人想划分角度（是的，将角度减半） - 尝试它并针对数字 PDF 绘图。​​

关于所有提案：我认为，将角度减半是行不通的。我认为 JohanC 的镜像/加/减是行不通的。我认为接受/拒绝会起作用。但是你不能安全地测试各种参数（mu、kappa）的所有可能性，看看它是如何进行的。即使有人会给你带来魔法盒，声称它会根据你想要的分布生成 RV，你也可以检查一下

您提出了一个很好的观点，我可以测试所有这些，而且我一直都是。我只是对这个地区有点陌生，并要求一些指导。感谢您的代码、帮助和时间。

【参考方案2】：

您可以通过 numpy 的过滤（theta=theta[(theta>=0)&(theta<=np.pi)]，缩短样本数组）丢弃所需范围之外的值。因此，您可以先增加生成样本的数量，然后过滤，然后获取所需大小的子数组。

或者您可以添加/减去 pi 以将它们全部放入该范围（通过theta = np.where(theta < 0, theta + np.pi, np.where(theta > np.pi, theta - np.pi, theta))）。正如@SeverinPappadeux 所指出的那样，这种改变分布并且可能是不希望的。

import matplotlib.pyplot as plt
from matplotlib.collections import LineCollection
import numpy as np
from scipy.stats import vonmises

mu = np.pi / 2
kappa = 4

orig_theta = vonmises.rvs(kappa, loc=mu, size=(10000))
fig, axes = plt.subplots(ncols=2, sharex=True, sharey=True, figsize=(12, 4))
for ax in axes:
    theta = orig_theta.copy()
    if ax == axes[0]:
        ax.set_title(f"$Von Mises, \\mu=mu:.2f, \\kappa=kappa$")
    else:
        theta = theta[(theta >= 0) & (theta <= np.pi)]
        print(len(theta))
        ax.set_title(f"$Von Mises, angles\\ filtered\\ (100 * len(theta) / (len(orig_theta)):.2f\\ \\%)$")
    segs = np.zeros((len(theta), 2, 2))
    segs[:, 1, 0] = np.cos(theta)
    segs[:, 1, 1] = np.sin(theta)
    line_segments = LineCollection(segs, linewidths=.1, colors='blue', alpha=0.5)
    ax.add_collection(line_segments)
    ax.autoscale()
    ax.set_aspect('equal')
plt.show()

【讨论】：

加/减是错误的，我相信。它不再遵循冯米塞斯的形状

是的，你是对的，尾巴会变得太肥。我错误地认为它们会以类似的方式分布。我会更新答案。

感谢您如此详细的回复！我希望有一种简单的方法可以改变方程式以绕半圈。我忘了提，输出必须始终是相同数量的元素。使用拒绝抽样是否合适？从 vonmises.rvs 请求一个值，然后根据它是否在距离所需平均值的 pi/2 范围内接受或拒绝它？

出于速度考虑，我真的会生成太多样本，然后再取一个子样本。 theta = vonmises.rvs(kappa, loc=mu, size=(15000)) 然后theta = theta[(theta >= 0) & (theta <= np.pi)][:10000]。