Introduction

To demonstrate the difference between frequentist and bayesian perspectives we will use a simple regression problem.

We will create $N$ points uniformly from $[0,1]$ and use them as input to $\sin(2\pi \mathbf{x})$. We add noise to each observation by adding a sample from a Gaussian, $\epsilon \sim \mathcal{N}(0,0.3)$. This creates our dependent random variable $t$.

$$\mathbf{t} = \sin(2 \pi \mathbf{x}) + \boldsymbol{\epsilon}$$

.

Notice that $\mathbf{t}$ is stochastic even though $\mathbf{x}$ is not. This is because $\mathbf{t}$ depends on the random variable $\boldsymbol{\epsilon}$. This is also called a nonstationary stochastic process. It is nonstationary due to the sinusoidal driving term. In these circumstances, we can calculate the analytic derivative $\frac{d\mathbf{t}}{d\mathbf{x}}$.

Given only the $N$ observations $\mathcal{D} = \lbrace \bf{x}, \bf{t} \rbrace$ our goal is to learn $\sin(2 \pi \mathbf{x})$. In other words we'd like to learn the true underlying data generating process, or find the signal through the noise. We should appreciate for a moment how strictly difficult this problem is. Given a discrete, finite, number of observations we wish to generalize to an infinite number of values. In general this is impossible to do exactly. As we shall soon see, the capacity of the model we choose $C$, the number of observations we have $N$, and the inference techniques we employ will all influence the accuracy of our approximation.

---

Let us begin with one of the simplest models, an $M^{\text{th}}$ order polynomial.

$$y(x,\mathbf{w}) = \sum_{i=0}^M w_i x^i$$

where the vector $\mathbf{w}$ contains the coefficients, or parameters of the model. Choosing a simple polynomial may at first seem naive, but remember that within some domain of convergence, a power series can approximate any function exactly. Though as we will see, using polynomials of higher and higer order does not immediately yeild better results as we might expect.

We have data and a model, how do we measure how accurately our model fits the data? There are many things we will need to consider when designing our cost function, or performance measure. For now however we use a very common technique which is to examine how far each predicted value $y(x,\mathbf{w})$ deviates from its observed value $t$.

Taking each deviation, squaring and summing the result gives us the aptly named, sum of squared errors cost function (SSE)

$$ J(\mathbf{w}) = \frac{1}{2N}\sum_{j=0}^N (t_j-y(x_j,\mathbf{w}))^2$$

where normalizing over the $N$ observations allows us to compare results for different sized datasets.

We now have all the peices of a typical machine learning problem, the experience, the task, and the performance measure. Our experience, E, is the data set $\mathcal{D} = \lbrace (x_i,t_i)\rbrace$. Our task, T, is to fit our model $y(x,\mathbf{w})$ to the dataset $\mathcal{D}$. Our performance, P, is the SSE given by $J(\mathbf{w})$. Now we need an algorithm, and according to Mitchell's definition, our algorithm will be a machine learning algorithm if its performance on the task improves with experience, or if our models' fit to the data, as measured by the SSE, improves with more observations $(x_i,t_i)$.

For this problem, the cost function is quadratic in $\mathbf{w}$ and therefore has a unique minimum. This can be found analytically, using statistics, calculus, and linear aglebra. All such derivations produce what is called the Normal Equation

$$\mathbf{w}^* = (X^TX)^{-1}X^T\mathbf{t}$$

where $X$ is the Vandermonde matrix for the $N$ points between $[0,1]$. Numerically this is a nightmare to solve, due to the inverse of the Vandermonde matrix, and is not used in practice. In practice we solve the linear system with factorization methods or iterative schemes.

---

Capacity and Overfitting

First let us examine the accuracy of our results as a function of the order of the polynomial.

Regularization

If we examine the magnitude of the parameters as a function of polynomial order, we see that they quickly become unwieldly.

Therefore we add a term to the cost function that penalizes the parameters for getting too large. Methods that employ such scale controlling terms are often referred to as shrinkage methods. In statistics it is called ridge regression and in machine learning literature it is referred to as weight decay.

$$ J(\mathbf{w}) = \frac{1}{2N}\sum_{j=0}^N (t_j-y(x_j,\mathbf{w}))^2 + \frac{\lambda}{2}\left\| \bf{w} \right\|^2$$

