Uncertainty Quantification is coming! Modern approaches to Bayesian Inference are leaving traditional predictive models in the dust. Get behind the wheel and be part of a new dawn of simple and interpretable probabilistic modelling!

> We will focus on giving you an intuition for Bayesian inference, allowing you to apply intricate methods to any problem.

**In this course we will cover:**

- The foundations of probability and uncertainty, including the frequentist and Bayesian interpretations.
- Interpretations of the components of Bayes' Theorem and how they interact.
- How to fit probabilistic models to data, to yield, not one result, but all the probable results!
- Formulating a Bayesian Inference problem from scratch and solving it using both analytical and numerical methods.

**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.

### Section 1: Intro to Bayes

In this first section, we introduce the foundations of probability and Bayes' Theorem.

- What is Bayes' Theorem and how to use it to compute conditional probabilities?
- How to apply Bayes' theorem to solve inference problems and compute posterior densities from expert knowledge and data.

We are going to start gently by introducing different perspectives on probability and uncertainty. We will then use those insights to formulate inference problems in the Bayesian formalism and solve them on a computer.

### Section 2: Bayesian Linear Regression and parameter Estimation

In section two, we will introduce the Bayesian framework for linear regression and non-linear parameter estimation.

- Why would you need Bayesian Uncertainty Quantification for something as simple as linear regression? An example of an ill-posed regression problem.
- How to formulate and solve a nonlinear parameter estimation problem using the Bayesian formalism to explore the full posterior of probable model parameters.

We are going to solve these problems using direct methods when possible and when not, we are going to harness the amazing power of the Metropolis-Hastings algorithm - a simple and elegant approach for exact Bayesian Inference!


If you want to deploy AI in the Wild with confidence, then you need Uncertainty Quantification (UQ). Fundamentally, we must be able to assess the reliability of our estimations. This course will introduce you to the many faces of Bayesian Inference - the gold standard of robust and explainable UQ - from the theoretical foundations to Bayesian parameter estimation.

Dr Mikkel Lykkegaard

<Callout type="success" heading={"📑 Learning Objectives"}>
- Have an overview of the different disciplines of Uncertainty Quantification.
- Understand the foundational difference between Frequentist and Bayesian statistics.
- Understand the different between objective and subjective undercertainty.
- Know Kolmogorov's Axioms of probability.

</Callout>

## An Brief Overview of UQ Methods

Uncertainty Quantification (UQ) can be many things, from simply doing model sensitivity analysis to solving statistical inverse problems. Here is a (by no means exhaustive) overview of some of the different mehtods used to perform UQ:

- **Forward Uncertainty Quantification** - propagating uncertainty through some operation.
  - (Quasi-) Monte Carlo Methods - using samples to propagate the uncertainty.
  - Spectral methods such as Polynomial Chaos Expansion - using clever polynomials to propagate the uncertainty.
- **Inverse Uncertainty Quantification or Parameter Estimation** - Estimating the posterior distribution of parameters given some data.
  - Markov Chain Monte Carlo - using clever algorithms to sample from the posterior distribution.
  - Variational Inference such ADVI - using a variational distribution to approximate the posterior distribution function.
- **Predictive Uncertainty Quantification** - Quantifying the uncertainty of predictions.
  - Gaussian Process Regression - using a distribution over functions to estimate the predictive uncertainty.
  - Engineered methods - e.g. NOMU neural networks.
  - Ad-Hoc methods - e.g. Monte Carlo Dropout.
  - Approximate methods - e.g. Laplace Approximation.

in this course, we are mainly going to focus on the second branch, name _inverse UQ_, in particular direct method for estimating the posterior distribution and additionally a bit of Markov Chain Monte Carlo (MCMC), since the implementation of the simplest algorithm, the Metropolis-Hastings algorithm, is straightforward.

However, what all of the above approaches have in common is that they operate with propbability distributions. Our first port of call is therefore a slightly deeper dive into the meaning and interpretations of probability and uncertainty.

## Bayesian vs Frequentist

Before we can get into the nuts and bolts of Uncertainty Quantification for Machine Learning, we have to take a few steps back. In fact, we have to to go back to the very beginning of things - the foundations of probability theory and statistics. Here, we find two fundamentally different perspectives and schools of thought, namely the _Frequentists_ and the _Bayesians_. Now, I'm not personally a partisan. I believe in using the right tool for the job. However, the differences between Frequentist and Bayesian approaches run deeper than simply using different formulae, and each school of thought comes with its own set of axioms and assumptions. To truly appreciate this perceived schism in statistical sciences, we have to briefly visit the underlying _epistemological_ differences.

### Objective vs Subjective Uncertainty

The main difference boils down to the different views of Frequentists and Bayesians on the nature of uncertainty. Broadly, probability and uncertainty can be viewed as either _objective_ or _subjective_:

- **Objective Uncertainty**: Quantifies the uncertainty due to _randomness of events_ (aleatoric uncertainty). More specifically, probability is interpreted as _the limit of the relative frequency of many trials_. From this perspective, we assume a perfect model, but we do not know the next outcome. An example of objective uncertainty would be the outcome of a perfect dice. A perfect dice has no imerfections, yet due to the inherent randomness of rolling it, we are uncertain about its next outcome. This is the Frequentist perspective.
- **Subjective Uncertainty**: Quantifies a degree of _belief_. At its most extreme, this type of uncertainty is not tied to any physical process at all, and we are quantifying our subjective belief. An example would be using probability to describe the guiltiness of an accusee in a trial. This is the Bayesian perspective.
- **Everything in between**. Here, we have some imperfect knowldge of the system we are trying to quantify the uncertainty of. An example would be a imperfect dice. We are not certain that all the sides have the same probability, and the goal is to attempt to quantify that uncertainty and make predictions based on our observations. In this case, both Frequentist and Bayesian methods are broadly used, but each method comes with its own set of assumptions and its unique ways of processing information.

For a more in-depth review of the differences between Frequentist and Bayesian statistics, I highly recommend Kristin Lennox's talk [All About That Bayes](https://www.youtube.com/watch?v=eDMGDhyDxuY). Also check out this brilliant [XKCD comic](https://xkcd.com/1132/).

## Kolmogorov Axioms

The reason I bring up the Kolmogorov axioms is not to bore you with with a load of (very fundamental) mathematics. It's because I want to make the point that while providing rigour to the foundations of probability, they also provide _freedom_ to your probabilistic models. That is to say, as long as the probability functions we use conform to these simple rules, we can use the entire machinery of probability theory to make inference.

Consider a probability space $(\Omega, F, P)$, where $\Omega$ is the _sample space_, that is the space of all possible outcomes, $F$ is the _event space_, that is the space of all relevant sets of possible outcomes and $P$ is the probability function, which measures the probability of an even $E \in F$. Then Kolmogorov's axiom are:

