Fitting functions to data - Fundamentals of Numerical Computation

In Polynomial interpolation we saw how a polynomial can be used to interpolate data—that is, derive a continuous function that evaluates to give a set of prescribed values. But interpolation may not be appropriate in many applications.

Example 3.1.1 (Interpolating temperature data)

Julia

MATLAB

Python

Example 3.1.1

Here are 5-year averages of the worldwide temperature anomaly as compared to the 1951–1980 average (source: NASA).

year = 1955:5:2000
temp = [ -0.0480, -0.0180, -0.0360, -0.0120, -0.0040,
       0.1180, 0.2100, 0.3320, 0.3340, 0.4560 ]
    
scatter(year, temp, label="data",
    xlabel="year", ylabel="anomaly (degrees C)", 
    legend=:bottomright)

A polynomial interpolant can be used to fit the data. Here we build one using a Vandermonde matrix. First, though, we express time as decades since 1950, as it improves the condition number of the matrix.

t = @. (year - 1950) / 10
n = length(t)
V = [ t[i]^j for i in 1:n, j in 0:n-1 ]
c = V \ temp

10-element Vector{Float64}:
  -14.114000001832462
   76.36173810552113
 -165.45597224550528
  191.96056669514388
 -133.27347224319684
   58.015577787494486
  -15.962888891734785
    2.6948063497166928
   -0.2546666667177082
    0.010311111113288083

The coefficients in vector c are used to create a polynomial. Then we create a function that evaluates the polynomial after changing the time variable as we did for the Vandermonde matrix.

using Polynomials, Plots
p = Polynomial(c)
f = yr -> p((yr - 1950) / 10)
plot!(f, 1955, 2000, label="interpolant")

As you can see, the interpolant does represent the data, in a sense. However it’s a crazy-looking curve for the application. Trying too hard to reproduce all the data exactly is known as overfitting.

Example 3.1.1

Here are 5-year averages of the worldwide temperature anomaly as compared to the 1951–1980 average (source: NASA).

year = (1955:5:2000)';
y = [ -0.0480; -0.0180; -0.0360; -0.0120; -0.0040;
    0.1180; 0.2100; 0.3320; 0.3340; 0.4560 ];
scatter(year, y), axis tight
xlabel('year')
ylabel(('anomaly ({\circ}C)'));

t = (year - 1950) / 10;  
n = length(t);
V = ones(n, 1);    % t^0
for j = 1:n-1
    V(:, j+1) = t .* V(:,j);
end
c = V \ y;    % solve for coefficients

We created the Vandermonde matrix columns in increasing-degree order. Thus, the coefficients in c also follow that ordering, which is the opposite of what MATLAB uses. We need to flip the coefficients before using them in polyval.

p = @(year) polyval(c(end:-1:1), (year - 1950) / 10);
hold on
fplot(p, [1955, 2000])    % plot the interpolating function

Example 3.1.1

Here are 5-year averages of the worldwide temperature anomaly as compared to the 1951–1980 average (source: NASA).

year = arange(1955,2005,5)
y = array([ -0.0480, -0.0180, -0.0360, -0.0120, -0.0040,
    0.1180, 0.2100, 0.3320, 0.3340, 0.4560 ])

fig, ax = subplots()
ax.scatter(year, y, color="k", label="data")
xlabel("year")
ylabel("anomaly (degrees C)")
title("World temperature anomaly");

t = (year - 1950) / 10
V = vander(t)
c = linalg.solve(V, y)
print(c)

[ 1.03111111e-02 -2.54666667e-01  2.69480635e+00 -1.59628889e+01
  5.80155778e+01 -1.33273472e+02  1.91960567e+02 -1.65455972e+02
  7.63617381e+01 -1.41140000e+01]

The coefficients in vector c are used to create a polynomial. Then we create a function that evaluates the polynomial after changing the time variable as we did for the Vandermonde matrix.

p = poly1d(c)    # convert to a polynomial
tt = linspace(1955, 2000, 500)
ax.plot(tt, p((tt - 1950) / 10), label="interpolant")
ax.legend();
fig

In many cases we can get better results by relaxing the interpolation requirement. In the polynomial case this allows us to lower the degree of the polynomial, which limits the number of local max and min points. Let $(t_i,y_i)$ for $i=1,\ldots,m$ be the given points. We will represent the data by the polynomial

