Noisy regression

If we want to minimize $F(a,b,c) \; = \; \sum_{j=1}^n\big(y_j - ax_j^2 - bx_j - c)^2$ then we require $\frac{\partial F}{\partial a} \; = \; \frac{\partial F}{\partial b} \; =\; \frac{\partial F}{\partial c} \; = \; 0$ and hence $\left(\begin{array}{ccc} E[x^4] & E[x^3] & E[x^2] \\ E[x^3] & E[x^2] & E[x] \\ E[x^2] & E[x] & 1 \end{array}\right)\left(\begin{array}{c} a \\ b \\ c \end{array}\right) \; = \; \left(\begin{array}{c} E[x^2y] \\ E[xy] \\ E[y] \end{array} \right)$ where $E[x^m] \; = \; \frac{1}{n}\sum_{j=1}^nx_j^m \hspace{2cm} 1 \le m \le 4$ is the average value of $x^m$ for the data, and $E[x^my] \; = \; \frac{1}{n}\sum_{j=1}^m x_j^m y_j \hspace{2cm} 0 \le m \le 2$ is the average value of $x^my$ for the data (note that $n$ is the number of data points - in this case, $n=10000$ ).

Calculating these expected values can be done easily in Excel, and simple Gaussian elimination solves the equation, giving $a = 1.003089322$ , $b = 3.132961855$ , $c = 9.88548454$ , obtaining $a+b+c = \boxed{14.02153572}$ .

It is well known that if we have 3 points we can define a parabola with the following matrix

$\left[\begin{matrix} y_0 \\ y_1 \\ y_2 \end{matrix}\right] = \left[\begin{matrix} x_0^2 && x_0 && 1 \\ x_1^2 && x_1 && 1 \\ x_2^2 && x_2 && 1 \\ \end{matrix}\right] \cdot \left[\begin{matrix} a \\ b \\ c \end{matrix}\right]$

or more commonly written as $A\vec x = \vec b$ where

$A = \left[\begin{matrix} x_0^2 && x_0 && 1 \\ x_1^2 && x_1 && 1 \\ x_2^2 && x_2 && 1 \\ \end{matrix}\right] \text{, } \vec x = \left[\begin{matrix} a \\ b \\ c \end{matrix}\right] \text{, } \vec b = \left[\begin{matrix} y_0 \\ y_1 \\ y_2 \end{matrix}\right]$

and we can calculate $\vec x$ by computing $\vec x = A^{-1} \vec b$ .

Well, we do have three points but we have even more than 3 points but, because the data is noisy, is not likely that any 2 triplets will fit to the same parabola, and if we are to include all points, our matrix $A$ is not square, so we can't invert it. It turns out that there is a notion of a pseudo-inverse that can be applied to $m \times n$ matrix and when applied it has the nice property of minimizing the squared error and it is defined by $(A^T A)^{-1} A^T$ so we can compute $\vec x$ by computing $\vec x = (A^T A)^{-1} A^T \vec b$ . Here is my python implementation of it with the plot of the data and the regression line:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19

from numpy.linalg import *
import matplotlib.pyplot as plt
import numpy as np

x, y = np.loadtxt("data.txt").T

A = np.array([x**2, x, x**0]).T
pseudo_inv = np.matmul(inv(np.matmul(A.T, A)), A.T)
regres = np.matmul(pseudo_inv, y)

f = np.poly1d(regres)
_y = f(x)

print(regres)
print(sum(regres))

plt.plot(x, y, 'ro', markersize=1)
plt.plot(x,_y, 'b-', markersize=1)
plt.show()