https://physikmix.lewinb.net for more!

3b.py

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Thu Apr 25 22:40:53 2019

@author: Lewin Bormann, Emily Leung, Kira Jürgens
"""

import matplotlib.pyplot as plt
import numpy.random as rd

import math

def sample_variance(sample):
    avg = sample.sum()/sample.size
    return sum([(x - avg)**2 for x in sample])

samplesize = 10000
sample = rd.exponential(scale=100, size=(samplesize))

avg = sample.sum()/sample.size
error = math.sqrt(1/(samplesize * (samplesize-1)) * sample_variance(sample))

print("Average: {} Error: {}".format(avg, error))

fig = plt.figure()
ax = fig.add_subplot(111)

ax.hist(sample, 100, zorder=0)
ax.grid(True,zorder=1)

sc = ax.scatter([avg], [0], zorder=2)
sc.set_color(((1,0,0,1)))

plt.show()

3c.py

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Thu Apr 25 22:40:53 2019

@author: lbo
"""

import matplotlib.pyplot as plt
import numpy.random as rd
import numpy as np

import math

def sample_variance_sum(sample):
    avg = sample.sum()/sample.size
    return sum([(x - avg)**2 for x in sample])

def stat_error(sample):
    return math.sqrt(1/(samplesize * (samplesize-1)) * sample_variance_sum(sample))

scale = 100

samplesizes, averages, errors, stddevs = [],[],[],[]

for exp in range(1, 10):
    samplesize = exp*100
    sample = rd.exponential(scale=scale, size=(samplesize))
    
    avg = sample.sum()/sample.size
    stddev = np.std(sample)
    error = stat_error(sample)

    samplesizes.append(samplesize)
    averages.append(avg)
    errors.append(error)
    stddevs.append(stddev)
    
    fig = plt.figure()
    ax = fig.add_subplot(111)
    ax.hist(sample, 100, zorder=0)
    ax.grid(True,zorder=1)
    ax.set_title('n = {}'.format(samplesize))
    sc = ax.scatter([avg], [0], zorder=2)
    sc.set_color('red')

    print('n: {} Average: {} Stddev: {} Error: {}'.format(samplesize, avg, stddev, error))

fig = plt.figure()
ax = fig.add_subplot(111, ylim=(0, 2 * scale))
ax.grid(True)
ax.set_ylabel("Avg/Stddev")
avgscatter = ax.plot(samplesizes, averages, 'bo-')[0]
stddevscatter = ax.plot(samplesizes, stddevs, 'go-')[0]

ax2 = fig.add_subplot(111, ylim=(0,max(errors)), sharex=ax, frameon=False)
ax2.yaxis.tick_right()
ax2.yaxis.set_label_position("right")
ax2.set_ylabel("Err")
errscatter = ax2.plot(samplesizes, errors, 'ro-')[0]

ax.legend((avgscatter, stddevscatter, errscatter), ("Avg", "Stddev", "Err"))

plt.show()