#Using Box Muller to generate pairs of independent standard Normals

import numpy as np
import matplotlib.pyplot as plt

N = 5000 # Size of sample
u1 = np.random.uniform(size = N)
u2 = np.random.uniform(size = N)

# Box Muller
z1 = np.multiply(np.sqrt(-2 * np.log(u1)),np.cos(2*pi*u2))
z2 = np.multiply(np.sqrt(-2 * np.log(u1)),np.sin(2*pi*u2))

# Checking the marginal distributions of z1,z2 
p1 = plt.figure(1)
plt.hist(z1, bins = 'auto')
plt.title('Histogram of z1')
#p1.show()

p2 = plt.figure(2)
plt.hist(z2, bins = 'auto')
plt.title('Histogram of z2')
#p2.show()

#Checking the joint distribution of z1, z2
p3 = plt.figure(3)
area = np.pi*3
plt.scatter(z1,z2, s = area, c = np.arange(N), alpha = 0.5)
plt.title('Scatter plot of joint distribution of z1, z2')
plt.show()

# Cholesky decomposition to create correlated Normals

rho = 0.4
x = z1
y = rho*z1 + np.sqrt(1 - rho**2) * z2

# Checking the marginal distributions of x,y 
p4 = plt.figure(4)
plt.hist(x, bins = 'auto')
plt.title('Histogram of x')

p5 = plt.figure(5)
plt.hist(y, bins = 'auto')
plt.title('Histogram of y')

#Checking the joint distribution of x,y
p6 = plt.figure(6)
area = np.pi*3
plt.scatter(x,y, s = area, c = np.arange(N), alpha = 0.5)
plt.title('Scatter plot of joint distribution of x, y')
plt.show()

# Calculate the probablity of x < 0, y > 0 and compare to theoretical result 

from math import *

count = 0
for i in range(N):
    if (x[i] <0 ) and (y[i] > 0) :
        count += 1

print('The Monte Carlo simulated probability is ', count/N)
theoretical_prob = 1/4 - 1/(2 * pi) * atan(rho/sqrt(1-rho**2))
print('The theoretical probability is', theoretical_prob )

The Monte Carlo simulated probability is  0.1876
The theoretical probability is 0.1845050597827727