In [1]:
%pylab inline
import pandas as pd
import scipy.optimize as op

Populating the interactive namespace from numpy and matplotlib
In [4]:
a = pd.read_csv("Vibrationbeam.csv")
a['L'] = a['L']/100
a['L_unc'] = 0.01
a['F_unc'] = 0.1
a
Out[4]:
L F L_unc F_unc
0 0.1558 102.0 0.01 0.1
1 0.1558 646.5 0.01 0.1
2 0.1558 103.8 0.01 0.1
3 0.1558 103.0 0.01 0.1
4 0.1558 104.0 0.01 0.1
5 0.1650 53.0 0.01 0.1
6 0.1650 91.0 0.01 0.1
7 0.1650 91.0 0.01 0.1
8 0.1650 91.0 0.01 0.1
9 0.1650 91.0 0.01 0.1
10 0.1745 81.5 0.01 0.1
11 0.1745 82.0 0.01 0.1
12 0.1745 82.0 0.01 0.1
13 0.1745 82.0 0.01 0.1
14 0.1745 81.5 0.01 0.1
15 0.1845 73.5 0.01 0.1
16 0.1845 73.5 0.01 0.1
17 0.1845 73.5 0.01 0.1
18 0.1845 73.5 0.01 0.1
19 0.1845 73.5 0.01 0.1
20 0.2162 53.5 0.01 0.1
21 0.2162 53.5 0.01 0.1
22 0.2162 53.5 0.01 0.1
23 0.2162 53.0 0.01 0.1
24 0.2162 53.5 0.01 0.1
25 0.2455 42.0 0.01 0.1
26 0.2455 42.0 0.01 0.1
27 0.2455 42.0 0.01 0.1
28 0.2455 42.0 0.01 0.1
29 0.2455 42.5 0.01 0.1
In [44]:
a = a.loc[delete(a.index,[1,5])]
a
Out[44]:
L F L_unc F_unc
0 0.1558 102.0 0.01 0.1
2 0.1558 103.8 0.01 0.1
3 0.1558 103.0 0.01 0.1
4 0.1558 104.0 0.01 0.1
6 0.1650 91.0 0.01 0.1
7 0.1650 91.0 0.01 0.1
8 0.1650 91.0 0.01 0.1
9 0.1650 91.0 0.01 0.1
10 0.1745 81.5 0.01 0.1
11 0.1745 82.0 0.01 0.1
12 0.1745 82.0 0.01 0.1
13 0.1745 82.0 0.01 0.1
14 0.1745 81.5 0.01 0.1
15 0.1845 73.5 0.01 0.1
16 0.1845 73.5 0.01 0.1
17 0.1845 73.5 0.01 0.1
18 0.1845 73.5 0.01 0.1
19 0.1845 73.5 0.01 0.1
20 0.2162 53.5 0.01 0.1
21 0.2162 53.5 0.01 0.1
22 0.2162 53.5 0.01 0.1
23 0.2162 53.0 0.01 0.1
24 0.2162 53.5 0.01 0.1
25 0.2455 42.0 0.01 0.1
26 0.2455 42.0 0.01 0.1
27 0.2455 42.0 0.01 0.1
28 0.2455 42.0 0.01 0.1
29 0.2455 42.5 0.01 0.1
In [45]:
plot(a['L'],a['F'],'o')
Out[45]:
[<matplotlib.lines.Line2D at 0x26899b154e0>]
In [46]:
errorbar(a['L'],a['F'],a['F_unc'],fmt='r.')
ylabel("frequency (Hz)")
xlabel("length (m)")
Out[46]:
<matplotlib.text.Text at 0x26899a8f9b0>
In [47]:
x = a['L']
y = a['F']
y_unc = a['F_unc']

def f(x,m,b):
    return m*x + b

# do the fit

# method 0: no uncertainty included
##ans = op.curve_fit(f,x,y,p0=[1.0,1.0])

# method 1: use uncertainty as relative weight
ans = op.curve_fit(f,x,y,p0=[1.0,1.0], sigma=y_unc)

# method 2 (preferred): use uncertainty as absolute sigma
ans = op.curve_fit(f,x,y,p0=[1.0,1.0], sigma=y_unc, absolute_sigma=True)


# interpret results of fit
m,b = ans[0]
m_unc,b_unc = sqrt(diag(ans[1]))

print("m = %.2f +/- %.2f"%(m,m_unc))
print("b = %.2f +/- %.2f"%(b,b_unc))

# plot results
errorbar(x,y,y_unc, fmt='r.', label="data")
xx = linspace(0.15,.26)
plot(xx,f(xx,m,b),label="fit")
legend(loc="upper right")
xlabel("length (m)")
ylabel("frequency (Hz)")


# plot residuals (obs - fit results)
figure()
errorbar(x, y-f(x,m,b), y_unc, fmt="r.")
hlines(0.1,0.15,.26)
xlabel("length (m)")
ylabel("residual (y - fit)")
ylim(-100,150)

m = -652.00 +/- 0.61
b = 197.96 +/- 0.12
Out[47]:
(-100, 150)
In [50]:
x = a['L']
y = a['F']
y_unc = a['F_unc']

def f(x,m,b):
    return m*x**b

# do the fit

# method 0: no uncertainty included
##ans = op.curve_fit(f,x,y,p0=[1.0,1.0])

# method 1: use uncertainty as relative weight
ans = op.curve_fit(f,x,y,p0=[1.0,1.0], sigma=y_unc)

# method 2 (preferred): use uncertainty as absolute sigma
ans = op.curve_fit(f,x,y,p0=[1.0,1.0], sigma=y_unc, absolute_sigma=True)


# interpret results of fit
m,b = ans[0]
m_unc,b_unc = sqrt(diag(ans[1]))

print("m = %.2e +/- %.2e"%(m,m_unc))
print("b = %.2e +/- %.2e"%(b,b_unc))

# plot results
errorbar(x,y,y_unc, fmt='r.', label="data")
xx = linspace(0.15,.26)
plot(xx,f(xx,m,b),label="fit")
legend(loc="upper right")
xlabel("length (m)")
ylabel("frequency (Hz)")


# plot residuals (obs - fit results)
figure()
errorbar(x, y-f(x,m,b), y_unc, fmt="r.")
hlines(0.1,0.15,.26)
xlabel("length (m)")
ylabel("residual (y - fit)")
ylim(-100,150)

m = 2.59e+00 +/- 9.21e-03
b = -1.98e+00 +/- 2.04e-03
Out[50]:
(-100, 150)
In [ ]: