Ok, linear models might not seem that exciting. People think at first that they're just straight lines...this is so wrong, I love them! They provide the basics you need to understand advanced topics like **deep neural networks** or **Gaussian Processes**. To be honest, they underpin everything in ML.

In this final section, we will work through the key concepts of fitting both Linear and Generalised Linear Models (GLMs). We will talk about feature selection - which includes **Predictive Power Score** (or PPS) and **Principle Component Analysis** (PCA). We'll also cover the important topic of over-fitting and how we address this using regularisation methods.

We'll be playing in ninja mode; so the gloves come off in this course to wrestle with the ultimate machine learning model: _the linear model_. Plenty of reasons for taking a deep dive here:

1. This is an intuitive tangible model on which to build fundamental concepts.
2. We can adapt the linear model to solve more complicated problems--hence the general linear model.
3. We can even develop an intuitive approach to classification problems with this model.
4. The linear model is the ultimate _explainable_ model.

We shall also review key ideas in data analysis such as Principle Component Analysis (PCA); under-fitting, over-fitting and Regularisation; and Correlation in data.

> We will keep our focus on useful insights and practical tips in order to become a full-blown ninja when using general linear models in your own projects.

**In this course we will cover:**

- An intuitive introduction to general linear models--learn the core foundations needed for building more advanced models.
- Building in key concepts of linear models to more advanced (general) linear models.
- A deep dive into the linear models that give you the core foundations needed for building more advanced models.
- Best practices and tips for building and assessing models applied to your data.

**We‘ll work on the basis that you’re pretty new to Python**, but have some basic understanding of fundamental programming concepts and can run code locally on your machine.

Linear models might not seem that exciting, but they provide the foundation you need to understand advanced topics like Deep Neural Networks or Gaussian Processes further down the line. In this course, we’ll work through key concepts and techniques, including Predictive Power Score (or PPS) and Principle Component Analysis (PCA). We'll also cover the important topic of over-fitting and how we address this using regularisation methods.

Prof Tim Dodwell

## Becoming a Ninja in General Linear Models : Seciton Overview

We are onto Section 3 - the final section of this course. Here we are going to look at linear models, or to be more precise general linear models.

They might not be the sexiest of machine learning models, but they provide the foundations for many more sophisticated algorithms which include deep neural networks and Gaussian Processes, which you will learn about throughout more of the courses in digiLab academy. They also give us the stage to talk about certain key principals in machine learning, like overfitting and regularisation in a simplier model setting.

In this section we will cover the following topics:

- Step back to **curve fitting** as a general process, with straight lines and polynomials in general.

- Talk about **General Linear Models**, how do we fit general linear model made up of general features and see how these can be very rich and expressive models.

- Understand the core concetps of **underfitting and overfitting**, and then how do we overcome these challenges in particular by using **regularisation**.

- Talk about approaches for discovering relationships in data by looking at techniques for **linear and nonlinear relationships**.

- Work through our first unsupervised learning approach for dimension reduction, **Principle Component Analysis**. Which fits with the view of linear models, since it is a linear transform of our input data.

- Finally we discover **Logistic Regression**, in particular how we can extend the linear models for regression tasks for classification problems.


Section Overview

<Callout type="success" heading={null}>
In this explainer we will give a general introduction to linear models.

- We will talk through what a linear model is?
- How it works for simple linear models, like fitting polynomials which we probably know about before.
- Then we look how we can build linear model, which are more expressive.
- We will show how linear models are trained, both for example feature problems and data sets, and for much larger problems.
- Finally we will touch on the interpretation of bias, which comes naturally out of linear models.

</Callout>

<Callout type="warning" heading={null}>
**Starter Question.** Is the following a **linear model**?

$$
f(x) = x + 2\sin(2x) + \frac{1}{2}\exp\left(-\frac{x^2}{\lambda}\right)
$$

You should be able to confidently answer this question, and give a reason for your answer by the end of this explainer!

</Callout>

So in this explainer, we are very much looking at supervised learning methods. For this we will start with introducing Generalised Linear Models (or GLMs) in the context of **regression**.

<Callout type="info" heading={null}>
	A **regression task** is to predict one or more continous target variables $y$,
	given a $D$-dimesional vector $x$ of input variables.
</Callout>

So why are GLMs super important?

1. The foundations for introducing the core concepts in machine learning.
2. They are widely used and (can be) highly effective.
3. If possible to build, they are often more explainable that other models, in because with our understand of "correlation" as a central concept.
4. They provide a stepping stone to understand more complex machine learning models, in particular Neural Networks and Gaussian Processes.

## Fitting Polynomials

So we are probably used to fitting straight lines through data. I remember doing this at school

$$
f(x) = w_0 + w_1x
$$

where $w_0$ has the interpretation of the intercept and $w_1$ is the gradient (or slope) of the curve.

We can extend this to high-dimensional inputs, so what we mean by this is that there are two input variables ${\bf x} = [x_1, x_2]^T$, and so

$$
f({\bf x}) = w_0 + w_1x_1 + w_2 x_2,
$$

is a linear model (of linear functions) for a two dimensional input. Equally we can have models which are higher order polynomials, where we fit quadratic, cubic or high order curves through the data, so in general we have

$$
f(x) = w_0 + w_1x + w_2x^2 + \ldots + w_Kx^K = w_0 + \sum_{k=1}^Kw_k x^k.
$$

So the key point here is that when we talk about linear models. It doesn't mean the functions themselves are linear (although they can be). So we can also go higher dimensions (more than one input) and higher order

$$
f({\bf x}) = w_0 + w_1x_1 + w_{2}x_2 + w_3x_1^2 + w_4x_1x_2 + w_5x_2^2
$$

for example, but can be arbitary complex

## A More "General" Linear Model

A linear model is defined by the linear combination of functions. These functions which are known as **basis functions** or in Machine Learning **features**, and can be any general nonlinear function.

In general we can write a linear model as follows

$$
f({\bf x}) = w_0\phi_0({\bf x}) + w_1\phi_1({\bf x}) + \ldots + w_N\phi_K({\bf x}) = \sum_{i=0}^K w_i\phi_i({\bf x}),
$$

where each of the $w_i$'s (or weights) are scalar numbers or to be fancy we write $w_i \in \mathbb R$. That is $w_i$ belongs to the set of real numbers. In machine learning you might also see (and will see me write most the time) such a linear model written in vector form as follows

$$
f({\bf x}) = {\bf w}^T{\boldsymbol \phi} \quad \text{where} \quad {\boldsymbol \phi} = [\phi_0({\bf x}), \phi_0({\bf x}), \ldots, \phi_K({\bf x})]^T.
$$

Maybe to labour the point a little, this is just a more general case, so if we think of our original model we have $f(x) = w_0 + w_1x$, then this is simple

$$
\phi_0(x) \equiv 1 \quad \text{and} \quad \phi_1(x) = x,
$$

and hence

$$
f(x) = {\bf w}^T{\boldsymbol \phi} = [w_0 \quad w_1]^T \begin{bmatrix} 1 \\ x\end{bmatrix}.
$$

<Callout type="info" heading={null}>
It is often convention to reserve the first feature function, i.e. $\phi_0({\bf x})$ to be the constant function. This means that $\phi_0({\bf x})\equiv 1$, and therefore

$$
f({\bf x}) = w_0 + \sum_{i=1}^K w_i\phi_i({\bf x}_i)
$$

The parameter $w_0$ therefore allows us to model any fixed offset in the data, and with that is often called the _bias_ parameter (different from _bias_ in a statistical sense).

</Callout>

The basis or feature functions themselves can be pre-defined. So for example for a single input parameter example the basis could be powers of $x$, i.e. $\phi_k(x) = x^k$. Other options include Mexican Hat Function, where basis functions are defined on overlapping regions (think finite element shape functions) such that

Other basis functions include

$$
\phi_j(x) = \exp\left(-\frac{(x - \mu_k)^2}{2\lambda^2}\right)
$$

where are determined by a general point $\mu_k$ and a length scale $\lambda$. These are often called **Gaussian basis functions**. Or another common option is the sigmoidal basis functions which are defined by

$$
\phi_j(x) = \sigma\left(\frac{x - \mu_j}{s}\right) \quad \text{where} \quad \sigma(x) = \frac{1}{1 + \exp(x)}.
$$

To see what these look like on a plot, let's have a look - polynomials on the left, Gaussians in the middle and sigmoid basis functions on the right.

<Image
	src={
		"/images/courses/foundations-in-machine-learning/the-foundations-of-ml-curve-fitting/img1.jpg"
	}
	alt={
		"digiLab Academy - General Linear Models - basis functions | academy.digilab.co.uk"
	}
	width={1200}
	height={303}
	caption={""}
/>

We will also see later in the course that the basis functions don't have to be well defined formulae as we have above. They can actually be discovered from the data. Whilst outside the scope of this explainer, but will be covered else where, unsupervised machine learning models can be used. For those interested in looking ahead consider reading about **Principle Component Analysis** (aka PCA) and **Autoencoders**. More on that later.

## Training / Fitting a Linear Model

Important concept now is that given a set of basis functions $\phi_i({\bf x})$, we now wish to find the weights ${\bf w}$ which minimise the error on the training data - so we need to introduce the idea of fit!

So a very natural way to define error is for $N$ pairs of training data $({\bf x}_j, y_j)$, we can define the Mean-square loss function for a given set of weights

$$
\mathcal L = \frac{1}{2N}\sum_{j=1}^N \left( f_{\bf w}({\bf x}_j) - y_j \right)^2 = \frac{1}{2N}\sum_{j=1}^N \left( {\bf w}^T\boldsymbol \phi({\bf x}_j) - y_j \right)^2
$$

Our objective is therefore to minimise this **loss function**.

A natural way to do this is to find the gradient of error with respect to the weight's $w_i$, and use this information to learn how to adjust (for nonlinear case) or pick the weights to minimise the error. The gradient basically tells us how much the error will increase/decrease if we increase/decrease the value of $w_i$. Optimisation strategies then seek to change the value of the weights to reduce the error.

$$
\frac{\partial \mathcal L}{\partial w_i} = \frac{1}{N}\sum_{j=1}^N\left({\bf w}^T{\boldsymbol \phi}({\bf x}_j) - y_j\right)\phi_i({\bf x}_j)
$$

When we have reached a minimum error, then the gradient with respect to each of the weights will be zero.

In general, for nonlinear model, we would use some gradient based optimisation strategy for finding this optimal set of weights which minimises the error. Here, due to the linearity of the model, the solution can be calculated directly in closed-form. We won't go into the detail here, but give the solution. Let us first define the matrix, which is a rectangulat matrix often refered to as the design matrix

