import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
sns.set(font_scale=1.0)
sns.set_style("white")
sns.set_style("ticks")
palette = sns.color_palette('deep')
df = pd.read_csv('notebooks/data/data100m.csv')
df.columns=['Year', 'Time']
N = df.shape[0]
max_year, min_year = df['Year'].values.max() , df['Year'].values.min()
x = (df['Year'].values.reshape(N,1) - min_year)/(max_year - min_year)
t = df['Time'].values.reshape(N,1)
X_func = lambda u : np.hstack([np.ones((u.shape[0],1)), u])
X = X_func(x)
w_hat = np.linalg.inv(X.T @ X) @ X.T @ t
loss_func = 1./N * (t - X @ w_hat).T @ (t - X @ w_hat)
xgrid = np.linspace(0, 1, 100).reshape(100,1)
Xg = X_func(xgrid)
time_grid = Xg @ w_hat
fig = plt.figure(figsize=(6,4))
plt.plot(df['Year'].values, df['Time'].values, 'o', color='crimson', label='Data')
plt.plot(xgrid*(max_year - min_year) + min_year, time_grid, '-', color='dodgerblue', label='Model')
plt.xlabel('Year')
plt.ylabel('Time (seconds)')
loss_title = r'Loss function, $\mathcal{L}=$'+str(np.around(float(loss_func), 5))+'; \t norm of $\hat{\mathbf{w}}$='+str(np.around(float(np.linalg.norm(w_hat,2)), 3))
plt.title(loss_title)
plt.legend()
plt.savefig('olympics_0.png', dpi=150, bbox_inches='tight')
plt.close()Lecture 6
Linear modelling
The three model levels
In this lecture, we will explore three distinct but related flavors of modelling.
Linear least squares model
Gaussian noise model (introducing the idea of likelihood)
Full Bayesian treatment (next time!)
Much of the exposition shown here is based on the first three chapters of Rogers and Girolami’s, A First Course in Machine Learning.
Linear least squares
- Consider the data shown in the plot below. It shows the winning times for the men’s 100 meter race at the Summer Olympics for many years.
Our goal will be to fit a model to this data. To begin, we will consider a linear model, i.e., t=f(x;w0,w1)=w0+w1x where x is the year and t is the winning time.
w0 and w1 are unknown model parameters that we need to ascertain.
Good sense would suggest that the best line passes as closely as possible through all the data points on the left.
Linear least squares
Defining a good model
One common strategy for defining this is based on the squared distance between the truth and the model. Thus, for a given year, ti, this is written as: Li=(ti−f(xi;w0,w1))2.
However, as we want a model that fits well across all the data, we may consider the average across the entire data set, i.e., all N data points. This is given by: L=N1i=1∑NLi=N1i=1∑N(ti−f(xi;w0,w1))2
Note that this loss function is always positive, and the lower it is the better! Finding optimal values for w0,w1 can be expressed as w0,w1argminN1i=1∑N(ti−f(xi;w0,w1))2
Note that other loss functions can be considered. A common example is the absolute loss, i.e., Li=∣ti−f(xi;w0,w1)∣
Linear least squares
Matrix-vector notation
It will be very useful to work with vectors and matrices. For convenience, we define: X=⎣⎡11⋮1x1x2⋮xN⎦⎤=⎣⎡x1Tx2T⋮xNT⎦⎤,t=⎣⎡t1t2⋮tN⎦⎤,w=[w0w1]
The loss function from the prior slide is equivalent to writing L=N1(t−Xw)T(t−Xw).
This can be expanded to yield L=N1(tT−(Xw)T)(t−Xw)=N1[tTt−2tTXw+wTXTXw]
Linear least squares
Minimizing the loss
As our objective is to minimize the loss, the obvious idea is to find out for which w, the derivative of the loss function, ∂L/∂w, goes to zero.
Note that in practice, we refer to these points as turning points as they may equally correspond to maxima, minima, or saddle points. A positive second derivative is a sure sign of a minima.
Prior to working out the derivatives, it will be useful to take note of the following identities on the left below.
| g(v) | ∂g/∂v |
|---|---|
| vTx | x |
| xTv | x |
| vTv | 2v |
| vTCv | 2Cv |
The derivative of the loss function is given by ∂w∂L=−N2XTt+N2XTXw
Setting the derivative to zero, we have XTXw=XTt⇒w^=(XTX)−1XTt where w^ represents the value of w that minimizes the loss.
Linear least squares
- For this result we set X=⎣⎡11⋮1x1x2⋮xN⎦⎤ and solve for w^=(XTX)−1XTt.
- Once these weights are obtained, we can extrapolate (blue line) over the years.
- Note the graph title shows the loss function value and the L2 norm, ∥w^∥2.
Linear least squares
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
sns.set(font_scale=1.0)
sns.set_style("white")
sns.set_style("ticks")
palette = sns.color_palette('deep')
df = pd.read_csv('notebooks/data/data100m.csv')
df.columns=['Year', 'Time']
N = df.shape[0]
max_year, min_year = df['Year'].values.max() , df['Year'].values.min()
x = (df['Year'].values.reshape(N,1) - min_year)/(max_year - min_year)
t = df['Time'].values.reshape(N,1)
X_func = lambda u : np.hstack([np.ones((u.shape[0],1)), u , u**2, u**3])
X = X_func(x)
w_hat = np.linalg.inv(X.T @ X) @ X.T @ t
loss_func = 1./N * (t - X @ w_hat).T @ (t - X @ w_hat)
xgrid = np.linspace(0, 1, 100).reshape(100,1)
Xg = X_func(xgrid)
time_grid = Xg @ w_hat
fig = plt.figure(figsize=(6,4))
plt.plot(df['Year'].values, df['Time'].values, 'o', color='crimson', label='Data')
plt.plot(xgrid*(max_year - min_year) + min_year, time_grid, '-', color='dodgerblue', label='Model')
plt.xlabel('Year')
plt.ylabel('Time (seconds)')
loss_title = r'Loss function, $\mathcal{L}=$'+str(np.around(float(loss_func), 5))+'; \t norm of $\hat{\mathbf{w}}$='+str(np.around(float(np.linalg.norm(w_hat,2)), 3))
plt.title(loss_title)
plt.legend()
plt.savefig('olympics_3.png', dpi=150, bbox_inches='tight')
plt.close()- For this result we set X=⎣⎡11⋮1x1x2⋮xNx12x22⋮xN2x13x23⋮xN3⎦⎤ and solve for w^=(XTX)−1XTt.
- Once these weights are obtained, we can extrapolate (blue line) over the years.
- Note the graph title shows the loss function value and the L2 norm, ∥w^∥2.
Linear least squares
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
sns.set(font_scale=1.0)
sns.set_style("white")
sns.set_style("ticks")
palette = sns.color_palette('deep')
df = pd.read_csv('notebooks/data/data100m.csv')
df.columns=['Year', 'Time']
N = df.shape[0]
max_year, min_year = df['Year'].values.max() , df['Year'].values.min()
x = (df['Year'].values.reshape(N,1) - min_year)/(max_year - min_year)
t = df['Time'].values.reshape(N,1)
X_func = lambda u : np.hstack([np.ones((u.shape[0],1)), u , u**2, u**3])
X = X_func(x)
w_hat = np.linalg.inv(X.T @ X) @ X.T @ t
loss_func = 1./N * (t - X @ w_hat).T @ (t - X @ w_hat)
xgrid = np.linspace(0, 1, 100).reshape(100,1)
Xg = X_func(xgrid)
time_grid = Xg @ w_hat
fig = plt.figure(figsize=(6,4))
plt.plot(df['Year'].values, df['Time'].values, 'o', color='crimson', label='Data')
plt.plot(xgrid*(max_year - min_year) + min_year, time_grid, '-', color='dodgerblue', label='Model')
plt.xlabel('Year')
plt.ylabel('Time (seconds)')
loss_title = r'Loss function, $\mathcal{L}=$'+str(np.around(float(loss_func), 5))+'; \t norm of $\hat{\mathbf{w}}$='+str(np.around(float(np.linalg.norm(w_hat,2)), 3))
plt.title(loss_title)
plt.legend()
plt.savefig('olympics_8.png', dpi=150, bbox_inches='tight')
plt.close()- For this result we set X=⎣⎡11⋮1x1x2⋮xNx12x22⋮xN2x13x23⋮xN3…………x18x28⋮xN8⎦⎤ and solve for w^=(XTX)−1XTt.
- Once these weights are obtained, we can extrapolate (blue line) over the years.
- Note the graph title shows the loss function value and the L2 norm, ∥w^∥2.
Linear least squares
With regularization
There is clearly a trade-off between the:
- Complexity of the model in terms of the number of weights, and
- the value of the loss function.
There is also the risk of over-fitting to the data. For instance, if we had only 9 data points, then the last model would have interpolated each point, at the risk of not being generalizable.
As we do not want our model to be too complex, there are two relatively simple recipes:
- Split the data into test and train (your homework!)
- Add a regularization term, i.e., L=L+λwTw where λ is a constant.
Maximum likelihood
The linear model from before is unable to capture each data point, and there are errors between the true data and the model.
Now we will consider a paradigm where these errors are explicitly modelled.
We consider a model of the form tj=f(xn;w)+ϵnϵn∼N(0,σ2),wherej∈[1,N](1)
Recall, we had previously learnt that adding a constant to a Gaussian random variable alters its mean. Thus, the random variable tj has a probability density function p(tj∣xj,w,σ2)=N(wTxj,σ2)
Carefully note the conditioning: the probability density function for tj depends on particular values of xj and w.
Maximum likelihood
Defining the likelihood
If we evaluate the linear model from before, and assume that in Equation 1 σ2=0.05, we would find that p(tj∣xj=[1,1980]T,w=[10.964,−1.31]T,σ2=0.05)=N(10.03,0.05) This quantity is known as the likelihood of the n-th data point.
Note that for a continuous random variable t, p(t) cannot be interpreted as a probability.
The height of the curve to the left tells us how likely it is that we observe a particular t for x=1980.
The most likely is B, followed by C and then A. Note the actual winning time is tn=10.25.
While we obviously cannot change the actual winning time, we can change w and σ2 to move the density to make it as high as possible at t=10.25.
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
from scipy.stats import multivariate_normal
# Get the data
df = pd.read_csv('notebooks/data/data100m.csv')
df.columns=['Year', 'Time']
N = df.shape[0]
max_year, min_year = df['Year'].values.max() , df['Year'].values.min()
x = (df['Year'].values.reshape(N,1) - min_year)/(max_year - min_year)
t = df['Time'].values.reshape(N,1)
X_func = lambda u : np.hstack([np.ones((u.shape[0],1)), u])
X = X_func(x)
w_hat = np.linalg.inv(X.T @ X) @ X.T @ t
# Specific year!
year_j = 1980
X_j = X_func(np.array( [ (year_j - min_year) / (max_year - min_year) ] ).reshape(1,1) )
time_j = float(X_j @ w_hat)
T_1980 = multivariate_normal(time_j, 0.05)
ti = np.linspace(9, 11, 100)
pt_x = T_1980.pdf(ti)
fig = plt.figure(figsize=(7,3))
plt.plot(ti, pt_x, '-', color='orangered', lw=3, label='From linear model')
plt.axvline(9.53, linestyle='-.', color='dodgerblue', label='A')
plt.axvline(10.08, linestyle='-', color='green', label='B')
plt.axvline(10.40, linestyle='--', color='navy', label='C')
plt.xlabel('Time (seconds)')
plt.ylabel(r'$p \left( t | x \right)$')
plt.title(r'For $x=1980$')
plt.legend()
plt.close()Maximum likelihood
Defining the likelihood
This idea of finding parameters that can maximize the likelihood is very important in machine learning.
However, in general, we are seldom interested in the likelihood of an isolated data point – we are interested in the likelihood across all the data.
This leads to the conditional distribution across all N data points p(t1,…,tN∣x1,…,xN,w,σ2)
If we assume the noise at each data point is independent, this conditional density can be factorized into N separate terms L=p(t∣X,w,σ2)=j=1∏Np(tj∣xj,w,σ2)=j=1∏NN(wTxn).(2)
Note that the tj values are not completely independent—times have clearly decreased over the years! They are conditionally independent. For a given value of w the tj are independent; otherwise they are not.
We will now maximize the likelihood (see Equation 2 ).
Maximum likelihood
Maximizing the logarithm of the likelihood
Plugging in the definition of a Gaussian probability density function into Equation 2 we arrive at L=j=1∏N2πσ21exp(−2σ21(tj−f(xj;w))2)
Taking the logarithm on both sides and simplifying: log(L)=j=1∑Nlog(2πσ21exp(−2σ21(tj−f(xj;w))2)) =j=1∑N(−21log(2π)−log(σ)−2σ21(tj−f(xj;w))2) =−2Nlog(2π)−Nlog(σ)−2σ21j=1∑N(tj−f(xj;w))2
Maximum likelihood
Maximizing the logarithm of the likelihood
Just as we did earlier, with the least squares solution, we can set the derivative of the logarithm of the loss function to be zero. ∂w∂log(L)=σ21j=1∑Nxj(tj−xjTw)=σ21j=1∑Nxjtj−xjxjTw≡0
Just as we did before, we can use matrix vector notation to write this out as
∂w∂log(L)=σ21(XTt−XTXw)=0
Solving this expression leads to
XTt−XTXw=0⇒w^=(XTX)−1XTt.Thus, the maximum likelihood solution for w is exactly the solution for the least squares problem! Minimizing the squared loss is equivalent to the maximum likelihood solution if the noise is assumed Gaussian.
Maximum likelihood
Maximizing the logarithm of the likelihood
What remains now is to compute the maximum likelihood estimate of the noise, σ. Assuming that w^=w we can write ∂σ∂log(L)=−σN+σ31j=1∑N(tj−xTw^)2≡0.
Rearranging, this yields σ2^=1/N∑j=1N(tj−xTw^).
This expression states that the variance is the averaged squared error, which intuitively makes sense. Re-writing this using matrix notation, we have σ2^=N1(t−Xw^)T(t−Xw^)=N1(tTt−2tTXw^+w^TXTXw^)
Plugging in w^=(XTX)−1XTt, we arrive at σ2^=N1(tTt−tTX(XTX)−1XTt)
Maximum likelihood
Visualizing the Gaussian noise model
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
sns.set(font_scale=1.0)
sns.set_style("white")
sns.set_style("ticks")
palette = sns.color_palette('deep')
df = pd.read_csv('notebooks/data/data100m.csv')
df.columns=['Year', 'Time']
N = df.shape[0]
max_year, min_year = df['Year'].values.max() , df['Year'].values.min()
x = (df['Year'].values.reshape(N,1) - min_year)/(max_year - min_year)
t = df['Time'].values.reshape(N,1)
X_func = lambda u : np.hstack([np.ones((u.shape[0],1)), u ])
X = X_func(x)
w_hat = np.linalg.inv(X.T @ X) @ X.T @ t
loss_func = 1./N * (t - X @ w_hat).T @ (t - X @ w_hat)
xgrid = np.linspace(0, 1, 100).reshape(100,1)
Xg = X_func(xgrid)
time_grid = Xg @ w_hat
xi = xgrid*(max_year - min_year) + min_year
xi = xi.flatten()
sigma_hat_squared = float( (1. / N) * (t.T @ t - t.T @ X @ w_hat) )
sigma_hat = np.sqrt(sigma_hat_squared)
yi = xi* 0 + sigma_hat
loss_func = -N/2 * np.log(np.pi * 2) - N * np.log(sigma_hat) - \
1./(2 * sigma_hat_squared) * np.sum((X @ w_hat - t)**2)
fig = plt.figure(figsize=(6,4))
a, = plt.plot(df['Year'].values, df['Time'].values, 'o', color='crimson', label='Data')
plt.plot(xi, time_grid, '-', color='r')
c = plt.fill_between(xi, time_grid.flatten()-yi, time_grid.flatten()+yi, color='red', alpha=0.2, label='Model')
plt.xlabel('Year')
plt.ylabel('Time (seconds)')
loss_title = r'Logarithm of loss function, $log \left( \mathcal{L} \right)=$'+str(np.around(float(loss_func), 5))
plt.title(loss_title)
plt.legend([a,c], ['Data', 'Model'])
plt.savefig('olympics_last.png', dpi=150, bbox_inches='tight')
plt.close()The graph on the left is the final model.
AE8803 | Gaussian Processes for Machine Learning
Lecture 6 Linear modelling