<jupyter_start><jupyter_text>Linear Modelling - Least Squares_Linear modelling_ is the process of learning a linear relationship between attributes and responces. Consider the following data set where $x_n$ denotes the nth olympic year and $t_n$ denotes the nth time of the mens 100m.<jupyter_code>%matplotlib inline

x = [1896, 1900, 1904, 1906, 1908, 1912, 1920, 1924, 1928, 1932, 1936, 1948, 1952, 1956, 1960, 1964, 
     1968, 1972, 1976, 1980, 1984, 1988, 1992, 1996, 2000, 2004, 2008]
t = [12.00, 11.00, 11.00, 11.20, 10.80, 10.80, 10.80, 10.60, 10.80, 10.30, 10.30, 10.30, 10.40, 10.50, 
     10.20, 10.00, 9.95, 10.14, 10.06, 10.25, 9.99, 9.92, 9.96, 9.84, 9.87, 9.85, 9.69]

import matplotlib.pyplot as plt
plt.plot(x, t, 'ro')
plt.show()<jupyter_output><empty_output><jupyter_text>Defining the modelWe can clearly see a statistical dependence between year and time (note this is not clausal).We will define are model as a function that maps the input attributes (olympic year) to the target values (winning time). Clearly this function will involves terms other than $x$ called _parameters_, these are the values that we learn from the data set. A function $f$ that maps input attributes $x_0, ..., x_n$ with parameters $w_0, ..., w_k$ is denoted as $f(x_0, ..., x_n; w_0, ... w_k)$.To pick which model we are going to use, we make the assumtion that the data can be modelled linearly, thus the function will be of the form.$$t = f(x; w_0, w_1) = w_0 + w_1 x$$  Loss functionTo find the _best_ model we must define what the _best_ means. In our case we a looking for the line which is as close to all data points as possible. A common way to measure this is the squared diffrence known as the squared loss function:$$\mathcal{L}_n(t_n, f(x_n; w_0, w_1)) = (t_n - f(x_n; w_0, w_1))^2$$The smaller the loss for year $n$, the closer the model at $x_n$ is to $t_n$. We need a low loss for all years, thus we find the average loss:$$\mathcal{L} = \frac{1}{N} \sum^{N}_{n = 1} L_n(t_n, f(x_n; w_0, w_1))$$Thus we will tune $w_0$ and $w_1$ such that the average loss is minimilized, expressed mathmatically:$$\operatorname*{argmin}_{w_0, w_1} = \frac{1}{N} \sum^N_{n=1} L_n(t_n, f(x_n; w_0, w_1))$$We have choosen this loss function because we can find values for $w_0$ and $w_1$ analytically however other loss functions could be used such as the absolute loss:$$\mathcal{L}_n = \left| t_n - f(x_n; w_0, w_1) \right|$$  Finding a solutionTo find the $\operatorname*{argmin}_{w_0, w_1} = \frac{1}{N} \sum^N_{n=1} L_n(t_n, f(x_n; w_0, w_1))$ we can differentiate to find the point in which the gradient is zero (the minima of the function).First we substiture the linear model into the expression for average loss:$$\begin{align}\mathcal{L} &= \frac{1}{N} \sum^{N}_{n = 1} L_n(t_n, f(x_n; w_0, w_1)) \\&= \frac{1}{N} \sum^{N}_{n = 1} (t_n - f(x_n; w_0, w_1))^2 \\&= \frac{1}{N} \sum^{N}_{n = 1} (t_n - (w_0 x_n + w_1))^2 \\&= \frac{1}{N} \sum^{N}_{n = 1} (w_1^2 x_n^2 + 2 w_1 x_n w_0 - 2 w_1 x_n t_n + w_0^2 - 2 w_0 t_n + t_n^2) \\&= \frac{1}{N} \sum^{N}_{n = 1} (w_1^2 x_n^2 + 2 w_1 x_n (w_0 - t_n) + w_0^2 - 2 w_0 t_n + t_n^2) \\\end{align}$$Differentiate with respect to $w_0$$$\begin{align}\frac{1}{N} \sum^{N}_{n = 1} (w_0^2 + 2 w_1 x_n w_0 - 2 w_0 t_n) && \text{Remove terms not including $w_0$} \\w_0^2 + 2 w_1 w_0 \frac{1}{N} \left( \sum^{N}_{n = 1} x_n \right) - 2 w_0 \frac{1}{N} \left( \sum^{N}_{n = 1} t_n \right) && \text{Rearrange} \\\frac{\partial{\mathcal{L}}}{\partial{w_0}} = 2 w_0 + 2 w_1 \frac{1}{N} \left( \sum^{N}_{n = 1} x_n \right) - \frac{2}{N} \left( \sum^{N}_{n = 1} t_n \right) && \text{Differentiate}\end{align}$$Differentiate with respect to $w_1$$$\begin{align}\frac{1}{N} \sum^{N}_{n = 1} (w_1^2 x_n^2 + 2 w_1 x_n w_0 - 2 w_1 x_n t_n)) && \text{Remove terms not including $w_1$} \\w_1^2\frac{1}{N} \left( \sum^{N}_{n = 1} x_n^2 \right) + 2 w_1 w_0 \frac{1}{N} \left( \sum^{N}_{n = 1} x_n \right) - 2 w_1 \frac{1}{N} \left( \sum^{N}_{n = 1} x_n t_n \right) && \text{Rearrange} \\\frac{\partial{\mathcal{L}}}{\partial{w_1}} = 2 w_1 \frac{1}{N} \left( \sum^{N}_{n = 1} x_n^2 \right) + 2 w_0 \frac{1}{N} \left( \sum^{N}_{n = 1} x_n \right) - \frac{2}{N} \left( \sum^{N}_{n = 1} x_n t_n \right) && \text{Differentiate}\end{align}$$Differentiate again, to check if the turning point is a minima$$\begin{align}\frac{\partial^2{\mathcal{L}}}{\partial{w_0}} &= 2 \\\frac{\partial^2{\mathcal{L}}}{\partial{w_1}} &= \frac{2}{N} \sum^N_{n=1}{x^2_n} \\\end{align}$$Set $\frac{\partial{\mathcal{L}}}{\partial{w_0}}$ to $0$$$\begin{align}2 w_0 + 2 w_1 \frac{1}{N} \left( \sum^{N}_{n = 1} x_n \right) - \frac{2}{N} \left( \sum^{N}_{n = 1} t_n \right) &= 0 \\2 w_0 &= \frac{2}{N} \left( \sum^{N}_{n = 1} t_n \right) - 2 w_1 \frac{1}{N} \left( \sum^{N}_{n = 1} x_n \right) \\w_0 &= \frac{1}{N} \left( \sum^{N}_{n = 1} t_n \right) - w_1 \frac{1}{N} \left( \sum^{N}_{n = 1} x_n \right) \\\widehat{w_0} &= \bar{t} - w_1 \bar{x} \\\end{align}$$Subing in the new expression for $w_0$ into $\frac{\partial{\mathcal{L}}}{\partial{w_1}}$ we get:$$\begin{align}\frac{\partial{\mathcal{L}}}{\partial{w_1}} &= 2 w_1 \frac{1}{N} \left( \sum^{N}_{n = 1} x_n^2 \right) + 2 \hat{w_0} \frac{1}{N} \left( \sum^{N}_{n = 1} x_n \right) - \frac{2}{N} \left( \sum^{N}_{n = 1} x_n t_n \right) \\&= 2 w_1 \frac{1}{N} \left( \sum^{N}_{n = 1} x_n^2 \right) + 2 (\bar{t} - w_1 \bar{x}) \frac{1}{N} \left( \sum^{N}_{n = 1} x_n \right) - \frac{2}{N} \left( \sum^{N}_{n = 1} x_n t_n \right) \\&= w_1 \frac{2}{N} \left( \sum^{N}_{n = 1} x_n^2 \right) + \bar{t} \frac{2}{N} \left( \sum^{N}_{n = 1} x_n \right) - w_1 \bar{x} \frac{2}{N} \left( \sum^{N}_{n = 1} x_n \right) - \frac{2}{N} \left( \sum^{N}_{n = 1} x_n t_n \right)\end{align}$$Simplify using $\bar{x} = \frac{1}{N} \sum^{N}_{n = 1} x_n$$$\frac{\partial{\mathcal{L}}}{\partial{w_1}} = w_1 \frac{2}{N} \left( \sum^{N}_{n = 1} x_n^2 \right) + 2 \bar{t} \bar{x} - 2 w_1 (\bar{x})^2 - \frac{2}{N} \left( \sum^{N}_{n = 1} x_n t_n \right)$$Set $\frac{\partial{\mathcal{L}}}{\partial{w_1}}$ to $0$$$\begin{align}w_1 \frac{2}{N} \left( \sum^{N}_{n = 1} x_n^2 \right) + 2 \bar{t} \bar{x} - 2 w_1 (\bar{x})^2 - \frac{2}{N} \left( \sum^{N}_{n = 1} x_n t_n \right) &= 0 \\w_1 \frac{2}{N} \left( \sum^{N}_{n = 1} x_n^2 \right) - 2 w_1 (\bar{x})^2 &= \frac{2}{N} \left( \sum^{N}_{n = 1} x_n t_n \right) - 2 \bar{t} \bar{x} \\w_1 \left[\frac{2}{N} \left( \sum^{N}_{n = 1} x_n^2 \right) - 2 (\bar{x})^2 \right] &= \frac{2}{N} \left( \sum^{N}_{n = 1} x_n t_n \right) - 2 \bar{t} \bar{x} \\w_1 &= \frac{\frac{1}{N} \left( \sum^{N}_{n = 1} x_n t_n \right) - \bar{t} \bar{x}}{\frac{1}{N} \left( \sum^{N}_{n = 1} x_n^2 \right) - (\bar{x})^2} \\\end{align}$$To simplify we will define a few more average quantitys, $\overline{x^2} = \frac{1}{N} \sum^{N}_{n = 1} x_n^2$ and $\overline{xt} = \frac{1}{N} \sum^{N}_{n = 1} x_n t_n$, thus we arrive at:$$\widehat{w_1} = \frac{\overline{xt} - \bar{x} \bar{t}}{\overline{x^2} - (\bar{x})^2}$$Thus all that we need to compute the best parameter values are:$$\begin{align}\widehat{w_0} &= \bar{t} - w_1 \bar{x}\\\widehat{w_1} &= \frac{\overline{xt} - \bar{x} \bar{t}}{\overline{x^2} - (\bar{x})^2}\\\end{align}$$  Least squares fit<jupyter_code>import numpy as np

