montecarlo.py

import math
import numpy as np

# integration interval
a = 0
b = 1

Width = b - a

# points to sample
nmax = 6942111

# the function we'd like to average on the interval
def f(x):
    return 4 / (1 + x ** 2)

# variables to keep track of the current average and number of points used
ave = 0
n = 0

# the loop that does the averaging
while nmax > n:
    x = np.random.uniform(a,b)
    ave = f(x) / (n + 1) + n * ave / (n + 1)
    n += 1
    print("The number of points sampled is now {:d}".format(n))
    print("The most recent random number used in the interval [{:f},{:f}] -> {:f}".format(a,b,x))
    print("The average value of the function on all the sampled points is {:f}".format(ave))
    print("The Monte-Carlo estimate for the integral is {:f}".format(Width*ave))
    print("==")