In [1]:
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline

1. A straight-line fit

In [2]:
xs, ys = np.loadtxt('ds1.txt')
plt.plot(xs, ys, 's');
In [3]:
# Part (a)
xbar = xs.mean()
ybar = ys.mean()
t0 = np.cov(xs,ys, ddof=1)
varx = t0[0,0]
covar = t0[0,1]
m = covar/varx
c = ybar - m*xbar
xline = np.linspace(-5.0, 5.0, 10)
yline = m*xline + c
print("Gradient = %.5f,\tIntercept = %.5f" %(m,c))
plt.plot(xs, ys, 's')
plt.plot(xline, yline, 'r-', linewidth=3);
Gradient = 1.47958,	Intercept = 0.45793
In [4]:
# Part (b)

2. A linear model

In [5]:
xs, ys = np.loadtxt('ds2.txt')
plt.plot(xs, ys, '.');
In [6]:
# (a)
N = len(xs)
m = 4 # the number of free parameters in a cubic
X = np.zeros(shape=(N, m))
for i in range(N):
    X[i,:] = np.array([xs[i]**k for k in range(m)])
In [7]:
# (b) The normal equations are in the form A x = b,
# with A = X^T X and b = X^T y and x the vector of parameters.
A = np.dot(X.transpose(), X)
b = np.dot(X.transpose(), ys)
par = np.linalg.solve(A, b) # solve the linear system to find the parameters
print("Cubic coefficients: ", par)

n = 100
xline = np.linspace(-4.1,4.1,n)
yline = np.zeros(n)
for i in range(m):
    yline += par[i]*xline**i
plt.plot(xs, ys, '.')
plt.plot(xline, yline, 'r-', lw=3)
Cubic coefficients:  [ 0.96909809  1.49844166  0.00643006 -0.19808062]
Out[7]:
[<matplotlib.lines.Line2D at 0x1133bc0f0>]
In [8]:
# (c)
from scipy.optimize import curve_fit
def f(x, a,b,c,d):
    return a + b*x + c*x**2 + d*x**3

par2, cov = curve_fit(f, xs, ys, p0=(1,1.5,0,-0.2))
print("From normal equations:     ", par)
print("From optimize.curve_fit(): ", par2)
From normal equations:      [ 0.96909809  1.49844166  0.00643006 -0.19808062]
From optimize.curve_fit():  [ 0.96909809  1.49844166  0.00643006 -0.19808062]

3. A non-linear model

In [9]:
xs, ys = np.loadtxt('ds3.txt')
plt.plot(xs, ys, '.');
In [10]:
def f(x, b0, b1, b2, b3):
    return b0 + b1*np.sin(b2*x + b3)

# guess the params by looking at the plot above.
par0 = [0.4, 0.4, 2, 0]
par, cov = curve_fit(f, xs, ys, par0)
print("Parameters : ", par)
n = 100
xline = np.linspace(-4.1,4.1,n)
yline = np.zeros(n)
for i in range(n):
    yline[i] = f(xline[i], par[0],par[1],par[2],par[3])
plt.plot(xs, ys, '.')
plt.plot(xline, yline, 'r-', lw=3);
Parameters :  [ 0.39605484  0.30269268  2.80655875  0.52262687]