★ Monte Carlo Simulation To Calculate PI ★



In [1]:

    
# Import modules
import time
import math
import numpy as np
import scipy
import matplotlib.pyplot as plt

Necesaary Function For Monte Carlo Simulation



In [2]:

    
def linear_congruential_generator(x, a, b, m):
    x = (a * x + b) % m
    u = x / m
    return u, x, a, b, m

def stdrand(x):
    return linear_congruential_generator(x, pow(7, 5), 0, pow(2, 31) - 1)[:2]

def halton(p, n):
    b = np.zeros(math.ceil(math.log(n + 1) / math.log(p)))
    u = np.zeros(n)
    for j in range(n):
        i = 0
        b[0] = b[0] + 1
        while b[i] > p - 1 + np.finfo(float).eps:
            b[i] = 0
            i += 1
            b[i] += 1 
        u[j] = 0
        for k in range(1, b.size + 1):
            u[j] = u[j] + b[k-1] * pow(p, -k)
    return u

Monte Carlo Simulation (with Minimal standard random number generator)



In [3]:

    
def monte_carlo_process_std(toss):
    x = time.time()
    hit = 0
    for i in range(toss):
        u1, x = stdrand(x)
        u2, x = stdrand(x)
        if pow(u1, 2) + pow(u2, 2) < 1.0:
            hit += 1
    return hit * 4.0 / toss

pi = monte_carlo_process_std(2000000)
print('pi = %.10f, err = %.10f' %(pi, abs(pi - np.pi)))









    



pi = 3.1411940000, err = 0.0003986536

Monte Carlo Simulation (with LCG where multiplier = 13, offset = 0 and modulus = 31)



In [4]:

    
def monte_carlo_process_customized(toss):
    x0 = time.time()
    args = (x0, 13, 0, 31)
    hit = 0
    for i in range(toss):
        u1, *args = linear_congruential_generator(*args)
        u2, *args = linear_congruential_generator(*args)
        if pow(u1, 2) + pow(u2, 2) < 1.0:
            hit += 1
    return hit * 4.0 / toss

pi = monte_carlo_process_customized(2000000)
print('pi = %.10f, err = %.10f' %(pi, abs(pi - np.pi)))









    



pi = 3.1255760000, err = 0.0160166536

Monte Carlo Simulation (with Quasi-random numbers)



In [5]:

    
def monte_carlo_process_quasi(toss):
    hit = 0
    px = halton(2, toss)
    py = halton(3, toss)
    for i in range(toss):
        u1 = px[i]
        u2 = py[i]
        if pow(u1, 2) + pow(u2, 2) < 1.0:
            hit += 1
    return hit * 4.0 / toss

pi = monte_carlo_process_quasi(2000000)
print('pi = %.10f, err = %.10f' %(pi, abs(pi - np.pi)))









    



pi = 3.1415720000, err = 0.0000206536