""" Poisson-verteilte Zufallszahlen """
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import poisson

mu = 2

x1 = np.arange(11)
x2 = np.linspace(0,10,200)

fig,ax = plt.subplots( 3, 1 )
ax[0].bar( x1, poisson.pmf( x1, mu ) )
ax[1].plot( x2, poisson.cdf( x2, mu ) )

""" Effizientes Verfahren fuer Poisson-Zufallszahlen:
Fuer nabhaengige exponentiell verteilte Ereignisse ist die 
Zahl der Ereignisse in einem festen Intervall poissonverteilt mit Parameter mu
1) Generiere exponentiell verteilte Zeiten t_i zwischen Ereignissen mit Mittelwert 1
2) Suche kleinstes I fuer das \sum_{i=1}^{I+1} t_i > mu
 """
npoints = 1000
randexp  = -np.log( np.random.rand( 2*mu*npoints ) )
randpois = np.empty( npoints )

rsum  = 0
i     = 0
i0    = 0
ipois = 0

while( ipois < npoints ):
    rsum += randexp[i]
    if rsum > mu:
        randpois[ ipois ] = i - i0 - 1
        i0 = i 
        rsum = 0
        ipois += 1
    i += 1

ax[2].hist( randpois, bins = x1-0.5 )

## plt.savefig( 'random_poisson.pdf' )
plt.show()