SciPy : 반원의 폰 미제스 분포?
반원으로 감싼 von-Mises 분포를 정의하는 가장 좋은 방법을 찾으려고합니다 (다른 농도에서 방향없는 선을 그리는 데 사용하고 있습니다). 현재 SciPy의 vonmises.rvs ()를 사용하고 있습니다. 본질적으로, 예를 들어 pi / 2의 평균 방향을 입력하고 분포가 pi / 2 이하로 잘 리기를 원합니다.
잘린 정규 분포를 사용할 수 있지만 von-mises의 래핑을 잃을 것입니다 (평균 방향 0을 원하는 경우).
섬유 방향 매핑을 연구하는 연구 논문에서이 작업을 수행 한 것을 보았지만 구현하는 방법을 알 수 없습니다 (파이썬에서). 나는 어디서부터 시작해야할지에 대해 약간 고민하고 있습니다.
내 von Mesis가 (numpy.vonmises에서) 다음과 같이 정의 된 경우 :
np.exp(kappa*np.cos(x-mu))/(2*np.pi*i0(kappa))
와:
mu, kappa = 0, 4.0
x = np.linspace(-np.pi, np.pi, num=51)
대신 반원을 감싸도록 변경하려면 어떻게해야합니까?
이 경험이있는 사람이 지침을 제공 할 수 있습니까?
답변
직접 수치 역 CDF 샘플링을하는 데 유용하며 제한된 도메인을 사용한 배포에 적합합니다. 다음은 코드 샘플, PDF 및 CDF 테이블 작성 및 역 CDF 방법을 사용한 샘플링입니다. 물론 최적화 및 벡터화 가능
코드, 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()
PDF, CDF 및 샘플링 히스토그램이있는 그래프
numpy의 필터링 ( theta=theta[(theta>=0)&(theta<=np.pi)], 샘플 배열 단축)을 통해 원하는 범위 밖의 값을 버릴 수 있습니다. 따라서 먼저 생성 된 샘플 수를 늘린 다음 필터링 한 다음 원하는 크기의 하위 배열을 가져올 수 있습니다.
또는 pi를 더하거나 빼서 모두 해당 범위에 넣을 수 있습니다 (를 통해 theta = np.where(theta < 0, theta + np.pi, np.where(theta > np.pi, theta - np.pi, theta))). @SeverinPappadeux에서 언급했듯이 이러한 배포는 변경되므로 바람직하지 않을 수 있습니다.
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()
