SciPy: distribusi von Mises dalam setengah lingkaran?
Saya mencoba mencari cara terbaik untuk mendefinisikan distribusi von-Mises yang dibungkus dalam setengah lingkaran (saya menggunakannya untuk menggambar garis tanpa arah pada konsentrasi yang berbeda). Saya saat ini menggunakan vonmises.rvs () SciPy. Pada dasarnya, saya ingin dapat memasukkan, katakanlah, orientasi rata-rata pi / 2 dan distribusinya terpotong tidak lebih dari pi / 2 di kedua sisi.
Saya dapat menggunakan distribusi normal yang terpotong, tetapi saya akan kehilangan pembungkus von-mises (katakanlah jika saya ingin orientasi rata-rata 0)
Saya telah melihat ini dilakukan di makalah penelitian yang melihat orientasi serat pemetaan, tetapi saya tidak tahu cara menerapkannya (dalam python). Saya agak bingung harus mulai dari mana.
Jika von Mesis saya didefinisikan sebagai (dari numpy.vonmises):
np.exp(kappa*np.cos(x-mu))/(2*np.pi*i0(kappa))
dengan:
mu, kappa = 0, 4.0
x = np.linspace(-np.pi, np.pi, num=51)
Bagaimana saya mengubahnya menjadi menggunakan pembungkus setengah lingkaran sebagai gantinya?
Adakah orang yang memiliki pengalaman dengan hal ini menawarkan beberapa panduan?
Jawaban
Berguna untuk memiliki pengambilan sampel CDF terbalik numerik langsung, seharusnya berfungsi baik untuk distribusi dengan domain terbatas. Berikut contoh kode, pembuatan tabel PDF dan CDF dan pengambilan sampel menggunakan metode CDF terbalik. Bisa dioptimalkan dan di-vektorisasi, tentunya
Kode, 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()
dan grafik dengan PDF, CDF dan histogram sampling

Anda dapat membuang nilai di luar rentang yang diinginkan melalui pemfilteran numpy ( theta=theta[(theta>=0)&(theta<=np.pi)]
, memperpendek larik sampel). Jadi, Anda dapat menambah jumlah sampel yang dihasilkan terlebih dahulu, lalu memfilter dan kemudian mengambil subarray dengan ukuran yang diinginkan.
Atau Anda bisa menambah / mengurangi pi untuk memasukkan semuanya ke dalam kisaran itu (melalui theta = np.where(theta < 0, theta + np.pi, np.where(theta > np.pi, theta - np.pi, theta))
). Seperti dicatat oleh @SeverinPappadeux, hal itu mengubah distribusi dan mungkin tidak diinginkan.
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()
