Simulace náhodné proměnné von Neumannovou zamítací metodou

import numpy as np
import matplotlib.pyplot as plt

rng=np.random.default_rng()
print(rng)

Generator(PCG64)

simulujeme náhodnou proměnnou s hustotou pravděpodobnosti
$$f(x)=\left (\frac{\sin(\pi x)}{\pi x}\right )^2$$

# simulované rozdělení
def sinc(x):
    return (np.sin(np.pi*x)/(np.pi*x))**2

jako majorizační funkci použijeme (nenormovaný) lorenzián
$$g(x)=\frac{w}{w^2+(x-x_0)^2}$$

#majorizační funkce
def major(x,x0,w):
    return w/(w**2+(x-x0)**2)

integrál majorizační funkce je
$$G(x)=\int_{-\infty}^x g(t)\, dt = \text{arctg}\left (\frac{x-x0}{w}\right )+\frac{\pi}{2}$$

# integrál majorizační funkce "něco jako distribuční funkce majorizační funkce"  
def G(x,x0,w):
    return np.atan((x-x0)/w)+np.pi/2

inverzní funkce k integrálu majorizační funkce je $$G^{-1}(r)=w\text{tg}\left (r-\frac{\pi}{2}\right ) + x_0$$

# inverzní funkce k integrálu majorizační funkce 
def Ginv(r,x0,w):
    return w*np.tan(r-np.pi/2)+x0

x0=0
w=1
fig,ax=plt.subplots()
xp=np.linspace(-5,5,1000)
plt.plot(xp,sinc(xp),label="f(x)")
plt.plot(xp,major(xp,x0,w),label="majorizační funkce g(x)")
ax.set_xlabel("x")
ax.set_ylabel("f(x)")
ax.legend()
plt.show()

fig,ax=plt.subplots()
plt.plot(xp,G(xp,x0,w),c='orange',label="integrál majorizační funkce G(x)")
ax.set_xlabel("x")
ax.set_ylabel("G(x)")
ax.legend()
plt.show()

#výběr a inicializace generátoru náhodných čísel
rng=np.random.default_rng()
print(rng)

Generator(PCG64)

x=[]
N=100000 #počet nasimulovaných hodnot
i_done=0 #kolik už máme nasimulovaných hodnot
while i_done<N:
    r1=rng.random()*np.pi
    x1=Ginv(r1,x0,w)
    r2=rng.random()*major(x1,x0,w)
    if r2<=sinc(x1): #akceptujeme, jinak odmítneme
        x=np.append(x,x1)
        i_done+=1

fig,ax=plt.subplots()
ax.hist(x,range=[-5,5],bins=100,density=True,label="simulovaný histogram")
ax.plot(xp,sinc(xp),c='red',label="teoretická pdf")
ax.set_xlabel("x")
ax.set_ylabel("f(x)")
ax.legend()
plt.show()