Bayesian inference

Difference between Bayesian and frequentist inference

The difference is in the definition of probability

Bayesian

Bayesian Approach to Statistics - Probability is a measure of uncertainty.
- Both parameters and data in a statistical model can contain uncertainty
- "The population mean is between 3.5 and 5 with probability 28%"

Frequentist

Frequentist Approach to Statistics - Probability is a proportion (frequency) of successes in an infinite number of repeated experiments
- This is how it was explained in the Probability - Fundamentals lecture from Eva's course
- Only data and estimators can be given a probabilistic statement, not parameters

Under a frequentist definition we cannot say for example that "The population mean is between 3.5 and 5 with probability 28%".

However, we can say "the estimator for population mean is between 3.5 and 5 in 28% of the experiments".

The Bayesian definition is inclusive of the frequentist definition, but not the other way around.

Example - fishes

This example requires an understanding of conditional distribution, marginalization and integration (calculus), which are essential in Bayesian stats

Background on these concepts

Conditional distribution was discussed as part of Categorical covariates. It's the distrubution of a random response variable $Y$ given that another variable $X$ is at a fixed position
Marginalization is basically the opposite of conditional distribution - the effect of a single variable on the response without the effect of another variable
Integration was discussed in Continuous Distributions lecture in Eva's class (area under a curve)

Let's go back to our fish population mean example. Assuming we have a population of individuals whose length ( $y_{i}$ ) follow a Gaussian distribution:

y_{i} | μ \sim N (μ, σ^{2})

This is our probability density function and observation model.

Assume that the variance $σ^{2}$ is known, so we are only trying to use statistical inference to find the population mean $μ$
As shown in Week 3 (maximum likelihood estimate), the likelihood function for a normal distribution is:

p (y | μ) = \prod_{i = 1}^{n} N (y_{i} | μ, σ^{2})

Therefore the maximum likelihood function is $\bar{μ} = max_{μ} p (y | μ)$ . This is a function of the parameters ( $μ$ ) when the data ( $y$ ) is fixed

Conditional Probability

The first plot shows how the probability density of $y$ is conditional on the value of $μ$ . Here are six discrete values of $μ$ from 0-5:
!lecture_notes_week5_1, p.3
Let's say that there is an equal chance of $μ$ being one of these six possible values ( $p (μ) = 1 / 6$ ). This value is the "prior probability" (as we learned in maths class as part of Bayes' Theorem). It's a known unconditional probability.

Joint probability density

From Bayes' Theorem, the joint probability density for $y_{i}$ and $μ$ is $p (y_{i}, μ) = p (y_{i} | μ) p (μ)$ .

This is a function with two parameters now. If we look at how they co-vary with the conditional distributions we will get