$$
\Phi =
\begin{bmatrix}
\phi_0({\bf x}_1) & \phi_1({\bf x}_1) & \ldots & \phi_K({\bf x}_1)\\
\phi_0({\bf x}_2) & \phi_1({\bf x}_2) & \ldots & \phi_K({\bf x}_2)\\
\vdots & \vdots & \ddots & \vdots \\
\phi_0({\bf x}_N) & \phi_1({\bf x}_N) & \ldots & \phi_K({\bf x}_N)
\end{bmatrix}
$$

Then the "optimal" weights for the model are given by

$$
{\bf w} = \left( \Phi^T \Phi\right)^{-1} \Phi^T{\bf y}
$$

We see that the calculation involves taking the inverse of a dense matrix i.e. $\left( \Phi^T \Phi\right)^{-1}$. This will be computationally expensive if the number of features $K$ and the number of data points $N$ is large. The time for calculation grows quadratically with $K$ and linearly in $N$. For large numbers of features ($K$) and therefore weights in such case we seek to use optimisation schemes like **gradient descent**.

<Callout type="info" heading={""}>

## Role of the bias parameter $w_0$

We will see in a work sample, how we can implement this and fit a linear model below.

Before we do this. It is a good opportunity to look at the role of $w_0$ the bias parameter, remembering that we conventionally choose $\phi_0({\bf x}) \equiv 1$. So if we look at what happens by setting the gradient of the error with respect to $w_0$ to zero

$$
\frac{\partial e}{\partial w_0} = \frac{1}{N}\sum_{j=1}^N\Big[{\bf w}^T{\boldsymbol \phi}({\bf x}_j) - y_j\Big]\phi_0({\bf x}_j) = 0
$$

So let us rewrite this a little noting that $\phi_0({\bf x}) \equiv 1$, then

$$
\frac{1}{N}\sum_{j=1}^N\left(w_0  + \sum_{k=1}^Kw_k\phi_k({\bf x}_j)- y_j\right) = w_0 - \underbrace{\frac{1}{N}\sum_{j=1}^Ny_j}_{\bar y} + \sum_{k=1}^Kw_k\underbrace{\frac{1}{N}\sum_{j=1}^N\phi_k({\bf x}_j)}_{\bar{\phi}_k}
$$

Ok, so all this equal zero, which means

$$
w_0 = \bar{y} - \sum_{k=1}^Kw_k\bar{\phi}_k
$$

What we notice is that the parameter $w_0$ makes a constant shift which corrects between the averages over the training set of the target variables and the weighted sum of averages of the basis functions evaluated at input points.

</Callout>


The Foundations of ML - Curve Fitting

In this walkthrough we will go through he basics of fitting linear models with sklearn using `LinearRegression()`. We will look at a simple polynomial model to get started, and focus on understanding how we can use `PolynomialFeature()` to build a general function on which to build a basic model.

### Adding Libraries

```python
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
```

### A Fishy Data set!

In this data set we are going to look at some fish data of 7 different species, where we will try and predict weight from simple measurements. Important stuff!

In total the data set contains 7 measurements of 159 fish.

```python
df = pd.read_csv('https://raw.githubusercontent.com/satishgunjal/datasets/master/Fish.csv')
df = df.drop(['Length1', 'Length2', 'Length3'], axis =1) # Can also use axis = 'columns'
df.sample(5) # Display random 5 records
```

<Image
	src={
		"/images/courses/foundations-in-machine-learning/regression-walkthrough/img1.jpg"
	}
	alt={"academy.digilab.co.uk"}
	width={400}
	height={196}
	caption={""}
/>

Ok Let us first look at the correlations in the data, for this we calculate and plot a correlation matrix.

```python
df.corr()

plt.rcParams["figure.figsize"] = (6,4) # Custom figure size in inches
sns.heatmap(df.corr(), annot =True)
plt.title('Correlation Matrix')
```

<Image
	src={
		"/images/courses/foundations-in-machine-learning/regression-walkthrough/img2.jpg"
	}
	alt={"academy.digilab.co.uk"}
	width={600}
	height={454}
	caption={""}
/>

```python
sns.pairplot(df, kind = 'scatter', hue = 'Species')
```

<Image
	src={
		"/images/courses/foundations-in-machine-learning/regression-walkthrough/img3.jpg"
	}
	alt={"academy.digilab.co.uk"}
	width={860}
	height={739}
	caption={""}
/>

If we want to predict weight, using height looks complex. There looks like two clear responses based on the fish type. Yet making predictions using width, looks like a much more unified relationship, much less depedent on fish type.

So let's do that.

```python
X = df['Width']
y = df['Weight']

X_train, X_test, y_train, y_test = train_test_split(X,y, test_size = 0.2)

plt.scatter(X_train,y_train, alpha = 0.4,color='lightcoral',label='Training Data')
plt.scatter(X_test,y_test, alpha = 0.8,color='lightblue',label='Testing Data')
plt.xlabel('Width')
plt.xlabel('Weight')
plt.legend()
plt.show()
```

<Image
	src={
		"/images/courses/foundations-in-machine-learning/regression-walkthrough/img4.jpg"
	}
	alt={"academy.digilab.co.uk"}
	width={1000}
	height={664}
	caption={""}
/>

There are clear outliers in this data set, yet since we have sufficient data their influence will be minimal. We could remove them, but we wont here.

### Building a Linear Model

Ok so the first simple model we are going to build is a polynomial model, which looks like this

$$
y = w_0 + w_1X + w_2X^2 = \begin{bmatrix} w_0, w_1, w_2 \end{bmatrix}^T\begin{bmatrix} 1 \\ X \\ X^2 \end{bmatrix} = {\bf w}^T \boldsymbol{\phi}
$$

So here we talk about $\boldsymbol{\phi}$ being the feature vectors, in this case they a polynomial features; really easy to do in ``sklearn`.

```python
from sklearn.linear_model import LinearRegression
from sklearn.preprocessing import PolynomialFeatures

polynomial_features= PolynomialFeatures(degree=2)

phi_train = polynomial_features.fit_transform(np.array(X_train).reshape(-1,1))

lr = LinearRegression().fit(phi_train, y_train)
```

Ok so now let us plot out the function and compare it to the data.

```python
X_all = np.linspace(1.0, 8.0, 1000).reshape(-1,1)
X_plot_poly = polynomial_features.fit_transform(X_all)
y_pred = lr.predict(X_plot_poly)
plt.scatter(X_train,y_train, alpha = 0.4,color='lightcoral',label='Training Data')
plt.scatter(X_test,y_test, alpha = 0.8,color='lightblue',label='Testing Data')
plt.plot(X_all,y_pred,'-k',alpha = 0.8,label='Linear Model')
plt.xlabel('Width')
plt.ylabel('Weight')
plt.legend()
plt.show()
```

<Image
	src={
		"/images/courses/foundations-in-machine-learning/regression-walkthrough/img5.jpg"
	}
	alt={"academy.digilab.co.uk"}
	width={1000}
	height={664}
	caption={""}
/>

It isn't a good model. Ok it fits on average ok, but we know that weight can't be negative. So we need to impose constraints.

`sklearn` gives us some options.

```python
lr_positive = LinearRegression(fit_intercept = False, positive = True).fit(phi_train, y_train)
y_pred_positive = lr_positive.predict(X_plot_poly)

