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$.
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,
