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

# generate "noisy" data
n = 20 # number of points
x = np.random.randn(n)
y = np.random.randn(1)*x*x + np.random.randn(1)*x + np.random.randn(1) # second-order polynomial
y = y + 0.1*np.random.randn(n) # add noise

# your least-squares estimation code here (2 lines)
X = np.stack( (np.square(x), x, np.ones((n))), axis=1)
u = inv(np.transpose(X)@X)@np.transpose(X)@y

# plot data and least-squares fit
xp = np.arange(-2,2,0.1)
yp = u[0]*xp*xp + u[1]*xp + u[2]
plt.plot(x,y,'bo') # original data
plt.plot( xp, yp, 'r-' ) # second-order polynomial fit
plt.show()