SciPy：冯·米塞斯分布在一个半圆上吗？

问题描述

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

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

任何有经验的人都可以提供一些指导吗？

解决方法

您可以通过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] = 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()

对直接数字逆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,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,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,density = True,color ='green',alpha = 0.7)
plt.show()

并带有PDF，CDF和采样直方图的图形

distribution distribution random random random scipy scipy scipy.stats