""" Gaussverteilte Zufallszahlen nach Box-Muller """
import numpy as np
import matplotlib.pyplot as plt

fig,ax = plt.subplots(3,1)
npoints = 10000

# 1. Variante: mit Auswertung von sin/cos
u1 =  np.random.rand( npoints )
u2 =  np.random.rand( npoints )
sqrtlog = np.sqrt( -2*np.log( u1 ) )
x = sqrtlog * np.cos( 2*np.pi*u2 )
y = sqrtlog * np.sin( 2*np.pi*u2 )

# beide Saetze von Zufallszahlen zusammenfuegen und als Histogram darstellen
gauss = np.concatenate( ( x, y ) )
ax[0].hist( gauss, 20 )

# 2. Variante: Polarmethode, ohne sin/cos
u1 =  np.random.rand( npoints )
u2 =  np.random.rand( npoints )

# zwei Zufallszahlen aus [-1,1[
v1 = 2*u1 - 1
v2 = 2*u2 - 1

# Verwerfungsmethode: akzeptieren, wenn in Einheitskreis
# Projektion auf neue Arrays mithilfe dieses Akzeptanzkriteriums
vsquare = v1**2 + v2**2
condition = np.logical_and( vsquare > 0, vsquare <= 1 )
v1 = v1[ condition ]
v2 = v2[ condition ]
vsquare = vsquare[ condition ]
v = np.sqrt( vsquare )

# alternative Berechnung von sin/cos
cos = v1/v
sin = v2/v
sqrtlog = np.sqrt( -2*np.log( vsquare ) )
x = sqrtlog * cos
y = sqrtlog * sin

# beide Saetze von Zufallszahlen zusammenfuegen und als Histogram darstellen
gauss = np.concatenate( ( x, y ) )
ax[1].hist( gauss, 20 )

# 3. Variante: korrelierte gaussverteilte Zufallszahlen mittels Cholesky-Zerlegung der inversen Kovarianzmatrix
mu = np.array( [ 2.0, -1.0 ] )
sigma_x = 2.0
sigma_y = 0.5
rho_xy  = -0.5
V = np.array( [ [ sigma_x**2,             rho_xy*sigma_x*sigma_y ],
                [ rho_xy*sigma_x*sigma_y, sigma_y**2 ] ] )
Vinv = np.linalg.inv( V )
L    = np.linalg.cholesky( Vinv )
Linv = np.linalg.inv( L )

# Vorbereitung der Daten: unkorrelierte gaussverteilte Zufallszahlen
mu0     = np.array( [ 0.0, 0.0 ] ) 
Vuncorr = np.eye( 2 )
x, y = np.random.multivariate_normal( mu0, Vuncorr, size=npoints ).T
xy = np.array( list( zip( x, y ) ) )

# Korrelierte Zufallszahlen ueber Cholesky-Zerlegung
xcorr = []
ycorr = []
for pair in xy:
    corr = np.dot( Linv, pair ) + mu
    xcorr.append( corr[0] )
    ycorr.append( corr[1] )

ax[2].scatter( xcorr, ycorr, s=1 )
ax[2].set_xlim( mu[0] - 5, mu[0] + 5 )
ax[2].set_ylim( mu[1] - 5, mu[1] + 5 )

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