t_bar = np.average(t)
x_bar = np.average(x)
xt_bar = np.average([x[i] * t[i] for i in range(len(x))])
x2_bar = np.average(list(map(lambda xi: xi**2, x)))

w1 = (xt_bar - x_bar * t_bar) / (x2_bar - x_bar ** 2)
w0 = t_bar - w1 * x_bar

print("w0 = {}, w1 = {}".format(w0, w1))
model_x = np.linspace(min(x), max(x), 100)
model_t = list(map(lambda xi: w0 + w1 * xi, model_x))

plt.plot(x, t, 'ro')
plt.plot(model_x, model_t)
plt.show()<jupyter_output><empty_output><jupyter_text>Extending to vectorsClearly their are lots of attributes which affect mens 100m times. We can extend are model by using vectors, consider the following:$$\mathbf{X} = \begin{bmatrix}1 & x_1\\1 & x_2\\\vdots & \vdots\\1 & x_n\end{bmatrix}\quad\mathbf{t} = \begin{bmatrix}t_1\\t_2\\\vdots\\t_n\end{bmatrix}\quad\mathbf{w} = \begin{bmatrix}w_0\\w_1\end{bmatrix}$$The loss can be written as:$$\mathcal{L} = \frac{1}{N}(\mathbf{t}-\mathbf{Xw})^T (\mathbf{t}-\mathbf{Xw})$$Which we can see by expanding out:$$\begin{align}\mathcal{L} &= \frac{1}{N} \times \left( \begin{bmatrix}t_1 - w_0 - w_1 x_1\\t_2 - w_0 - w_1 x_2\\ \vdots \\ t_n - w_0 - w_1 x_n \end{bmatrix} \right)^T \times \left( \begin{bmatrix}t_1 - w_0 - w_1 x_1\\t_2 - w_0 - w_1 x_2\\ \vdots \\ t_n - w_0 - w_1 x_n \end{bmatrix} \right)\\&= \frac{1}{N} \left[ (t_1 - w_0 + w_1 x_1)^2 + ... + (t_n - w_0 + w_1 x_n)^2 \right] \\&= \frac{1}{N} \sum^{N}_{n = 1} L_n(t_n, f(x_n; w_0, w_1))\end{align}$$As before we differentiate with respect to $\mathbf{w}$ (the parameters in vector form), to do this we can use the following rules:| $f(\mathbf{w})$ | $\frac{\partial f}{\partial \mathbf{w}}$ || :-------------: |:-------------:|| $\mathbf{w}^T \mathbf{x}$ | $\mathbf{x}$ || $\mathbf{x}^T \mathbf{w}$ | $\mathbf{x}$ || $\mathbf{w}^T \mathbf{w}$ | $2 \mathbf{w}$ || $\mathbf{w}^T \mathbf{Cw}$ | $2 \mathbf{Cw}$ |Thus we get the following:$$\frac{\partial \mathcal{L}}{\partial \mathbf{w}} = \frac{2}{N} \mathbf{X}^T \mathbf{Xw} - \frac{2}{N} \mathbf{X}^T \mathbf{t} $$To find $\widehat{\mathbf{w}}$ (the optimum parameters) we set $\frac{\partial \mathcal{L}}{\partial \mathbf{w}}$ and use inverses since division isnt defined for matricies.$$\begin{align}\frac{2}{N} \mathbf{X}^T \mathbf{Xw} - \frac{2}{N} \mathbf{X}^T \mathbf{t} &= 0 \\\frac{2}{N} \mathbf{X}^T \mathbf{Xw} &= \frac{2}{N} \mathbf{X}^T \mathbf{t} \\\mathbf{X}^T \mathbf{Xw} &= \mathbf{X}^T \mathbf{t} \\\mathbf{I} \mathbf{w} &= (\mathbf{X}^T \mathbf{X})^{-1} \mathbf{X}^T \mathbf{t} \\\widehat{\mathbf{w}} &= (\mathbf{X}^T \mathbf{X})^{-1} \mathbf{X}^T \mathbf{t} \\\end{align}$$Now we can define linear models with any number of parameters.  Non-linear responcesSo far the model takes the form $f(x; \mathbf{w}) = w_0 + w_1 x$, however by extending matrix $\mathbf{X}$ we get a polynomial function of any order:$$\mathbf{X} = \begin{bmatrix}x_1^0 & x_1^1 & \cdots & x_1^K \\x_2^0 & x_2^1 & \cdots & x_2^K \\\vdots & \vdots & \ddots & \vdots \\x_N^0 & x_N^1 & \cdots & x_N^K \\\end{bmatrix}$$Once more we can extend this past polynomial function with $K$ functions of $x$, $h_k(X)$:$$\mathbf{X} = \begin{bmatrix}h_1(x_1) & h_2(x_1) & \cdots & h_K(x_1) \\h_1(x_2) & h_2(x_2) & \cdots & h_K(x_2) \\\vdots & \vdots & \ddots & \vdots \\h_1(x_N) & h_2(x_N) & \cdots & h_K(x_N) \\\end{bmatrix}$$<jupyter_code>from scipy.optimize import curve_fit

