ex1.py

import random
import math


class my_random_vector:
    def __init__(self,N_dim,min_value = -5, max_value = 5):
        self.vec = [min_value + (max_value-min_value)*random.random() for i in range(N_dim)]

    def norm_squared(self):
        nn = 0
        for x in self.vec:
            nn+= x**2
        return nn

def pdf(vec):
    return math.exp(-0.5*vec.norm_squared())


def integrate_normal(n_dim, n_tests, n_sigma = 5):
    n_tests = int(n_tests)
    integral = 0
    for test in range(n_tests):
        integral += pdf(my_random_vector(n_dim,-n_sigma,n_sigma))
    return 1./n_tests * integral * (2*n_sigma)**n_dim



for n_dim in [1,2,5,15]:
    mc = integrate_normal(n_dim,1e5,3)
    truth = math.sqrt((2*math.pi)**n_dim)
    print("N_dim = {}".format(n_dim))
    print("Truth = {}, Monte-Carlo = {}".format(truth,mc))
    print("Difference : {:.2f}%".format(100*(truth-mc)/truth))