In [8]:
plt.plot(year[keep],temp[keep], marker='.',linewidth=0)
Out[8]:
[<matplotlib.lines.Line2D at 0x73c52c22acf0>]
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
09.01
np.logspace(-20,-10,11)
np.array(10*[2.])
#
2*np.ones(10)
np.zeros(5, dtype=np.float32)
09.02
x=np.random.poisson(100,10)
print (x)
print (x[1:] - x[:-1])
[ 94 107 96 95 91 91 105 95 105 101] [ 13 -11 -1 -4 0 14 -10 10 -4]
09.03
year, temp = np.loadtxt('data/munich_temperatures_average_with_bad_data.txt',
unpack=True)
temp_red=temp[(temp>-90) & (temp<90)]
print (len(temp)-len(temp_red),'elements were removed')
60 elements were removed
10.01
year, temp = np.loadtxt('data/munich_temperatures_average_with_bad_data.txt',
unpack=True)
keep=(temp>-90) & (temp<90)
plt.plot(year[keep],temp[keep], marker='.',linewidth=0)
[<matplotlib.lines.Line2D at 0x73c52c22acf0>]
plt.plot(year[keep]%1,temp[keep], marker='.',linewidth=0)
[<matplotlib.lines.Line2D at 0x73c52c276240>]
10.02
# define normalized Gaussian with mean 'mu' and sigma 'sig'
def gaussian(x,mu,sig):
return (1/(np.sqrt(2*np.pi)*sig))*np.exp(-(x-mu)**2/(2*sig**2))
# Gaussian with mean 10 and sigma sqrt(10)
x=np.linspace(-5,30,100)
d=np.random.normal(7,np.sqrt(7),10000)
_=plt.hist(d,bins=np.linspace(-5,30,36))
plt.plot(x,10000*gaussian(x,7,np.sqrt(7)),color='red',linewidth=3)
[<matplotlib.lines.Line2D at 0x73c52c13fe00>]
# I plot here the sample from the Poisson distribution together with the Gaussian
x=np.linspace(-5,30,100)
d=np.random.poisson(7,10000)
_=plt.hist(d,bins=np.linspace(-5,30,36))
plt.plot(x,10000*gaussian(x,7,np.sqrt(7)),color='red',linewidth=3)
[<matplotlib.lines.Line2D at 0x73c52c1fd460>]
10.03
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
x=np.arange(-10,10,0.1)
y=np.arange(-10,10,0.1)
# Note we calculate the meshgrid first here
# this makes the computations below
# more readable
xx,yy = np.meshgrid(x,y)
r1sq=xx**2+(yy-5)**2 # work out radii
r2sq=xx**2+(yy+5)**2
b1x,b1y = - (yy-5)/r1sq, xx/r1sq # calculate magnetic field vectors
b2x,b2y = + (yy+5)/r2sq, -xx/r2sq
plt.streamplot(x,y,b1x+b2x,b1y+b2y,color='k')
<matplotlib.streamplot.StreamplotSet at 0x73c52c0c88c0>
Solution for the coil
a=5 #radius of the coil
n=90 #elements on the coil
x=np.arange(-10,10,0.05)
y=np.arange(-10,10,0.05)
dphi=2*np.pi/n
xx,yy=np.meshgrid(x,y) #generate coordinate grid
bx=0*xx
by=0*yy
for phi in np.arange(0,2*np.pi,dphi):
cos=np.cos(phi)
#separation coil element-observer
rp=np.sqrt((xx-a*cos)**2+a**2*(1.0-cos**2)+yy**2)
#magnetic field components (cross product)
bx+=-a*cos*yy*dphi/rp**3
by+=-a*(a-cos*xx)*dphi/rp**3
plt.streamplot(x,y,bx,by,color='k')
<matplotlib.streamplot.StreamplotSet at 0x73c4a26adaf0>
# Example "Helmholtz coils"
a=5 #radius of the coil
n=90 #elements on the coil
x=np.arange(-10,10,0.05)
y=np.arange(-10,10,0.05)
dphi=2*np.pi/n
xx,yy=np.meshgrid(x,y) #generate coordinate grid
bx=0*xx
by=0*yy
for phi in np.arange(0,2*np.pi,dphi):
cos=np.cos(phi)
#separation coil element-observer
rp=np.sqrt((xx-a*cos)**2+a**2*(1.0-cos**2)+(yy-a/2)**2)
rp2=np.sqrt((xx-a*cos)**2+a**2*(1.0-cos**2)+(yy+a/2)**2)
#magnetic field components (cross product)
bx+=-a*cos*(yy-a/2)*dphi/rp**3
by+=-a*(a-cos*xx)*dphi/rp**3
bx+=-a*cos*(yy+a/2)*dphi/rp2**3
by+=-a*(a-cos*xx)*dphi/rp2**3
plt.streamplot(x,y,bx,by,color='k')
<matplotlib.streamplot.StreamplotSet at 0x73c4a26f0b60>
11.1
import glob
import os
for file in glob.glob('data/*'):
print (str(int(os.path.getsize(file)/1024+0.5))+'kB', file)
# Defined constant(s)
G = 6.67384e-11 # SI units
# Define values to sample from
mean_m_1 = 40.e4
sigma_m_1 = 0.05e4
mean_m_2 = 30.e4
sigma_m_2 = 0.1e4
mean_r = 3.2
sigma_r = 0.01
Standard error propagation
# Compute mean and error of force
mean_f = G * mean_m_1 * mean_m_2 / mean_r ** 2
sigma_f = mean_f * np.sqrt((sigma_m_1 / mean_m_1) ** 2
+ (sigma_m_2 / mean_m_2) ** 2
+ 4. * (sigma_r / mean_r) ** 2)
print(mean_f, sigma_f)
0.7820906249999999 0.005625407272639265
Generate random samples for m_1, m_2 and r and calculate the force
N = 1000000
m_1 = np.random.normal(mean_m_1, sigma_m_1, N)
m_2 = np.random.normal(mean_m_2, sigma_m_2, N)
r = np.random.normal(mean_r, sigma_r, N)
# calculate force
F = G * m_1 * m_2 / r ** 2
print(np.mean(F), np.std(F))
0.7821141212409272 0.005630767952733316
Plot the Monte-Carlo histogram using the option density=True and the Gaussian from the error propagation
# Define range of output values for plotting
xmin = 0.75
xmax = 0.82
# use Gaussian definition from above!
x = np.linspace(xmin, xmax, 1000)
y = gaussian(x, mean_f, sigma_f)
plt.hist(F, bins=50, range=[xmin, xmax], density=True)
plt.plot(x, y, color='red', lw=3)
plt.xlabel("Force [N]")
plt.ylabel("Relative Probability")
plt.xlim(xmin, xmax)
(0.75, 0.82)
Repeat this with larger error bars:
# Define values to sample from
mean_m_1 = 40.e4
sigma_m_1 = 2.e4
mean_m_2 = 30.e4
sigma_m_2 = 10.e4
mean_r = 3.2
sigma_r = 1.0
# Define range of output values for plotting
xmin = -1.
xmax = 5.
# Define number of samples
N = 1000000
## STANDARD ERROR PROPAGATION
# Compute mean and error of force
mean_f = G * mean_m_1 * mean_m_2 / mean_r ** 2
sigma_f = mean_f * np.sqrt((sigma_m_1 / mean_m_2) ** 2
+ (sigma_m_2 / mean_m_2) ** 2
+ 2. * (sigma_r / mean_r) ** 2)
# Define Gaussian function
x = np.linspace(xmin, xmax, 1000)
y = gaussian(x, mean_f, sigma_f)
## MONTE-CARLO ERROR PROPAGATION
# Sample from initial values
m_1 = np.random.normal(mean_m_1, sigma_m_1, N)
m_2 = np.random.normal(mean_m_2, sigma_m_2, N)
r = np.random.normal(mean_r, sigma_r, N)
# Compute final values
F = G * m_1 * m_2 / r ** 2
## PLOTTING
plt.hist(F, bins=50, range=[xmin, xmax], density=True)
plt.plot(x, y, color='red', lw=3)
plt.xlabel("Force [N]")
plt.ylabel("Relative Probability")
plt.xlim(xmin, xmax)
(-1.0, 5.0)