plt.scatter(X_train,y_train, alpha = 0.4,color='lightcoral',label='Training Data')
plt.scatter(X_test,y_test, alpha = 0.8,color='lightblue',label='Testing Data')
plt.plot(X_all,y_pred,'-k',alpha = 0.8,label='Linear Model')
plt.plot(X_all,y_pred_positive,'-g',alpha = 0.8,label='Linear Model w. +ve constraint')
plt.xlabel('Width')
plt.ylabel('Weight')
plt.legend()
plt.show()
```

<Image
	src={
		"/images/courses/foundations-in-machine-learning/regression-walkthrough/img6.jpg"
	}
	alt={"academy.digilab.co.uk"}
	width={550}
	height={372}
	caption={""}
/>

Our model now at least make physically meaning predictions, although we see the performance of the model is reduce, particular at smalle width values.

How might we do better?

- Get rid of outliers?
- Try a different form of linear model?
- Apply a transform of the data first? If so, then what might you try?


Regression Walkthrough

<Callout type="success" heading={null}>
In this explainer we will talk about the key machine learning principles of under-fitting and over-fitting in the context of linear models. However, they apply generally to many types of machine learning models, the best known examples are Neural Networks.

Summary of the key learning outcomes are:

- Understand the importance and connection between over and underfitting and generalisation.
- We will demo what over and underfitting is with a simple polynomial model.
- We will described the general strategies we have for overcoming over-fitting.

</Callout>

For me Machine Learning is all about "Generalisation". To recap what we mean by generalisation. It is the ability of a model to extend to giving good predictions (or answers) on cases outside it's general training set. In general for me it is what we recognise as "intelligence" in an algorithm.

From this viewpoint, generalisation should be a (or some cases the) primary metric we would look at when evaluating a model. So when do models not generalise well?

So Let's think of a really simple example. Suppose we have an underlying process (a damped oscillator for example) which generates data as follows

$$
y_j = \exp\left(-\frac{1}{8}x_j^2\right)\sin \left(4\pi x_j\right) + \varepsilon
$$

where $\varepsilon \sim \mathcal N(0, 0.1)$ to simulate noisy data. Generating data $N + M$ evenly spaced on $x \in [0,1]$. We are going to do a "test/train split", so use $N$ for training (green dots), whilst keeping $M$ back for testing (red dots).

Ok, so now let us try and fit a polynomial model to this of the form

$$
f_K(x) = w_0 + w_1x + w_2x^2 +\ldots + w_Kw^K = \sum_{k=0}^Kw_0x^k.
$$

Let us make $K$ (the order of the polynomial) large, so $K=25$.

<Image
	src={
		"/images/courses/foundations-in-machine-learning/underfitting-and-overfitting/img1.jpg"
	}
	alt={
		"digiLab Academy - General Linear Models - basis functions | academy.digilab.co.uk"
	}
	width={550}
	height={421}
	caption={""}
/>

So what are we looking at here. The true funciton is shown as a solid black line, whilst the polynomial fit is shown in blue. We see that the fit to training points (shown in green) is very good, however the accuracy at unseen points (red points) is poor.

In this case we see that the model does not generalise well, to unseen data. It isn't a good model. This is a simple example of overfitting. So overfitting is characterised by a low training error, and a high validation error. Training curves with this sort of shape

<Image
	src={
		"/images/courses/foundations-in-machine-learning/underfitting-and-overfitting/img2.jpg"
	}
	alt={
		"digiLab Academy - General Linear Models - basis functions | academy.digilab.co.uk"
	}
	width={375}
	height={249}
	caption={""}
/>

The opposite to overfitting is underfitting. This is the case where the representative power of the model is not sufficient to describe the features in the data. So a much lower order polynomial will demonstrate this. So for example lets look at a third order polynomial

<Image
	src={
		"/images/courses/foundations-in-machine-learning/underfitting-and-overfitting/img3.jpg"
	}
	alt={
		"digiLab Academy - General Linear Models - basis functions | academy.digilab.co.uk"
	}
	width={570}
	height={415}
	caption={""}
/>

Here we see that a cubic (third order polynomial) is insufficiently expressive to describe the oscillations in the data. Here we see that the polynomial captures the general decaying trend of the data, but not the refined oscillations. This model is described as underfitting the data. Underfitting is characterised by both a high residual training error (i.e. it can be reduced beyond some point) and a high validation error.

Ok so how to build model which do not overfit or underfit, i.e. how do we build models which are "just right"!

In a polynomial model, there is a nice way to control this over and underfitting. We optimise over the order of the polynomial.
So let's look at the a plot of training and validation error against order of the polynomial.

So let us build models from order $k=5$ all the way up to $40$, and plot the validation error against the order polynomial.

<Image
	src={
		"/images/courses/foundations-in-machine-learning/underfitting-and-overfitting/img4.jpg"
	}
	alt={
		"digiLab Academy - General Linear Models - basis functions | academy.digilab.co.uk"
	}
	width={560}
	height={431}
	caption={""}
/>

We see that models of order $K=10$ to $K=15$ minimises the validation error, anything in this range gives a "good" model! Here we do not look at the training error, after the model order exceeds the amount of data the model fits the training data perfectly.

Ok so we can play with the order of a polynomial, but in general how are the ways in which we can control overfitting:

1. **Control the number of parameters in the model**. So this is the example that we have just looked. We can control the order of the polynomial for a simple polynomial model, but in general for machine learning we can control the number of parameters we have. So for example in a neural network this would be the number of nodes or layers we have, in a random forest the number of leaves we have. The more parameters we have the more likely we are going to overfit. In the next section will see general strategies for constraining the number of parameters in a model, rather than simple playing with the maximum order. There is a more sophisticated approach for control the number of parameters, which strategically drops parameters which offer little benefit to fit, this process is called **regularisation** which we talk about in a following explainer.

2. **Increase the amount of data**. Overfitting is a function of the amount of data. which is dependent on the number of parameters we have to play with, but also on the amount of data we have in relation to the number of parameters. Therefore it is more likely in cases where we have only limited data, and highly parameterised models. Limited data is a function of a large number of inputs, so what is enough data is often a complex questions. But in general collecting or using more data will help with overfitting.

3. **Dropout**. So we aren't going to cover this alot here, but a technique widely used in training of larger machine learning models like neural networks, which always have the tendency to be overfit to some degree, is called dropout. Dropout is a strategy used at training time, where a percentage of parameters are randomly set to zero at each training set. What does this do? It removes the dependence on a particular or all terms, and the general approach is one that is both expressive but robust in it's representation. This is often done alongside batch training or optimisation, where the loss or fit is only access over a subset of the data. More on this in our deep learning courses.

## Conclusion

Overfitting and Underfitting should be a central consideration when building any machine learning model. It is particularly important in cases where you have high dimensional (many) inputs and / or smaller datasets. A good model balances fitting and complexity, with the ultimate aim of ensuring good generalisation.

Whilst there are various strategies to control overfitting, the widely used approaches are more data collection, regularisation and dropout. We will look at regularisation in the next walk through.


Underfitting and Overfitting

## Controlling ML Models - Regularisation in Practise

<Callout type="success" heading={null}>
	In this walkthrough we are going to look at regularisation, a really central
	approach to prevent overfitting in machine learning models. Here you learn how
	to build a linear model, and apply regularisation methods (e.g. LASSO). We will
	clearly see the benefits of this approach.
</Callout>

### Adding Libaries

```python
import numpy as np
import pandas as pd

from sklearn.linear_model import LinearRegression
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import PolynomialFeatures

import matplotlib.pyplot as plt
```

### Loading the dataset, scaling and transformation into features.

After loading various libraries, we now collecting the Boston housing data set.

In this example we will consider a single input and a single output for our regression challenge.

$x$ - the input is LSTAT which is a measure of the number of low income families

$y$ - the ouput is the median (local) house price over the data set

```python
url = "https://raw.githubusercontent.com/jbrownlee/Datasets/master/housing.csv"

df = pd.read_csv(url, header = None)

data = df.values

X = data[0:50, 12].reshape(-1, 1) # Single input data - LSTAT
y = data[0:50, 13].reshape(-1, 1) # Target Variable - Median House Price
```

Nice and easy here. Read data in upload it into a dataframe.

We now can get input and output data in good shape. The call `reshape(-1,1)` makes sure that the array vector and so the shape is of the array is `(n,1)` rather just `(n,)`.

We then move on to a sklearn pipeline. Here we first apply a scaling `MinMaxScaler()` to scale this input between $0$ and $1$. Instead of looking at a single feature, we then expand the input representation to a polynomial feature space of order 10, i.e.

$$
\boldsymbol{\phi} = \Big[ 1, x, x^2, \ldots, x^9, x^{10} \Big]
$$

```python
from sklearn.pipeline import Pipeline
from sklearn.preprocessing import MinMaxScaler

pipe = Pipeline(
    [
        ("minmax", MinMaxScaler()),
        ("feature", PolynomialFeatures(degree=10)),
    ]
)

X_poly = pipe.fit_transform(X)
```

### Test / train split of the data

We then do the normal test train / split of the data set.

```python
X_train, X_test, y_train, y_test = train_test_split(X_poly, y, test_size=0.3, random_state=123)
```

### Looking at the data

```python
fig, ax = plt.subplots(1, 1)

ax.plot(X_train[:,1], y_train, 'ob', alpha=0.4)
ax.plot(X_test[:,1], y_test, 'og', alpha=0.4)

plt.ylabel('Output - Media House Value')
plt.xlabel('Input - LSTAT (scaled)')

plt.show()
```

<Image
	src={
		"/images/courses/foundations-in-machine-learning/controlling-ml-models-regularisation-in-practise/img1.jpg"
	}
	alt={"academy.digilab.co.uk"}
	width={560}
	height={430}
	caption={""}
/>

## Fitting a linear model

Ok so now it is time to fit the linear model.

```python
lr = LinearRegression().fit(X_train, y_train)
```

Then let us plot out how to see how we did over the range of the data set.

```python
X_plot = np.linspace(0.0, 1.0, 100).reshape(-1,1)
X_plot_poly = PolynomialFeatures(degree=10).fit_transform(X_plot)

y_plot = lr.predict(X_plot_poly)

fig, ax = plt.subplots(1, 1)

ax.plot(X_train[:,1], y_train, 'ob', alpha=0.4, label='Training Data')
ax.plot(X_test[:,1], y_test, 'og', alpha=0.4, label='Testing Data')
ax.plot(X_plot, y_plot, '-',color='lightcoral', label='Polynomial Fit')

plt.ylabel('Output - Media House Value - Thousand Dollars')
plt.xlabel('Input - LSTAT (scaled) ')

plt.ylim([0.0, 50.])
plt.xlim([0.0, 1.])

plt.show()
```

<Image
	src={
		"/images/courses/foundations-in-machine-learning/controlling-ml-models-regularisation-in-practise/img2.jpg"
	}
	alt={"academy.digilab.co.uk"}
	width={570}
	height={436}
	caption={""}
/>

We see this doesn't look good, particularly to the righthand side between the values of 0.7 of 1.0, the predicted model does not generalise well.

The model is too expressive, the model is clearly overfitting the data.

### Applying Regularisation

Let us follow what we discussed in the previous explainer on regularsiation, and use a LASSO model to fit to the data.

We aren't sure what the best regularisation strength $\alpha$ (which is $= 1/lambda$ from the notes) should be. We therefore fit lots of models over a range of values between $0.001$ and $0.2$. Here we use a simple loop but we could of used `sklearn`'s grid search functions

```python
from sklearn.linear_model import Lasso

alpha = np.linspace(0.001, 0.2, 1000)

training = []
testing = []

for a in alpha:
    lass = Lasso(alpha=a).fit(X_train, y_train)
    training.append(lass.score(X_train, y_train))
    testing.append(lass.score(X_test, y_test))
```

We can now plot out the training and testing scores

```python
plt.plot(alpha, training, '-b', label='Training')
plt.plot(alpha, testing, '-g', label='Testing')

plt.ylabel('Testing and Training')
plt.xlabel('Regularisation Strength')


plt.legend()
plt.show()
```

<Image
	src={
		"/images/courses/foundations-in-machine-learning/controlling-ml-models-regularisation-in-practise/img3.jpg"
	}
	alt={"academy.digilab.co.uk"}
	width={570}
	height={428}
	caption={""}
/>

<Callout type="info" heading={null}>
The function ```score()``` returns the coefficient of determination of the prediction.

The coefficient of determination is given by

$$
R^2 = \Big( 1 - \frac{u}{v} \Big)
$$

where $u$ is the residual sum of squares ((y_true - y_pred)** 2).sum() and $v$ is the total sum of squares ```((y_true - y_true.mean()) ** 2).sum()```.

The best possible score is 1.0 and it can be negative (because the model can be arbitrarily worse). A constant model that always predicts the expected value of y, disregarding the input features, would get a score of $0.0$.

</Callout>

Ok let us pick the regularisation strength where we do best over our testing data. So about a value of $0.02$

```python
lass = Lasso(alpha=0.025).fit(X_train, y_train)

y_plot = lass.predict(X_plot_poly)

fig, ax = plt.subplots(1, 1)

ax.plot(X_train[:,1], y_train, 'ob', alpha=0.4, label='Training Data')
ax.plot(X_test[:,1], y_test, 'og', alpha=0.4, label='Testing Data')
ax.plot(X_plot, y_plot, '-',color='lightcoral', label='Polynomial Fit')

plt.ylabel('Output - Media House Value - Thousand Dollars')
plt.xlabel('Input - LSTAT (scaled)')

plt.ylim([0.0, 50])

