
# Bayesian Methods

Data Science Africa, Nyeri, Kenya

### What is Machine Learning?

$\text{data} + \text{model} \xrightarrow{\text{compute}} \text{prediction}$

• data : observations, could be actively or passively acquired (meta-data).
• model : assumptions, based on previous experience (other data! transfer learning etc), or beliefs about the regularities of the universe. Inductive bias.
• prediction : an action to be taken or a categorization or a quality score.

### What is Machine Learning?

$\text{data} + \text{model} \xrightarrow{\text{compute}} \text{prediction}$

• To combine data with a model need:
• a prediction function $$\mappingFunction(\cdot)$$ includes our beliefs about the regularities of the universe
• an objective function $$\errorFunction(\cdot)$$ defines the cost of misprediction.

### Olympic Marathon Data

 Gold medal times for Olympic Marathon since 1896. Marathons before 1924 didn’t have a standardised distance. Present results using pace per km. In 1904 Marathon was badly organised leading to very slow times. Image from Wikimedia Commons http://bit.ly/16kMKHQ

### Regression: Linear Releationship

$\dataScalar_i = m \inputScalar_i + c$

• $$\dataScalar_i$$ : winning pace.

• $$\inputScalar_i$$ : year of Olympics.

• $$m$$ : rate of improvement over time.

• $$c$$ : winning time at year 0.

### $$\dataScalar = m\inputScalar + c$$

point 1: $$\inputScalar = 1$$, $$\dataScalar=3$$ $3 = m + c$

point 2: $$\inputScalar = 3$$, $$\dataScalar=1$$ $1 = 3m + c$

point 3: $$\inputScalar = 2$$, $$\dataScalar=2.5$$

$2.5 = 2m + c$

### $$\dataScalar = m\inputScalar + c + \noiseScalar$$

point 1: $$\inputScalar = 1$$, $$\dataScalar=3$$ $3 = m + c + \noiseScalar_1$

point 2: $$\inputScalar = 3$$, $$\dataScalar=1$$ $1 = 3m + c + \noiseScalar_2$

point 3: $$\inputScalar = 2$$, $$\dataScalar=2.5$$ $2.5 = 2m + c + \noiseScalar_3$

### A Probabilistic Process

Set the mean of Gaussian to be a function. $p \left(\dataScalar_i|\inputScalar_i\right)=\frac{1}{\sqrt{2\pi\dataStd^2}}\exp \left(-\frac{\left(\dataScalar_i-\mappingFunction\left(\inputScalar_i\right)\right)^{2}}{2\dataStd^2}\right).$

This gives us a ‘noisy function’.

This is known as a stochastic process.

### The Gaussian Density

• Perhaps the most common probability density.

\begin{align} p(\dataScalar| \meanScalar, \dataStd^2) & = \frac{1}{\sqrt{2\pi\dataStd^2}}\exp\left(-\frac{(\dataScalar - \meanScalar)^2}{2\dataStd^2}\right)\\& \buildrel\triangle\over = \gaussianDist{\dataScalar}{\meanScalar}{\dataStd^2} \end{align}

### Gaussian Density

The Gaussian PDF with $${\meanScalar}=1.7$$ and variance $${\dataStd}^2=0.0225$$. Mean shown as cyan line. It could represent the heights of a population of students.

### Gaussian Density

$\gaussianDist{\dataScalar}{\meanScalar}{\dataStd^2} = \frac{1}{\sqrt{2\pi\dataStd^2}} \exp\left(-\frac{(\dataScalar-\meanScalar)^2}{2\dataStd^2}\right)$

$$\dataStd^2$$ is the variance of the density and $$\meanScalar$$ is the mean.

### Sum of Gaussians

Sum of Gaussian variables is also Gaussian.

$\dataScalar_i \sim \gaussianSamp{\meanScalar_i}{\sigma_i^2}$

And the sum is distributed as

$\sum_{i=1}^{\numData} \dataScalar_i \sim \gaussianSamp{\sum_{i=1}^\numData \meanScalar_i}{\sum_{i=1}^\numData \sigma_i^2}$

(Aside: As sum increases, sum of non-Gaussian, finite variance variables is also Gaussian because of central limit theorem.)

### Scaling a Gaussian

Scaling a Gaussian leads to a Gaussian.

$\dataScalar \sim \gaussianSamp{\meanScalar}{\sigma^2}$

And the scaled variable is distributed as