Dataset Size and Overfitting

Now we examine how our training and test error changes as the number of observations increase.

Probability and Bayesian Primer

Let our dataset be given by $\mathcal{D} = \lbrace \mathbf{x}, \mathbf{t} \rbrace$ and our model by $y(x,\mathbf{w})$. Our model is parameterized by $\mathbf{w}$ where $\dim{\mathbf{w}}=L$.

The joint distribution tells us everything about the data and is given by

$$p(w, \mathcal{D})$$

The product rule for probability theory allows us to express the joint distribution in terms of conditional and marginal distributions and is given by

$$p(w, \mathcal{D}) = p(w \lvert \mathcal{D})p(\mathcal{D}) = p(\mathcal{D} \lvert w)p(w)$$

The sum rule for probability theory allows us to eliminate certain variables and find what are called marginal distributions. It is given by

$$p(w) = \int p(w, \mathcal{D}) d\mathcal{D}$$$$p(\mathcal{D}) = \int p(w,\mathcal{D}) dw$$

for continuous random variables, and $$p(w) = \sum_{\mathcal{D}} p(w, \mathcal{D})$$ $$p(\mathcal{D}) = \sum_{w} p(w,\mathcal{D})$$ for discrete random variables.

Dividing by $p(\mathcal{D})$ in the product rule above we arrive at Bayes Formula.

$$p(w \lvert \mathcal{D}) = \frac{p(\mathcal{D}\lvert w)p(w)}{p(\mathcal{D})}$$

where

$$p(\mathcal{D}\lvert w) \hspace{5pt} \text{is called the Likelihood}$$$$p(w) \hspace{5pt} \text{is called the Prior}$$$$p(\mathcal{D}) \hspace{5pt} \text{is called the Marginal}$$$$p(w \lvert \mathcal{D}) \hspace{5pt} \text{is called the Posterior}$$

Another common way of expressing Bayes Formula involves combining both the product and sum rules

$$p(w \lvert \mathcal{D}) = \frac{p(\mathcal{D}\lvert w)p(w)}{\sum_{w}p(\mathcal{D}\lvert w)p(w)}$$

Probabilistic Perspective

Let us now cast the regression problem in a stochastic framework. Again our goal is to predict the target value $t$ given an input $x$ based on a set of traning data $\mathbf{x} = (x_1,\dots,x_n)^T$ and and their corresponding target values $\mathbf{t} = (t_1,\dots,t_n)^T$. We know there is noise in the data, there always is, and until now we have not taken it into consideration. Now however we will directly consider the uncertainty in our predictions. We will assume that given an input $x$, our target value $t$ is given by a Gaussian distribution with a mean $y(x,\mathbf{w})$ and a standard deviation $\sigma$. As per our previous formulation the mean $y(x,\mathbf{w})$ is our prediction, however unlike our previous formulation here we explicitly consider our uncertainty via the standard deviation $\sigma$.

$$p(t \lvert x,\mathbf{w},\sigma) = \mathcal{N}(t\lvert y(x,\mathbf{w}),\sigma)$$

Now if we assume that our training data $\lbrace \mathbf{x}, \mathbf{t} \rbrace$ are independent and identically distributed (iid), then we can express our likelihood function as

$$p(\mathbf{t} \lvert \mathbf{x},\mathbf{w},\sigma) = \prod_{n=1}^N \mathcal{N}(t_n\lvert y(x_n,\mathbf{w}),\sigma)$$

We notice a few things here. Firstly in casting this problem stochastically we have increased the number of parameters from $\dim{w}$ to $\dim{w}+1$ for the extra $\sigma$ term. By assuming $\sigma$ to be constant for all values of $x$ we're assuming the noise is what is called, homoscedastic. If instead we believed the noise were to change for different values of $x$, then we could allow the standard deviation to be a function of the input $\sigma(x)$, in which case we call the noise, heteroscedastic.

This new formulation changes our task T. Before we were fitting an $M^{th}$ order polynomial to the data, here however we wish to fit parameters $\theta = \lbrace{\mathbf{w}, \sigma}\rbrace$. Our goal this time will be to find the parameters $\theta^*$ that maximize the probability $p(\mathbf{t} \lvert \mathbf{x}, \mathbf{w}, \sigma)$. This goal, maximizing the likelihood, is the quintessential frequentist example. Where before we explicity chose the performance measure P to be the (SSE), here we will analytically optimize the likelihood and use the result as our performance measure.