plt.show()
```

<Image
	src={
		"/images/courses/foundations-in-machine-learning/controlling-ml-models-regularisation-in-practise/img4.jpg"
	}
	alt={"academy.digilab.co.uk"}
	width={560}
	height={436}
	caption={""}
/>

This looks alot better. No oscillations that try to model the noise in the data, generalises well particularly between input values of $0.7$ and $1$ where the other model performed very badly.

We can now look to see what the final coefficients of this model were

```python
print(lass.intercept)
print(lass.coef_)
```

```
array([ 31.27 ])
array([  0.        , -33.77058217,   0.        ,  17.24305027,
         0.        ,   0.        ,  -0.        ,  -0.        ,
        -0.        ,  -0.        ,  -0.        ])
```

So this gives us a nice simple model as the prediction

$$
f(x) = 31.27 - 33.77x + 17.24x^3
$$

We go through this result in more detail in the explainer.

<Callout type="success" heading={null}>
Now you fitted a LASSO model why not try to fit a Ridge Regression model and compare the outputs of the functions you fit.

Hint you need to import the following

```python
from sklearn.linear_model import Ridge
```

which can be used in the same way as `Lasso`

</Callout>


Controlling ML Models - Regularisation in Practise

<Callout type="success" heading={null}>
**Correlation** is a statistical measure of the dependencey, casual or not, between two variables. In this explainer we consider four key points:

- The Pearson's Correlation Coefficient.
- The Spearman's (or Rank) Correlation Coefficient.
- How to show and interpret correlation plots in many variable problems.
- Understand the key limitations in expressing general correlations within real world data.

</Callout>

Whilst in general, correlation may mean any type association between variables; it typical refers to the how linearly dependent two variables are. When people talk about the correlation coefficient between variables, they typically mean the **Pearson's correlation coefficient** (we will talk about others below though). The formulate for Pearon's correlations coefficient is

$$
\rho_{X,Y} = \text{Corr}(X,Y) = \frac{\text{Cov}(X,Y)}{\sigma_X\sigma_Y} = \frac{\mathbb E [(X - \mu_X)(Y - \mu_Y)]}{\sigma_X\sigma_Y}.
$$

When this comes down to calculating it in practice, given two sets of variables $X = \{x_0, x_1, \ldots, x_N \}$ and $Y = \{y_0, y_1, \ldots, y_N \}$, Pearson's correlation can be computed

$$
\tilde{\rho}_{X,Y} = \frac{\sum_{j=0}^N(x_j - \bar{x})(y_j - \bar{y})}{\sqrt{\sum_{j=0}^N(x_i - \bar{x})^2\sum_{j=0}^N(y_i - \bar{y})^2}}
$$

The correlation coefficient which ranges between $-1$ and $+1$.

- The correlation coefficient is $+1$ if two variables have a perfect positive / increasing linear relationship.

- The correlation coefficient is $-1$ if two variables have a perfect negative / decreasing linear relationship.

- If the variables are independent, Pearson's correlation coefficient is $0$, but the converse is not true. This is because the correlation coefficient detects only linear dependencies between two variables. Here are ?? examples, where the correlation is $0$, but there is clearly a nonlinear relationship between variables in each case.

<Callout type="info" heading={null}>
	In practise correlation coefficients will be estimated by a finitie set of
	samples. Therefore the means $\mu_X$ and $\mu_Y$, and the expectations $\mathbb
	E[\cdot]$ will be replaced by sample means, and standard deviations computed
	using unbiased estimates. For more information on the later - have a look at -
	[unbiased estimates of standard deviation of
	wiki](https://en.wikipedia.org/wiki/Unbiased_estimation_of_standard_deviation)
</Callout>

## Rank Correlation Coefficients

Spearman's rank coefficient, is a different measure of the dependency between two variables. The measure is nonparametric measure, this means that it does not assume any parameterised distribution (i.e. normal) or it does not impose a particular type of relationship between variables.

Spearman's rank coefficient is calculated based on the ranks of the values rather than the raw data. The coefficient can range from -1 to 1, with -1 indicating a strong negative relationship, 1 indicating a strong positive relationship, and 0 indicating no relationship.

To calculate Spearman's rank coefficient, the ranks of the values in each variable are determined, and the difference between the ranks is calculated for each pair of values. The differences are then squared, and the sum of the squares is divided by the number of pairs. The resulting value is then used to calculate the coefficient using the following formula:

$$
\rho = 1 - \frac{6 \sum_{j=1}^n (r_X^{(j)} - r_Y^{(j)})^2}{ n(n^2 - 1)}
$$

where $n$ is the sample size, and $r_X^{(j)}$ and $r_Y^{(j)}$, is the ranking of the $j^{th}$ component of the datasets $X$ and $Y$.

Spearman's rank coefficient is often used when the assumptions of Pearson's correlation coefficient, which is based on linear relationships, are not met.

## Predictive Body Fat Example

For this example we use a real world data set - called the 'Body Fat Prediction Dataset'. It can be found to download from the kaggle website - [here](https://www.kaggle.com/datasets/fedesoriano/body-fat-prediction-dataset)!

<Callout type="info" heading={null}>
	The **Body Fat Prediction Dataset** is a data set of measurements of 252 men's
	body fat, along with $15$ other measurements, e.g. neck, chest and wait
	circumference. It is a really good data set to illustrate multiple regression
	techniques. We use it in a number of our examples.
</Callout>

The downloaded csv (`bodyfat.csv`) is then loaded into a `pandas` dataframe. Columns for 252 measurements of `Abdomen` and `BodyFat` can be convert to lists $x$ and $y$ respectively as follows

```Python
from scipy.stats import pearsonr
import matplotlib.pyplot as plt
import seaborn as sns
import pandas as pd
import numpy as np

from sklearn.linear_model import LinearRegression

df = pd.read_csv('data/bodyfat.csv')

# Call the scipy function for pearson correlation
corr, pVal = pearsonr(x= df['Abdomen'], y=df['BodyFat'])

# Also make a linear fit using sklearn's function Linear Regression

x = df['Abdomen'].tolist()
y = df['BodyFat'].tolist()

model = LinearRegression().fit(np.array(x).reshape(-1, 1), np.array(y).reshape(-1, 1))

# Plot a Scatter & Linear of Best Fit

plt.scatter(x, y, marker = 'o', c = 'g', alpha=0.3)

plt.plot(np.linspace(min(x),max(x),2), model.coef_[0]*x_grid + model.intercept_[0], '-b', alpha=0.3, label = 'Linear Fit')

# Lots of plotting options can be added - not not included here

plt.show()
```

So the output of this segment of code, with a few extra plotting commands to make it look pretty is.

<Image
	src={
		"/images/courses/foundations-in-machine-learning/exploring-simple-relationships-correlations/img1.jpg"
	}
	alt={"academy.digilab.co.uk"}
	width={580}
	height={458}
	caption={""}
/>

So here for this dataset, _Person Correlation Coefficient_ is $0.81$ (p-value is 0.0, see below) and coefficient of the linear fit $0.63$. This makes it clear that the gradient of best linear fit is not the same as correlation.

<Callout type="info" heading={null}>
	The `scipy` function `personr()`, not only returns a correlation, but also a
	p-value. This means that under the hood, the function performs a test, with the
	null hypothesis that the distributions underpinning the samples are actually
	uncorrelated and normally distributed. The p-value gives the probability of an
	uncorrelated system producing the datasets that have a sample correlation as
	"strong" (either positive or negative) as the one compuited from this dataset.
</Callout>

```Python
from scipy.stats import spearmanr
corr, pVal = spearmanr(x, y)
```

## Correlation Plot

To plot a correlation plot in Python, you can use the seaborn library. Here is an example of how you might create a correlation plot using the `seaborn.heatmap()` function:

```Python
import seaborn as sns
import pandas as pd

sns.set(rc={'figure.figsize':(11.7,8.27)})

# Calculate the correlation matrix
corr = df.corr()

# Plot the correlation matrix as a heatmap
sns.heatmap(corr, xticklabels=corr.columns, yticklabels=corr.columns, cmap='RdBu', annot=True)
```

This will create a heatmap that shows the correlation between the variables x and y. The heatmap uses color to indicate the strength and direction of the correlation, with blue indicating a strong negative relationship and red indicating a strong positive relationship.

<Image
	src={
		"/images/courses/foundations-in-machine-learning/exploring-simple-relationships-correlations/img2.jpg"
	}
	alt={"academy.digilab.co.uk"}
	width={950}
	height={738}
	caption={""}
/>

You can customize the appearance of the plot by adjusting the arguments passed to the heatmap function. For example, you can change the color map, adjust the font size, or add a title to the plot.

## Conclusions : Good and Bad

Correlation is a useful tool to understanding linear dependencies in data. If variables are linearly dependent, it provides perfect information for building simple predictive models, or doing engineer dimension reduction (i.e. elimination of dependent variables).

<Callout type="warning" heading={null}>
	Pearson's correlation coefficients are based on the assumption of **linear
	relationship**, it clearly will fall over where there are nonlinear
	correlations between two variables.
</Callout>

Below is a series of engineered examples which demonstrate relationships, each of which have zero Pearson's correlation coefficient.

<Image
	src={
		"/images/courses/foundations-in-machine-learning/exploring-simple-relationships-correlations/img3.jpg"
	}
	alt={"academy.digilab.co.uk"}
	width={1200}
	height={189}
	caption={""}
/>

For such nonlinear relationships, other forms of nonlinear correlation metrics are required.

This leads us to a discussion of Predictive Power Score (PPS). PPS is an asymmetric, data-type-agnostic score that can detect linear or non-linear relationships between two variables. Which we will look at.

Want to read more have a look at this article ["Are Correlations any guide to predictive value?"](https://www.jstor.org/stable/2985494). This is not a new challenge this article was published by Journal Royal Statistical Society in 1956 by Dr Robert Ferber!


Exploring Simple Relationships - Correlations

<Callout type="success" heading={null}>
The **Predictive Power Score** (PPS) can be used to detect linear and non-linear relationships between two variables. The score ranges from 0 to 1, from no predictive power to perfect predictive power. It can be used as a powerful alternative to (linear) correlation - and over comes some of the key issues the linearity assumptions which underpin correlation analysis.

In this explainer, we demonstrate how PPS can be used:

1. How it is a PPS used as an early diagnositic tool in exploring the relationships between variables in data?
2. How does PPS address the key limitations of assuming linearity in standard correlations measures (see . link).
3. How it can be used in model design and inform feature engineering or model reduction in practise.

</Callout>

## Explainer through a toy problem.

As the name suggests predictive power score gives a score (between 0 and 1) which indicates how informative one variable ($x$) is in predicting the value of another $y$.

A key difference between PPS and correlations is that this measure is in general non-symmetric. What do we mean by this?

Well let's consider the following relation between two parameters $x$ and $y$,

$$
y = x^2 + \epsilon.
$$

Here $\epsilon$ is a noise parameter defined by a parameterised Gaussian e.g. $\mathcal{N}(0, \sigma^2)$. The plot below shows the data.

<Image
	src={
		"/images/courses/foundations-in-machine-learning/finding-nonlinear-relations-predictive-power-score/img1.jpg"
	}
	alt={"academy.digilab.co.uk"}
	width={600}
	height={436}
	caption={""}
/>

In this simple example the $x$ has a strong predictive power in terms of $y$. It can exploit the underlying quadratic relationship. However, if I give you a value of $y$ the associated value of $x$ is ill-posed. What do I mean by this? For this simple quadratic, there are two underpinning values which are good predictions. The variable $y$ has poor predictive power when trying predicting $x$.

The PPS has be developed in an open python library called `ppscore` this library can be found [https://github.com/8080labs/ppscore](https://github.com/8080labs/ppscore), it's use is nice and simple.

To install we start with install the package.

```Python
pip install -U ppscore
```

Then from here we can do this example as follows

```Python
import pandas as pd
import numpy as np
import ppscore as pps

df = pd.DataFrame()
df["x"] = np.random.uniform(-1, 1, 1000)
df["y"] = df["x"] * df["x"] + np.random.normal(0.0, 0.05,  1000)

pps.matrix(df)

import seaborn as sns
matrix_df = pps.matrix(df)[['x', 'y', 'ppscore']].pivot(columns='x', index='y', values='ppscore')
sns.heatmap(matrix_df, vmin=0, vmax=1, cmap="Blues", linewidths=0.5, annot=True)
```

The last part is the key output but which generates our matrix plot, which we can identify as the output from PPS, this looks like

<Image
	src={
		"/images/courses/foundations-in-machine-learning/finding-nonlinear-relations-predictive-power-score/img2.jpg"
	}
	alt={"academy.digilab.co.uk"}
	width={535}
	height={438}
	caption={""}
/>

So let's go through this output. We see that the diagonal is and always will be 1. A variable predicts itself perfectly (not surpise here). The bottom left hand corner says that if I know $x$ then I can predict $y$ with confidence. The $0.76$ means it is not perfect, and this is driven by the noise. Finally as expected the predictive power of $y$ predict $x$ is $0$ this is because the problem that there are two values of $x$ each $y$. Simple as that.

## How does it work under the hood?

Under the hood, the PPS algorithm uses a machine learning algorithm called **XGBoost** (eXtreme Gradient Boosting) to learn the mapping from the feature variable to the target variable. It therefore considers pair-wise predictive capabilities from one variable to another.

This idea of pairwise predictions is important. It therefore misses complex dependents which involve the interaction of multiple features in a non-linear way.

<Callout type="info" heading={null}>
**XGBoost is a powerful and widely used gradient boosting framework** that uses decision trees as base models. Gradient boosting is an ensemble learning technique that combines multiple weak learners (in this case, decision trees) to create a strong predictor.

The **XGBoost** algorithm is particularly well-suited for PPS because it can handle both numerical and categorical features, it can handle missing values, and it can capture complex relationships between variables.

</Callout>

## Applications of the PPS and the PPS matrix

After we learned about the advantages of the PPS, let’s see where we can use the PPS in the wild.

1. **Find patterns in the data:** The PPS finds every relationship that the correlation finds — and more. Thus, you can use the PPS matrix as an alternative to the correlation matrix to detect and understand linear or nonlinear patterns in your data. This is possible across data types using a single score that always ranges from 0 to 1.

2. **Feature selection:** In addition to your usual feature selection mechanism, you can use the predictive power score to find good predictors for your target column. Also, you can eliminate features that just add random noise. Those features sometimes still score high in feature importance metrics. In addition, you can eliminate features that can be predicted by other features because they don’t add new information. Besides, you can identify pairs of mutually predictive features in the PPS matrix — this includes strongly correlated features but will also detect non-linear relationships.

3. **Detect information leakage:** Use the PPS matrix to detect information leakage between variables — even if the information leakage is mediated via other variables.

4. **Data Normalization:** Find entity structures in the data via interpreting the PPS matrix as a directed graph. This might be surprising when the data contains latent structures that were previously unknown.

## Back to the Body fat Prediction Dataset

Just to go back to the example we looked at in the correlation explainer, around measurements from 252 males which provided weight, age, hip, chest, adohem measurements alongside body fat measurements. Here is the correlation heatmap we had in the "correlation" explainer

<Image
	src={
		"/images/courses/foundations-in-machine-learning/finding-nonlinear-relations-predictive-power-score/img3.jpg"
	}
	alt={"academy.digilab.co.uk"}
	width={950}
	height={738}
	caption={""}
/>

As we see that correlation give measures on how well two variables are linearly dependent. Correlation is symmetric, and gives limited information if relationships are nonlinear, and with that how well one variable can predict another.

We apply `ppscore` to this dataset.

```Python
import ppscore as pps
import seaborn as sns
matrix_df = pps.matrix(df)[['x', 'y', 'ppscore']].pivot(columns='x', index='y', values='ppscore')
sns.heatmap(matrix_df, vmin=0, vmax=1, cmap="Blues", linewidths=0.5, annot=True)
```

Here is the result

<Image
	src={
		"/images/courses/foundations-in-machine-learning/finding-nonlinear-relations-predictive-power-score/img4.jpg"
	}
	alt={"academy.digilab.co.uk"}
	width={960}
	height={764}
	caption={""}
/>

- Clear that if I want to predict Body Fat I look along the row associated with body fat. We see that two easy measurements we could use (using a tape measure) are measuring Abdomen ($x_1$) and Chest ($x_2$).

### Application 1. Spotting Data Leakage

First let's look at the "Density", we see that density has a very strong predict power, and this is symmetric. By this we mean that BodyFat and Density are both equally good predictors of each other.

On the face of it, this makes it looks like a perfect candidate for an input feature when building a machine learning model. But this would actually be a dishonest model. This is because BodyFat is actually a derived quantity from Density, and therefore in calculating one your have the other. This defeats the point of needing a simple predictive model for BodyFat, rather than taking expensive measurements.

As a general rule, inputs which have such strong (symmetric) predictive power score indicate "something is going on". Here best practise is to explore what is the relationship between these variables.

This is an example of "data leakage", where you are inplanting direct information into the inputs about the outputs without knowing about it.

<Callout type="warning" heading={null}>
	We have additional examples of this throughout the course. If you want a great
	example of data leakage have a look at [this twitter
	feed](https://twitter.com/svpino/status/1592140348905517056) for some.
</Callout>

### Application 2. Predictive Power Score as part of your Machine Learning Workflow

Using the PPS score in the example above it clear the two easily measureable features I can use (Abdomen and Chest) to build a predictive model for Body Fat. PPS can therefore be used in an early exploration phase, which can inform appropriate model selection.

Whilst Machine Learning models can work with large numbers of inputs (and outputs), the required amount of data will increase with the number of input features. It is therefore prudent to select those only which are most informative with respect out the outputs. This is where predictive power comes in.

Let's consider

$$
\hat{y} = f({\bf x}) = w_1 x_1 + w_2 x_2 = {\bf w}^T{\bf x},
$$

where $x_1$ represents Abdomen and $x_2$ Chest measurements.

<Callout type="info" heading={null}>
**Challenge for you.** Can you write a bit of code to fit this Linear Model?

You can find the data set you need here . . . [body fat link](https://www.kaggle.com/datasets/fedesoriano/body-fat-prediction-dataset)

</Callout>


Finding Nonlinear Relations - Predictive Power Score

:::success
**📂 Resources**
Download the resources for this lesson <Link href="/resources/general-linear-models-for-machine-learning/correlation-and-pps-walkthrough/Resources.zip" download="Resources.zip">here.</Link>
:::

In this walkthrough we walk through methods for finding linear correlations and exploring nonlinear relationships with in new data sets. Look at a range of reall datasets as well as synthetics examples to demonstrate the points.

### Adding libaries

```python
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import seaborn as sns
```

## The first data set - what is the body fat?

In this first data set we look at just over 250 measurements of males, for which we have lots of different body measurements as well as the results of a body fat measurement. Can we build a good predictive model between the measurements (which are easy to do at home) with a percentage of body fat (more complex to measure).

As a starting point, can we we understanding the correlations in our data, which will guide us in building a good linear model using our data.

Let us first look at the data

```python
df = pd.read_csv('bodyfat.csv')

df.sample(5)
```

<Image
	src={
		"/images/courses/foundations-in-machine-learning/correlation-and-pps-walkthrough/img1.jpg"
	}
	alt={"academy.digilab.co.uk"}
	width={1000}
	height={163}
	caption={""}
/>

First we are going to flyby and show how we calculate the Pearson and Spearman correlations in our dataset. For this we use the well known package `scipy.stats`, which makes this easy.

```python
from scipy.stats import pearsonr, spearmanr

corr_p, pVal_p = pearsonr(x=df['Abdomen'], y=df['BodyFat'])

corr_s, pVal_s = spearmanr(df['Abdomen'], df['BodyFat'])

print(pearsonr(x=df['Abdomen'], y=df['BodyFat']))
print(spearmanr(df['Abdomen'], df['BodyFat']))
```

Here not only do we spit out correlations but we get a p - value. A p-value is about how confident we are in a correlation, it gives the probability that a random dataset would display the same correlations. Here we can see very strong correlations, where we are confident that they are inherent in the data we are observing not simple random from sampling.

### Plotting a Correlation Matrix

`pandas` provides a really nice function to be able to calculate all the correlations across variables, combined with this we can use plotting functions in `seaborn` to make really nice visuals of correlations matrices.

```python
sns.set(rc={'figure.figsize':(11.7,8.27)})

corr = df.corr()

sns.heatmap(corr, xticklabels=corr.columns, yticklabels=corr.columns, cmap='RdBu', annot = True)
```

<Image
	src={
		"/images/courses/foundations-in-machine-learning/correlation-and-pps-walkthrough/img2.jpg"
	}
	alt={"academy.digilab.co.uk"}
	width={950}
	height={738}
	caption={""}
/>

### Use the correlations to build a Linear Model

Here we our discovery of correlations can be used to help build a simple model. We see a strong correlation between Abdomen measurement and Body Fat.

```python
from sklearn.linear_model import LinearRegression

x = df['Abdomen'].tolist()
y = df['BodyFat'].tolist()

linearModel = LinearRegression().fit(np.array(x).reshape(-1,1), np.array(y).reshape(-1,1))

x_grid = np.linspace(np.min(x), np.max(x), 2)

plt.scatter(x, y, marker='o', c = 'g', alpha = 0.4)

plt.plot(x_grid, linearModel.coef_[0]*x_grid + linearModel.intercept_[0], '-b', alpha = 0.4)

