Codes for Crude Monte Carlo Integration Approach and Importance Sampling for a given set of function

from math import exp

from random import random

import numpy as np

import matplotlib.pyplot as plt

from scipy.stats import norm

from IPython.display import clear_output

Q1.

Monte Carlo Integration Approach

We define a function fh which is equal to f(x) × h(x) to calculate an accurate we know E[h(x)] = integration of fh, which is the primary solution of our question!

#Function for Question 1 and it's inverse(not used)

def fh1(x):
    
    ret = (exp(-2 * abs(x-5)))
    
    return ret

fh = np.vectorize(fh1)

def invfh1(x1):
    
    num = -(math.log(x)/2)
    
    if num < 0:
        
        ret = 5 + num
    else:
        ret = 5 - num
    
    return ret

invfh = np.vectorize(invfh1)

#Visualizing the function

X = np.linspace(0, 10, 1000)

plt.plot(X, fh(X))

plt.show()

Code of the function to calculate the crude Monte Carlo estimation of the given function

x1 = Start value of the integrand

x2 = End Value of the integrand

func = function to be used (default = fh)

n = number of iterations (default = 10000)

def mcintegrate(x1, x2, func = fh, n = 10000):
     
    ylist = []
    
    for i in range(n):
        
        x = x1 + (x2 - x1) * random()
        
        y = func(x)
        
        ylist.append(y)
            
    return(np.mean(ylist)*(x2-x1), np.var(ylist)*(x2-x1) 

We iterate the function 100 times to get a much more accurate view.

mc = []
mcvar = []

for i in range(0,100):
    
    mc.append(mcintegrate(0,10)[0])
    mcvar.append(mcintegrate(0,10)[1])
    
answer = np.mean(mc)
var = np.var(mcvar)

print(" The approximate value of Integrand is: " + str(answer))
    The approximate value of Integrand is: 0.9976833402743391
print(" The mean of the var of the 100 tests : " + str(var))
    The mean of the var of the 100 tests : 0.00014382726284691412

Q1.(ii) Importance sampling approach. Use N(5,1) as a candidate probability distribution, i.e., g(x).

Code to calculate the value of integration using Importance sampling approach

x1 = Start value of the integrand

x2 = End Value of the integrand

func = function to be used (default = fh)

n = number of iterations (default = 10000)

def isintegrateq1(x1, x2, n = 10000, func = fh):
    
    running_total = []
    
    for i in range(n):
        
        rand_norm = np.random.normal(loc = 5, scale = 1)
        
        running_total.append(func(rand_norm)/norm.pdf(rand_norm, loc = 5))
    
    approximation = np.mean(running_total)
    variance = np.var(running_total)
    
    return approximation, variance

IS1 = []
IS1var = []
for i in range(0,100):
    
    IS1.append(isintegrateq1(0,10)[0])
    IS1var.append(isintegrateq1(0,10)[1])
    
answer = np.mean(IS1)

var = np.mean(IS1var)
print(" The approximate value of Integrand is: " + str(answer))
 The approximate value of Integrand is: 0.9997532248288528
print(" The mean of the var of the 100 tests : " + str(var))
 The mean of the var of the 100 tests : 0.3572089251079059

Both the solutions are really, really close to the actual answer = 1 (Checked on Wolfram)

Q1. (iii) Compare the answers of (i) and (ii) for different sample sizes.

We get a more accurate answer in the Importance Sampling as it is closer to 1 than the answer in Crude Monte Carlo integration.

Q2.

We define a function f which is equal to f(x) = 0.5×exp(-|x|) to calculate an accurate view of the integrand with different moments

#Function for Question 2 with differing moments

def f1(x, m):
    
    ret = 0.5 * exp(-1 * abs(x)) * x**(m)
    
    return ret

f = np.vectorize(f1)

#Visualizing the function

X = np.linspace(-10, 10, 1000)

plt.plot(X, f(X, 0))

plt.show()

Code to calculate the value of integration using the Importance Sampling approach for the given function

n = number of iterations

m = moment of the function

def isintegrateq2(n, m):
    
    running_total = []
    
    for i in range(n):
        
        rand_norm = np.random.normal(loc = 0, scale = 4)
        
        running_total.append(f(rand_norm, m)/(norm.pdf(rand_norm, loc = 0, scale = 4)))
                
    approximation = np.mean(running_total)
    variance = np.var(running_total)
    
    return approximation, variance

# run simulation just to check once

n = 10000

approx = isintegrateq2(n, 2)

approx

We iterate the function 100 times to get a much more accurate view.

We get a list of means and variances of the 100 samples.

meanlist = [ [] for _ in range(3) ]

varlist = [ [] for _ in range(3) ]

for z in range(0,3):
    
    n = [10000, 1000, 100]
    
    for i in range(0,5):
    
        IS2 = []
        IS2var = []

        for j in range(0,100):

            IS2.append(isintegrateq2(n[z], i+1)[0])
            IS2var.append(isintegrateq2(n[z], i+1)[1])
        
        answer = np.mean(IS2)

        var = np.mean(IS2var)
        
        meanlist[z].append(answer)
        varlist[z].append(var)

Creating a Pandas DataFrame to visualize better

from pandas import DataFrame

df_mean = DataFrame(meanlist,columns=['M = 1', 'M = 2', 'M = 3', 'M = 4', 'M = 5'])

df_var = DataFrame(varlist,columns=['M = 1', 'M = 2', 'M = 3', 'M = 4', 'M = 5'])

df_mean.index = ['N = 10000', 'N = 1000', 'N = 100']

df_var.index = ['N = 10000', 'N = 1000', 'N = 100']

print(df_mean)

print(df_var)

I have got all the Solutions answered

Somewhat a bad practice though for the triple loop! O(n^3) This is not my homework! I seek out to find such questions to toughen up my Coding and Maths skills!

This is the conclusion of the Article!