PHYS 391 - Day 6

Part 1 - SDOM demonstration

Generation parameters

import math
import random
mean = 1.2
var = 0.2

Generate 25 random numbers, then make a histogram of this data to see what we have.

%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
import statistics as stat

# Set bins explicitly
binlo = 0.6
binhi = 1.8
binwidth = 0.1
binvec = np.arange(binlo, binhi, binwidth)
binvec2 = np.arange(binlo, binhi, binwidth/5.)
meanvec = []

Repeat this as many times as desired

N=25
rvec = [random.gauss(mean, var) for i in range(N)]
plt.hist(rvec, bins=binvec)
plt.xlabel('Voltage (V)')
plt.ylabel('Readings')
plt.show()

meanvec.append(stat.mean(rvec))
plt.hist(meanvec, bins=binvec2)
plt.xlabel('Mean Voltage (V)')
plt.ylabel('Trials')
plt.show()

print('N: {0:3d}'.format(len(meanvec)))
print('SDOM:          {0:.3f}'.format(round(stat.stdev(meanvec), 3)))
print('sigma/root(N): {0:.3f}'.format(var/math.sqrt(N)))
print('diff:          {0:.2f}%'.format(100*(round(stat.stdev(meanvec), 3)/(var/math.sqrt(N))-1)))

N:  50
SDOM:          0.041
sigma/root(N): 0.040
diff:          2.50%

Part 2 - Flat parent distribution

Try again with a flat (uniform) distribution

meanvec = []

N = 25
rvec = [random.uniform(mean-1.7*var, mean+1.7*var) for i in range(N)]
plt.hist(rvec, bins=binvec)
plt.xlabel('Voltage (V)')
plt.ylabel('Readings')
plt.show()

meanvec.append(stat.mean(rvec))
plt.hist(meanvec, bins=binvec2)
plt.xlabel('Mean Voltage (V)')
plt.ylabel('Trials')
plt.show()

print('N: {0:3d}'.format(len(meanvec)))
print('SDOM: {0:.3f}'.format(round(stat.stdev(meanvec), 3)))
print('sigma/root(N): {0:.3f}'.format(var/math.sqrt(N)))
print('diff: {0:.2f}%'.format(100*(round(stat.stdev(meanvec), 3)/(var/math.sqrt(N))-1)))

N:  40
SDOM: 0.042
sigma/root(N): 0.040
diff: 5.00%

Part 3 - Small N

Try again with a smaller number of events

N = 5
meanvec = []
ssdvec = []

for trials in range(1000):
    rvec = [random.uniform(mean-1.7*var, mean+1.7*var) for i in range(N)]
    meanvec.append(stat.mean(rvec))
    ssdvec.append(stat.stdev(rvec))

plt.hist(rvec, bins=binvec)
plt.xlabel('Voltage (V)')
plt.ylabel('Readings')
plt.show()

plt.hist(meanvec, bins=binvec2)
plt.xlabel('Mean Voltage (V)')
plt.ylabel('Trials')    
plt.show()

plt.hist(ssdvec, bins=np.arange(0.0, 0.4, 0.01))
plt.xlabel('Sample SD (V)')
plt.ylabel('Trials')    
plt.show()


print('N: {0:3d}'.format(len(meanvec)))
print('mean(SD): {0:.3f}'.format(stat.mean(ssdvec)))
print('SDOM: {0:.3f}'.format(round(stat.stdev(meanvec), 3)))
print('sigma/root(N): {0:.3f}'.format(var/math.sqrt(N)))
print('diff: {0:.2f}%'.format(100*(round(stat.stdev(meanvec), 3)/(var/math.sqrt(N))-1)))    

N: 1000
mean(SD): 0.185
SDOM: 0.092
sigma/root(N): 0.089
diff: 2.86%