Google is often the best tool to find things in the python and matplotlib documentation if you need to find how to do something, it’s what I do.
Students sometimes misunderstood the method, opting for scatter plots to represent the data instead of histograms.
Some students filled the sample with the values $f(x)$ instead of filling it with the $x$ values.
There were a few submissions with inappropriate maximum, which leads to inefficiencies if very large, or incorrect results if too small.
Not all plots showed a reference curve for the histogram.
Somes student used too few sample points.
Generally well solved.
Some students did not do a loglog
plot, or failed to take the absolute value of the error.
Some students plotted the number of dimensions on the x-axis and a different curve for each number of points.
Many students did not use a list/array for the different dimensions and instead wrote them each out. You can create lists/arrays that you can then call in loops, for example an array with the volumes in the different dimensions.
def genSample(npts):
sample = []
# find the maximum
xs = numpy.linspace(0,10,1000)
fs = f(xs)
maxf = max(fs) + 0.01 # small margin of safety
# create the sample
while len(sample) < npts:
x = 10 * random.random()
r = random.random()
if r * maxf < f(x) :
sample.append(x)
return numpy.array(sample)
def integrate(npoints, dim):
# the random numbers
rs = numpy.random.uniform(-1, 1, size=(npoints,dim))
radii = numpy.sum( rs * rs, axis = 1)
#print(radii)
number_in = 0
for radius in radii:
if radius < 1:
number_in += 1
return number_in/npoints * 2**dim
here are example plots for the assignment:
sample = genSample(10**5)
plt.hist(sample, bins=40, alpha=0.5, density=True);
plt.plot(xs,fs,'k--')
plt.xlabel('x')
plt.ylabel('Sampling density');
plt.title("Histogram of the generated sample")
ns = [2**ii for ii in range(4,20)]
dimensions = range(2,7)
targets = {
2: numpy.pi,
3: 4/3*numpy.pi,
4: 1/2*numpy.pi**2,
5: 8/15*numpy.pi**2,
6: 1/6*numpy.pi**3,
}
for dim in dimensions:
vals = [ abs(targets[dim] - integrate(ii, dim)) for ii in ns]
plt.plot(ns, vals, label="D={}".format(dim))
plt.plot(ns,[ 1/nn**(1/2) for nn in ns],'k-',linewidth=3,label="$1/\sqrt{N}$")
plt.loglog()
plt.legend()
plt.title("Error scaling for Monte Carlo integration")
plt.xlabel("Number of points")
plt.ylabel("Monte Carlo integration error");