$\mappingScalar \dataScalar \sim \gaussianSamp{\mappingScalar\meanScalar}{\mappingScalar^2 \sigma^2}.$

### A Probabilistic Process

Set the mean of Gaussian to be a function.

$p\left(\dataScalar_i|\inputScalar_i\right)=\frac{1}{\sqrt{2\pi\dataStd^2}}\exp\left(-\frac{\left(\dataScalar_i-f\left(\inputScalar_i\right)\right)^{2}}{2\dataStd^2}\right).$

This gives us a ‘noisy function’.

This is known as a stochastic process.

### Height as a Function of Weight

In the standard Gaussian, parametized by mean and variance.

Make the mean a linear function of an input.

This leads to a regression model. \begin{align*} \dataScalar_i=&\mappingFunction\left(\inputScalar_i\right)+\noiseScalar_i,\\ \noiseScalar_i \sim & \gaussianSamp{0}{\dataStd^2}. \end{align*}

Assume $$\dataScalar_i$$ is height and $$\inputScalar_i$$ is weight.

### Data Point Likelihood

Likelihood of an individual data point $p\left(\dataScalar_i|\inputScalar_i,m,c\right)=\frac{1}{\sqrt{2\pi \dataStd^2}}\exp\left(-\frac{\left(\dataScalar_i-m\inputScalar_i-c\right)^{2}}{2\dataStd^2}\right).$ Parameters are gradient, $$m$$, offset, $$c$$ of the function and noise variance $$\dataStd^2$$.

### Data Set Likelihood

If the noise, $$\epsilon_i$$ is sampled independently for each data point. Each data point is independent (given $$m$$ and $$c$$). For independent variables: $p(\dataVector) = \prod_{i=1}^\numData p(\dataScalar_i)$ $p(\dataVector|\inputVector, m, c) = \prod_{i=1}^\numData p(\dataScalar_i|\inputScalar_i, m, c)$

### For Gaussian

i.i.d. assumption $p(\dataVector|\inputVector, m, c) = \prod_{i=1}^\numData \frac{1}{\sqrt{2\pi \dataStd^2}}\exp \left(-\frac{\left(\dataScalar_i- m\inputScalar_i-c\right)^{2}}{2\dataStd^2}\right).$ $p(\dataVector|\inputVector, m, c) = \frac{1}{\left(2\pi \dataStd^2\right)^{\frac{\numData}{2}}}\exp\left(-\frac{\sum_{i=1}^\numData\left(\dataScalar_i-m\inputScalar_i-c\right)^{2}}{2\dataStd^2}\right).$

### Log Likelihood Function

• Normally work with the log likelihood: $L(m,c,\dataStd^{2})=-\frac{\numData}{2}\log 2\pi -\frac{\numData}{2}\log \dataStd^2 -\sum_{i=1}^{\numData}\frac{\left(\dataScalar_i-m\inputScalar_i-c\right)^{2}}{2\dataStd^2}.$

### Consistency of Maximum Likelihood

• If data was really generated according to probability we specified.
• Correct parameters will be recovered in limit as $$\numData \rightarrow \infty$$.
• This can be proven through sample based approximations (law of large numbers) of “KL divergences”.
• Mainstay of classical statistics.

### Probabilistic Interpretation of the Error Function

• Probabilistic Interpretation for Error Function is Negative Log Likelihood.
• Minimizing error function is equivalent to maximizing log likelihood.
• Maximizing log likelihood is equivalent to maximizing the likelihood because $$\log$$ is monotonic.
• Probabilistic interpretation: Minimizing error function is equivalent to maximum likelihood with respect to parameters.

### Error Function

• Negative log likelihood is the error function leading to an error function $\errorFunction(m,c,\dataStd^{2})=\frac{\numData}{2}\log \dataStd^2+\frac{1}{2\dataStd^2}\sum _{i=1}^{\numData}\left(\dataScalar_i-m\inputScalar_i-c\right)^{2}.$
• Learning proceeds by minimizing this error function for the data set provided.

### Connection: Sum of Squares Error

• Ignoring terms which don’t depend on $$m$$ and $$c$$ gives $\errorFunction(m, c) \propto \sum_{i=1}^\numData (\dataScalar_i - \mappingFunction(\inputScalar_i))^2$ where $$\mappingFunction(\inputScalar_i) = m\inputScalar_i + c$$.
• This is known as the sum of squares error function.
• Commonly used and is closely associated with the Gaussian likelihood.