1. The probability of an event is a positive, real number: $P(E) \in \mathbb{R}, P(E) \geq 0 \quad \forall E \in F$.
2. The probability that at least one event from the sample sapce will occur is 1: $P(\Omega) = 1$.
3. The probability of the union of disjoint (non-overlapping) events is equal to the sum of their individual probabilities: $P(\cup_{i=1}^{\infty}E_i) = \sum_{i=1}^{\infty}P(E_i)$

Moreover, these axiom has various consequences that we are not going into detail with in this lesson, but you can look them up in any textbook about probability.


Introduction to Uncertainty Quantification

<Callout type="success" heading={"📑 Learning Objectives"}>
- Identify the various components of Bayes' formula: the posterior, prior, likelihood and evidence.
- Use Bayes' formula to compute simple conditional probabilities for events.

</Callout>

## Bayes' Theorem

Here is a fun fact: the inventor of Bayes' Theorem (Thomas Bayes) was not himself a Bayesian, as we understand it today. He simply saw the Bayesian Theorem as a way to calculate conditional probablities. So let's first have a look at Bayes' Theorem, the way it was intended by its creator.

### For events

For events $A$ and $B$, we can write Bayes' Theorem as:

$$
P(A|B) = \frac{P(A)P(B|A)}{P(B)}
$$

where

- $P(A)$ is the _prior_ probability of $A$.
- $P(B|A)$ is the _likelihood_ of $B$ given $A$.
- $P(B)$ is _marginal_ probability or _evidence_, and,
- $P(A|B)$ is the _posterior_ probability of $A$ given $B$.

The marginal likelihood can be expanded by summing up over all possible events ${A_i} \in \mathcal A$:

$$
P(B) = \sum_i P(A_i)P(B|A_i)
$$

I appreciate that there is a lot to unpack here, and it doesn't get any simpler as we move from events to continuous random variables, and from conditional probabilities to Bayesian Inference. So let's take a moment to dive into an example of how to use Bayes' Theorem for a real-world example of conditional event estimation.

#### Example

Suppose that you have tested positive for COVID-19 using a lateral flow test, and you want to calculate the posterior probability that you actually have the disease. Let's introduce some notation first. The test can either be positive or negative, so we will call the event that the test is positive $t_+$, and, conversely, that it is negative $t_-$. Similarly, you can either be COVID-positive or COVID-negative, and we'll call those events $c_+$ and $c_-$, respectively:

- $t_+$ : Positive test
- $t_-$ : Negative test
- $c_+$ : COVID-positive
- $c_-$ : COVID-negative

If we encode this into Bayes' Theorem, we would like to find the posterior probability of having COVID-19, given that you just tested positive $P(c_+|t_+)$:

$$
P(c_+|t_+) = \frac{P(c_+)P(t_+|c_+)}{P(t_+)}
$$

Let's unpack this. $P(c_+)$ is the prior probability of being COVID-positive. If we have no more information, let's assume that 10% of the population has COVID-19, so that $P(c_+)=0.1$.

The likelihood is a bit more involved, and we will have to look into the actual _sensitivity_ and _specificity_ of the lateral flow test. We'll have to take those numbers from somewhere, and this [news article](https://www.ox.ac.uk/news/2020-11-11-oxford-university-and-phe-confirm-lateral-flow-tests-show-high-specificity-and-are) tells us that they are 76.8% and 99.7%, respectively. Another word for the sensitivity is the _true positive rate_, or, in other words, $P(t_+|c_+) = 0.768$. The specificity is also called the _true negative rate_, i.e. $P(t_-|c_-) = 0.997$.

The _marginal likelihood_ expresses the likelihood of observing the data, considering all possible values of the variable of interest, in this case your COVID status. This can be written in the following way:

$$
{P(t_+)} = P(c_+)P(t_+|c_+) + P(c_-)P(t_+|c_-)
$$

To calculate this, we need the two final ingredients, namely the prior probability of being COVID-negative $P(c_-)$ and the false positive rate of the lateral flow test $P(t_+|c_-)$. These can easily be obtained from the numbers we already have:

$$
\begin{aligned}
P(c_-) = 1 - P(c_+) = 1 - 0.1 &= 0.9 \\
P(t_+|c_-) = 1 - P(t_-|c_-) = 1 - 0.997 &= 0.003
\end{aligned}
$$

Putting it all together, we get:

$$
P(c_+|t_+) = \frac{P(c_+)P(t_+|c_+)}{P(c_+)P(t_+|c_+) + P(c_-)P(t_+|c_-)} = \frac{0.1 \cdot 0.768}{0.1 \cdot 0.768 + 0.9 \cdot 0.003} = 0.966
$$

or, 96.6%. Let's flip the calculations around and see what is the probability of being _negative_, given that you get a negative test $P(c_-|t_-)$:

$$
P(c_-|t_-) = \frac{P(c_-)P(t_-|c_-)}{P(c_+)P(t_-|c_+) + P(c_-)P(t_-|c_-)} = \frac{0.9 \cdot 0.997}{0.1 \cdot 0.232 + 0.9 \cdot 0.997} = 0.974
$$

or, 97.9%. This means that if you test negative using a lateral flow testing kit, you are almost certainly negative, and you can get on with your life.


The Importance of Being Bayes: Introducing Bayes' Theorem

<Callout type="success" heading={"📑 Learning Objectives"}>
- Identify various Bayesian schools of thought, objective, subjective and strongly subjective.
- Understand the different roles of the prior.
- Define the likelihood function and understand how it differs from a probability distribution.
- Appreciate the challenge of estimating the evidence / marginal likelihood.

</Callout>

## Unpacking Bayes' Theorem

We will now move on to the Bayesian interpretation of Bayes' Theorem. As mentioned in the previous session, Thomas Bayes himself was not a Bayesian, but simply saw Bayes' Theorem as a way to calculate conditional probabilities. So what does it mean to be Bayesian? The short answer is that it involves interpreting probability as an expression of _belief_, but the long answer is that it varies.

### The Different Flavours of Bayesian

The different schools of Baysian statistics can broadly be separated into the _objective_ Bayesians, the _subjective_ Bayesians and the _strongly subjective_ Baysieans. The difference between each school, boil down to their view of the _prior_ (well get more into that later).