plt.show()
```

<Image
	src={
		"/images/courses/foundations-in-machine-learning/correlation-and-pps-walkthrough/img3.jpg"
	}
	alt={"academy.digilab.co.uk"}
	width={1000}
	height={752}
	caption={""}
/>

## When do correlations not work!

Ok so let's look at an example which has zero correlation! Here is a model

$$
y = x^2 + \epsilon
$$

where $\epsilon$ is noise. Let's look at this data.

```python
np.random.seed(123)

N = 500

X = -1. + 2.0 * np.random.uniform(size=(N,))
y = X**2 + np.random.normal(0.0, 0.05, size=X.shape)

plt.scatter(X,y, marker='o', c='g', alpha = 0.2)
plt.show()
```

Ok, let's test these correlations. We there is nothing there but there is clearly a relationship!

```python
print(spearmanr(X,y))
print(pearsonr(X,y))
```

where the output is

```python
SpearmanrResult(correlation=-0.008570914283657133, pvalue=0.8483882540669414)
PearsonRResult(statistic=0.018690496801745102, pvalue=0.6767350220071778)
```

This data has no linear correlation. We have talked about this in the explainer.

### Predictive Power Score - PP score

So here we apply predictive power score - which we have talked about in explainers. If you haven't install it then you can do this using `pip install -U ppscore`.

To use PP score we need to stuff our data back in a dataframe.

```python
import ppscore as pps

df = pd.DataFrame()
df["X"] = X
df["y"] = y
```

We can now simple call `ppscore`

```
pps.matrix(df)
```

So look at this in a nice way we can then plot the results of the output as well. So let's do that.

```python
matrix_df = pps.matrix(df)[['x', 'y', 'ppscore']].pivot(columns='x', index='y', values='ppscore')
sns.heatmap(matrix_df, vmin=0, vmax=1, cmap="Blues", linewidths=0.5, annot=True)
```

<Image
	src={
		"/images/courses/foundations-in-machine-learning/correlation-and-pps-walkthrough/img4.jpg"
	}
	alt={"academy.digilab.co.uk"}
	width={800}
	height={645}
	caption={""}
/>

So here, we first say the clear thing, the diagonal are 1. This is because if you have a variable say $x$ then you can perfectly predict itself, clear.

Here we see that given $y$ we can not predict $x$ this is because for any $y$ in our data there is a both a positive and negative value of $x$ which would make a possible output. So in summary the predict power of $y$ for variable $x$ is zero.

Alternatively, we see that the predictive power of $x$ to predict $y$. Here we see a pp score of 0.78.

### Applying it to the Body Fat dataset

So let's do that

```python
df = pd.read_csv('bodyfat.csv')


sns.set(rc={'figure.figsize':(11.7,8.27)})

matrix_df = pps.matrix(df)[['x', 'y', 'ppscore']].pivot(columns='x', index='y', values='ppscore')
sns.heatmap(matrix_df, vmin=0, vmax=1, cmap="Blues", linewidths=0.5, annot=True)
```

We see that pp-score pulls out a much simplier picture for building a model than correlations. For body fat we see that the strongest predictor s are Abdomen and Chest. We also see that the Density has a huge predictive power, this is actually false since you need density to calculate body fat, so they are not independent variables.


Correlation & PPS Walkthrough

<Callout type="success" heading={null}>
In this explainer we look at how we can "generalise" linear models to work for classification problems. This is called **logistic regression**. At the end of this explainer your should:

- Understand the key idea of logistic regression.
- Understandit's make up (e.g. linear predictor + link function).
- Understand how logisitic models are trained (e.g. loss function + gradient descent).
- Know how logistic regression can be extended to multiple classes.
- Know the issues of overfitting, and the solutions to this (regularisation).

</Callout>

Logistic regression is a statistical model which is used for classification problems. Logisitic regression estimates the probabilility $p$ that an event will occur. Hence the output of the model is between 0 and 1.

So we have a supervised learning problem, with our normal data set

$$
\{({\bf x}_0, y_0), ({\bf x}_1, y_1), \ldots, ({\bf x}_N, y_N) \}.
$$

Where ${\bf x}_i$ are our inputs and $y_i$ are our target variables, in this case a binary label of either $A$ or $B$.

Instead of fitting a straight line or hyper-plane (or any linear model), the logistic regression model uses the logistic function to squeeze the output of a linear equation between $0$ and $1$ so it represents a probability of giving label $A$.

The logistic function is defined as:

$$
p_A({\bf x}) = \frac{1}{1 + \exp(-f({\bf x}))} \quad \text{where} \quad f({\bf x}) = {\bf w}^T{\boldsymbol \phi}.
$$

Probably easier if we sketch this function out, so let's do that

<Image
	src={
		"/images/courses/foundations-in-machine-learning/from-regression-to-classification-logistic-regression/img1.jpg"
	}
	alt={"academy.digilab.co.uk"}
	width={1200}
	height={609}
	caption={""}
/>

The pointe here is that the input to the sigmoid or logit function is the output of a linear model itself $(f{\bf x}) = {\bf w}^T{\boldsymbol \phi}$. As we will see this gives great flexibility.

To give an interpretation of this, we can rearrange the equations in terms of $f({\bf x})$, the linear model, so we get

$$
f({\bf x}) = \ln\left( \frac{p_A}{1-p_A}\right) = \ln\left( \frac{p_A}{p_B}\right)= w_0 + w_1\phi_1({\bf x}) + \ldots + w_M\phi_M({\bf x}) = {\bf w}^T{\boldsymbol \phi}.
$$

So $p_A/(1-p_A) = p_A/p_B$ is the 'odds', the ratio of label A and not A (i.e. B). So the interpretation of a logistic model is building a linear model for the 'log odds'.

Before moving on how we train logistic regression model. Let us sketch some toy examples.

So here we look a simple linear model $f({\bf x}) = w_0 + w_1x$, this function is then push through the sigmoid function, generating the probability of $p_A$ against input values $x$. In this example we decision boundary is at $x = 1/2$.

<Image
	src={
		"/images/courses/foundations-in-machine-learning/from-regression-to-classification-logistic-regression/img2.jpg"
	}
	alt={"academy.digilab.co.uk"}
	width={1200}
	height={457}
	caption={""}
/>

As discussed, a linear model doesn't mean a linear function. So with logistic regression models there is the possibility of great flexibility. Here is an example where $f$ is a quadratic function, resulting in a more complex decision boundary

<Image
	src={
		"/images/courses/foundations-in-machine-learning/from-regression-to-classification-logistic-regression/img3.jpg"
	}
	alt={"academy.digilab.co.uk"}
	width={1200}
	height={501}
	caption={""}
/>

<Callout type="info" heading={null}>
Before looking at training a logistic model, we note that logistic regression is an general extension of ordinary linear models. There are two ingredients

- Linear Predictor just as in a standard linear model $f = {\bf w}^T{\boldsymbol \phi}$.
- A **link function**, in our case the sigmoid function, but could be more general. The link function provides the relationship between the linear prediction and the mean of the distribution function to be modelled.

</Callout>

## Training a Logistic Model

For a Logistic Regression problem we can use a categorical cross-entropy loss, which is given by

$$
\mathcal L = - \frac{1}{N}\sum_{j=1}^N \Big[y_{j}\log(p_{j}) + (1 - y_{j})\log(1 - p_{j}) \Big].
$$

The optimal weights can then be found via gradient based optimisation scheme (e.g. steepest descent, Newton or Quasi-Newton).

<Callout type="info" heading={null}>
The gradients can be calculated using chain rule, we don't do the full calculation here but

$$
\frac{d \mathcal L}{d {\bf w}} = \frac{d \mathcal L}{d f} \cdot \frac{d f}{d {\bf w}} = \frac{d \mathcal L}{d f} \cdot\boldsymbol{\phi}({\bf x})
$$

We note that models are linear with respect to their weights, so if we differentiate $f$ with respect to ${\bf w}$ we such get the basis functions ${\boldsymbol \phi}({\bf x})$.

This leaves us to calculate $\partial L / \partial f$. The loss is just the loss $\ell_j$ for each samples added up. So we have

$$
\begin{align}
\frac{d \mathcal L}{d {\bf w}} &= \frac{d \mathcal L}{d f} \cdot\boldsymbol{\phi}({\bf x}) = \Big[\frac{1}{N}\sum_{j=1}^N \frac{\partial \ell_j}{\partial f}\Big]\cdot\boldsymbol{\phi}({\bf x}) \\&= \frac{1}{N}\sum*{j=1}^N \frac{\partial}{\partial f}\Big[ y*{j}\log(p*{j}) + (1 - y*{j})\log(1 - p*{j})\Big]\cdot\boldsymbol{\phi}({\bf x}) \\
&= \frac{1}{N}\sum*{j=1}^N \Big[y_j - p_j\Big] \cdot\boldsymbol{\phi}({\bf x})
\end{align}
$$

The last part set requires a bit of manipulation, but remember that

$$
p_j = \frac{1}{1 + \exp(-f({\bf x}_j))}
$$

</Callout>

## Importance of Regularisation

Logistic regression is a convex optimization problem (the likelihood function is concave), and it's known to not have a finite solution when it can fully separate the data, so the loss function can only reach its lowest value asymptomatically as the weights tend to $\pm \infty$.

Data is fully seperable if we have limited data relative to the flexibility in the model. Which defines the classic trade off between over and underfitting.

This has the effect of tightening decision boundaries around each data point when the data is separable, and the linear model is sufficiently expressive, with this asymptotically overfitting on the training set.

<Image
	src={
		"/images/courses/foundations-in-machine-learning/from-regression-to-classification-logistic-regression/img4.jpg"
	}
	alt={"academy.digilab.co.uk"}
	width={900}
	height={607}
	caption={""}
/>

Without regularization, the asymptotic nature of logistic regression would keep driving the loss towards 0 in high dimensions. There are therefore two well used strageties.

- $L^2$ (or Ridge) Regularisation, where an additional term $\lambda \|{\bf w}\|^2$ is added to the loss.

- Early Stopping.

We deal with regularisation as a seperate topic. For now `sklearn` automatically applies ridge regularisation ('l2') by default with $\lambda$ set to 1. So you now know what this this does, and understand that $\lambda$ the regularisation parameter is actually a hyperparameter which you should also optimise over during training.

<Image
	src={
		"/images/courses/foundations-in-machine-learning/from-regression-to-classification-logistic-regression/img5.jpg"
	}
	alt={"academy.digilab.co.uk"}
	width={1200}
	height={893}
	caption={""}
/>

Here is a snapshot from the class documentation for logistic regression. Whilst there is no need to write your own code for doing it, since `sklearn`'s implementation is good, it is important you understand the meaning of the default assumptions. Fitting a good logistic model, will often require tuning of the regularisation parameter. Here you will see they allow a user to set $C = 1/\lambda$ the inverse of the regularisation strength.

## Multi-class Classification

In what we have discussed so far we have considered binary classification. This is where there are a choice of only two classes, labels or outcomes. In general classification can map to any number of classes, which is referred to as multi-class classification.

Multi-class classification can be achieve via a simple extension of binary classification, described above, by following an approach called **One-vs-all** or **One-vs_rest**.

<Image
	src={
		"/images/courses/foundations-in-machine-learning/from-regression-to-classification-logistic-regression/img6.jpg"
	}
	alt={"academy.digilab.co.uk"}
	width={800}
	height={851}
	caption={""}
/>

Suppose we now have $M$ classes, logistic regression is classifiered for each class $p_i$, from $i = 1$ to $M$.

We can build a model which predicts the probability of each of the classes seperate, which is a binary classifier. I.e. is it class $i$ or is it not class $i$.

The classifier with the highest probability wins

$$
\text{argmax}_i \Big[p_i({\bf x})\Big]
$$

## Conclusion

Logistic regression is a extension of linear models to a classification problem. It is part of a broader class of models which are called "Generalised Linear Models". These are the composition of two models

- a linear predictor (here we define as $f = {\bf w}^T{\boldsymbol \phi}$)
- link function - in this case the logit or sigmoid function.

This second function plays the role of squashing the function between $0$ and $1$, transforming the output of the model to probability, which can be used for classification.

Regularisation is an important part of classification models, and when you use packages there are often chosen default parameters which play a central role in the outputs or quality of the results.

Finally, binary classification (two classes / labels) can be extended easily to multi-class classification using the so called "one-vs-all" strategy. This is nothing more than building multiple binary classifers, one for each class.


From Regression to Classification - Logistic Regression

## True of False - Logistic Regression in the Wild

<Callout type="success" heading={null}>
	In this walk through we are going to look at logistic regression, an example of
	generalised linear model which enables us to achieve classification in a
	regression framework.
</Callout>

### Import Libraries

```python
import numpy as np
import pandas as pd
```

### Pima Indians Diabetes Dataset

Here we load a [dataset from kaggle](https://www.kaggle.com/datasets/uciml/pima-indians-diabetes-database/download?datasetVersionNumber=1)

This is a data set with data about women, various measurements and the outcome of weather they had diabetes or not. We are going to build a classification model using logistic regression to try and predict the result of these tests.

```python
df  = pd.read_csv('diabetes.csv', sep=',')
df.head()
```

<Image
	src={
		"/images/courses/foundations-in-machine-learning/logistic-regression-in-wild/img1.jpg"
	}
	alt={"academy.digilab.co.uk"}
	width={1000}
	height={172}
	caption={""}
/>

We can see there are 9 coloumns, and in total 768 entries, we can see this by calling the following

```python
df.info()
```

<Image
	src={
		"/images/courses/foundations-in-machine-learning/logistic-regression-in-wild/img2.jpg"
	}
	alt={"academy.digilab.co.uk"}
	width={600}
	height={451}
	caption={""}
/>

Ok, to let us pull out some feature (or input) coloumns we are interested in. and then assign `y` to the outcome!

```python
feature_cols = ['Pregnancies', 'Glucose', 'BloodPressure', 'SkinThickness','Insulin','BMI','DiabetesPedigreeFunction', 'Age']
X = df[feature_cols]
y = df.Outcome
```

As normal we will split out data into training data and testing data

```python
from sklearn.model_selection import train_test_split