### Reminder

• Two functions involved:
• Prediction function: $$\mappingFunction(\inputScalar_i)$$
• Error, or Objective function: $$\errorFunction(m, c)$$
• Error function depends on parameters through prediction function.

### Mathematical Interpretation

• What is the mathematical interpretation?
• There is a cost function.
• It expresses mismatch between your prediction and reality. $\errorFunction(m, c)=\sum_{i=1}^\numData \left(\dataScalar_i - m\inputScalar_i-c\right)^2$
• This is known as the sum of squares error.

### Learning is Optimization

• Learning is minimization of the cost function.
• At the minima the gradient is zero.
• Coordinate ascent, find gradient in each coordinate and set to zero. $\frac{\text{d}\errorFunction(c)}{\text{d}c} = -2\sum_{i=1}^\numData \left(\dataScalar_i- m \inputScalar_i - c \right)$ $0 = -2\sum_{i=1}^\numData\left(\dataScalar_i- m\inputScalar_i - c \right)$

### Learning is Optimization

• Fixed point equations $0 = -2\sum_{i=1}^\numData \dataScalar_i +2\sum_{i=1}^\numData m \inputScalar_i +2n c$ $c = \frac{\sum_{i=1}^\numData \left(\dataScalar_i - m\inputScalar_i\right)}{\numData}$

### Learning is Optimization

• Learning is minimization of the cost function.
• At the minima the gradient is zero.
• Coordinate ascent, find gradient in each coordinate and set to zero. $\frac{\text{d}\errorFunction(m)}{\text{d}m} = -2\sum_{i=1}^\numData \inputScalar_i\left(\dataScalar_i- m \inputScalar_i - c \right)$ $0 = -2\sum_{i=1}^\numData \inputScalar_i \left(\dataScalar_i-m \inputScalar_i - c \right)$

### Learning is Optimization

• Fixed point equations $0 = -2\sum_{i=1}^\numData \inputScalar_i\dataScalar_i+2\sum_{i=1}^\numData m \inputScalar_i^2+2\sum_{i=1}^\numData c\inputScalar_i$ $m = \frac{\sum_{i=1}^\numData \left(\dataScalar_i -c\right)\inputScalar_i}{\sum_{i=1}^\numData\inputScalar_i^2}$

$m^* = \frac{\sum_{i=1}^\numData (\dataScalar_i - c)\inputScalar_i}{\sum_{i=1}^\numData \inputScalar_i^2}$

### Fixed Point Updates

Worked example. \begin{aligned} c^{*}=&\frac{\sum _{i=1}^{\numData}\left(\dataScalar_i-m^{*}\inputScalar_i\right)}{\numData},\\ m^{*}=&\frac{\sum _{i=1}^{\numData}\inputScalar_i\left(\dataScalar_i-c^{*}\right)}{\sum _{i=1}^{\numData}\inputScalar_i^{2}},\\ \left.\dataStd^2\right.^{*}=&\frac{\sum _{i=1}^{\numData}\left(\dataScalar_i-m^{*}\inputScalar_i-c^{*}\right)^{2}}{\numData} \end{aligned}

### Important Concepts Not Covered

• Other optimization methods:
• Second order methods, conjugate gradient, quasi-Newton and Newton.
• Effective heuristics such as momentum.
• Local vs global solutions.

• Section 1.1-1.2 of Rogers and Girolami (2011) for fitting linear models.
• Section 1.2.5 of Bishop (2006) up to equation 1.65.

### Multi-dimensional Inputs

• Multivariate functions involve more than one input.
• Height might be a function of weight and gender.
• There could be other contributory factors.
• Place these factors in a feature vector $$\inputVector_i$$.
• Linear function is now defined as $\mappingFunction(\inputVector_i) = \sum_{j=1}^p w_j \inputScalar_{i, j} + c$

### Vector Notation

• Write in vector notation, $\mappingFunction(\inputVector_i) = \mappingVector^\top \inputVector_i + c$
• Can absorb $$c$$ into $$\mappingVector$$ by assuming extra input $$\inputScalar_0$$ which is always 1. $\mappingFunction(\inputVector_i) = \mappingVector^\top \inputVector_i$

### Objective Functions and Regression

• Classification: map feature to class label.
• Regression: map feature to real value our prediction function is

$\mappingFunction(\inputScalar_i) = m\inputScalar_i + c$

• Need an algorithm to fit it.

• Least squares: minimize an error.

