# COMP90051 Workshop 2
## Part B: Bayesian inference

***

In this part of the workshop, we'll develop some intuition for priors and posteriors, which are crucial to Bayesian inference.


```python
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import bernoulli, beta
```

### 1. A lucky find

On the way to class, you discover an unusual coin on the ground.



As a dedicated student in statistical ML, you're interested in determining whether the coin is biased. 
More specifically, you want to estimate the probability $\theta$ that the coin will land heads-up when you toss it.

You can use the function below to simulate a coin toss: it returns `1` for heads and `0` for tails.


```python
def toss_coin():
    if bernoulli.rvs(p = (int.from_bytes("coin".encode(), 'little') % 10000)/10000):
        return 1
    return 0
```

### 2. Prior belief
Before you even toss the coin, you notice that the heads side appears to have more mass. 
Thus, your _prior belief_ is that $\theta$ is slightly biased away from 0.5 towards 0—i.e. you expect tails are more likely.

To quantify this prior belief, we assume that the prior distribution for $\theta$ is $\mathrm{Beta}(a,b)$, for some choice of the hyperparameters $a, b > 0$. 
(See [link](https://en.wikipedia.org/wiki/Beta_distribution) for info about the Beta distribution.)
The prior probability density function for $\theta$ is therefore given by:

$$ p(\theta) = \frac{1}{B(a,b)} \theta^{a-1} (1 - \theta)^{b-1} $$
where $B(a,b)$ is a special function called the _Beta function_.

Select appropriate values for $a$ and $b$ by looking at the plot of $p(\theta)$ below.


```python
a = 2
b = 1
theta = np.linspace(0, 1, 1000)
plt.plot(theta, beta.pdf(theta, a, b), 'b-', lw=1)
plt.xlabel(r'$\theta$')
plt.ylabel(r'$p(\theta)$')
plt.show()
```

### 3. Posterior updates
Now toss the coin once and denote the outcome by $x_1$.


```python
x1 = toss_coin()
print(x1)
```

    1


We can update our belief about $\theta$, based on this new evidence $x_1$.
To do this we apply Bayes' rule to compute the posterior for $\theta$:
$$ p(\theta | x_1) = \frac{p(x_1 | \theta) \, p(\theta)}{p(x_1)} \propto p(x_1 | \theta) \, p(\theta)$$
where $p(\theta)$ is the prior given above and 
$$ p(x_1 | \theta) = \theta^{x_1} (1 - \theta)^{1 - x_1} $$
is the likelihood.

***
**Exercise:** Show (on paper) that
$$ p(\theta | x_1) \propto \theta^{x_1 + a - 1} (1 - \theta)^{(1 - x_1) + b - 1} $$
which implies that $\theta | x_1 \sim \mathrm{Beta}[x_1 + a, (1 - x_1) + b]$.
***

Toss the coin a second time, denoting the outcome by $x_2$.


```python
x2 = toss_coin()
print(x2)
```

    1


Again, we want to update our belief about $\theta$ based on the new information $x_2$. 
We take the previous posterior $p(\theta|x_1)$ as the new prior and apply Bayes' rule:
$$ p(\theta | x_1, x_2) \propto p(x_2 | \theta) p(\theta | x_1)$$
\[Note: We assume the tosses are independent, otherwise the likelihood for $x_2$ would depend on $x_1$.\] 
This gives $\theta | x_1, x_2 \sim \mathrm{Beta}[x_1 + x_2 + a, (2 - x_1 - x_2) + b]$.

***
**Exercise:** Show that for $n$ coin tosses, the posterior is $\theta | x_1, \ldots, x_n \sim \operatorname{Beta}[n_H + a, n - n_H + b]$ where $n_H = \sum_{i = 1}^{n} x_i$ is the number of heads observed.
***

### 4. MAP estimator and MLE estimator

The posterior $\theta|x_1, \ldots, x_n$ contains all the information we know about $\theta$ after observing $n$ coin tosses.
One way of obtaining a point estimate of $\theta$ from the posterior, is to take the value with the maximum a posteriori probability (MAP):
$$
\begin{align}
    \hat{\theta}_{MAP} &= \arg \max_{\theta} p(\theta|x_1, \ldots, x_n) \\
        & = \frac{n_H + a - 1}{n + a + b - 2}
\end{align}
$$

In general, the MAP estimator gives a different result to the maximum likelihood estimator (MLE) for $\theta$:
$$
\begin{align}
    \hat{\theta}_{MLE} &=\arg \max_{\theta} p(x_1, \ldots, x_n|\theta) \\
        & = \frac{n_H}{n}
\end{align}
$$

***
**Exercise:** Derive the above results for $\hat{\theta}_{MAP}$ and  $\hat{\theta}_{MLE}$.
***

### 5. Convergence of the estimates

Let's now toss the coin an additional 48 times (so that $n = 50$), recording $\hat{\theta}_{MLE}$ and $\hat{\theta}_{MAP}$ after each toss.


```python
extra_tosses = 48
num_tosses = 2 + extra_tosses
num_heads = 0
theta_map = np.zeros(num_tosses)
theta_mle = np.zeros(num_tosses)
for i in range(0, num_tosses):
    if i == 0: 
        num_heads += x1 
    elif i == 1:
        num_heads += x2
    else:
        num_heads += toss_coin()
    theta_map[i] = ... # fill in
    theta_mle[i] = ... # fill in
```

We plot the results below.


```python
plt.plot(theta_map, label = "MAP")
plt.plot(theta_mle, label = "MLE")
plt.xlabel('Number of draws')
plt.ylabel(r'$\hat{\theta}$')
plt.legend()
plt.show()
```

**Questions:** 

1. Is the coin biased?
1. Do the MAP and MLE estimates converge to the same value for $\theta$?
1. What happens if you set $a = 1; b = 1$?
1. How does the posterior distribution for $\theta$ compare to the prior plotted above? (Use the code block below to plot the posterior.)


```python
plt.plot(theta, beta.pdf(theta, a + num_heads, b + num_tosses - num_heads), 'b-', lw=1)
plt.xlabel(r'$\theta$')
plt.ylabel(r'$p(\theta|x_1, \ldots, x_n)$')
plt.show()
```


```python

```

e276cf23ecfa081eee816a5a7830922c485ae9cf
ipynb
Jupyter Notebook
Statistical_Machine_Learning/workshop/worksheet02b.ipynb
infinityglow/Unimelb-CS-Subjects
07bdb49fd4c50035b7f2e80ca218ac2b620098e4
MIT
2022-02-14T16:31:07.000Z
2022-02-14T16:31:07.000Z
Statistical_Machine_Learning/workshop/worksheet02b.ipynb
hidara2000/Unimelb-CS-Subjects
07bdb49fd4c50035b7f2e80ca218ac2b620098e4
MIT
Statistical_Machine_Learning/workshop/worksheet02b.ipynb
hidara2000/Unimelb-CS-Subjects
07bdb49fd4c50035b7f2e80ca218ac2b620098e4
MIT
2021-06-14T11:59:13.000Z
2021-06-14T11:59:13.000Z
Qwen/Qwen-72B
1. YES
2. YES