Assignment 6

General feedback

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.

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.

Somes student used too few sample points.

Monte Carlo integration

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.

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