Assignment 6

General feedback

Google/GenAI 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.

Distributions

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.

Some student used too few sample points.

Monte Carlo integration

Generally well solved.

Code failing because they had changed the “import numpy” to “import numpy as np”.

Some cases of using variables that were not defined? This usually happens when you had defined the variables and then removed them. Upon the kernel refreshing, the code then fails.

Some students did not do a loglog plot, or failed to take the absolute value of the error, or did not know how to show the 1/sqrt(N) relationship

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.

Autograded tasks

Distribution

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)

Monte Carlo Integration

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

Plot examples

here are example plots for the assignment:

Distribution

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")

Monte Carlo Integration

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");