# 3.5. Least-squares Fitting to Data#

**References:**

Chapter 4

*Least Squares*of [Sauer, 2022], Sections 1 and 2.Section 8.1

*Discrete Least Squares Approximation*of [Burden*et al.*, 2016].

```
import numpy as np
from matplotlib.pyplot import figure, plot, title, legend, grid, loglog
from numpy.random import random
from numericalMethods import solveLinearSystem
from numericalMethods import evaluatePolynomial
```

```
import numericalMethods as nm
```

We have seen that when trying to fit a curve to a large collection of data points, fitting a single polynomial to all of them can be a bad approach. This is even more so if the data itself is inaccurate, due for example to measurement error.

Thus an important approach is to find a function of some simple form that is close to the given points but not necesarily fitting them exactly: given \(N\) points

we seek a function \(f(x)\) so that the errors at each point,

are “small” overall, in some sense.

Two important choices for the function \(f(x)\) are

(a) polynomials of low degree, and

(b) periodic sinusidal functions;

we will start with the simplest case of fitting a straight line.

## 3.5.1. Measuring “goodness of fit”: several options#

The first decision to be made is how to measure the overall error in the fit, since the error is now a vector of values \(e = \{e_i\}\), not a single number. Two approaches are widely used:

*Min-Max*: minimize the maximum of the absolute errors at each point, \(\displaystyle \|e\|_{max} \text{ or } \|e\|_\infty, = \max_{1 \leq i \leq n} |e_i|\)*Least Squares*: Minimize the sum of the squares of the errors, \(\displaystyle \sum_1^n e_i^2\)

## 3.5.2. What doesn’t work#

Another seemingly natural approach is:

Minimize the sum of the absolute errors, \(\displaystyle \|e\|_1 = \sum_1^n |e_i|\)

but this often fails completely. In the following example, all three lines minimize this measure of error, along with infinitely many others: any line that passes below half of the points and above the other half.

```
xdata = [1., 2., 3., 4.]
ydata = [0., 3., 2., 5.]
figure(figsize=[12,6])
plot(xdata, ydata, 'b*', label="Data")
xplot = np.linspace(0.,5.)
ylow = xplot -0.5
yhigh = xplot + 0.5
yflat = 2.5*np.ones_like(xplot)
plot(xplot, ylow, label="low")
plot(xplot, yhigh, label="high")
plot(xplot, yflat, label="flat")
legend(loc="best");
```

The Min-Max method is important and useful, but computationally difficult. One hint is the presence of absolute values in the formula, which get in the way of using calculus to get equations for the minimum.

Thus the easiest and most common approach is **Least Squares**, or equivalently, minimizing the root-mean-square error, which is just the Euclidean length \(\|e\|_2\) of the error vector \(e\).
That “geometrical” interpretation of the goal can be useful.
So we start with that.

## 3.5.3. Linear least squares#

The simplest approach is to seek the straight line \(y = f(x) = c_0 + c_1x\) that minimizes the total square sum error,

Note well that the unknowns here are just the two values \(c_0\) and \(c_1\), and \(E\) is s fairly simple polynomial function of them. The minimum error must occur at a critical point of this function, where both partial derivatives are zero:

These are just simultaneous linear equations, which is the secret of why the least squares approach is so much easier than any alternative. The equations are:

where of course \(\sum_i 1\) is just \(N\).

It will help later to introduce the notation

so that the equations are

with

(0-based indexing wins this time.)

```
def linefit(x, y):
m0 = len(x)
m1 = sum(x)
m2 = sum(x**2)
M = np.array([[m0, m1], [m1, m2]])
p = np.array([sum(y), sum(x*y)])
c = solveLinearSystem(M, p)
return c
```

```
N = 10
x = np.linspace(-1, 1, N)
# Emulate a straight line with measurement errors:
# random(N) gives N values uniformly distributed in the range [0,1],
# and so with mean 0.5.
# Thus subtracting 1/2 simulates more symmetric "errors", of mean zero.
yline = 2*x + 3
y = yline + (random(N) - 0.5)
figure(figsize=[12,6])
plot(x, yline, 'g', label="The original line")
plot(x, y, '*', label="Data")
c = linefit(x, y)
print("The coefficients are", c)
xplot = np.linspace(-1, 1, 100)
plot(xplot, evaluatePolynomial(xplot, c), 'r', label="Linear least squares fit")
legend(loc="best");
```

```
The coefficients are [3.1081582 1.82658183]
```

## 3.5.4. Least squares fiting to higher degree polynomials#

The method above extends to finding a polynomial

that gives the best least squares fit to data

in that the coefficients \(c_k\) given the minimum of

Note that when \(N=n+1\), the solution is the interpolating polynomial, with error zero.

The necessary conditions for a minimum are that all \(n+1\) partial derivatives of \(E\) are zero:

This gives