y \approx f(t) = c_1 + c_2t + \cdots + c_{n-1} t^{n-2} + c_n t^{n-1},

(3.1.1)

with $n<m$ . Just as in (2.1.3), we can express a vector of $f$ -values by a matrix-vector multiplication. In other words, we seek an approximation

\begin{bmatrix} y_1 \\ y_2 \\ y_3 \\ \vdots \\ y_m \end{bmatrix} \approx \begin{bmatrix} f(t_1) \\ f(t_2) \\ f(t_3) \\ \vdots \\ f(t_m) \end{bmatrix} = \begin{bmatrix} 1 & t_1 & \cdots & t_1^{n-1} \\ 1 & t_2 & \cdots & t_2^{n-1} \\ 1 & t_3 & \cdots & t_3^{n-1} \\ \vdots & \vdots & & \vdots \\ 1 & t_m & \cdots & t_m^{n-1} \\ \end{bmatrix} \begin{bmatrix} c_1 \\ c_2 \\ \vdots \\ c_n \end{bmatrix} = \mathbf{V} \mathbf{c}.

(3.1.2)

Note that $\mathbf{V}$ has the same structure as the Vandermonde matrix in (2.1.3) but is $m\times n$ , thus taller than it is wide. It’s impossible in general to satisfy $m$ conditions with $n<m$ variables, and we say the system is overdetermined. Rather than solving the system exactly, we have to find the best approximation. Below we specify precisely what is meant by this, but first we note that Julia uses the same backslash notation to solve the problem in both the square and overdetermined cases.

Example 3.1.2 (Fitting temperature data)

Julia

MATLAB

Python

Example 3.1.2

Here are the 5-year temperature averages again.

year = 1955:5:2000
t = @. (year - 1950) / 10
temp = [ -0.0480, -0.0180, -0.0360, -0.0120, -0.0040,
          0.1180, 0.2100, 0.3320, 0.3340, 0.4560 ]

10-element Vector{Float64}:
 -0.048
 -0.018
 -0.036
 -0.012
 -0.004
  0.118
  0.21
  0.332
  0.334
  0.456

The standard best-fit line results from using a linear polynomial that meets the least-squares criterion.

V = [ t.^0 t ]    # Vandermonde-ish matrix
@show size(V)
c = V \ temp
p = Polynomial(c)

f = yr -> p((yr - 1955) / 10)
scatter(year, temp, label="data",
    xlabel="year", ylabel="anomaly (degrees C)", leg=:bottomright)
plot!(f, 1955, 2000, label="linear fit")

If we use a global cubic polynomial, the points are fit more closely.

V = [ t[i]^j for i in 1:length(t), j in 0:3 ]   
@show size(V);

size(V) = (10, 4)

Now we solve the new least-squares problem to redefine the fitting polynomial.

p = Polynomial( V \ temp )
plot!(f, 1955, 2000, label="cubic fit")

If we were to continue increasing the degree of the polynomial, the residual at the data points would get smaller, but overfitting would increase.

Example 3.1.2

Here are the 5-year temperature averages again.

year = (1955:5:2000)';
y = [ -0.0480; -0.0180; -0.0360; -0.0120; -0.0040;
    0.1180; 0.2100; 0.3320; 0.3340; 0.4560 ];

The standard best-fit line results from using a linear polynomial that meets the least-squares criterion.

t = (year - 1955) / 10;    % better matrix conditioning later
V = [ t.^0 t ];            % Vandermonde-ish matrix
size(V)

c = V \ y;
f = @(year) polyval(c(end:-1:1), (year - 1955) / 10);

clf
scatter(year, y), axis tight
xlabel('year'), ylabel('anomaly ({\circ}C)')
hold on
fplot(f, [1955, 2000])

If we use a global cubic polynomial, the points are fit more closely.

V = [t.^0, t.^1, t.^2, t.^3];    % Vandermonde-ish matrix  
size(V)

Now we solve the new least-squares problem to redefine the fitting polynomial.

c = V \ y;
f = @(year) polyval(c(end:-1:1), (year - 1955) / 10);
fplot(f, [1955, 2000]) 
legend('data', 'linear', 'cubic', 'Location', 'northwest');

