Simulace matematického kyvadla metodou inverzní funkce

Nasimulujeme histogram poloh (x-ových sořadnic kyvadla)
tak, že budeme generovat náhodnou proměnnou s husotou pravděpodobnosti
$f(x) = \frac{1}{\pi\sqrt{x_0^2-x^2}},$ kde $x_0$ je maximální výchlka kyvadla

Použijeme metodu inverzní funkce. Distribuční funkce odpovídající hustotě pravděpodobnosti $f(x)$ je $F(x)=\frac{1}{2}+\frac{1}{\pi}\text{arcsin}\left (\frac{x}{x_0}\right ).$
Inverzní funkce k distribuční funkci je $F^{-1}(x) = x_0\sin\left (\pi (r-\frac{1}{2})\right )$, kde $r$ je náhodná proměnná s rovnoměrným rozdělením na intervalu $(0,1)$, tj. $r\in \text{U(0,1)}$

import numpy as np
import matplotlib.pyplot as plt

ndata=100000 #počet simulovaných hodnot
nbins=100 #počet binů
x0=1 #max. výchylka
r=np.random.random(ndata) #vygenerujeme náhodné proměnné s rovnoměrným rozdělením U(0,1)
x=x0*np.sin(np.pi*(r-0.5)) #transformujeme náhodnou proměnnou r tak abychom dostali náhodnou proměnnou s požadovaným rozdělením f(x)

plt.hist(x,nbins,density=True,label="data")
#pro srovnání vykreslíme také teoretickou hustotu pravděpodobnosti f(x)
epsilon=0.005 #ćtvrtina binu
#epsilon je tam protože pro x=x0 f(x) diverguje a dělili bychom nulou 
x_teor=np.linspace(-x0+epsilon,x0-epsilon,1000)
y_teor=1/(np.pi*np.sqrt(x0**2-x_teor**2))
plt.plot(x_teor,y_teor,c='red',label="teoretická f(x)")
plt.xlabel("x")
plt.ylabel("f(x)")
plt.legend()
plt.show()