SciPy: distribución de von Mises en un semicírculo?

Aug 17 2020

Estoy tratando de encontrar la mejor manera de definir una distribución de von-Mises envuelta en un semicírculo (la estoy usando para dibujar líneas sin dirección en diferentes concentraciones). Actualmente estoy usando vonmises.rvs () de SciPy. Esencialmente, quiero poder poner, digamos, una orientación media de pi / 2 y tener la distribución truncada a no más de pi / 2 en cualquier lado.

Podría usar una distribución normal truncada, pero perderé la envoltura de von-mises (digamos si quiero una orientación media de 0)

He visto esto en artículos de investigación que analizan el mapeo de orientaciones de fibra, pero no puedo encontrar la manera de implementarlo (en Python). Estoy un poco atascado sobre por dónde empezar.

Si mi von Mesis se define como (de 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)

¿Cómo lo alteraría para usar una envoltura alrededor de un semicírculo en su lugar?

¿Alguien con algo de experiencia con esto podría ofrecer alguna orientación?

Respuestas

1 SeverinPappadeux Aug 17 2020 at 17:57

Es útil tener un muestreo CDF inverso numérico directo, debería funcionar muy bien para la distribución con dominio acotado. Aquí hay un ejemplo de código, la creación de tablas PDF y CDF y el muestreo utilizando el método CDF inverso. Podría optimizarse y vectorizarse, por supuesto

Código, 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()

y gráfico con PDF, CDF e histograma de muestreo

1 JohanC Aug 17 2020 at 14:45

Puede descartar los valores fuera del rango deseado a través del filtrado de numpy ( theta=theta[(theta>=0)&(theta<=np.pi)]acortando la matriz de muestras). Por lo tanto, primero puede incrementar el número de muestras generadas, luego filtrar y luego tomar un subarreglo del tamaño deseado.

O puede sumar / restar pi para ponerlos todos en ese rango (vía theta = np.where(theta < 0, theta + np.pi, np.where(theta > np.pi, theta - np.pi, theta))). Como señaló @SeverinPappadeux, tales cambios cambian la distribución y probablemente no sean deseados.

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