SciPy: von Mises กระจายครึ่งวงกลม?

Aug 17 2020

ฉันพยายามหาวิธีที่ดีที่สุดในการกำหนดการแจกแจงแบบฟอน - มิเซสที่ห่อด้วยครึ่งวงกลม (ฉันใช้มันเพื่อวาดเส้นไร้ทิศทางที่ความเข้มข้นต่างกัน) ฉันกำลังใช้ vonmises.rvs () ของ SciPy โดยพื้นฐานแล้วฉันต้องการที่จะใส่, พูด, การวางแนวค่าเฉลี่ยของ pi / 2 และให้การกระจายถูกตัดให้ไม่เกิน pi / 2 ทั้งสองข้าง

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)

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.

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?

คำตอบ

1 SeverinPappadeux Aug 17 2020 at 17:57

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

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()

and graph with PDF, CDF and sampling histogram

1 JohanC Aug 17 2020 at 14:45

You could discard the values outside the desired range via numpy's filtering (theta=theta[(theta>=0)&(theta<=np.pi)], shortening the array of samples). So, you could first increment the number of generated samples, then filter and then take a subarray of the desired size.

Or you could add/subtract pi to put them all into that range (via theta = np.where(theta < 0, theta + np.pi, np.where(theta > np.pi, theta - np.pi, theta))). As noted by @SeverinPappadeux such changes the distribution and is probably not desired.

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()