If we had a continuous (non-discrete) values of $μ$ (unlike how its's described in the lecture notes at this part in the example)we have something like:

A bit confusing but I think the important part is that the joint probability density is affected by both variables, and these both contribute to the marginal probability (below) ## Marginalization Now, based on marginalization (same as Law of total probability) we can sum up all the discrete possibilities to get the "marginal distribution" - the probability regardless of which $\mu$ value is being used: $$p(y_i)=p(y_i|\mu=0)p(\mu=0)+p(y_i|\mu=1)p(\mu=1)+...+p(y_i|\mu=5)p(\mu=5)$$ $$p(y_i)=\sum_{k=0}^5p(y_i|\mu=k)p(\mu=k)$$ Visualization of this: !lecture_notes_week5_1, p.4

Of course, our example discrete prior probability values are not very applicable to real life, especially to our fish example. In theory, their mean could be any positive real number (infinitely many possible values). Therefore, instead of summing up discrete values, we should instead use an integral from $- \infty$ to $\infty$ (or whatever range of values is plausible for your data)

p (y_{i}) = \int_{- \infty}^{\infty} p (y_{i} | μ) p (μ) d μ

Bayes theorem

As discussed in math class we can use Bayes' Theorem to describe this:

p (μ | y) = \frac{p (y | μ) p (μ)}{p (y)}

Anatomy

$p (μ | y)$ - Posterior probability - the probability of getting this mean value ( $μ$ ) given this data ( $y$ ).
$p (y | μ)$ - Likelihood function - the likelihood of getting these value of $y$ given $μ$ -in this example this is the Gaussian distribution $p (y | μ) \sim N (μ, σ^{2})$
$p (μ)$ - Prior probability density - the likelihood of having this value of $μ$ in any scenario.
$p (y) = \int p (y | μ) p (μ) d μ$ - Marginal probability density or marginal likelihood of model. We called this "total probability" in Eva's course, these are the same thing.

Applied

If we sampled one individual from the population which has $y = 1$ , and the prior probability is known (there is an equal (uniform) chance of $y$ being any number between -2 and 7 and no chance of anything else), we can get the probability of the individual having this value. This is something we cannot do in the frequentist domain.

In R:

# possible values for mu
mu = seq(-2,7,by=0.1)
# Uniform probability for values of mu
Pr.mu = rep(1/length(mu),length(mu)) 
# e.g. equal probability of 1/91 for being 91 discrete values.

# Likelihood function for mu
y1 = 1
likelihood = dnorm(y1, mean = mu, sd = sqrt(sigma2))

## Bayes theorem
posterior.prob = Pr.mu * likelihood / sum(Pr.mu * likelihood)

# Plot the posterior distribution
plot(mu, posterior.prob)

!lecture_notes_week5_1, p.7
This result is a probability density at each of these 90 discrete values from -2 to 7 at 0.1 intervals. Since this is a probability density, we can do statistical inference like getting the mean and variance.

The Posterior mean, which is similar to the maximum likelihood estimate:

E [μ | y] = \sum_{k = 1}^{| μ |} μ_{k} p (μ_{k} | y)

In R, we would calculate this with:

posterior.mean = sum(mu*posterior.prob) # posterior.prob is a vector with 91 values
print(posterior.mean)
## [1] 1

We can also calculate the posterior probability of the population mean being above zero:

Pr (μ > 0 | y) = \sum_{μ_{k} > 0} p (μ_{k} | y)

posterior.prob.over.zero = sum(posterior.prob * (mu>0) )
print(posterior.prob.over.zero)
## [1] 0.9833605

With integrals

We are still pretty far from actually applying this - these uniform distributions are not realistic. To get something more realistic, we need to integrate.

More realisitc Posterior mean with integrals:

E [μ | y] = \int μ p (μ | y) d μ

In frequentist statistics, the most complex operation to do is finding maxima and minima - AKA derivatives. For Bayesian statistics, estimating parameters requires the use of integrals, which are very difficult to calculate by hand or with a computer (uses much processing power). This is why Bayesian statistics has only started to become more useful after the 1990s when more powerful computers became available.

Bayesian Computation

In our example, we used 91 points as an approximation for integration. If we are using many more points and adding more parameters this kind of calculation becomes terrible - for example if we had a sequence of 100 plausible values and ten parameters, we would need to calculate $100^{10}$ posterior probability densities.

The most common way to do these tough calculations is the Markov chain Monte Carlo methods. Used to sample random realizations from the posterior distribution, then use this posterior sample to calculate the inference summaries.

Using some R techniques which aren't actually important for the lesson (just use for visualizing) we can visualize a histogram of the posterior probability density based on the random samples. Now, using these samples we can calculate the "inferential summaries" such as the Posterior mean with this formula:

E [μ | y] \approx \frac{1}{n_{s a m p l e}} \sum_{k = 1}^{n_{s a m p l e}} μ_{s a m p l e, k}

Where $μ_{s a m p l e, k}$ is the $k$ th sample from the posterior distribution - essentially we're just taking the population mean from our sample.

Bayesian applications, the next lecture, goes into more detail about libraries used for Bayesian analysis.

Prior density

In order for Bayesian stats to work, the "Prior distribution" needs to be known
There are debates about whether this is appropriate for objective science
However, for basic applications, we don't need to worry about it too much
- The default prior probabilities for fixed effects, variance, and Random Effects are well justified and available as default options in Bayesian software like R