<jupyter_start><jupyter_text>MVA - MCMC                                             TP2 - Importance Sampling                                        <NAME><jupyter_code># Import the libraries : 
import numpy as np 
%matplotlib inline 
import matplotlib.pyplot as plt 
# Import chi2 distribution 
from scipy.stats import chi
# Import the Normal distribution 
from scipy.stats import norm<jupyter_output><empty_output><jupyter_text>Exercice 3 : Importance sampling   3.A - Poor Importance Sampling :  We want the expectation of the function : $$ f(x)  = 2 \sin\left( \frac{\pi}{1.5} x \right)  \mathbb{1}_{\mathbb{R}^+ } (x) $$  Function f :<jupyter_code>def fx(x) :
    # Function to define the function to plot
    return np.multiply(2*np.sin((np.pi * x)/1.5),x >0) 
df = 1.65
mu = 0.8
std = np.sqrt(1.5)
x = np.linspace(0,5, 1000)
plt.plot(x, chi.pdf(x, df),'r-',label = "p")
plt.plot(x, norm.pdf(x, loc = mu,scale = std ),'b-',label = "q")
plt.plot(x,fx(x)*chi.pdf(x, df),'g-',label = "f * p")
plt.legend()
plt.title(r"$\mu = 0.8$")

df = 1.65
mu = 6
std = np.sqrt(0.5)
x = np.linspace(0,9, 1000)
plt.plot(x, chi.pdf(x, df),'r-',label = "p")
plt.plot(x, norm.pdf(x, loc = mu,scale = std ),'b-',label = "q")
#plt.plot(x,fx(x)*chi.pdf(x, df),'g-',label = "f * p")
plt.legend()
plt.title(r"$\mu = 6$")<jupyter_output><empty_output><jupyter_text>Function to sample from p and q :<jupyter_code>## Sample from a Chi
def samplep(df,N):
    return chi.rvs(df,size = N)
## Sample from a gaussian with keeping only positive values 
def sampleq(mu,var,N): 
    l = []
    while len(l) != N : 
        x = norm.rvs(loc=mu,scale=var)
        if x > 0 : 
            l.append(x)
    return np.array(l)<jupyter_output><empty_output><jupyter_text>1. Implementation of the important sampling method :<jupyter_code>def true_sampling(df, N) : 
    ## Sampling from p : Chi distribution 
    p  = samplep(df,N) 
    ## Evaluate f on p 
    f = fx(p)
    ## Evalute the weights 
    wp  = chi.pdf(p,df) 
    return np.divide(np.sum( np.multiply(f ,wp)) , N ) 
def important_sampling(df,mu,var,N) : 
    ## Sampling from q : normal distribution
    q  = sampleq(mu,var,N) 
    ## Evaluation of f in q 
    f  = fx(q) 
    ## Compute the weights on q : 
    wq = norm.pdf(q,loc = mu, scale = var) 
    wp = chi.pdf(q,df) 
    w  = np.divide(wp,wq) 
    ## I normalize with N as I used the true probabilities :
    w_norm = np.divide(w,N)
    return  np.sum( np.multiply(f , w_norm ))<jupyter_output><empty_output><jupyter_text>2. Comparison of the estimates with respect to N :   Set the parameters :<jupyter_code># Parameter of the Chi distribution 
df = 1.65 
# Parameter of the Gaussian distribution 
mu = 0.8
std = np.sqrt(1.5)
# list of the sample size 
N = [10,100,10**3,10**4,10**5]<jupyter_output><empty_output><jupyter_text>True sampling with different sample size :<jupyter_code>l = [] 
for n in N: 
    l.append(true_sampling(df,n))
print l<jupyter_output>[0.41322920127714691, 0.48997762951558327, 0.53601275231833823, 0.5322356841946122, 0.53967274294581036]<jupyter_text>Important sampling with different sample size :<jupyter_code>l = [] 
for n in N: 
    l.append(important_sampling(df,mu,np.sqrt(1.5),n) /2 )
print l<jupyter_output>[0.69454607665977197, 0.64194455727773236, 0.55421345157103108, 0.52314512477821995, 0.51736142066785551]<jupyter_text>We converge with important sampling to the true estimate.  3. With shifted mean<jupyter_code># Update of the parameter of the Gaussian distribution 
mu = 6
std = np.sqrt(0.5)
l = [] 
for n in N: 
    l.append(important_sampling(df,mu,std,n) /2 )
print l<jupyter_output>[-2.6942207065158731e-05, 3.127669595481433e-05, 0.010787366723203991, 0.0023308467978117565, 0.0038807259904314537]
{'fraction_non_alphanumeric': 0.144962486602358, 'fraction_numerical': 0.09297963558413719, 'mean_word_length': 3.5747549019607843, 'pattern_counts': {'":': 0, '<': 31, '<?xml version=': 0, '>': 33, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '8', 'dir': 'jupyter-structured-clean-dedup', 'id': '76002', 'license': 'MIT', 'n_tokens_mistral': 1705, 'n_tokens_neox': 1465, 'n_words': 350, 'path': 'TP2/tp2_exo3.ipynb', 'repo_name': 'abyoussef/MCMC', 'size': '3748'}