$\errorFunction(m, c) = \sum_{i=1}^\numData (\dataScalar_i * \mappingFunction(\inputScalar_i))^2$

### Regression

• Create an artifical data set.

We now need to decide on a true value for $$m$$ and a true value for $$c$$ to use for generating the data.

We can use these values to create our artificial data. The formula $\dataScalar_i = m\inputScalar_i + c$ is translated to code as follows:

### Plot of Data

We can now plot the artifical data we’ve created.

These points lie exactly on a straight line, that’s not very realistic, let’s corrupt them with a bit of Gaussian ‘noise’.

### Contour Plot of Error Function

• Visualise the error function surface, create vectors of values.

• create a grid of values to evaluate the error function in 2D.

• compute the error function at each combination of $$c$$ and $$m$$.

### Contour Plot of Error

• We can now make a contour plot.

### Steepest Descent

• Minimize the sum of squares error function.
• One way of doing that is gradient descent.
• Initialize with a guess for $$m$$ and $$c$$
• update that guess by subtracting a portion of the gradient from the guess.
• Like walking down a hill in the steepest direction of the hill to get to the bottom.

### Algorithm

• We start with a guess for $$m$$ and $$c$$.

• Now we need to compute the gradient of the error function, firstly with respect to $$c$$,

$\frac{\text{d}\errorFunction(m, c)}{\text{d} c} = -2\sum_{i=1}^\numData (\dataScalar_i - m\inputScalar_i - c)$

• This is computed in python as follows

### Deriving the Gradient

To see how the gradient was derived, first note that the $$c$$ appears in every term in the sum. So we are just differentiating $$(\dataScalar_i - m\inputScalar_i - c)^2$$ for each term in the sum. The gradient of this term with respect to $$c$$ is simply the gradient of the outer quadratic, multiplied by the gradient with respect to $$c$$ of the part inside the quadratic. The gradient of a quadratic is two times the argument of the quadratic, and the gradient of the inside linear term is just minus one. This is true for all terms in the sum, so we are left with the sum in the gradient.

The gradient with respect tom $$m$$ is similar, but now the gradient of the quadratic’s argument is $$-\inputScalar_i$$ so the gradient with respect to $$m$$ is

$\frac{\text{d}\errorFunction(m, c)}{\text{d} m} = -2\sum_{i=1}^\numData \inputScalar_i(\dataScalar_i - m\inputScalar_i - c)$

which can be implemented in python (numpy) as

### Update Equations

• Now we have gradients with respect to $$m$$ and $$c$$.
• Can update our inital guesses for $$m$$ and $$c$$ using the gradient.
• We don’t want to just subtract the gradient from $$m$$ and $$c$$,
• We need to take a small step in the gradient direction.
• Otherwise we might overshoot the minimum.
• We want to follow the gradient to get to the minimum, the gradient changes all the time.

### Update Equations

• The step size has already been introduced, it’s again known as the learning rate and is denoted by $$\learnRate$$. $c_\text{new}\leftarrow c_{\text{old}} - \learnRate \frac{\text{d}\errorFunction(m, c)}{\text{d}c}$

• gives us an update for our estimate of $$c$$ (which in the code we’ve been calling c_star to represent a common way of writing a parameter estimate, $$c^*$$) and $m_\text{new} \leftarrow m_{\text{old}} - \learnRate \frac{\text{d}\errorFunction(m, c)}{\text{d}m}$
• Giving us an update for $$m$$.

### Update Code

• These updates can be coded as
• Fit model by descending gradient.

### Stochastic Gradient Descent

• If $$\numData$$ is small, gradient descent is fine.
• But sometimes (e.g. on the internet $$\numData$$ could be a billion.
• Stochastic gradient descent is more similar to perceptron.
• Look at gradient of one data point at a time rather than summing across all data points)
• This gives a stochastic estimate of gradient.

### Stochastic Gradient Descent

• The real gradient with respect to $$m$$ is given by

$\frac{\text{d}\errorFunction(m, c)}{\text{d} m} = -2\sum_{i=1}^\numData \inputScalar_i(\dataScalar_i - m\inputScalar_i - c)$

but it has $$\numData$$ terms in the sum. Substituting in the gradient we can see that the full update is of the form

