SciPy：von Mises分布在半个圆上？ [英] SciPy: von Mises distribution on a half circle?

问题描述

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

I'm trying to figure out the best way to define a von-Mises distribution wrapped on a half-circle (I'm using it to draw directionless lines at different concentrations). I'm currently using SciPy's vonmises.rvs(). Essentially, I want to be able to put in, say, a mean orientation of pi/2 and have the distribution truncated to no more than pi/2 either side.

我可以使用a截断正态分布，但是我将失去von-mises的包裹（例如，如果我想要平均方向为0）

I could use a truncated normal distribution, but I will lose the wrapping of the von-mises (say if I want a mean orientation of 0)

我已经在研究映射光纤的研究论文中看到了这一点方向，但我不知道如何实现它（在python中）。我有点从哪里开始。

I've seen this done in research papers looking at mapping fibre orientations, but I can't figure out how to implement it (in python). I'm a bit stuck on where to start.

如果我的冯·梅西斯被定义为（摘自numpy.vonmises）：

If my von Mesis is defined as (from numpy.vonmises):

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

with：

mu, kappa = 0, 4.0

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

我如何更改它以使用环绕

How would I alter it to use a wrap around a half-circle instead?

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

Could anyone with some experience with this offer some guidance?

推荐答案

对直接数字逆CDF采样很有用，它对于带界域分布非常有用。这是代码示例，构建PDF和CDF表并使用反CDF方法进行采样。当然可以进行优化和向量化

Is is useful to have direct numerical inverse CDF sampling, it should work great for distribution with bounded domain. Here is code sample, building PDF and CDF tables and sampling using inverse CDF method. Could be optimized and vectorized, of course

代码，Python 3.8，x64 Windows 10

Code, 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和采样直方图的图形

and graph with PDF, CDF and sampling histogram

这篇关于SciPy：von Mises分布在半个圆上？的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！