or with the notation \(m_k = \sum_i x_i^k\), \(p_k = \sum_i x_i^k y_i\) introduced above,

That is, the equations are again \(Mc = p\), but now with

(Alternative geometrical derivation)

In the next section Least-squares Fitting to Data: Appendix on The Geometrical Approach, another way to derive this result is given, using geometry and linear algebra instead of calculus.

```
def fitPolynomialLeastSquares(x, y, n):
"""Compute the coeffients c_i of the polynomial of degree n that give the best least squares fit to data (x[i], y[i]).
"""
N = len(x)
m = np.zeros(2*n+1)
for k in range(2*n+1):
m[k] = sum(x**k)
M = np.zeros([n+1,n+1])
for i in range(n+1):
for j in range(n+1):
M[i, j] = m[i+j]
p = np.zeros(n+1)
for k in range(n+1):
p[k] = sum(x**k * y)
c = solveLinearSystem(M, p)
return c
```

```
N = 10
n = 3
xdata = np.linspace(0, np.pi/2, N)
ydata = np.sin(xdata)
figure(figsize=[14,5])
plot(xdata, ydata, 'b*', label="sin(x) data")
xplot = np.linspace(-0.5, np.pi/2 + 0.5, 100)
plot(xplot, np.sin(xplot), 'b', label="sin(x) curve")
c = fitPolynomialLeastSquares(xdata, ydata, n)
print("The coefficients are", c)
plot(xplot, evaluatePolynomial(xplot, c), 'r', label="Cubic least squares fit")
legend(loc="best");
```

```
The coefficients are [-0.00108728 1.02381733 -0.06859178 -0.11328717]
```

The discrepancy between the original function \(\sin(x)\) and this cubic is:

```
figure(figsize=[14,5])
plot(xplot, np.sin(xplot) - evaluatePolynomial(xplot, c))
title("errors fitting at N = 10 points")
grid(True);
```

What if we fit at more points?

```
N = 50
xdata = np.linspace(0, np.pi/2, N)
ydata = np.sin(xdata)
figure(figsize=[14,5])
plot(xdata, ydata, 'b.', label="sin(x) data")
plot(xplot, np.sin(xplot), 'b', label="sin(x) curve")
c = fitPolynomialLeastSquares(xdata, ydata, n)
print("The coefficients are", c)
plot(xplot, evaluatePolynomial(xplot, c), 'r', label="Cubic least squares fit")
legend(loc="best")
grid(True);
```

```
The coefficients are [-0.00200789 1.02647568 -0.06970463 -0.1137182 ]
```

```
figure(figsize=[14,5])
plot(xplot, np.sin(xplot) - evaluatePolynomial(xplot, c))
title("errors fitting at N = 50 points")
grid(True);
```

Not much changes!

This hints at another use of least squares fitting: fitting a simpler curve (like a cubic) to a function (like \(\sin(x)\)), rather than to discrete data.

## 3.5.5. Nonlinear fitting: power-law relationships#

When data \((x_i, y_i)\) is inherently positive, it is often natural to seek an approximate power law relationship

That is, one seeks the power \(p\) and scale factor \(c\) that minimizes error in some sense.

When the magnitudes and the data \(y_i\) vary greatly, it is often appropriate to look at the relative errors

and this can be shown to be very close to looking at the absolute errors of the logarithms

Introducing the new variables \(X_i = \ln(x_i)\), \(Y_i = \ln(Y_i)\) and \(C = \ln(c)\), this becomes the familiar problem of finding a linear approxation of the data \(Y_i\) by \(C + p X_i\).

## 3.5.6. A simulation#

```
cexact = 2.0
pexact = 1.5
x = np.logspace(0.01, 2.0, 10)
xplot = np.logspace(0.01, 2.0) # For graphs later
yexact = cexact * x**pexact
y = yexact * (1.0 + (random(len(yexact))- 0.5)/2)
```

```
figure(figsize=[12,6])
plot(x, yexact, '.', label='exact')
plot(x, y, '*', label='noisy')
legend()
grid(True);
```

```
figure(figsize=[12,6])
loglog(x, yexact, '.', label='exact')
loglog(x, y, '*', label='noisy')
legend()
grid(True)
```

```
X = np.log(x)
Y = np.log(y)
Cp = fitPolynomialLeastSquares(X, Y, 1)
C = np.exp(Cp[0])
p = Cp[1]
print(f"{C=}, {p=}")
```

```
C=2.0413344277069996, p=1.474126908463272
```

```
figure(figsize=[12,6])
plot(x, yexact, '.', label='exact')
plot(x, y, '*', label='noisy')
plot(xplot, C * xplot**p)
legend()
grid(True);
```

```
figure(figsize=[12,6])
loglog(x, yexact, '.', label='exact')
loglog(x, y, '*', label='noisy')
loglog(xplot, C * xplot**p)
legend()
grid(True);
```