SciPy: distribuzione di von Mises su un semicerchio?

Question

Sto cercando di capire il modo migliore per definire una distribuzione di von Mises avvolta su un semicerchio (lo sto usando per disegnare linee senza direzione a diverse concentrazioni). Attualmente sto usando vonmises.rvs () di SciPy. Essenzialmente, voglio essere in grado di inserire, diciamo, un orientamento medio di pi / 2 e avere la distribuzione troncata a non più di pi / 2 su entrambi i lati.

Potrei usare una distribuzione normale troncata, ma perderò il wrapping dei von-mises (diciamo se voglio un orientamento medio di 0)

Ho visto questo fatto in documenti di ricerca che esaminano la mappatura degli orientamenti delle fibre, ma non riesco a capire come implementarlo (in Python). Sono un po 'bloccato su da dove cominciare.

Se il mio von Mesis è definito come (da numpy.vonmises):

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

con:

mu, kappa = 0, 4.0

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

Come lo altererei per utilizzare invece un avvolgimento attorno a un semicerchio?

Qualcuno con una certa esperienza con questo potrebbe offrire una guida?

Answer 1

È utile avere un campionamento CDF inverso numerico diretto, dovrebbe funzionare alla grande per la distribuzione con dominio limitato. Ecco un esempio di codice, creazione di tabelle PDF e CDF e campionamento utilizzando il metodo CDF inverso. Potrebbe essere ottimizzato e vettorializzato, ovviamente

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

e grafico con PDF, CDF e istogramma di campionamento

Answer 2

È possibile scartare i valori al di fuori dell'intervallo desiderato tramite il filtro di numpy ( theta=theta[(theta>=0)&(theta<=np.pi)], accorciando l'array di campioni). Quindi, potresti prima incrementare il numero di campioni generati, quindi filtrare e quindi prendere un sottoarray della dimensione desiderata.

Oppure puoi aggiungere / sottrarre pi greco per metterli tutti in quell'intervallo (tramite theta = np.where(theta < 0, theta + np.pi, np.where(theta > np.pi, theta - np.pi, theta))). Come notato da @SeverinPappadeux tale modifica la distribuzione e probabilmente non è desiderata.

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