SciPy: Yarım daire üzerinde von Mises dağılımı?
Yarım daireye sarılmış bir von-Mises dağılımını tanımlamanın en iyi yolunu bulmaya çalışıyorum (farklı konsantrasyonlarda yönsüz çizgiler çizmek için kullanıyorum). Şu anda SciPy's vonmises.rvs () kullanıyorum. Esasen, diyelim ki ortalama bir pi / 2 oryantasyonu koyabilmek ve dağılımın her iki tarafta pi / 2'den fazla olmayacak şekilde kesilmesini istiyorum.
Kesilmiş bir normal dağılım kullanabilirim, ancak von miseslerin sarılmasını kaybedeceğim (ortalama 0 yönelim istiyorsam diyelim)
Bunun, fiber yönelimlerini haritalandırmaya bakan araştırma makalelerinde yapıldığını gördüm, ancak nasıl uygulanacağını bulamıyorum (python'da). Nereden başlayacağım konusunda biraz takılıp kaldım.
Benim von Mesis'im şu şekilde tanımlanmışsa (numpy.vonmises'den):
np.exp(kappa*np.cos(x-mu))/(2*np.pi*i0(kappa))
ile:
mu, kappa = 0, 4.0
x = np.linspace(-np.pi, np.pi, num=51)
Bunun yerine yarım dairenin etrafına sarmak için nasıl değiştirebilirim?
Bu konuda biraz deneyime sahip olan herhangi biri size rehberlik edebilir mi?
Yanıtlar
Doğrudan sayısal ters CDF örneklemesine sahip olmak yararlıdır, sınırlı alanla dağıtım için harika çalışmalıdır. İşte kod örneği, PDF ve CDF tabloları oluşturma ve ters CDF yöntemini kullanarak örnekleme. Elbette optimize edilebilir ve vektörleştirilebilir
Kod, 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()
ve PDF, CDF ve örnekleme histogramı içeren grafik

Numpy'nin filtrelemesi yoluyla istenen aralığın dışındaki değerleri atabilirsiniz theta=theta[(theta>=0)&(theta<=np.pi)]
(örnek dizisini kısaltarak). Böylece, önce üretilen örneklerin sayısını artırabilir, ardından filtreleyebilir ve ardından istenen boyutta bir alt dizi alabilirsiniz.
Ya da hepsini bu aralığa (aracılığıyla theta = np.where(theta < 0, theta + np.pi, np.where(theta > np.pi, theta - np.pi, theta))
) koymak için pi ekleyebilir / çıkarabilirsiniz . @SeverinPappadeux tarafından belirtildiği gibi, bu tür dağıtımı değiştirir ve muhtemelen istenmez.
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()
