Low-Dimensional Fittings¶
We use the classical least-square method to fit low degree polynomials to the data.
See e.g. Golub and Van Loan (1996) and Trefethen and Bau III (1997) for details on the least-square method.
Importing the libraries¶
import csv
import numpy as np
import matplotlib.pyplot as plt
Loading the data¶
This time we use the csv library to read the data from file and use numpy to work with the data.
Loading with csv¶
We first load the data into a python list with csv.
file = open('water.csv',"r")
water_csv = list(csv.reader(file, delimiter=","))
water_csv
Header and datapoints¶
We separate the header from the datapoints and have the datapoints put into a numpy array as follows.
header = dict([(water_csv[0][i],water_csv[1][i]) for i in range(3)])
print(header)
datapoints = np.array(water_csv[2:]).astype("float")
print(datapoints)
Linear approximation¶
Given a set of temperature and density data in the form $(T_j, \rho_j)$, we look for a linear relation $\rho^{(1)}(T) = c + m T$ which is the "best fit" for the data.
One way to approach this problem is to interpret the best fit as minimizing the sum of squares of the residuals. The residual for each $j$ measurement is
$$ r_j = \rho_j - \rho^{(1)}(T_j), $$and the sum of squares of the residuals is
$$ \min_{c,m\in \mathbb{R}} \sum_j (\rho_j - \rho^{(1)}(T))^2. $$Matrix form¶
This can be written in matrix form as
$$ \displaystyle \min_{\mathbf{u}\in \mathbb{R}^2} \|A\mathbf{u} - \mathbf{f}\|_2^2, $$where $\|\cdot\|_2$ is the Euclidian norm of a vector and
$$ A = \left[ \begin{matrix} T_1 & 1 \\ \vdots & 1 \\ T_n & 1 \end{matrix}\right], \qquad \mathbf{u} = \left( \begin{matrix} m \\ c \end{matrix}\right), \qquad \mathbf{f} = \left( \begin{matrix} \rho_1 \\ \vdots \\ \rho_n \end{matrix} \right). $$The matrix $A$ is a simple Vandermonde type matrix obtained from the temperature data, $\mathbf{u}$ is the unknown vector with the desired coefficients for the approximation, and $\mathbf{f}$ is the vector with the density measurements.
The Vandermonde matrix¶
We use the numpy function numpy.vstack() to build the Vandermonde matrix $A$.
T = datapoints[:,0]
A1 = np.vstack([T,np.ones(len(T))]).T
print(A1)
The density measurements¶
The density data is the second column of the data array.
f = datapoints[:,1]
print(f)
Solution of the least-square problem for the linear approximation¶
We use the numpy function numpy.linagl.lstsq() to solve the least-square problem.
m, d = np.linalg.lstsq(A1, f, rcond=None)[0]
print(m,d)
Visualizing the result¶
Now we plot the linear approximation along with the data to visualize the quality of the approximation
plt.figure(figsize=(10,5))
plt.plot(T, f, 'o', label='Data', color='tab:blue')
plt.plot(T, m*T + d, 'b', label='Linear approximation', color='tab:red')
plt.title('Plot of the data and of the linear approximation', fontsize=14)
plt.xlabel(header['temp'], fontsize=12)
plt.ylabel(header['density'], fontsize=12)
plt.legend()
plt.show()
Second-degree approximation¶
One can see from the plot above that the first-degree approximation doesn't seem to be a very good approximation. We look, now, for a second-degree approximation.
The least-square problem¶
For the second-degree approximation, we look for a second-degree polynomial $\rho^{(2)}(T) = aT^2 + bT + c$ that best approximates the data in the sense of minimizing the sum of the square of the residuals
$$ r_j = \rho_j - \rho^{(2)}(T_j) $$Matrix form¶
This can be written in matrix form as
$$ \displaystyle \min_{\mathbf{u}\in \mathbb{R}^2} \|A\mathbf{u} - \mathbf{f}\|_2^2, $$where $\mathbf{f}$ is as before but the Vandermonde matrix and the vector of unknowns take the form
$$ A = \left[ \begin{matrix} T_1^2 & T_1 & 1 \\ \vdots & 1 \\ T_n^2 & T_n & 1 \end{matrix}\right], \qquad \mathbf{u} = \left( \begin{matrix} a \\ b \\ c \end{matrix}\right). $$Vandermonde matrix¶
In this case, the Vandermonde matrix is
A2 = np.vstack([T**2, T,np.ones(len(T))]).T
print(A2)
Solution¶
a, b, c = np.linalg.lstsq(A2, f, rcond=None)[0]
print(a,b,c)
Visualizing the result¶
plt.figure(figsize=(10,5))
plt.plot(T, f, 'o', label='Data', color='tab:blue')
plt.plot(T, a*T**2 + b*T + c, 'r', label='Quadratic approximation', color='tab:green')
plt.title('Plot of the data and of the quadratic approximation', fontsize=14)
plt.xlabel(header['temp'], fontsize=12)
plt.ylabel(header['viscosity'], fontsize=12)
plt.legend()
plt.show()
This seems much better
Comparing the two approximations¶
Visual comparison¶
Visually, the second-degree approximation is way better.
plt.figure(figsize=(10,5))
plt.plot(T, f, 'o', label='Data', color='tab:blue')
plt.plot(T, m*T + d, 'b', label='Linear approximation', color='tab:red')
plt.plot(T, a*T**2 + b*T + c, 'r', label='Quadratic approximation', color='tab:green')
plt.title('Plot of the data along with the approximations', fontsize=14)
plt.xlabel(header['temp'], fontsize=12)
plt.ylabel(header['viscosity'], fontsize=12)
plt.legend()
plt.show()
Comparing the quadratic error¶
The quadratic error is the sum of the square of the residual errors for each measurement, i.e.
$$ \Delta = \sum_j r_j^2 = \|A\mathbf{u} - \mathbf{f}\|_2^2, $$for the best approximation obtained. But this error is not normalized by the number of measurements. This is achieved with the mean quadratic error:
$$ E = \frac{1}{N} \sum_j r_j^2. $$N = len(f)
print(f'Number of measurements: {N}\n')
print(f'Quadratic error for the linear approximation: {np.linalg.lstsq(A1, f, rcond=None)[1][0]:.2e}')
print(f'Quadratic error for the quadratic approximation: {np.linalg.lstsq(A2, f, rcond=None)[1][0]:.2e}\n')
print(f'Mean quadratic error for the linear approximation: {np.linalg.lstsq(A1, f, rcond=None)[1][0]/N:.2e}')
print(f'Mean quadratic error for the quadratic approximation: {np.linalg.lstsq(A2, f, rcond=None)[1][0]/N:.2e}')