Monte Carlo integral approximation - importance sampling Question.

% matplotlib inline
from sympy import *
from scipy import random
import numpy as np
import matplotlib.pyplot as plt

a = 0
b = 1        
M = 100000   
xrand = np.zeros(M)
yrand = np.zeros(M)
zrand = np.zeros(M)

for i in range(len(xrand)):
    xrand[i] = random.uniform(a,b)
    yrand[i] = random.uniform(a,b)
    zrand[i] = random.uniform(a,b)

import random

N = 900
def monte_carlo(the_list) :
    if len(the_list) > N :
        return the_list
    else :
        x = random.random()
        y = random.random()
        z = random.random()
        if p(y)/p(x) > z :
            the_list.append(y)    
            return monte_carlo(the_list)
        else :
            the_list.append(x)   
                # Recursive
            return monte_carlo(the_list)

new_list = []
print(monte_carlo(new_list)) 

from sympy import *
x = Symbol('x')
y = Symbol('y')
z = Symbol('z')
a = Symbol('a')

def f(x):
    return 3 * x **2
def p(x):
    return 2 * x 
def g(x) :
    return f(x)/p(x)

Integral(f(x),x).doit()
pprint(Integral(f(x),x).doit()) 
pprint(Integral(f(x), (x,0,1)).doit()) # accurate integration value
pprint(Integral(p(x), (x,0,1)).doit()) # Normalization check

# Attempt to Solve Function Static Equations with Parcel Spanning Method

% % matplotlib inline
from sympy import *
from scipy import random
import numpy as np
import matplotlib.pyplot as plt

a = 0
b = 1 # limit of integration 0 to pi function integration
M = 100000 #Number of splits
xrand = np.zeros(M)

for i in range(len(xrand)):
    xrand[i] = random.uniform(a,b)

integral = 0.0
for i in range(M):
    integral += f(xrand[i])

answer = (b-a)/float(M)*integral
error = abs(1-answer)
print("The integral from 0 to 1 of f(x): ", answer)
print("error: ", error)


import random

# Try importance sampling.

'''
 Dear Yupto.

 The basic idea is to obtain the integral(f(x))dx (from 0 to 1) as an approximation process using random numbers.
 f(x) = (f(x)/p(x)) * p(x) = g(x) * p(x) where p(x) is called weight function or probability distribution function (pdf).

 At this time, if you extract N x values that satisfy a specific determination condition (conditional expression in the code below) and obtain the average value of g(x), 
 " Mean value of g(x) = the original static integral value of f(x)" is known to be valid. The larger the N value, the smaller the difference between the two values.

 (For your information, assume that f(x) = 3x^2, p(x) = 2x and g(x) = 3x^2/2x.)
 (So the final value to get close is the static value of 1.)
''' 

# If you give an empty list, if you satisfy a specific expression, you will return up to N random values in the list.
N = 900
def monte_carlo(the_list) :
    if len(the_list) > N :
        return the_list
    else :
        x = random.random()
        y = random.random()
        z = random.random()
        If p(y)/p(x) > z : # p(X) is a first-order polynomial in this case, so y/x is fine.
            The_list.append(y) # Re-run this function itself
            return monte_carlo(the_list)
        else :
            Put the_list.append(x) #x on the list,
                # Recursive
            return monte_carlo(the_list)

new_list = []
#print(monte_carlo(new_list)) 

I = 0
for i in range(N):
    I += g(monte_carlo(new_list)[i])/ N
print(I)