SciPy: वॉन मिल्स एक आधे सर्कल पर वितरण करता है?

Aug 17 2020

मैं एक अर्ध-चक्र पर लिपटे वॉन-मिज़ वितरण को परिभाषित करने का सबसे अच्छा तरीका जानने की कोशिश कर रहा हूं (मैं इसका उपयोग विभिन्न सांद्रता में दिशाहीन रेखाएं खींचने के लिए कर रहा हूं)। मैं वर्तमान में SciPy के vonmises.rvs () का उपयोग कर रहा हूं। अनिवार्य रूप से, मैं पीआई / 2 का एक मतलब अभिविन्यास, कहना, सक्षम करना चाहता हूं और वितरण को पीआई / 2 दोनों तरफ से अधिक नहीं करना है।

मैं एक काटे गए सामान्य वितरण का उपयोग कर सकता हूं, लेकिन मैं वॉन-मीज़ की रैपिंग खो दूंगा (यदि मुझे 0 का माध्य अभिविन्यास चाहिए)

मैंने इसे शोध पत्रों में फाइबर मैपिंग के माध्यम से देखते हुए देखा है, लेकिन मैं इसे (अजगर में) कैसे लागू कर सकता हूं, इसका पता नहीं लगा सकता। मैं थोड़ा सा फंस गया हूं कि कहां से शुरू करूं।

यदि मेरे वॉन मेसिस को (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)

मैं इसके बजाय एक आधे-चक्र के चारों ओर एक लपेट का उपयोग करने के लिए इसे कैसे बदलूंगा?

इस अनुभव के साथ किसी को भी कुछ मार्गदर्शन दे सकता है?

जवाब

1 SeverinPappadeux Aug 17 2020 at 17:57

प्रत्यक्ष संख्यात्मक उलटा सीडीएफ नमूना लेने के लिए उपयोगी है, यह बाध्य डोमेन के साथ वितरण के लिए बहुत अच्छा काम करना चाहिए। यहाँ कोड नमूना, पीडीएफ और सीडीएफ तालिकाओं का निर्माण और उलटा सीडीएफ विधि का उपयोग करके नमूना है। अनुकूलित और सदिश किया जा सकता है, निश्चित रूप से

कोड, पायथन 3.8, x64 विंडोज 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()

और पीडीएफ, सीडीएफ और नमूना हिस्टोग्राम के साथ ग्राफ

1 JohanC Aug 17 2020 at 14:45

आप सुन्न के फ़िल्टरिंग के माध्यम से वांछित सीमा के बाहर के मूल्यों को छोड़ सकते हैं ( theta=theta[(theta>=0)&(theta<=np.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()