SciPy: distribusi von Mises dalam setengah lingkaran?

Aug 17 2020

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

1 SeverinPappadeux Aug 17 2020 at 17:57

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

1 JohanC Aug 17 2020 at 14:45

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