The log function is monotonic, therefore, the extrema of $\ln f(x)$ equal the extrema of $f(x)$. Considering this and the fact that the logarithm of exponential distributions are analytically convenient, we will optimize the log of the likelihood.

$$ \ln p(\mathbf{t} \lvert \mathbf{x}, \mathbf{w}, \sigma) = -\frac{1}{2\sigma} \sum_{n=1}^N \lbrace y(x_n,\mathbf{w}) - t_n \rbrace^2 + \frac{N}{2} \ln \frac{1}{\sigma}- \frac{N}{2} \ln(2\pi)$$

Let us maximize the log likelihood with respect to $\mathbf{w}$. We notice two things, firstly that the second two terms are independent of $\mathbf{w}$, and secondly that maximizing the likelihood is equivalent to minimizing the negative log likelihood.

$$ \mathbf{w}_{ML} = \max_\mathbf{w} \hspace{4pt} \ln p(\mathbf{t} \lvert \mathbf{x}, \mathbf{w}, \sigma) = \min_{\mathbf{w}} \frac{1}{2} \sum_{n=1}^N \lbrace y(x_n,\mathbf{w}) - t_n \rbrace^2 = \min \text{SSE}$$

Thus the SSE has arisen as a consequence of maximizing the likelihood under the assumption of Gaussian noise. The optimal parameters found with such a technique are aptly called the maximum likelihood solutions (ML). With the maximum likelihood solution for $\mathbf{w}$ found we can also find the maximum likelihood solution for the standard deviation.

$$\sigma_{ML} = \frac{1}{N} \sum_{n=1}^N \lbrace y(x_n, \mathbf{w})-t_n \rbrace^2$$

With the optimal parameters found we can define our predictive distribution

$$p(t \lvert x, \mathbf{w}_{ML}, \sigma_{ML}) = \mathcal{N}(t \lvert y(x, \mathbf{w}_{ML}), \sigma_{ML})$$

A Step Towards Bayes

Let us first recall the distribution used in the previous example.

$$p(\bf{t} \lvert \bf{x},\bf{w},\sigma)$$

This is the likelihood function. It poses the problem as "what is the probability of the data, given a certain set of parameters". In the previous example we found the set of parameters that maximized this distribution. In other words, we're finding the parameters that maximize our data ON AVERAGE. This is the underlying distinction between Bayes and frequentist perspectives. Let's go deeper...

$$ \text{Frequentist} \hspace {110pt} \text{Bayesian}$$$$p(\bf{t} \lvert \bf{x},\bf{w},\sigma) \hspace{70pt} p(\bf{w} \lvert \bf{x}, \bf{t}, \sigma) \propto p(\bf{t} \lvert \bf{x},\bf{w},\sigma)p(\bf{w})$$$$ $$

On the left we have the frequentist perspective, or the perspective of frequencies of random repeatable events. In the frequentist perspective, the data is considered the random variable and it aims at finding solutions considering our uncertainty of observed data. The parameters are considered fixed and their value is determined by some form of 'estimator'. Error bars or uncertainty in the parameters are found by considering the distribution of possible datasets $\mathcal{D}$.

In a frequentist perspective the data is the random variable.

On the right we have the Bayesian viewpoint. From this perspective there is only one single data set $\mathcal{D}$, the one that is actually observed. The uncertainty in the parameters is expressed via the probability distribution describing them.

In a Bayesian perspective the parameters are the random variable.

Let us examine the curve fitting problem with a prior over the parameters. For ease of exposition we will assume a simple Gaussian form

$$p(\bf{w} \lvert \alpha) = \mathcal{N}(\bf{w} \lvert \bf{0}, \alpha \bf{I}) = \bigg ( \frac{1}{2\pi \alpha} \bigg)^{(M+1)/2}\exp\bigg\{-\frac{1}{2\alpha}\bf{w}^T\bf{w}\bigg\}$$

Using Bayes theorem above we can express the posterior in terms of the likelihood and the prior. We interpret their product as a distribution describing the probability of $\bf{w}$ given the observed data $\{\bf{x},\bf{t}\}$. If we find the $\bf{w}$ that maximize this posterior distribution we are finding the so-called maximum posterior solutions (MAP).