Random number generation from arbitrary PDF

Numerical determination of CDF

Given some PDF

\[f: \Omega \rightarrow \mathbb{R}\]

which may or may not be normalized, we can determine the CDF using the function nbi_stat.eval_cdf(). This function integrates the PDF from a lower to an upper bound and returns the cumulative integral as an array. This array can late be fed to nbi_stat.sample_pdf() to draw random numbers from the PDF.

Sampling a PDF

The function nbi_stat.sample_pdf() takes uniform random variables as input and produces random numbers according to the passed evaluated CDF.

Let’s make an example. As our PDF we take a single resonance peak over a quadratic background

\[f(E,a_1,a_2,a_3,A,\Gamma,E_0) = a_1+a_2x+a_3x^2+A\frac{\Gamma/(2\pi)}{(E-E_0)^2+(\Gamma/2)^2}\]

with

\[\begin{split}a_1 &= -0.3\\ a_2 &= 73\\ a_3 &= -18\\ A &= 33.8\\ \Gamma &= 0.138\\ E_0 &= 0.9667\end{split}\]
>>> import nbi_stat
>>>
>>> def pdf(E,a1=-.3,a2=73,a3=-18,A=33.8,Gamma=0.138,E0=0.9667):
>>>     return a1+a2*E+a3*E**2+A*(Gamma/(2*np.pi))/((E-E0)**2+(Gamma/2)**2)
>>>
>>> x = np.linspace(0,3,31)
>>> cdf = nbi_stat.eval_cdf(pdf,x)
>>> o = nbi_stat.sample_pdf(np.random.random(200),x,cdf)
>>>
>>> h, b, w, e = nbi_stat.histogram(o,x)
>>>
>>> plt.errorbar(b,h,e,w/2,".",label="Samples")
>>> plt.plot(x,pdf(x),label="PDF")
>>> plt.legend()
>>> plt.xlabel(r"$E$")

(Source code, png, hires.png, pdf)

_images/random-1.png