import numpy as np
import matplotlib.pyplot as plt

# Question 2: we indeed observe that the difference decreases as 1/sqrt{N}
Nsamples = np.logspace(2, 6, 5, dtype=int)
Integral = np.zeros_like(Nsamples, dtype=float)
nrepet = 20 # number of times we repeat the MC sampling for averaging
for n, N in enumerate(Nsamples):
    for k in range(nrepet):
        samples = np.random.random(N) # we generate N samples according to a uniform distribution in [0,1]
        Integral[n] += np.mean(1. / samples ** (1. / 3.) + samples / 10.)
    Integral[n] /= nrepet
print(Integral)
print('Exact value =', 31./20.)

plt.figure()
plt.loglog(Nsamples, np.absolute(Integral - 31. / 20.), '-o')
plt.loglog(Nsamples, 1. / np.sqrt(Nsamples), '--')
plt.xlabel('N')
plt.ylabel('|I - 31/20|')
plt.grid()
plt.show()

# Question 4: x=y^{3/2} with y drawn from a uniform distribution between 0 and 1 gives you x drawn from p(x).
def sampling_p(N):
    y = np.random.random(N)
    return y ** 1.5
    
# Question 5: Be careful that you have to divide by p(x) to estimate the integral. You can do this with a pen and a paper before implementation.
# We still have the rror which decreases as 1/sqrt(N), but the error is smaller than for the uniform distribution.
Nsamples = np.logspace(2, 6, 5, dtype=int)
Integral_2 = np.zeros_like(Nsamples, dtype=float)
nrepet = 20 # number of times we repeat the MC sampling for averaging
for n, N in enumerate(Nsamples):
    for k in range(nrepet):
        samples = sampling_p(N) # we generate N samples according to a uniform distribution in [0,1]
        Integral_2[n] += np.mean(1.5 + 3. * samples ** (4. / 3.) / 20.)
    Integral_2[n] /= nrepet
print(Integral_2)

plt.figure()
plt.loglog(Nsamples, np.absolute(Integral - 31. / 20.), '-o', label='uniform')
plt.loglog(Nsamples, np.absolute(Integral_2 - 31. / 20.), '-o', label='p(x)')
plt.loglog(Nsamples, 1. / np.sqrt(Nsamples), '--')
plt.legend()
plt.xlabel('N')
plt.ylabel('|I - 31/20|')
plt.grid()
plt.show()