If we were to continue increasing the degree of the polynomial, the residual at the data points would get smaller, but overfitting would increase.

Example 3.1.2

Here are the 5-year temperature averages again.

year = arange(1955, 2005, 5)
y = array([-0.0480, -0.0180, -0.0360, -0.0120, -0.0040,
    0.1180, 0.2100, 0.3320, 0.3340, 0.4560])
t = (year - 1950) / 10

The standard best-fit line results from using a linear polynomial that meets the least-squares criterion.

V = array([ [t[i], 1] for i in range(t.size) ])    # Vandermonde-ish matrix
print(V.shape)

(10, 2)

from numpy.linalg import lstsq
c, res, rank, sv = lstsq(V, y)
p = poly1d(c)
f = lambda year: p((year - 1950) / 10)


```{code-cell}
fig, ax = subplots()
ax.scatter(year, y, color="k", label="data")
yr = linspace(1955, 2000, 500)
ax.plot(yr, f(yr), label="linear fit")

xlabel("year")
ylabel("anomaly (degrees C)")
title("World temperature anomaly");
ax.legend();

If we use a global cubic polynomial, the points are fit more closely.

V = array([ [t[i]**3,t[i]**2,t[i],1] for i in range(t.size) ])    # Vandermonde-ish matrix
print(V.shape)

(10, 4)

Now we solve the new least-squares problem to redefine the fitting polynomial.

c, res, rank, sv = lstsq(V, y, rcond=None)
yr = linspace(1955, 2000, 500)
ax.plot(yr, f(yr), label="cubic fit")
fig

If we were to continue increasing the degree of the polynomial, the residual at the data points would get smaller, but overfitting would increase.

3.1.1The least-squares formulation¶

In the most general terms, our fitting functions take the form

f(t) = c_1 f_1(t) + \cdots + c_n f_n(t)

(3.1.3)

where $f_1,\ldots,f_n$ are all known functions with no undetermined parameters. This leaves only $c_1,\ldots,c_n$ to be determined. The essential feature of a linear least-squares problem is that the fit depends only linearly on the unknown parameters. For instance, a function of the form $f(t)=c_1 + c_2 e^{c_3 t}$ is not of this type.

At each observation $(t_i,y_i)$ , we define a residual, $y_i - f(t_i)$ . A sensible formulation of the fitting criterion is to minimize

R(c_1,\ldots,c_n) = \sum_{i=1}^m\, [ y_i - f(t_i) ]^2,

(3.1.4)

over all possible choices of parameters $c_1,\ldots,c_n$ . We can apply linear algebra to write the problem in the form $R=\mathbf{r}^T \mathbf{r}$ , where

\mathbf{r} = \begin{bmatrix} y_1 \\ y_2 \\ \vdots \\y_{m-1} \\ y_m \end{bmatrix} - \begin{bmatrix} f_1(t_1) & f_2(t_1) & \cdots & f_n(t_1) \\[1mm] f_1(t_2) & f_2(t_2) & \cdots & f_n(t_2) \\[1mm] & \vdots \\ f_1(t_{m-1}) & f_2(t_{m-1}) & \cdots & f_n(t_{m-1}) \\[1mm] f_1(t_m) & f_2(t_m) & \cdots & f_n(t_m) \\[1mm] \end{bmatrix} \begin{bmatrix} c_1 \\ c_2 \\ \vdots \\ c_n \end{bmatrix}.

(3.1.5)

Recalling that $\mathbf{r}^T\mathbf{r}=\| \mathbf{r} \|_2^2$ , and renaming the variables to standardize the statement, we arrive at the general linear least-squares problem.

While we could choose to minimize in any vector norm, the 2-norm is the most common and convenient choice. For the rest of this chapter we exclusively use the 2-norm. In the edge case $m=n$ for a nonsingular $\mathbf{A}$ , the definitions of the linear least-squares and linear systems problems coincide: the solution of $\mathbf{A}\mathbf{x}=\mathbf{b}$ implies $\mathbf{r}=\boldsymbol{0}$ , which is a global minimum of $\| \mathbf{r} \|_2^2 \ge 0$ .

3.1.2Change of variables¶

The most familiar and common case of a polynomial least-squares fit is the straight line, $f(t) = c_1 + c_2 t$ . Certain other fit functions can be transformed into this situation. For example, suppose we want to fit data using $g(t)= a_1 e^{a_2 t}$ . Then

\log y \approx \log g(t) = (\log a_1) + a_2 t = c_1 + c_2 t.

(3.1.7)

While the fit of the $y_i$ to $g(t)$ is nonlinearly dependent on fitting parameters, the fit of $\log y$ to a straight line is a linear problem. Similarly, the power-law relationship $y\approx f(t)=a_1 t^{a_2}$ is equivalent to

\log y \approx (\log a_1) + a_2 (\log t).

(3.1.8)

Thus, the variable $z=\log y$ can be fit linearly in terms of the variable $s=\log t$ . In practice these two cases—exponential fit and power law—are easily detected by using log-linear or log-log plots, respectively.

Example 3.1.3 (Fitting a power law)

Finding numerical approximations to π has fascinated people for millennia. One famous formula is

\displaystyle \frac{\pi^2}{6} = 1 + \frac{1}{2^2} + \frac{1}{3^2} + \cdots.

(3.1.9)

Say $s_k$ is the sum of the first $k$ terms of the series above, and $p_k = \sqrt{6s_k}$ . Here is a fancy way to compute these sequences in a compact code.

Julia

MATLAB

Python

Example 3.1.3

a = [1/k^2 for k=1:100] 
s = cumsum(a)        # cumulative summation
p = @. sqrt(6*s)

scatter(1:100, p;
    title="Sequence convergence",
    xlabel=L"k",  ylabel=L"p_k")

This graph suggests that maybe $p_k\to \pi$ , but it’s far from clear how close the sequence gets. It’s more informative to plot the sequence of errors, $\epsilon_k= |\pi-p_k|$ . By plotting the error sequence on a log-log scale, we can see a nearly linear relationship.

ϵ = @. abs(π - p)    # error sequence
scatter(1:100, ϵ;
    title="Convergence of errors",
    xaxis=(:log10,L"k"),  yaxis=(:log10,"error"))

The straight line on the log-log scale suggests a power-law relationship where $\epsilon_k\approx a k^b$ , or $\log \epsilon_k \approx b (\log k) + \log a$ .

k = 1:100
V = [ k.^0 log.(k) ]     # fitting matrix
c = V \ log.(ϵ)          # coefficients of linear fit

2-element Vector{Float64}:
 -0.18237524972829994
 -0.9674103233127929

In terms of the parameters $a$ and $b$ used above, we have

a, b = exp(c[1]), c[2];
@show b;

b = -0.9674103233127929

It’s tempting to conjecture that the slope $b\to -1$ asymptotically. Here is how the numerical fit compares to the original convergence curve.

plot!(k, a * k.^b, l=:dash, label="power-law fit")

Example 3.1.3

k = (1:100)';
a = 1./k.^2;      % sequence
s = cumsum(a);    % cumulative summation
p = sqrt(6*s);
clf
plot(k, p, 'o-')
xlabel('k'), ylabel('p_k')
title('Sequence converging to \pi')

ep = abs(pi - p);    % error sequence
loglog(k, ep, 'o')
title('Convergence')
xlabel('k'), ylabel('|p_k - \pi|'), axis tight

The straight line on the log-log scale suggests a power-law relationship where $\epsilon_k\approx a k^b$ , or $\log \epsilon_k \approx b (\log k) + \log a$ .

V = [ k.^0, log(k) ];    % fitting matrix
c = V \ log(ep)          % coefficients of linear fit

In terms of the parameters $a$ and $b$ used above, we have

a = exp(c(1)),  b = c(2)

It’s tempting to conjecture that the slope $b\to -1$ asymptotically. Here is how the numerical fit compares to the original convergence curve.

hold on
loglog(k, a * k.^b)
legend('sequence', 'power-law fit');

Example 3.1.3

a = array([1 / (k+1)**2 for k in range(100)])
s = cumsum(a)        # cumulative summation
p = sqrt(6*s)

plot(range(100), p, "o")
xlabel("$k$") 
ylabel("$p_k$") 
title("Sequence convergence");

ep = abs(pi - p)    # error sequence
loglog(range(100), ep, "o")
xlabel("$k$") 
ylabel("error") 
title("Sequence convergence");

The straight line on the log-log scale suggests a power-law relationship where $\epsilon_k\approx a k^b$ , or $\log \epsilon_k \approx b (\log k) + \log a$ .

V = array([ [1, log(k+1)] for k in range(100) ])     # fitting matrix
c = lstsq(V, log(ep), rcond=None)[0]           # coefficients of linear fit
print(c)

[-0.18237525 -0.96741032]

In terms of the parameters $a$ and $b$ used above, we have

a, b = exp(c[0]), c[1]
print(f"b: {b:.3f}")

b: -0.967

It’s tempting to conjecture that the slope $b\to -1$ asymptotically. Here is how the numerical fit compares to the original convergence curve.

loglog(range(100), ep, "o", label="sequence")
k = arange(1,100)
plot(k, a*k**b, "--", label="power fit")
xlabel("$k$");  ylabel("error"); 
legend(); title("Sequence convergence");

3.1.3Exercises¶

Exercise 3.1.2

⌨ Here are counts of the U.S. population in millions from the census performed every ten years, beginning with 1790 and ending with 2010.

3.929, 5.308, 7.240, 9.638, 12.87, 17.07, 23.19, 31.44, 39.82, 50.19, 62.95, 76.21,
92.22, 106.0, 122.8, 132.2, 150.7, 179.3, 203.3, 226.5, 248.7, 281.4, 308.7

(a) Find a best-fitting cubic polynomial for the data. Plot the data as points superimposed on a (smooth) graph of the cubic over the full range of time. Label the axes. What does the fit predict for the population in the years 2000, 2010, and 2020?

(b) Look up the actual U.S. population in 2000, 2010, and 2020 and compare to the predictions of part (a).

Exercise 3.1.3

⌨ The following are weekly box office earnings (in dollars) in the U.S. for the 2012 film The Hunger Games. (Source: boxofficemojo.com.)

189_932_838, 79_406_327, 46_230_374, 26_830_921, 18_804_290,
13_822_248, 7_474_688, 6_129_424, 4_377_675, 3_764_963, 2_426_574,
1_713_298, 1_426_102, 1_031_985, 694_947, 518_242, 460_578, 317_909

(Note that Julia lets you use _ where you would normally put a comma in a long number.) Fit these values to a function of the form $y(t)\approx a e^{b t}$ . Plot the data together with the fit using standard linear scales on the axes, and then plot them again using a log scale on the vertical axis.

Exercise 3.1.4

⌨ In this problem you are trying to find an approximation to the periodic function $g(t)=e^{\sin(t-1)}$ over one period, $0 < t \le 2\pi$ . As data, define

t_i = \frac{2\pi i}{60}, \; y_i = g(t_i), \quad i=1,\ldots,60.

(3.1.10)

(a) Find the coefficients of the least-squares fit

y(t) \approx c_1 + c_2t + \cdots + c_7 t^6.

(3.1.11)

Superimpose a plot of the data values as points with a curve showing the fit.

(b) Find the coefficients of the least-squares fit

y \approx d_1 + d_2\cos(t) + d_3\sin(t) + d_4\cos(2t) + d_5\sin(2t).

(3.1.12)

Unlike part (a), this fitting function is itself periodic. Superimpose a plot of the data values as points with a curve showing the fit.

Exercise 3.1.5

⌨ Define the following data in Julia.

t = 0:.5:10
y = tanh.(t)

(a) Fit the data to a cubic polynomial. Plot the data together with the polynomial fit over the interval $0 \le t \le 10$ .

(b) Fit the data to the function $c_1 + c_2z + c_3z^2 + c_4z^3$ , where $z=t^2/(1+t^2)$ . Plot the data together with the fit. What feature of $z$ makes this fit much better than the original cubic?

Exercise 3.1.7

⌨ Kepler found that the orbital period τ of a planet depends on its mean distance $R$ from the sun according to $\tau=c R^{\alpha}$ for a simple rational number α. Perform a linear least-squares fit from the following table in order to determine the most likely simple rational value of α.

Planet	Distance from sun in Mkm	Orbital period in days
Mercury	57.59	87.99
Venus	108.11	224.7
Earth	149.57	365.26
Mars	227.84	686.98
Jupiter	778.14	4332.4
Saturn	1427	10759
Uranus	2870.3	30684
Neptune	4499.9	60188

Preface

3. Overdetermined linear systems

Preface

The normal equations