$m_\text{new} \leftarrow m_\text{old} + 2\learnRate \left[\inputScalar_1 (\dataScalar_1 - m_\text{old}\inputScalar_1 - c_\text{old}) + (\inputScalar_2 (\dataScalar_2 - m_\text{old}\inputScalar_2 - c_\text{old}) + \dots + (\inputScalar_n (\dataScalar_n - m_\text{old}\inputScalar_n - c_\text{old})\right]$

This could be split up into lots of individual updates $m_1 \leftarrow m_\text{old} + 2\learnRate \left[\inputScalar_1 (\dataScalar_1 - m_\text{old}\inputScalar_1 - c_\text{old})\right]$ $m_2 \leftarrow m_1 + 2\learnRate \left[\inputScalar_2 (\dataScalar_2 - m_\text{old}\inputScalar_2 - c_\text{old})\right]$ $m_3 \leftarrow m_2 + 2\learnRate \left[\dots\right]$ $m_n \leftarrow m_{n-1} + 2\learnRate \left[\inputScalar_n (\dataScalar_n - m_\text{old}\inputScalar_n - c_\text{old})\right]$

which would lead to the same final update.

### Updating $$c$$ and $$m$$

• In the sum we don’t $$m$$ and $$c$$ we use for computing the gradient term at each update.
• In stochastic gradient descent we do change them.
• This means it’s not quite the same as steepest desceint.
• But we can present each data point in a random order, like we did for the perceptron.
• This makes the algorithm suitable for large scale web use (recently this domain is know as ‘Big Data’) and algorithms like this are widely used by Google, Microsoft, Amazon, Twitter and Facebook.

### Stochastic Gradient Descent

• Or more accurate, since the data is normally presented in a random order we just can write $m_\text{new} = m_\text{old} + 2\learnRate\left[\inputScalar_i (\dataScalar_i - m_\text{old}\inputScalar_i - c_\text{old})\right]$

### SGD for Linear Regression

Putting it all together in an algorithm, we can do stochastic gradient descent for our regression data.

### Reflection on Linear Regression and Supervised Learning

1. What effect does the learning rate have in the optimization? What’s the effect of making it too small, what’s the effect of making it too big? Do you get the same result for both stochastic and steepest gradient descent?

2. The stochastic gradient descent doesn’t help very much for such a small data set. It’s real advantage comes when there are many, you’ll see this in the lab.

### Log Likelihood for Multivariate Regression

The likelihood of a single data point is

$p\left(\dataScalar_i|\inputScalar_i\right)=\frac{1}{\sqrt{2\pi\dataStd^2}}\exp\left(-\frac{\left(\dataScalar_i-\mappingVector^{\top}\inputVector_i\right)^{2}}{2\dataStd^2}\right).$

Leading to a log likelihood for the data set of

$L(\mappingVector,\dataStd^2)= -\frac{\numData}{2}\log \dataStd^2-\frac{\numData}{2}\log 2\pi -\frac{\sum_{i=1}^{\numData}\left(\dataScalar_i-\mappingVector^{\top}\inputVector_i\right)^{2}}{2\dataStd^2}.$

### Error Function

And a corresponding error function of $\errorFunction(\mappingVector,\dataStd^2)=\frac{\numData}{2}\log\dataStd^2 + \frac{\sum_{i=1}^{\numData}\left(\dataScalar_i-\mappingVector^{\top}\inputVector_i\right)^{2}}{2\dataStd^2}.$

### Expand the Brackets

\begin{align*} \errorFunction(\mappingVector,\dataStd^2) = & \frac{\numData}{2}\log \dataStd^2 + \frac{1}{2\dataStd^2}\sum _{i=1}^{\numData}\dataScalar_i^{2}-\frac{1}{\dataStd^2}\sum _{i=1}^{\numData}\dataScalar_i\mappingVector^{\top}\inputVector_i\\&+\frac{1}{2\dataStd^2}\sum _{i=1}^{\numData}\mappingVector^{\top}\inputVector_i\inputVector_i^{\top}\mappingVector +\text{const}.\\ = & \frac{\numData}{2}\log \dataStd^2 + \frac{1}{2\dataStd^2}\sum _{i=1}^{\numData}\dataScalar_i^{2}-\frac{1}{\dataStd^2} \mappingVector^\top\sum_{i=1}^{\numData}\inputVector_i\dataScalar_i\\&+\frac{1}{2\dataStd^2} \mappingVector^{\top}\left[\sum _{i=1}^{\numData}\inputVector_i\inputVector_i^{\top}\right]\mappingVector +\text{const}. \end{align*}

### Multivariate Derivatives

• We will need some multivariate calculus.
• For now some simple multivariate differentiation: $\frac{\text{d}{\mathbf{a}^{\top}}{\mappingVector}}{\text{d}\mappingVector}=\mathbf{a}$ and $\frac{\mappingVector^{\top}\mathbf{A}\mappingVector}{\text{d}\mappingVector}=\left(\mathbf{A}+\mathbf{A}^{\top}\right)\mappingVector$ or if $$\mathbf{A}$$ is symmetric (i.e. $$\mathbf{A}=\mathbf{A}^{\top}$$) $\frac{\text{d}\mappingVector^{\top}\mathbf{A}\mappingVector}{\text{d}\mappingVector}=2\mathbf{A}\mappingVector.$

### Differentiate the Objective

Differentiating with respect to the vector $$\mappingVector$$ we obtain $\frac{\partial L\left(\mappingVector,\dataStd^2 \right)}{\partial \mappingVector}=\frac{1}{\dataStd^2} \sum _{i=1}^{\numData}\inputVector_i \dataScalar_i-\frac{1}{\dataStd^2} \left[\sum _{i=1}^{\numData}\inputVector_i\inputVector_i^{\top}\right]\mappingVector$ Leading to $\mappingVector^{*}=\left[\sum _{i=1}^{\numData}\inputVector_i\inputVector_i^{\top}\right]^{-1}\sum _{i=1}^{\numData}\inputVector_i\dataScalar_i,$

### Differentiate the Objective

Rewrite in matrix notation: $\sum_{i=1}^{\numData}\inputVector_i\inputVector_i^\top = \inputMatrix^\top \inputMatrix$ $\sum_{i=1}^{\numData}\inputVector_i\dataScalar_i = \inputMatrix^\top \dataVector$

### Update Equations

• Update for $$\mappingVector^{*}$$. $\mappingVector^{*} = \left(\inputMatrix^\top \inputMatrix\right)^{-1} \inputMatrix^\top \dataVector$
• The equation for $$\left.\dataStd^2\right.^{*}$$ may also be found $\left.\dataStd^2\right.^{{*}}=\frac{\sum_{i=1}^{\numData}\left(\dataScalar_i-\left.\mappingVector^{*}\right.^{\top}\inputVector_i\right)^{2}}{\numData}.$

### Solving the Multivariate System

• Section 1.3 of Rogers and Girolami (2011) for Matrix & Vector Review.

### Basis Functions

• Basis functions can be global. E.g. quadratic basis: $\basisVector = [1, \inputScalar, \inputScalar^2]$

\begin{align*} \basisFunc_1(\inputScalar) = 1, \\ \basisFunc_2(\inputScalar) = x, \\ \basisFunc_3(\inputScalar) = \inputScalar^2. \end{align*}

$\basisVector(\inputScalar) = \begin{bmatrix} 1\\ x \\ \inputScalar^2\end{bmatrix}.$

### Matrix Valued Function

$\basisMatrix(\inputVector) = \begin{bmatrix} 1 & \inputScalar_1 & \inputScalar_1^2 \\ 1 & \inputScalar_2 & \inputScalar_2^2\\ \vdots & \vdots & \vdots \\ 1 & \inputScalar_n & \inputScalar_n^2 \end{bmatrix}$

### Functions Derived from Quadratic Basis

$\mappingFunction(\inputScalar) = {\color{cyan}\mappingScalar_0} + {\color{green}\mappingScalar_1 \inputScalar} + {\color{yellow}\mappingScalar_2 \inputScalar^2}$

### Underdetermined System

• What about two unknowns and one observation? $\dataScalar_1 = m\inputScalar_1 + c$

Can compute $$m$$ given $$c$$. $m = \frac{\dataScalar_1 - c}{\inputScalar}$

### Alan Turing

Alan Turing, in 1946 he was only 11 minutes slower than the winner of the 1948 games. Would he have won a hypothetical games held in 1946? Source: Alan Turing Internet Scrapbook.

### Probability Winning Olympics?

• He was a formidable Marathon runner.
• In 1946 he ran a time 2 hours 46 minutes.
• That’s a pace of 3.95 min/km.
• What is the probability he would have won an Olympics if one had been held in 1946?

### Prior Distribution

• Bayesian inference requires a prior on the parameters.

• The prior represents your belief before you see the data of the likely value of the parameters.

• For linear regression, consider a Gaussian prior on the intercept: $c \sim \gaussianSamp{0}{\alpha_1}$

### Posterior Distribution

• Posterior distribution is found by combining the prior with the likelihood.
• Posterior distribution is your belief after you see the data of the likely value of the parameters.
• The posterior is found through Bayes’ Rule $p(c|\dataScalar) = \frac{p(\dataScalar|c)p(c)}{p(\dataScalar)}$

### Stages to Derivation of the Posterior

• Multiply likelihood by prior
• they are “exponentiated quadratics”, the answer is always also an exponentiated quadratic because $$\exp(a^2)\exp(b^2) = \exp(a^2 + b^2)$$.
• Complete the square to get the resulting density in the form of a Gaussian.
• Recognise the mean and (co)variance of the Gaussian. This is the estimate of the posterior.

### Main Trick

$p(c) = \frac{1}{\sqrt{2\pi\alpha_1}} \exp\left(-\frac{1}{2\alpha_1}c^2\right)$ $p(\dataVector|\inputVector, c, m, \dataStd^2) = \frac{1}{\left(2\pi\dataStd^2\right)^{\frac{\numData}{2}}} \exp\left(-\frac{1}{2\dataStd^2}\sum_{i=1}^\numData(\dataScalar_i - m\inputScalar_i - c)^2\right)$

$p(c| \dataVector, \inputVector, m, \dataStd^2) = \frac{p(\dataVector|\inputVector, c, m, \dataStd^2)p(c)}{p(\dataVector|\inputVector, m, \dataStd^2)}$

$p(c| \dataVector, \inputVector, m, \dataStd^2) = \frac{p(\dataVector|\inputVector, c, m, \dataStd^2)p(c)}{\int p(\dataVector|\inputVector, c, m, \dataStd^2)p(c) \text{d} c}$

$p(c| \dataVector, \inputVector, m, \dataStd^2) \propto p(\dataVector|\inputVector, c, m, \dataStd^2)p(c)$

\begin{aligned} \log p(c | \dataVector, \inputVector, m, \dataStd^2) =&-\frac{1}{2\dataStd^2} \sum_{i=1}^\numData(\dataScalar_i-c - m\inputScalar_i)^2-\frac{1}{2\alpha_1} c^2 + \text{const}\\ = &-\frac{1}{2\dataStd^2}\sum_{i=1}^\numData(\dataScalar_i-m\inputScalar_i)^2 -\left(\frac{\numData}{2\dataStd^2} + \frac{1}{2\alpha_1}\right)c^2\\ & + c\frac{\sum_{i=1}^\numData(\dataScalar_i-m\inputScalar_i)}{\dataStd^2}, \end{aligned}

complete the square of the quadratic form to obtain $\log p(c | \dataVector, \inputVector, m, \dataStd^2) = -\frac{1}{2\tau^2}(c - \mu)^2 +\text{const},$ where $$\tau^2 = \left(\numData\dataStd^{-2} +\alpha_1^{-1}\right)^{-1}$$ and $$\mu = \frac{\tau^2}{\dataStd^2} \sum_{i=1}^\numData(\dataScalar_i-m\inputScalar_i)$$.

### Two Dimensional Gaussian

• Consider height, $$h/m$$ and weight, $$w/kg$$.
• Could sample height from a distribution: $p(h) \sim \gaussianSamp{1.7}{0.0225}$
• And similarly weight: $p(w) \sim \gaussianSamp{75}{36}$

### Height and Weight Models

Gaussian distributions for height and weight.

### Independence Assumption

• This assumes height and weight are independent. $p(h, w) = p(h)p(w)$

• In reality they are dependent (body mass index) $$= \frac{w}{h^2}$$.

### Independent Gaussians

$p(w, h) = p(w)p(h)$

### Independent Gaussians

$p(w, h) = \frac{1}{\sqrt{2\pi \dataStd_1^2}\sqrt{2\pi\dataStd_2^2}} \exp\left(-\frac{1}{2}\left(\frac{(w-\meanScalar_1)^2}{\dataStd_1^2} + \frac{(h-\meanScalar_2)^2}{\dataStd_2^2}\right)\right)$

### Independent Gaussians

$p(w, h) = \frac{1}{\sqrt{2\pi\dataStd_1^22\pi\dataStd_2^2}} \exp\left(-\frac{1}{2}\left(\begin{bmatrix}w \\ h\end{bmatrix} - \begin{bmatrix}\meanScalar_1 \\ \meanScalar_2\end{bmatrix}\right)^\top\begin{bmatrix}\dataStd_1^2& 0\\0&\dataStd_2^2\end{bmatrix}^{-1}\left(\begin{bmatrix}w \\ h\end{bmatrix} - \begin{bmatrix}\meanScalar_1 \\ \meanScalar_2\end{bmatrix}\right)\right)$

### Independent Gaussians

$p(\dataVector) = \frac{1}{\det{2\pi \mathbf{D}}^{\frac{1}{2}}} \exp\left(-\frac{1}{2}(\dataVector - \meanVector)^\top\mathbf{D}^{-1}(\dataVector - \meanVector)\right)$

### Correlated Gaussian

Form correlated from original by rotating the data space using matrix $$\rotationMatrix$$.

$p(\dataVector) = \frac{1}{\det{2\pi\mathbf{D}}^{\frac{1}{2}}} \exp\left(-\frac{1}{2}(\dataVector - \meanVector)^\top\mathbf{D}^{-1}(\dataVector - \meanVector)\right)$

### Correlated Gaussian

Form correlated from original by rotating the data space using matrix $$\rotationMatrix$$.

$p(\dataVector) = \frac{1}{\det{2\pi\mathbf{D}}^{\frac{1}{2}}} \exp\left(-\frac{1}{2}(\rotationMatrix^\top\dataVector - \rotationMatrix^\top\meanVector)^\top\mathbf{D}^{-1}(\rotationMatrix^\top\dataVector - \rotationMatrix^\top\meanVector)\right)$

### Correlated Gaussian

Form correlated from original by rotating the data space using matrix $$\rotationMatrix$$.

$p(\dataVector) = \frac{1}{\det{2\pi\mathbf{D}}^{\frac{1}{2}}} \exp\left(-\frac{1}{2}(\dataVector - \meanVector)^\top\rotationMatrix\mathbf{D}^{-1}\rotationMatrix^\top(\dataVector - \meanVector)\right)$ this gives a covariance matrix: $\covarianceMatrix^{-1} = \rotationMatrix \mathbf{D}^{-1} \rotationMatrix^\top$

### Correlated Gaussian

Form correlated from original by rotating the data space using matrix $$\rotationMatrix$$.

$p(\dataVector) = \frac{1}{\det{2\pi\covarianceMatrix}^{\frac{1}{2}}} \exp\left(-\frac{1}{2}(\dataVector - \meanVector)^\top\covarianceMatrix^{-1} (\dataVector - \meanVector)\right)$ this gives a covariance matrix: $\covarianceMatrix = \rotationMatrix \mathbf{D} \rotationMatrix^\top$

### Sampling the Prior

• Always useful to perform a ‘sanity check’ and sample from the prior before observing the data.
• Since $$\dataVector = \basisMatrix \mappingVector + \noiseVector$$ just need to sample $\mappingVector \sim \gaussianSamp{0}{\alpha\eye}$ $\noiseVector \sim \gaussianSamp{\zerosVector}{\dataStd^2}$ with $$\alpha=1$$ and $$\dataStd^2 = 0.01$$.

### Computing the Posterior

$p(\mappingVector | \dataVector, \inputVector, \dataStd^2) = \gaussianDist{\mappingVector}{\meanVector_\mappingScalar}{\covarianceMatrix_\mappingScalar}$ with $\covarianceMatrix_\mappingScalar = \left(\dataStd^{-2}\basisMatrix^\top \basisMatrix + \alpha^{-1}\eye\right)^{-1}$ and $\meanVector_\mappingScalar = \covarianceMatrix_\mappingScalar \dataStd^{-2}\basisMatrix^\top \dataVector$

### Marginal Likelihood

• The marginal likelihood can also be computed, it has the form: $p(\dataVector|\inputMatrix, \dataStd^2, \alpha) = \frac{1}{(2\pi)^\frac{n}{2}\left|\kernelMatrix\right|^\frac{1}{2}} \exp\left(-\frac{1}{2} \dataVector^\top \kernelMatrix^{-1} \dataVector\right)$ where $$\kernelMatrix = \alpha \basisMatrix\basisMatrix^\top + \dataStd^2 \eye$$.

• So it is a zero mean $$\numData$$-dimensional Gaussian with covariance matrix $$\kernelMatrix$$.

### References

Bishop, C.M., 2006. Pattern recognition and machine learning. springer.

Rogers, S., Girolami, M., 2011. A first course in machine learning. CRC Press.