Problem

Compute the following 10-dimensional integral via Monte Carlo Rejection sampling method,

\[I = \int_{x_1 = 0}^{x_1 = 1} dx_1 \cdots \int_{x_{10} = 0}^{x_{10} = 1} dx_{10} \bigg(\sum_{i=1}^{i=10} ~ x_i ~ \bigg) ~,\]

Ensure the accuracy of your integration result is better than $0.00001$ by progressively comparing the integral with actual integration result $I = 155 / 6$.

Solution

Python

Here is an example Python script that computes the above integral: main.py,

import numpy as np

ndim = 10
xmin = 0
xmax = 1

getFunc = lambda point: np.sum(point)**2
maxHeightFunc = getFunc(ndim*[max(xmin,xmax)])
boxVolume = maxHeightFunc
numAcceptedRejected = 0
numAccepted = 0
truth = 155/6.
integral = []

while True:

    numAcceptedRejected += 1

    if maxHeightFunc * np.random.rand() <= getFunc(np.random.rand(ndim)): numAccepted += 1
    integral.append(boxVolume * numAccepted / numAcceptedRejected)

    delta = np.abs(1 - integral[-1] / truth)
    #if numAcceptedRejected%100 == 0: print(numAccepted, numAcceptedRejected, delta, integral[-1])
    if numAcceptedRejected > 1 and delta <= 0.00001:
        #print(np.abs((integral[-1] - integral[-2]) / integral[-1]))
        break

print("integral = {}, relative accuracy = {}".format(integral[-1],np.double(delta)))

import seaborn as sns
import matplotlib.pyplot as plt
sns.set()
plt.figure()
ax = plt.plot   ( range(1,len(integral)+1)
                , integral
                )
plt.xscale("log")
plt.title("integral = {},\n relative accuracy = {}".format(integral[-1],np.double(delta)))
plt.xlabel("log10( Number of Monte Carlo Rejection Sampling Steps )")
plt.ylabel("Integral Result")

Here is a plot of the integral progression,

Comments