# -*- coding: utf-8 -*-
"""
Created on Mon Nov 13 17:14:54 2017

@author: florent

In this text, for X=N(0,1), we compute P(X>a) in 3 different ways:
- we give the exact value of the probability thanks to the function scipy.stats.norm.sf
- we give a standard Monte Carlo estimate out of n independant simulations of X
- we give an importance sampling estimate out of n independant simulations of X'=N(theta,1) for theta=a
"""

import numpy as np
import scipy.stats

a=6.
n=int(1e7)

print("\n"+"Exact value : ")

real_p=scipy.stats.norm.sf(a)

print(str(real_p))

print("\nStandard Monte Carlo : ")

X=np.random.randn(n)

pest=np.mean(X>a)

sigma = np.sqrt(pest*(1-pest))

print(str(pest)+" +/- "+str(1.96*sigma/np.sqrt(n)))

print("\nIS Monte Carlo : ")

theta=a+.095

X+=theta

Ech=(X>a)*(np.exp(-theta*X+.5*theta**2))

pest=np.mean(Ech)

sigma=np.std(Ech)

er=1.96*sigma/np.sqrt(n)

print(str(pest)+" +/- "+str(er))

print("\nExact value in CI : ")

print(pest-er<real_p<pest+er)