X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.25, random_state= 16)
```

Now we are in a good position to fit our logistic model, and then use the train model to predict our testing data

```python
from sklearn.linear_model import LogisticRegression

logreg = LogisticRegression(solver='liblinear')

logreg.fit(X_train, y_train)

y_pred = logreg.predict(X_test)
```

So how did we do? Let us plot a confusion matrix to see how the classification did!

```python
import matplotlib.pyplot as plt
import seaborn as sns

from sklearn import metrics

cnf_matrix = metrics.confusion_matrix(y_test, y_pred)

class_names=[0,1] # name  of classes
fig, ax = plt.subplots()
tick_marks = np.arange(len(class_names))
plt.xticks(tick_marks, class_names)
plt.yticks(tick_marks, class_names)
# create heatmap
sns.heatmap(pd.DataFrame(cnf_matrix), annot=True, cmap="YlGnBu" ,fmt='g')
ax.xaxis.set_label_position("top")
plt.tight_layout()
plt.title('Confusion matrix', y=1.1)
plt.ylabel('Actual label')
plt.xlabel('Predicted label')
```

<Image
	src={
		"/images/courses/foundations-in-machine-learning/logistic-regression-in-wild/img3.jpg"
	}
	alt={"academy.digilab.co.uk"}
	width={620}
	height={535}
	caption={""}
/>

Next we can look at a ROC_Curve as well

```python
y_pred_proba = logreg.predict_proba(X_test)[::,1]

fpr, tpr, _ = metrics.roc_curve(y_test,  y_pred_proba)