- **Objective Bayesians**: Believe that the prior should be non-informative, and designed in such a way that it avoids biasing the result. A lot of effort has been invested into this branch of Bayesian statistics, since un-biasedness is attractive to most statisticians. However, it is a rather complicated and theoretical topic, and we are not going to cover it detail in this course. It is often achieved throught the application of the [Jeffreys Prior](https://en.wikipedia.org/wiki/Jeffreys_prior), which we will look into in the next session.
- **Subjective Bayesians**: Believe that priors should be constructed on the bases of _expert elicitation_. That is, the prior should incorporate the current body of knowledge about the parameter of the prior, as described by experts in the field. This allows for encoding prior _beliefs_ into the problems and to exploit previous literature on the topic. A good example is the probability of the sun coming up tomorrow, as explained in this [XKCD comic](https://xkcd.com/1132/). We are mainly going to take this approach in this course.
- **Strongly Subjective Bayesians**: This school of thought is representative of the idea of _probability as belief_ taken to the extreme. Since we can only ever express probability as subjective belief, everything goes, as long as you are explicit about which prior you are using, and that you can qualify and defend your choice of prior. Then your audience can decide for themselves whether they buy into your analysis, since they are fully informed about your prior beliefs.

To put this into perspective, let's briefly revisit Bayes' Theorem for events, and then move on to expanding to continuous random variables.

### For events

As explained in the previous session, for events $A$ and $B$, we can write Bayes' Theorem as:

$$
P(A|B) = \frac{P(A)P(B|A)}{P(B)}
$$

where

- $P(A)$ is the _prior_ probability of $A$.
- $P(B|A)$ is the _likelihood_ of $B$ given $A$.
- $P(B)$ is _marginal_ probability or _evidence_, and,
- $P(A|B)$ is the _posterior_ probability of $A$ given $B$.

The marginal likelihood can be expanded by summing up over all possible events ${A_i} \in \mathcal A$:

$$
P(B) = \sum_i P(A_i)P(B|A_i)
$$

However, Bayes' Theorem can also be employed to continous random variables. Let's have a look at that.

### For continuous random variables

Suppose you have some sensor data $d$, and you want to perform Bayesian calibration of some model parameters $\theta$ using those sensor data. We can then find the posterior density of the model parameters conditional on the sensor data:

$$
p(\theta|d) = \frac{p(\theta)p(d|\theta)}{p(d)}
$$

The definitions are the same as above, but since the parameters $\theta$ are now continuous, the marginal likelihood becomes an integral:

$$
p(d) = \int p(\theta)p(d|\theta) \: d\theta
$$

This rarely has an analytical sulution and is often intractible to compute. Workarounds for the challenge exist, including Markov Chain Monte Carlo (MCMC) and Variational Inference (VI), but those are both highly advanced methods that each warrant multiple lectures of their own. A more straightfoward worksaround is the use of _conjugate priors_, which we will have a brief look at in the next session.

## The meaning of it all

### Prior

As mentioned in the introduction to this session, the prior means different things to different types of Bayesian. However, the most common perspective is that it expresses a prior _belief_ about the parameters. There are different ways to construct such a prior.

- First, it could be extracted from the **previous body of literature** on the topic. For example, if you wanted to test whether actually vaccines work, it would be sensible to take into account the immense body of historical evidence that they do, in fact, work. This (fairly conclusive) historical evidence could then be encoded as a prior, so that it would require a lot of new data to change that conclusion.
- Second, priors can be distilled through **expert elicitation**. In this case, domain experts are be interrogated on their beliefs about model parameters, and theor knowlegde is then converted into a prior distribution. If we wanted to calculate the probablity that the sun is going to rise tomorrow, we might ask an atrophysicist their opinions. They would look at their models, maybe check a table of the positions of heavenly bodies, and provide their prior belief as to whether the sun will rise tomorrow.
- Third, the prior could express **physical constraints** of the model parameters. Maybe we are trying the infer the posterior distribution of cyclists passing our house each day. We know that such a parameter must a positive integer, and we can encode that knowledge into our prior.

In any case, the prior is nevessary for Bayes' Theorem to do its job, and it is not always obvious how to define it. In the next session I will show you how the choice of prior affects the outcome, and give some practical guidance on the choice of prior.

### Likelihood

The likelihood is something that, broadly, measures the misfit between the model and the observations, as a function of the model parameters. If the we ignore the prior, and simply maximise the likelihood function, we get the Maximum Likelihood Estimate, which is equivalent of the Frequentist parameter estimate.

Importantly, the likelihood looks a lot like a probability distribution, but **it is not** a probability distribution, since it does not integrate to 1. Let's see why.

Suppose you are flipping a coin a bunch of times and counting the heads. That process is called a Bernoulli process, and is modelled using a [_binomial_ distribution](https://en.wikipedia.org/wiki/Binomial_distribution):

$$
k \sim \mathcal B(n,h)
$$

here, the random variable $k$, the number of heads, follows a binomial distribution with parameters $n$, the number of trials, and $h$, the probability of getting a heads. This distribution has support $k = \{0, \dots, n\}$.

Suppose now we have completed an experiment and have a measurement of $k$ heads after $n$ trials. We would now like to compute our posterior probability distribution $h$, i.e. the probability of getting a heads. If we fix $n$ (the number of trials) we get:

$$
p(h|k) = \frac{p(h) \mathcal B(k|h)}{p(k)}
$$

Since $\mathcal B(k|h)$ is now seen as a function of $h$, it is no longer a valid probability distribution. The binomial dsitribution is constructed such that the probability mass of any possible outcome $k$ sums up to 1. When we flip it around and make it a function of $h$, it is not even a discrete function anymore, but a continuous function with domain $h \in [0,1]$.

This is a bit of an intricate detail, and it surprises many people a first. It certainly surprised me. However, getting an intuition of what a likelihood function does and how it works is very important for the remainder of this course. I will give an example in the next session where we will also look depper into the coin flipping experiment.

### Evidence

The evidence is broadly a normalisation constant which makes sure that the posterior is a valid probability distribution. The deeper meaning of it is not very relevant for this course, but it is important to appreciate that it can be hard to actually calculate in practice, since it is an integral over all model parameters:

$$
p(d) = \int p(\theta)p(d|\theta) \: d\theta
$$

This makes Baysian Inference particularly challenging in high dimensions. Different workarounds to this challenge exist, namely:

- **Conjugate priors**: Using appropriate priors, the posterior will have a closed-form solution, and we can avoid calculating the evidence alltogether. This is fast, exact and works for a variety of problems. However, most problems cannot be solved this way.
- **Variational Inference**: Here we use a _variational_ distribution that does have a closed form solution to approximate the posterior. This is a very flexible approach, but it is only approximate and does introduce some bias in the posterior.
- **Markov chain Monte Carlo**: This is an exact method for drawing samples from the posterior. However, it can be extremely computiationally costly, and should only be used when no other method is adequate.

### Posterior

The posterior is the actual target of all this shenanigans. It expresses our beliefs about our model parameters, given that we have observed some data, conditioned on our prior beliefs about the model parameters.

I am not going to dwell further on this. Let's instead have a look at an example.


Unpacking Bayes' Theorem: Prior, Likelihood and Posterior

:::success
**📂 Resources**
Download the resources for this lesson <Link href="/resources/uncertainty-quantification-for-machine-learning/prior-and-posterior-coin-flipping-walkthrough/Resources.zip" download="Resources.zip">here.</Link>
:::

<Callout type="success" heading={"📑 Learning Objectives"}>
- Using Bayes' Theorem to solve a coin flipping parameter estimation problem.
- Identify different prior distributions and use them.
- Use simple numerical methods to find the evidence in Bayes' Theorem.

</Callout>

## Using Bayes' Theorem to perform inference on a coin flipping example.

We have had a look at the interpretation of the prior, likelihood and posterior. Now let's have a look at how they might be used in practice. We start by defining a function that returns the prior distribution, likelihood function and posterior distribution for a given number of successes $k$ out of $n$ trials, with several different options for the prior distribution, namely the Jeffreys prior, the conjugate (Beta) prior, a normal prior and the uniform prior.

```python
def bernoulli_bayes(k, n, prior='jeffreys', prior_parameters=None):

    # define the Jeffreys prior for a binomial likelihood.
    if prior == 'jeffreys':
        prior = lambda p: 1 / np.sqrt(p * (1-p))

    # define the conjugate (Beta) prior for a binomial likelihood.
    elif prior == 'beta':
        # extract the prior parameters if they were given.
        if prior_parameters is not None:
            a, b = prior_parameters
        # otherwise set the default a = b = 0.5, corresponding to the Jeffreys prior.
        else:
            a, b = 0.5, 0.5
        prior = lambda p: beta(a=a, b=b).pdf(p)

    # define a normal prior.
    elif prior == 'normal':
        # extract the prior parameters if they were given.
        if prior_parameters is not None:
            loc, scale = prior_parameters
        # otherwise set it to a reasonably narrow normal centered at 0.5 (unbiased coin).
        else:
            loc, scale = 0.5, 0.1
        prior = lambda p: norm(loc=loc, scale=scale).pdf(p)

    # define a uniform prior.
    elif prior == 'uniform':
        prior = lambda p: uniform().pdf(p)

    # define the likelihood. this must be a function of the distribution parameter p,
    # rather than the support of the distribution.
    likelihood = lambda p: binom(n=n, p=p).pmf(k)

    # set the evaluation range (some priors are undefined at 0 and 1.
    range = (1e-3, 1 - 1e-3)
    p = np.linspace(*range, 100)

    # evaluate the prior and likelihood.
    prior_eval = prior(p)
    likelihood_eval = np.array([likelihood(_p) for _p in p])

    # if we use the conjugate prior (Beta), the posterior has a closed-form solution.
    if prior == 'beta':
        posterior = beta(a=a+k, b=b+(n-k))
        posterior_eval = posterior.pdf(p)

    # otherwise, we are using quadrature to evaluate the evidence.
    else:
        evidence = quad(lambda _p: prior(_p)*likelihood(_p), 0, 1)[0]
        posterior_eval = prior_eval*likelihood_eval / evidence

    # return prior, likelihood and posterior.
    return prior_eval, likelihood_eval, posterior_eval
```

Let's also define a plotting function that can visualise our different distributions.

```python
def plot_distributions(prior_eval, likelihood_eval, posterior_eval):

    # set the range
    range = (1e-3, 1 - 1e-3)
    p = np.linspace(*range, 100)

    fig, ax = plt.subplots(figsize=(16,4), ncols=3, nrows=1)

    # plot the prior.
    ax[0].set_title('Prior')
    ax[0].plot(p, prior_eval)
    ax[0].set_xlabel(r'$\gamma$')
    ax[0].set_ylabel(r'$\rho$')

    # plot the likelihood.
    ax[1].set_title('Likelihood')
    ax[1].plot(p, likelihood_eval)
    ax[1].set_xlabel(r'$\gamma$')
    ax[1].set_ylabel(r'$\rho$')

    # plot the posterior.
    ax[2].set_title('Posterior')
    ax[2].plot(p, posterior_eval)
    ax[2].set_xlabel(r'$\gamma$')
    ax[2].set_ylabel(r'$\rho$')

    plt.show()
```

Let's try and use it to make inference for an example where we got $k=2$ successes out of $n=3$ trials:

```python
successes = 2
trials = 3

for prior in ['jeffreys', 'beta', 'normal', 'uniform']:
    print('{} prior'.format(prior))
    plot_distributions(*bernoulli_bayes(successes, trials, prior))
```

<Image
	src={
		"/images/courses/uncertainty-quantification-for-machine-learning/prior-and-posterior-coin-flipping-walkthrough/posteriors_small_sample.png"
	}
	alt={"Equator split"}
	width={877}
	height={1120}
/>

As you can see, the posteriors for the Jeffreys prior and the Beta prior with $a=b=0.5$ are identical, since the priors are identical up to a constant. For the informed (Normal) prior, the posterior is near-identical to the likelihood, since for such a small dataset with $n=3$ trials, there is simply not enough evidence to suggest that the coin is any different than what we expected. For the Uniform prior, the posterior is identical to the likelihood function up to a constant, since the likelihood function is not a proper probability distribution.

Now let's try and see what happens when we observe more data. Let's preserve the same ratio between successes and trials, so that $k=66$ and $n=100$.

```python
successes = 66
trials = 100

for prior in ['jeffreys', 'beta', 'normal', 'uniform']:
    print('{} prior'.format(prior))
    plot_distributions(*bernoulli_bayes(successes, trials, prior))
```

<Image
	src={
		"/images/courses/uncertainty-quantification-for-machine-learning/prior-and-posterior-coin-flipping-walkthrough/posteriors_large_sample.png"
	}
	alt={"Equator split"}
	width={879}
	height={1121}
/>

With a larger number of trials, we see that the posterior contracts around (what we now think is) the true bias of the coin, namely $p=0.66$, for all the four different prior distributions. Since all the posterior distributions are more or less identical, it is clear that the prior does not have much influence on the posterior when more data is available.


Prior and Posterior - Coin Flipping Walkthrough

<Callout type="success" heading={"📑 Learning Objectives"}>
- Identify ill-conditioned systems of equations.
- Use a ridge regressor to regularise the linear least squares solution.
- Find the Bayesian posterior distribution for a linear model with a Gaussian prior and likelihood.

</Callout>

## Linear Least Squares Breakdown

### A Motivating Example

Consider the system of equations

$$
\bf Ax = b
$$

with

$$
{\bf A} = \begin{bmatrix}0.16 & 0.10 \\ 0.17 & 0.11 \\ 2.02 & 1.29\end{bmatrix}
\quad \text{and} \quad
{\bf b} = \begin{bmatrix}0.27 \\ 0.25 \\ 3.33\end{bmatrix}
$$

This is a minimal example from my favourite book on inverse problems, Hansen (2010). While this is a contrived example, it illustrates an very common issue with least squares estimates.

The least squares estimator for this equation is

$$
\hat{\bf x} = ({\bf A}^T {\bf A})^{-1} {\bf A}^T {\bf b} = \begin{bmatrix} 7.021 \\ -8.40 \end{bmatrix}.
$$

Now that looks fine at first glance, but what if I told you that the right hand side $b$ was noisy?

Our vector of measurements $b$ had been perturbed by some noise $\varepsilon$, i.e. ${\bf b} = {\bf b}_{true} + \varepsilon$, with

$$
\varepsilon = \begin{bmatrix}0.01 \\ -0.03 \\ 0.02\end{bmatrix}.
$$

This is not a huge amount of noise, but the least squares solution to ${\bf Ax} = {\bf b}_{true}$ is in fact ${\bf \hat{x}}_{true} = \begin{bmatrix} 1 \\ 1 \end{bmatrix}$! We can modify the original equation a little to take the noise into account:

$$
\bf Ax + \varepsilon = b
$$

with $\varepsilon \sim \mathcal N(0, \sigma)$, that is the signal $\bf b$ includes additive white Gaussian noise.

In this example, I knew exactly how much noise was added, so I could recover the true coefficients, but in reality you will not know the exact measurement errors. If you did, they wouldn't be errors! In the best case scenario, you will know the approximate variance of the noise.

How did this precariuous situation come to be? Well, as it turns out, $\bf A$ is what is commonly referred to as _ill-conditioned_, meaning that (certain) large pertubations in the input $\bf x$ do not change the output $\bf b$ much. And conversely, small pertubations in the measurements can lead to large pertubatons in the estimated coefficients.

### Regularisation

A common approach to alleviate the above problem is _ridge regression_, also sometimes referred to as _Tikhonov regularisation_. The idea here is to add a shift to the diagonal of the moment matrix ${\bf A}^T {\bf A}$ before inverting it, leading to the ridge estimator

$$
\hat{\bf x} = ({\bf A}^T {\bf A} + \lambda {\bf I})^{-1} {\bf A}^T {\bf b}
$$

where $\lambda$ is the ridge parameter, which determines the "stength" of regularisation. For example, for $\lambda = 0.001$, we get $\hat{\bf x}_{\lambda=0.001} = \begin{bmatrix} 1.20 \\ 0.70 \end{bmatrix}$, which is much closer to the true parameters than the naïve least squares estimator.

This is equivalent of solving the regularised least squares problem

$$
\hat{\bf x} = \underset{\bf x}{\text{min}} \; \lVert{\bf b} - {\bf Ax}\rVert_2^2 + \lambda \lVert{\bf x}\rVert_2^2.
$$

That is all fine and well, but it leaves an open question of _how to choose_ $\lambda$? In practice, this is often done heuristically, and according to certain well established criteria, as one can study in much more detail in e.g. Hansen (2010). However, these methods still produce a single estimate, rather than a distribution of possible estimates. This is why we will now give this problem the Bayesian treatment.

### A Bayesian Approach

To keep this section simple, we are going to make a few relatively strong assumptions. First, we will assume that the noise comes from a zero-mean Gaussian with known variance $\sigma_\varepsilon^2$, $\varepsilon \sim \mathcal N({\bf 0}, \sigma_\varepsilon^2 {\bf I})$. That allows is to write the likelihood as $p({\bf b}|{\bf x}) = \mathcal N({\bf b} - {\bf Ax}, \sigma_\varepsilon^2 {\bf I})$. Furthermore, we assume that the prior distribution of model parameters is also Gaussian, $p({\bf x}) = \mathcal N({\bf \mu}_0, {\bf \Sigma}_0)$. That alows us to write the _posterior_ distribution as (you guessed it!) a Gaussian:

$$
p({\bf x} | {\bf b}) = \mathcal N(\mu_t, {\bf \Sigma}\_t)
$$

with

$$
\mu*t = {\bf \Sigma}\_t \:({\bf \Sigma}\_0^{-1} \mu_0 + \sigma*\varepsilon^{-2} {\bf A}^T {\bf b});\; \text{and}\\
{\bf \Sigma}_t^{-1} = {\bf \Sigma}\_0^{-1} + \sigma_\varepsilon^{-2} {\bf A}^T{\bf A}
$$

Please refer to e.g. Bishop (2006) for more details. For example, if we for the problem outlined above choose $p({\bf x}) = \mathcal N({\bf 0}, {\bf I})$, and $\sigma_\varepsilon^2 = 0.01$,
we get

$$
{\bf \mu_t} = \begin{bmatrix}1.17 \\ 0.74 \end{bmatrix}
\quad \text{and} \quad
{\bf \Sigma}\_t = \begin{bmatrix}0.29 & -0.45 \\ -0.45 & 0.71\end{bmatrix}.
$$

### References

Hansen, Per Christian. Discrete Inverse Problems: Insight and Algorithms. Fundamentals of Algorithms. Philadelphia: Society for Industrial and Applied Mathematics, 2010.

Bishop, Christopher M. Pattern Recognition and Machine Learning. Information Science and Statistics. New York: Springer, 2006.


Linear Least Squares Breakdown and Bayesian Linear Regression

:::success
**📂 Resources**
Download the resources for this lesson <Link href="/resources/uncertainty-quantification-for-machine-learning/bayesian-linear-regression-a-coding-walkthrough/Resources.zip" download="Resources.zip">here.</Link>
:::

<Callout type="success" heading={"📑 Learning Objectives"}>
- Use ridge regression to regularise an ill-conditioned system of equations.
- Compute the Bayesian posterior distribution analytically for a linear model with a Gaussian prior and likelihood.
- Code a simple Markov Chain Monte Carlo (MCMC) algorithm and draw samples from the posterior distribution.

</Callout>

## Baysian Linear Regression in Practice: A Coding Walkthrough.

Let's define an ill-conditioned linear least squares problem and find and estimate for the regression coefficients using ordinary least squares.

```python
# ill-conditioned least squares problem, see Hansen (2010)
A = np.array([[0.16, 0.10],
              [0.17, 0.11],
              [2.02, 1.29]])
b = np.array([0.27, 0.25, 3.33])

# find the ordinary least-squares solution
beta = np.linalg.lstsq(A, b, rcond=None)[0]

# print the solution
print(beta)
[ 7.00888731 -8.39566299]
```

This is looks fine. It is _a_ solution, but it turns out to be an outlier in the Bayesian posterior. The first thing we will do to try and find a better solution is to regularise the least squares problem using ridge regression.

```python
# set the range of regularisation parameters that we test.
lambdas = np.logspace(-4,-3,101)

# set up a list to hold all the solutions.
betas = []

# iterate through the regularisation parameters.
for lamb in lambdas:
    # get the ridge estimate for the current regularisation parameters.
    beta = np.linalg.multi_dot((np.linalg.inv(np.dot(A.T, A) + lamb*np.eye(2)), A.T, b))
    # append it to the list of estimates.
    betas.append(beta)

# convert the list of estimates to a NumPy array.
betas = np.array(betas)

# plot the regularisation parameter against the norm of the ridge estimates.
plt.plot(lambdas, np.linalg.norm(betas, axis=1))
plt.show()
```

<Image
	src={
		"/images/courses/uncertainty-quantification-for-machine-learning/bayesian-linear-regression-a-coding-walkthrough/ridge_decay.png"
	}
	alt={"Equator split"}
	width={558}
	height={418}
/>

Here we see the norm of the regression coefficients decaying as we turn up the intensity of the regularisation. Let's print out the regularised solution corresponding to $\lambda = 0.001$.

```python
# print the beta corresponding to lambda = 0.001
print(betas[-1])
[1.19815118 0.70322534]
```

This is much closer to the "true" solution of $x = [1,1]^\intercal$, but it is still a point estimate. Let's give the problem the Bayesian treatment.

```python
# import the multivariate normal for the MCMC
from scipy.stats import multivariate_normal

# set the prior mean
mu_prior = np.zeros(2)
# set the prior covariance
sigma_prior = 1
cov_prior = sigma_prior**2*np.eye(2)

# define the prior.
prior = multivariate_normal(mean=mu_prior,
                            cov=cov_prior)

# set the standard deviation of the observation noise.
sigma_data = 0.1
# set the covariance of the observation noise
cov_data = sigma_data**2*np.eye(3)

# define the likelihood
# NOTE: this is not strictly correct since the mean should be the model output
# and the support should be the observations, and here they are swapped. it is
# just slightly more convenient like this and it doesn't matter for the Gaussian
# since that is symmetric with respect to the mean and the support.
likelihood = multivariate_normal(mean=b,
                                 cov=cov_data)

# compute the posterior mean and covariance
# see Bishop's Pattern Recognition and Machine learning
cov_post = np.linalg.inv(np.linalg.inv(cov_prior) + sigma_data**(-2)*np.dot(A.T, A))
mu_post = np.dot(cov_post, np.dot(np.linalg.inv(cov_prior), mu_prior) + sigma_data**(-2)*np.dot(A.T, b))

# define the posterior as a multivariate normal distribution
posterior = multivariate_normal(mean=mu_post, cov=cov_post)

# print the posterior mean
print(mu_post)
[1.17108637 0.74162625]

# print the posterior covariance
print(cov_post)
[[ 0.29074825 -0.45261217]
 [-0.45261217  0.71048368]]
```

Here, clearly the posterior mean is very close to the ridge estimate. This is not a coincidence and there is a subtle connection between the two. However, we also get an estimate for the posterior covariance, which is something that was missing in the ridge estimate. We get the entire distribution of possible model parameters, rather than a point estimate. Let's plot the posterior and see how it looks like.

```python
# set up a grid to evaluate the posterior distribution.
x = np.linspace(-1, 4,1000)
y = np.linspace(-2, 3,1000)
X, Y = np.meshgrid(x,y)
grid = np.c_[X.reshape(-1,1), Y.reshape(-1,1)]

# evaluate the posterior on the grid
posterior_test = posterior.pdf(grid).reshape(X.shape)

# plot the posterior
plt.figure(figsize=(6,6))
plt.imshow(posterior_test, origin='lower', extent=(x.min(), x.max(), y.min(), y.max()))
plt.show()
```

<Image
	src={
		"/images/courses/uncertainty-quantification-for-machine-learning/bayesian-linear-regression-a-coding-walkthrough/posterior_distribution.png"
	}
	alt={"Equator split"}
	width={516}
	height={511}
/>
Here we can see that the "true" solution $x = [1,1]^\intercal$ is covered by the
posterior distribution and even has a fairly high density. The probable solutions
are distributed on a line, however, and there are clearly many possible parameters
that can reproduce the data. This is because the moment matrix $A^T A$ is near-singular.
However, using the Bayesian approach we were able to identify all the possible solutions
with very little effort. Let's finally set up an MCMC sampler and draw some samples
from the posterior.

```python
# set up a Markov Chain Monte Carlo sampler to sample the posterior.

# set the preconditioned Crank-Nicolson (pCN) scaling
beta = 0.1
# set the number of samples
n_samples = 12000

# draw the initial sample
samples = [prior.rvs()]

# iterate through the number of samples
for i in range(n_samples):

    # generate a pCN proposal
    proposal = np.sqrt(1 - beta**2)*samples[-1] + beta*prior.rvs()

    # compute the acceptance probability
    alpha = np.exp(likelihood.logpdf(np.dot(A, proposal)) - \
                   likelihood.logpdf(np.dot(A, samples[-1])))

    # draw a random number and accept/reject
    if np.random.uniform() < alpha:
        samples.append(proposal)
    else:
        samples.append(samples[-1])

# convert the samples (minus the burnin of 2000 samples) to a NumPy array
samples = np.array(samples[2000:])

# plot the postterior samples
plt.figure(figsize=(6,6))
plt.scatter(samples[:,0], samples[:,1], alpha=0.01, c='k')
plt.xlim(x.min(), x.max())
plt.ylim(y.min(), y.max())
plt.show()
```

<Image
	src={
		"/images/courses/uncertainty-quantification-for-machine-learning/bayesian-linear-regression-a-coding-walkthrough/posterior_samples.png"
	}
	alt={"Equator split"}
	width={519}
	height={511}
/>

These samples indeed look like they are from the same distribution as the analytical solution, even though they were generated in a completely different way. And they are exactly from that distribution. But why would we use MCMC when we can just compute the exact posterior? Most often we are no so fortunate that there exists a analytical posterior, but MCMC works anyway. We will have a closer look at that in the next lesson. Here are some good ressources on MCMC for the curious of you:

### References

Brooks, S. (Ed.). (2011). Handbook for Markov chain Monte Carlo. Taylor & Francis.
Liang, F., Liu, C., & Carroll, R. J. (2010). Advanced Markov chain Monte Carlo methods: Learning from past samples. Wiley.
Liu, J. S. (2008). Monte Carlo strategies in scientific computing (2. ed). Springer.


Bayesian Linear Regression: A Coding Walkthrough

<Callout type="success" heading={"📑 Learning Objectives"}>
- Formulate the general data-generating model for Bayesian parameter estimation.
- Choose a prior and likelihood function for your Bayesian parameter estimation problem.
- Identify algorithms for Bayesian Inference, namely Markoc Chain Monte Carlo (MCMC) and Variational Inference (VI)

</Callout>

## Bayesian Parameter Estimation

### The data-generating model

Let's first revisit the model from the previous lesson, namely

$$
\bf Ax + \varepsilon = b
$$

This is written in the usual ordinary linear least squares notation, but in the Bayesian formalism it usually takes the form

$$
d_{obs} = \mathcal F(\theta) + \varepsilon
$$

where $d_{obs}$ is a vector of observations ($\bf b$ in the previous notation), $\theta$ is a vector of model parameters ($\bf x$ in the previous notation) for which we seek the posterior distribution, and $\mathcal F$ is the so-called _forward operator_ which describes how the model maps parameters $\theta$ to noise-free observations $d_{obs} - \varepsilon$. This operator can be anything, including a linear map as before, e.g. $\mathcal F(\theta) = \bf A \theta$.

### Another look at the likelihood and the prior

When the noise $\varepsilon$ is Gaussian, we can write the likelihood as $p(d_{obs}|\theta) = \mathcal N(d_{obs} - \mathcal F(\theta), \sigma_\varepsilon^2 {\bf I})$, or, more generally, $p(d_{obs}|\theta) = \mathcal N(d_{obs} - \mathcal F(\theta), \Sigma_\varepsilon)$, which is equivalent of

$$
p(d_{obs}|\theta) = (2\pi)^{-d/2} |\Sigma_{\varepsilon}|^{-1/2} \exp \left( -\frac{1}{2}(d_{obs} - \mathcal F(\theta))^\intercal \: \Sigma_{\varepsilon}^{-1} \: (d_{obs} - \mathcal F(\theta))\right)
$$

where $\Sigma_\varepsilon$ is the covariance matrix of the Gaussian noise. However, the noise is also not always Gaussian. For example, if we are trying to model something countable (number of cars passing a certain location, number of animals in a nature reserve, etc.) we may want to use e.g. a Poisson likelihood instead.

After correctly identifying the noise model and hence the likelihood function, the task is now to choose a prior for the model parameters $\theta$, i.e. $p(\theta)$ to complete Bayes' Theorem:

$$
p(\theta|d*{obs}) = \frac{p(\theta) p(d*{obs}|\theta)}{p(d\_{obs})}
$$

As we saw in the previous session, when the likelihood is Gaussian, the conjugate prior is also Gaussian. However, the conjugate prior may not always be the most adequate representation of our prior knowledge. Instead of choosing the conjugate prior for convenience, we have to put on our modelling glasses and make an informed choice for our prior.

### Bringing out the big guns

The Bayesian formalism allows for using any prior you like with any noise model you like, and even allows for using hierachical models with hyperpriors on the prior parameters. However, the convenient closed-form solution that we encountered in the previous session is unfortunately a rarity. If we want to exploit the full flexibility of the Bayesian formalism, we have to introduce some bigger guns - **inference algorithms**.

As mentioned in the very first lesson, there are two principal methods for Bayesian parameter estimation, namely Markov Chain Monte Carlo (MCMC) and Variational Inference (VI). We had a brief look at MCMC in the previous walkthrough, and we are going to expand on that in the next session, where we will model COVID cases. MCMC is asymptotically exact and also easy to implement, and we are going to stick to MCMC for this course for these reasons. VI is usually much faster than MCMC but it requires a little more background to understand. Moreover, it only provides an approximation to the posterior which is biased by the choice of variational distribution.


Bayesian Parameter Estimation

:::success
**📂 Resources**
Download the resources for this lesson <Link href="/resources/uncertainty-quantification-for-machine-learning/bayesian-parameter-estimation-walkthrough-with-covid-example/Resources.zip" download="Resources.zip">here.</Link>
:::

<Callout type="success" heading={"📑 Learning Objectives"}>
- Use ordinary least squares regression on a log-transformed variable to get an exponential model.
- Define the forward operator of a Baysian parameter estimation problem.
- Set up Bayesian prior distributions for discrete and continuous variables and set up a likelihood function.
- Write a simple MCMC sampler that samples from the posterior distribution.

</Callout>

## Bayesian Parameter Estimation: Walkthrough with COVID Example.

In this example, we are going to model COVID cases in the UK in the period 1-15 March 2020, just after the first cases were detected in the country. Lets first load the data and do a few transformations.

```python
# read the data from the CSV and convert the date into a datetime
data = pd.read_csv('covid.csv')
data['date'] = pd.to_datetime(data['date'])

# create a new variable, which is days since the first measurement
data['day'] = (data['date'] - data['date'].iloc[0]).dt.days
```

We are now going to perform ordinary least-squares (OLS) regression on the log-transformed output variable, the number of cases, i.e. $\log y = \beta_0 + \beta_1 \: t$. We are going to transform the OLS coefficients back into the original space so that $y(t) = e^{\beta_0} e^{\beta_1 \: t} = y_0 e^{r \: t}$. Then we will plot the model along with the datapoints.

```python
# find the least squares estimate for the log-transformed output variable
lr = np.linalg.lstsq(np.c_[np.ones(len(data)), data['day']], np.log(data['cases']), rcond=None)
print(lr)
(array([3.59852482, 0.26155652]), array([0.19268671]), 2, array([32.02976831,  2.02334923]))

# convert the coefficients back into the original space
y0 = np.exp(lr[0][0])
r = lr[0][1]

# compute the response of the model according to the least squares estimate
t = np.linspace(0,14)
y = y0 * np.exp(r*t)

# plot the model response and and data
data.plot.scatter('day', 'cases', c='r')
plt.plot(t, y)
plt.show()
```

<Image
	src={
		"/images/courses/uncertainty-quantification-for-machine-learning/bayesian-parameter-estimation-walkthrough-with-covid-example/ols_and_data.png"
	}
	alt={"Equator split"}
	width={580}
	height={432}
/>

That looks pretty good. However, this is just the OLS estimate and it is, as such, a point estimate. Not all the datapoints are well-described by the model. But using the Bayesian formalism that we looked at in the previous lesson, we can instead find the posterior distribution of the model parameters, and use that to the evaluate all the possible scenarios. We will first extract the data from the dataframe, then we will define the forward model $\mathcal F(\theta)$ and then we will set up some priors for our parameters $\theta = (y_0, r)$.

```python
# extract the data as NumPy arrays
x_data = data.day.values
y_data = data.cases.values

# set up the model
def model(y0, r):
    return (y0*np.exp(r*x_data)).astype(int)

# define the priors
prior_y0 = stats.poisson(36)              # poisson prior for the initial value
prior_r = stats.lognorm(s=0.5, scale=0.2) # lognorm prior for the rate

# plot the priors
fig, ax = plt.subplots(figsize=(12,4), nrows=1, ncols=2)

x = np.arange(0,60)
p = prior_y0.pmf(x)
ax[0].set_xlabel('$y_0$')
ax[0].set_ylabel('$p(y_0)$')
ax[0].scatter(x, p)

x = np.linspace(0, 1)
p = prior_r.pdf(x)
ax[1].set_xlabel('$r$')
ax[1].set_ylabel('$p(r)$')
ax[1].plot(x, p)

plt.show()
```

<Image
	src={
		"/images/courses/uncertainty-quantification-for-machine-learning/bayesian-parameter-estimation-walkthrough-with-covid-example/priors.png"
	}
	alt={"Equator split"}
	width={1014}
	height={372}
/>

These are the priors, and they encode the soft knowledge we have about the parameters before conditioning on the data. Now let's try and draw some samples from the prior, pass them through our forward model and plot the result against the data.

```python
# plot some draws from the prior and run them through the model
for i in range(1000):
    y0 = prior_y0.rvs()
    r = prior_r.rvs()
    y = model(y0, r)
    plt.plot(x_data, y, c='k', alpha=0.01)

# plot the data along with the prior draws
plt.scatter(x_data, y_data, c='r')
plt.ylim(0,2000)
plt.show()
```

<Image
	src={
		"/images/courses/uncertainty-quantification-for-machine-learning/bayesian-parameter-estimation-walkthrough-with-covid-example/prior_samples.png"
	}
	alt={"Equator split"}
	width={560}
	height={418}
/>

That looks pretty good. It is clear that all the samples are using the exponential model, but they are more dispersed than the data would permit. We could maybe allow for more flexibility in the $y_0$ prior, but let's keep it like this for simplicity. Now let's define the joint prior and the likelihood function.

```python
# define the joint prior
def prior(x):
    return sum([prior_y0.logpmf(x[0]), prior_r.logpdf(x[1])])

# define the likelihood
def likelihood(x):
    # we are using a negative binomial likelihood since it is a bit more flexible than
    # the poisson distribution. It allows us to also model the variance independent of the
    # mean. it is typically used as a kind of overdispersed Poisson.
    return stats.nbinom.logpmf(k=0, n=np.abs(model(x[0], x[1]) - y_data)+1, p=0.99).sum()
```

Let's also plot the likelihood function for a range of inputs.

```python
# plot the likelihood function
x = np.arange(-300, 300)
p = stats.nbinom.pmf(k=0, n=np.abs(x)+1, p=0.99)
plt.scatter(x, p, s=1)
plt.show()

```

<Image
	src={
		"/images/courses/uncertainty-quantification-for-machine-learning/bayesian-parameter-estimation-walkthrough-with-covid-example/likelihood.png"
	}
	alt={"Equator split"}
	width={547}
	height={413}
/>

What this shows us is that the likelihood is highest for $\epsilon = 0$, that is zero miscounts. It then decays slowly as we increase number of (absolute) miscounts. It is still pretty likely to have e.g. 100 miscounts on any given day of the timeseries that we are modelling. Lets' plot the joint prior and the likelihood function evaluated on the actual data.

```python
# evaluate the joint prior and likelihood for our problem.
y0_range = np.arange(10,60)
r_range = np.linspace(0.1, 0.5)
y0_grid, r_grid = np.meshgrid(y0_range, r_range)

prior_eval = np.zeros(y0_grid.shape)
likelihood_eval = np.zeros(y0_grid.shape)

for i in range(r_range.size):
    for j in range(y0_range.size):
        prior_eval[i,j] = prior([y0_grid[i,j], r_grid[i,j]])
        likelihood_eval[i,j] = likelihood([y0_grid[i,j], r_grid[i,j]])

# plot the joint prior and the likelihood function.
fig, ax = plt.subplots(figsize=(12,4), nrows=1, ncols=2)

ax[0].set_title('Prior')
y0_plot = ax[0].imshow(np.exp(prior_eval), extent=(10, 60, 0.1, 0.35), origin='lower', aspect='auto')
ax[0].set_xlabel('$y_0$')
ax[0].set_ylabel('$r$')
plt.colorbar(y0_plot, ax=ax[0])

ax[1].set_title('Likelihood')
r_plot = ax[1].imshow(np.exp(likelihood_eval), extent=(10, 60, 0.1, 0.35), origin='lower', aspect='auto')
ax[1].set_xlabel('$y_0$')
ax[1].set_ylabel('$r$')
plt.colorbar(r_plot, ax=ax[1])

plt.show()
```

<Image
	src={
		"/images/courses/uncertainty-quantification-for-machine-learning/bayesian-parameter-estimation-walkthrough-with-covid-example/prior_and_likelihood.png"
	}
	alt={"Equator split"}
	width={1011}
	height={394}
/>

Here we see that the joint prior is pretty diffuse and permits most values across this range of parameters. The likelihood, however, is very concentrated along a narrow ridge in this parameter space. Let's set up and MCMC sampler and draw some samples from the posterior distribution. We are then going to plot the posterior samples.

```python
# set up a proposal
def proposal(x):

    # set the step mechanism for each parameter
    y0_dx = int(stats.norm(0, 3).rvs()) # zero-centered normal cast as int for y0
    r_dx = stats.norm(0, 0.01).rvs()    # zero-centered normal for r

    # create a vector of the steps
    dx = np.array([y0_dx, r_dx])

    # compute the proposal
    p = x + dx

    # return the proposal
    return p

# set an empty list for the samples
samples = []

# draw the first MCMC state from the prior
y0_0 = prior_y0.rvs()
r_0 = prior_r.rvs()
samples.append(np.array([y0_0, r_0]))

# iterate through the MCMC
for i in tqdm(range(12000)):

    # create a proposal from the current state
    p = proposal(samples[-1])

    # comoute the acceptance probability
    alpha = np.exp(prior(p) + likelihood(p) - prior(samples[-1]) - likelihood(samples[-1]))

    # accept or reject the proposal
    if np.random.uniform() < alpha:
        samples.append(p)
    else:
        samples.append(samples[-1])

# create numpy array and remove burnin
samples = np.array(samples)[2000:,:]

# plot the posterior samples
plt.scatter(samples[:,0], samples[:,1], c='k', alpha=0.01)
plt.show()
```

<Image
	src={
		"/images/courses/uncertainty-quantification-for-machine-learning/bayesian-parameter-estimation-walkthrough-with-covid-example/posterior_samples.png"
	}
	alt={"Equator split"}
	width={556}
	height={413}
/>

As we can see, the posterior samples are fairly close to the likelihood that we plotted earlier. This is not so surprising since the prior is relatively diffuse and so the likelihood function dominates the posterior. Here, we can also clearly see the bands of $y_0$ on the x-axis, which occur because we model that variable as an integer. Finally, let's run some posterior samples throught the model and do some forecasting.

```python
# create a new time vector to do some forecasting
t = np.linspace(0, 20)

# draw 1000 samples from the posterior
for i in range(1000):
    idx = np.random.randint(0, samples.shape[0])
    y = samples[i,0]*np.exp(samples[i,1]*t)

    # plot each posterior sample
    plt.plot(t, y, c='k', alpha=0.01, zorder=0)

# plot the data along with the posterior samples
plt.scatter(x_data, y_data, c='r')
plt.savefig('../images/forecasting.png', bbox_inches='tight')
```

<Image
	src={
		"/images/courses/uncertainty-quantification-for-machine-learning/bayesian-parameter-estimation-walkthrough-with-covid-example/forecasting.png"
	}
	alt={"Equator split"}
	width={560}
	height={413}
/>

Here, we have cranked the time-horizon up to 20, so we predict for 6 days into the future. Since we are using a probabilistic (Bayesian) model, we also get the uncertainty of the forecast. This tells us that the uncertainty increases the futher we try to predict, and after only 6 days into the future, there could be anything between 4000 and 8000 cases. This is a really important takeaway, since it allows for decision makers such as hospitals to make preparations for the worst-case scenario. Using the Bayesian formalism, obtaining such worst-case estimates is straightforward and relatively painless.


An Introduction to Bayesian Uncertainty Quantification

Course Description

Section 1: Intro to Bayes

Section 2: Bayesian Linear Regression and parameter Estimation

Lessons

1. Introduction to Uncertainty Quantification

2. The Importance of Being Bayes: Introducing Bayes' Theorem

3. Unpacking Bayes' Theorem: Prior, Likelihood and Posterior

4. Prior and Posterior - Coin Flipping Walkthrough

5. Linear Least Squares Breakdown and Bayesian Linear Regression

6. Bayesian Linear Regression: A Coding Walkthrough

7. Bayesian Parameter Estimation

8. Bayesian Parameter Estimation: Walkthrough with COVID Example