Kernel > Restart & Clear output
before adding any changes to git!In this tutorial, we make our lives a bit easier and use some pre-defined fitting functions.
# Numerical operations:
import numpy as np
# Plotting library:
import matplotlib.pyplot as plt
# Minimizing/fitting library:
import scipy.optimize
np.array
by typing
np.array?
Fitting a polynomial to a point cloud.
For this we use the polyfit
function of numpy (https://docs.scipy.org/doc/numpy/reference/generated/numpy.polyfit.html).
Let's start by defining some points and their x and y projections:
points = np.array([(1, 1), (2, 4), (3, 1), (9, 3), (10, 10)])
x = points[:,0]
y = points[:,1]
Plot them:
Fitting a polynomial of degree deg
is easy as calling:
coeffs = np.polyfit(x, y, deg=3)
The vector coeffs
now contains the coefficients of the polynomial as a vector of length deg+1
:
coeffs
It's easy to also get the corresponding function:
f = np.poly1d(coeffs)
Let's finally have a look at the fit:
# Some points where we evaluate our new function
x_new = np.linspace(x[0], x[-1], 50)
y_new = f(x_new)
# Plot the datapoints
plt.plot(x, y, 'ko', label="data")
# Plot our fitted polynomial
plt.plot(x_new, y_new, 'r-', label="fit")
# Add legend etc
plt.legend()
plt.xlim([x[0]-1, x[-1] + 1 ])
plt.show()
x_exercise3 = np.linspace(0, 1, 20)
y_exercise3 = x_exercise3 + np.random.random(len(x_exercise3))
plt.plot(x_exercise3, y_exercise3, ".")
Again let's generate some points:
x = np.asarray(range(10))
y = 0.3*x + np.asarray([0,1,2,3,4,5,4,3,2,1])
Here, we define our own function that we want to fit to our data points:
def gaus(x, norm, mean, sigma):
# Note that this function takes a whole vector x of data!
return norm * np.exp(-(x-mean)**2/(2*sigma**2))
This function has 3 free parameters, norm
, mean
, sigma
, which we will now fit using scipy.optimize.curve_fit
(https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.curve_fit.html)
popt, pcov = scipy.optimize.curve_fit(gaus, x, y)
The variable popt
holds the parameters:
popt
The other variable, pcov
holds the covariance matrix of the fitted values.
This gets relevant once you want to give errors on your fitted quantities, but we'll ignore this for now.
# Plot points
plt.plot(x, y,'ko',label='data')
# Plot our gaussian
# Define some points on the xaxis
x_fine = np.linspace(min(x), max(x), 200)
y_fine = gaus(x_fine, *popt)
# The * unpacks the values from popt and uses them as
# parameters.
plt.plot(x_fine, y_fine, 'r-', label='fit')
# ...
plt.legend()
plt.show()
f(x) = x + b + gaus(x, norm, mean, sigma)
Hint:
Create a new function mymodel
which takes the parameters x
, as well as all fitted parameters, i.e. b
, norm
, mean
, sigma
. In the definition you can also use the function gaus
from above!
def mymodel(x, b, norm, mean, sigma):
return x + b + gaus(x, norm, mean, sigma)
b2
line(x, a1, b1, a2, c)
x
is a vector! Thus if x<c: ...
won't work. Rather, take a look at the function np.where
.
def line(x, a1, b1, a2, c):
return np.where(x<c, a1*x+b1, a2*x+a1*c+b1-a2*c)