import numpy as np
import scipy.stats 

from matplotlib import pyplot as plt
from IPython.core.pylabtools import figsize
%matplotlib inline

# We want to estimate the rare event 
#   N(0,1)>gamma

# crude monte carlo has high variance
gamma = 3.
nb_samples = 10000
X = np.random.normal(size=nb_samples)
(X>gamma).sum()/nb_samples

TODO: as exercise

def likelihood_ratio(X, mu_u=0, mu_v=gamma):
    assert mu_u==0
    return np.exp(mu_v**2/2) * np.exp(-mu_v*X)

#### importance sampling - low variance
mu_v=gamma
X = np.random.normal(loc=mu_v, size=nb_samples)
X_ = X[X>gamma]
#### imprtance sampling estimator has lower variance
likelihood_ratio(X_).sum() / nb_samples

likelihood_ratio(X_)