# Original linear model
def f1(x, w0, w1):
    return w0 + np.multiply(w1, x)

# Second order polinomial
def f2(x, w0, w1, w2):
    return w0 + np.multiply(w1, x) + np.multiply(w2, np.power(x, 2))

# 5-th order polinomial
def f3(x, w0, w1, w2, w3, w4, w5):
    return w0 + np.multiply(w1, x) + np.multiply(w2, np.power(x, 2)) + np.multiply(w3, np.power(x, 3)) \
         + np.multiply(w4, np.power(x, 4)) + np.multiply(w5, np.power(x, 5))

# With sin
def f4(x, w0, w1, w2, a, b):
    return w0 + np.multiply(w1, x) + np.multiply(w2, np.power(x, 2)) + np.multiply(w1, np.sin((x - a)/b))

functions = [
    (f1, 'r', "First order polynomial (linear)"), 
    (f2, 'g', "Second order polynomial"), 
    (f3, 'y', "5th order polynomial"), 
    (f4, 'b', "Function with sin term")
]

for (f, color, title) in functions:
    popt, pcov = curve_fit(f, x, t)
    plt.plot(x, t, 'ro')
    plt.plot(model_x, f(model_x, *popt), color)
    plt.title(title)
    plt.show()<jupyter_output><empty_output><jupyter_text>Generalization and overfittingThe question becomes, what is the best model. We want a model that _generalize_ beyond the training data to make future predictions. Clearly the 5th order polynomial has a lower loss that lower order polynomials, however it doesnt make good future predictions since it pays to much attention to the training data, and is _overfitted_.<jupyter_code>plt.plot(x, t, 'ro')