auc = np.around(metrics.roc_auc_score(y_test, y_pred_proba), decimals = 2)
plt.plot(fpr,tpr,label="data 1, auc="+str(auc))
plt.plot([0.0, 1.0], [0.0, 1.0], '-k', alpha=0.4, label="Random Classifier")
plt.xlim([-0.01, 1.0])
plt.ylim([0.0, 1.01])
plt.xlabel("False Positive Rate")
plt.ylabel("True Positive Rate")
plt.legend(loc=4)
plt.show()
```

<Image
	src={
		"/images/courses/foundations-in-machine-learning/logistic-regression-in-wild/img4.jpg"
	}
	alt={"academy.digilab.co.uk"}
	width={1200}
	height={884}
	caption={""}
/>


Logistic Regression in Wild

<Callout type="success" heading={null}>
In this explainer we look out our first unsupervised learning algorithm which is used for dimension reduction, called Principal Component Analysis (PCA). By the end of this sections you will

1. Understand the key principles of PCA, what it is and why would it be used.
2. Understand the key concepts of what is actually calculated when PCA is run.
3. Introduce the key concepts of explained variance (ratio) for choosing the amount of model reduction.
4. (Optional) Understand why the principal components are eigenvectors of the covariance matrix of your data set.
5. Walkthrough an example of PCA on a real world data set.

</Callout>

**Principal Component Analysis** (PCA) is a linear dimesion reduction method. The method is unsupervised, and so we just have a set of data

$$
{\bf X} := \{ {\bf x}_1, {\bf x}_2, \ldots, {\bf x}_N \}
$$

and no labels. Each of the samples are vector representations of the input, so that ${\bf x}_j \in \mathbb R^D$, this means the sample is represented by D-numbers, or put another way is D-dimensional input / sample space. Often input data can be very high dimensional, i.e $D$ is very large.

<Callout type="info" heading={null}>
	**A good example** of this, which often encounters are images. Take for an
	example a picture made up of 640 by 480 pixels, in which each pixel is
	represented by 3 color values - Red, Green, Blue (RGB). This gives a total
	number of dimensions of $640 \times 480 \times 3 = 921,600$.
</Callout>
High dimensional data is often challenging to analyse, interpret, visualise and store.
Whilst data might be represented in a high-dimensional way, the data is often much
more structured and correlation than this.

**Dimension Reduction** methods seek to exploit this structure and correlations to give a compact, lower dimensional representation of the data, without losing information.

<Image
	src={
		"/images/courses/foundations-in-machine-learning/looking-through-the-right-lense-principle-component-analysis/img1.jpg"
	}
	alt={"academy.digilab.co.uk"}
	width={600}
	height={461}
	caption={""}
/>

We are interested in finding a compact representation of the data ${\bf z}_j \in \mathbb R^M$ where $M$ is much less $N$ the original dimension of ${\bf x}_j$. This dimension reduction is achieved by a matrix operation ${\bf V}^T$, such that

$$
{\bf z}_j = {\bf V}^T{\bf x}_j
$$

where ${\bf V} := [ {\bf v}_1, \ldots, {\bf v}_M  ] \in \mathbb R^{D \times M}$, and the the reconstruction given by the matrix operation

$$
\tilde{{\bf x}}_j = {\bf V}{\bf z}_j.
$$

To visualise these matrix operations we have

<Image
	src={
		"/images/courses/foundations-in-machine-learning/looking-through-the-right-lense-principle-component-analysis/img2.jpg"
	}
	alt={"academy.digilab.co.uk"}
	width={1200}
	height={445}
	caption={""}
/>

PCA is a linear dimension reduction method, since components of the new representation ${\bf z}$ are linear combinations of full representation ${\bf x}$.
The question is now how do we pick the columns ${\bf v}_j$. We choose two:

- Choose directions ${\bf v}_j$ so capture maximum variance of the original data set $X$.
- Choose coloumns ${\bf v}_j$ and ${\bf v}_i$ so they are ortogonal ${\bf v}_j \cdot {\bf v}_i = 0$.

To visualise what this looks like on a simple data set, we see two principal components. The principal component pointing the direction of the greatest variance, and the second component orthogonal (i.e. at 90 degrees to the first).

<Image
	src={
		"/images/courses/foundations-in-machine-learning/looking-through-the-right-lense-principle-component-analysis/img3.jpg"
	}
	alt={"academy.digilab.co.uk"}
	width={600}
	height={496}
	caption={""}
/>

The principal components ${\bf v}_j$ are actually eigenvectors of the covariance matrix

$$
{\bf S} = \frac{1}{N}\sum_{j=1}^N{\bf x}_{j}{\bf x}_j^{T},
$$

i.e. they are solutions of the equation

$$
{\bf S}{\bf v} = \lambda {\bf v}.
$$

The matrix ${\bf S}$ is a symmetric matrix, with $D$ eigenvectors and associated eigenvalues, such that

$$
\lambda_1 > \lambda_2 > \ldots > \lambda_D.
$$

The eigenvalues give us a value of the contribution towards the total variance an eigenvector makes, this is also called the explained variance

$$
\alpha_i = \frac{\lambda_i}{\sum_{i=1}^D \lambda_i}.
$$

Plotting $\alpha_i$ against eigenvalues, shows a natural cutoff of how many principal modes should be include $M$ in the reduce order model, without losing to much information.

To show a conceptal plot of what this looks like, see below.

<Image
	src={
		"/images/courses/foundations-in-machine-learning/looking-through-the-right-lense-principle-component-analysis/img4.jpg"
	}
	alt={"academy.digilab.co.uk"}
	width={875}
	height={591}
	caption={""}
/>

The plot shows individual contributions of the principal components to explained variance (as a percentage). The orange curve shows the culmative effect. So what we mean by this is that if you include the first three components you will include 90% of the variance from the original data set.

We will show a practical example of in the python example walkthrough for PCA.

## PCA and General Linear Models

In linear models, we extended the concept away from simple linear functions e.g. $f(x) = w_0 + w_1x$, to more general representations of linear combinations of general functions, so that

$$
f(x) = w_0 + w_1\phi_1(x) + \ldots w_2\phi_M(x) = \sum_{i=1}^M w_i\phi_i(x)
$$

where we call $\phi_i(x)$ the **basis functions**.

So an interest point comes that how do you choose the basis functions. You could engineer them from experience. So, for example, if you were modelling a system that was dependent on a tide for example, you could easily put in some sensible periodic functions perhaps.

In general you might not know what to choose, and for this the principal components become a natural choice. So this means we take

$$
\phi_i({\bf x}) = v_i^T{\bf x} = z_i
$$

so this would be a sum of just linear features, but could be generally extended to higher order polynomials by taking

$$
\phi_{(k-1)M + i}({\bf x}) = z_i^k.
$$

## Concluding Remarks

PCA is an unsupervised learning algorithm. It is a parameteric model, since it imposes a linear project / transformation of the data only.

It is very useful method which is widely used in the data preparation, input representation and explorations stages of a machine learning workflow. As an AI practioner it tells me how correlated my inputs / features are that I am going to use.

We see that in many machine learning workflows that algorithms are either no unsupervised or supervised, but combine the good bits of both! For example, we have show how we can use PCA to inform the basis functions we use in a General Linear Models.

Bottomline, PCA is a very useful tool to have in your toolbox. Implementation is packages like `sklearn` mean that implementations are just a few lines of code.

## Bonus Information

<Callout type="info" heading={null}>
**A Bit of Maths - Not required for using PCA.** There are packages that implementing PCA, so you don't worry about how to calculate the eigenvectors and values of ${\bf S}$, if you don't want to. This becomes more complex when data sets or dimensions are large.
For those interested in the question "Why are the principal components the eigenvectors of ${\bf S}$?". . . Here it is.
</Callout>
We wish to maximise variance of the first reduce coordinate $z_1$ of ${\bf z}$ which is

$$
V_1 = \frac{1}{N}\sum_{j=1}^N z_{1j}^2
$$

is maximised. Noting that $z_{1j} = {\bf v}_1^T{\bf x}_{j}$. This is to project of sample ${\bf x}_{j}$ into the direction ${\bf v}_1$. So we can rewrite the variance

$$
V_1 = \frac{1}{N}\sum_{j=1}^N({\bf v}_1^T{\bf x}_{j})^2 = \frac{1}{N}\sum_{j=1}^N {\bf v}_1^T{\bf x}_{j}{\bf x}_j^{T}{\bf v}_1 =  {\bf v}_1^T\left(\frac{1}{N}\sum_{j=1}^N{\bf x}_{j}{\bf x}_j^{T}\right){\bf v}_1 = {\bf v}_1^T{\bf S}{\bf v}_1.
$$

We wish find ${\bf v}_1$ which maximises this variance, but we note that we can just make the magnitude of ${\bf v}_1$ arbitarily large, so this problem is ill-posed. We therefore have to solve a constrained maximisation problem, such that

$$
\max_{{\bf v}_1} \left[{\bf v}_1^T{\bf S}{\bf v}_1 \right]\quad \text{such that} \quad \|{\bf v}_1\|^2 = {\bf v}_1^T {\bf v}_1 = 1.
$$

This then becomes a constrained optimisation problem. So we write down the Lagrangian

$$
\mathcal L = {\bf v}_1^T{\bf S}{\bf v}_1 + \lambda_1 \left(1 - {\bf v}_1^T {\bf v}_1 \right).
$$

We can then find maximisers of $\lambda$ by differentiating with respect to ${\bf v}_1$ and $\lambda_1$, so that

$$
\frac{\partial\mathcal L}{\partial{\bf v}} = 2{\bf v}_1^T{\bf S} - 2\lambda_1{\bf v}_1^T \quad \text{and} \quad \frac{\partial\mathcal L}{\partial\lambda} = 1 - {\bf v}_1^T {\bf v}_1
$$

Setting these partial derivatives to zero which gives us our eigenproblem

$$
{\bf S}{\bf v}_1 = \lambda_1 {\bf v}_1 \quad \text{and} \quad {\bf v}_1^T{\bf v}_1 = 1.
$$

subsequent eigenvalues will then give us maximisers of $\mathcal L$ orthgonal to $v_{1}$. Hence we get the order sequence of eigenvectors (principal components).

We don't show it hear but is possible to show that maximising the variance of the projected data is the same as minimising the reconstruction error over the data set, i.e.

$$
\mathcal J = \frac{1}{N}\sum_{j=1}^N\|{\bf x}_j - \tilde{{\bf x}}_j\|^2 = \frac{1}{N}\sum_{j=1}^N\|{\bf x}_j - {\bf V}{\bf z}_j\|^2
$$


Looking through the Right Lens - Principle Component Analysis

In this quick explainer we are walkthrough our first simple example of Principal Component Analysis.

We will see how to compute the principal components, and then how to transform our data to be represented in this new space, and then build a linear model in this new feature space.

## Add libaries

```python
import matplotlib.pyplot as plt
import seaborn as sns
import pandas as pd
import numpy as np
```

The data we will use in this example is some red wine quality data. Lots of different measurements, and then measure of their quality out of 10 as the target value.

In total we have data for 1599 wines. So let's load that data set and pull out some random samples

```python
df = pd.read_csv('winequality-red.csv', sep=";")
df.sample(5)
```

<Image
	src={
		"/images/courses/foundations-in-machine-learning/looking-in-the-right-direction-exploring-pca/img1.jpg"
	}
	alt={"academy.digilab.co.uk"}
	width={1000}
	height={235}
	caption={""}
/>

So now we pull out our target variable quality, and then do the standard train / test split.

```python
from sklearn.model_selection import train_test_split

target = 'quality'
X = df.drop(target,axis=1)
y = df[target]

X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42)
```

## Re-scaling your data

Again like most data science work flows rescaling data is an important step. If we don't do this variables which happen to have higher values (because of what they measure) will bias the results. We therefore rescale all input parameters between a range of 0 and 1.

```python
from sklearn.preprocessing import scale

X_train_scaled = scale(X_train)
X_test_scaled = scale(X_test)
```

## Calculate Principal Components

Calling PCA in sklearn is nice and easy.

```python
from sklearn.decomposition import PCA

pca = PCA() # Default n_components = min(n_samples, n_features)
X_train_pc = pca.fit_transform(X_train_scaled)
X_test_pc = pca.transform(X_test_scaled)
```

Once we have fitted the pca transformation, behind the scenes sklearn will have computed the explained_variance ratio. We can then also calculate the accumulated explained variance using `np.cumsum()`.

```python
explained_var = pca.explained_variance_ratio_

accum_var = np.cumsum(explained_var)

Once calculated we can make a nice two y-axis plot.

fig, ax1 = plt.subplots()

ax2 = ax1.twinx()

ax1.plot(np.arange(0, explained_var.shape[0]), explained_var, '-ob', alpha = 0.5)

ax2.plot(np.arange(0, explained_var.shape[0]), accum_var, '-o',c='lightcoral', alpha = 0.5)

ax1.set_xlabel('Principal Component')
ax1.set_ylabel('Explained Variance', color='b')
ax2.set_ylabel('Accumulated Variance', color='lightcoral')
```

<Image
	src={
		"/images/courses/foundations-in-machine-learning/looking-in-the-right-direction-exploring-pca/img2.jpg"
	}
	alt={"academy.digilab.co.uk"}
	width={1000}
	height={695}
	caption={""}
/>

## Looking at High Dimension data :)

```python
fig, ax = plt.subplots(figsize=(6, 4))

unique_label = np.unique(y)
num_labels = len(unique_labels)

plt.scatter(X_train_pc[:,0], X_train_pc[:,1], c=y_train, alpha = 0.7, cmap='coolwarm')

cbar = plt.colorbar()
cbar.solids.set_edgecolor("face")


ax.legend()
ax.set_xlabel('Principal Component 0')
ax.set_ylabel('Principal Component 1')
```

<Image
	src={
		"/images/courses/foundations-in-machine-learning/looking-in-the-right-direction-exploring-pca/img3.jpg"
	}
	alt={"academy.digilab.co.uk"}
	width={1000}
	height={750}
	caption={""}
/>

## Build a Predictive Model

So here we are going to build a predictive model using Gradient Boosted Classifier. These are generally the best simplest multi-class classifiers to use out the box. They are my go to.

```python
from sklearn.ensemble import GradientBoostingClassifier

perf = []
num_pca = []

for i in np.arange(1,11):

    clf = GradientBoostingClassifier(n_estimators=2000, learning_rate=1.0, max_depth=2, random_state=0).fit(X_train_pc[:,:i], y_train)

    num_pca.append(i)
    perf.append(clf.score(X_test_pc[:,:i], y_test))
```

Let's pull out another double axis plot.

```python
fig, ax1 = plt.subplots()

ax2 = ax1.twinx()


ax1.plot(num_pca, perf, '-ob', alpha = 0.5)

ax2.plot(num_pca, (1.0 - perf/perf[-1])*100., '-o',c='lightcoral', alpha = 0.5)


ax1.set_xlabel('Principal Component')
ax1.set_ylabel('Score', color='b')
ax2.set_ylabel('Loss in Performance', color='lightcoral')
```

<Image
	src={
		"/images/courses/foundations-in-machine-learning/looking-in-the-right-direction-exploring-pca/img4.jpg"
	}
	alt={"academy.digilab.co.uk"}
	width={1000}
	height={714}
	caption={""}
/>

"Score" is a harsh metric for performance for the algorithm. So probably more important look at the relative metrics. We see that even with just 2 principal components we can obtain model with a score with 10% of the full input model.


Using General Linear Models for Machine Learning

Course Description

Lessons

1. Section Overview

2. The Foundations of ML - Curve Fitting

3. Regression Walkthrough

4. Underfitting and Overfitting

5. Controlling ML Models - Regularisation in Practise

6. Exploring Simple Relationships - Correlations

7. Finding Nonlinear Relations - Predictive Power Score

8. Correlation & PPS Walkthrough

9. From Regression to Classification - Logistic Regression

10. Logistic Regression in Wild

11. Looking through the Right Lens - Principle Component Analysis

12. Looking in the Right Direction - Exploring PCA