import numpy as np
import time
import matplotlib.pyplot as plt
%matplotlib inline

plt.rcParams['figure.figsize'] = 10, 10

# Linear Congruential random number generator
prev = 0

# From: https://en.wikipedia.org/wiki/Linear_congruential_generator#Parameters_in_common_use
a = 1103515245
b = 12345
c = 2 ** 31

def random_seed(x):
    global prev
    prev = (a * int(x) + b) % c

# generate number between 0 - 1
def random_gen():
    global prev
    x = (a * prev + b) % c
    prev = x
    return x / c

N = 10000

# seed
random_seed(time.time())

# generate two vectors of random uniform numbers
x1 = np.array([random_gen() for i in range(N)])
x2 = np.array([random_gen() for i in range(N)])

plt.scatter(x1, x2)
plt.xlabel("x1")
plt.ylabel("x2")
plt.show()

# apply transformation
y1 = np.sqrt(-2 * np.log(x1)) * np.sin(x2)
y2 = np.sqrt(-2 * np.log(x1)) * np.cos(x2)

plt.scatter(y1, y2)
plt.xlabel("y1")
plt.ylabel("y2")
plt.show()

plt.scatter(y1 * x1, y2 * x1)
plt.xlabel("y1 * x1")
plt.ylabel("y2 * x1")
plt.show()

C_1 = np.array([np.mean(y1), np.mean(y2)])
C_1

array([0.57667055, 1.05876345])

C_2 = [np.mean(y1 ** 2), np.mean(y2 ** 2)] - C_1 ** 2
C_2

array([0.22210667, 0.34833728])

C_3 = [np.mean(y1 ** 3), np.mean(y2 ** 3)] - (3 * C_1 * C_2) - (C_1 ** 3)
C_3

array([0.12587177, 0.16467505])