# training range
popt, pcov = curve_fit(f3, x, t)
plt.plot(model_x, f3(model_x, *popt), 'y')

# prediction range
pred_x = np.linspace(max(x), max(x) + 40, 40)
popt, pcov = curve_fit(f3, x, t)
plt.plot(pred_x, f3(pred_x, *popt), 'y--')

plt.title("5th order polynomial")
plt.show()<jupyter_output><empty_output><jupyter_text>Determining the optimum model complexity such that it generalizes well without overfitting is refered to as the bias-variance tradeoff and is very challenging.  Validation dataA common way to overcome this problem is splitting the dataset into training data and validation data. Several models are trained with the training data, then the loss is computed for each against the validation data. Models that generalize well will have a lower loss against the validation data, where as models that are overfitted will have higher losses.<jupyter_code>training_x = x[:-10]
training_t = t[:-10]

validation_x = x[-10:]
validation_t = t[-10:]

for (f, color, title) in functions:
    popt, pcov = curve_fit(f, training_x, training_t)
    plt.plot(training_x, training_t, 'ro')
    plt.plot(validation_x, validation_t, 'rs')
    plt.plot(model_x, f(model_x, *popt), color)
    plt.title(title)
    plt.show()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.21490166414523448, 'fraction_numerical': 0.049697428139183056, 'mean_word_length': 3.938737392603661, 'pattern_counts': {'":': 0, '<': 21, '<?xml version=': 0, '>': 21, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '5', 'dir': 'jupyter-structured-clean-dedup', 'id': '255460', 'license': 'MIT', 'n_tokens_mistral': 5943, 'n_tokens_neox': 5486, 'n_words': 1546, 'path': 'Linear Modelling - Least Squares.ipynb', 'repo_name': 'bongo227/machine-learning-notes', 'size': '13220'}