What is Bayesian inference?

A reminder on conditional probabilities

$\Pr(A \mid B)$: Probability of A given B
The ordering matters: $\Pr(A \mid B)$ is not the same as $\Pr(B \mid A)$.
$\Pr(A \mid B) = \displaystyle{\frac{\Pr(A \cap B)}{\Pr(B)}}$

Screening for vampirism

The chance of the test being positive given you are a vampire is $\Pr(+|\text{vampire}) = 0.90$ (sensitivity).
The chance of a negative test given you are mortal is $\Pr(-|\text{mortal}) = 0.95$ (specificity).

What is the question?

From the perspective of the test: Given a person is a vampire, what is the probability that the test is positive? $\Pr(+|\text{vampire}) = 0.90$.
From the perspective of a person: Given that the test is positive, what is the probability that this person is a vampire? $\Pr(\text{vampire}|+) = \; ?$
Assume that vampires are rare, and represent only $0.1\%$ of the population. This means that $\Pr(\text{vampire}) = 0.001$.

What is the answer? Bayes’ theorem to the rescue!

$\Pr(\text{vampire}|+) = \displaystyle{\frac{\Pr(\text{vampire and } +)}{\Pr(+)}}$
$\Pr(\text{vampire and } +) = \Pr(\text{vampire}) \; \Pr(+ | \text{vampire}) = 0.0009$
$\Pr(+) = 0.0009 + 0.04995 = 0.05085$
$\Pr(\text{vampire}|+) = 0.0009/0.05085 = 0.02$

\[\Pr(\text{vampire}|+)= \displaystyle{\frac{ \Pr(+|\text{vampire}) \; \Pr(\text{vampire})}{\Pr(+)}}\]

let’s have an example

Screening for vampirism

Suppose the diagnostic test has the same sensitivity and specificity but vampirism is more common: $10\%$ of the population is vampire.
What is the probability that a person is a vampire, given that the test is positive?

Solution

The probability that a person is a vampire, given that the test is positive

$\Pr(+ | \text{vampire}) = 0.9$
$\Pr(- | \text{mortal}) = 0.95$
$\Pr(\text{vampire}) = 0.1$

\[\begin{aligned} \Pr(+) &= \Pr(+ | \text{vampire}) \Pr(\text{vampire}) + \Pr(+ | \text{mortal}) \Pr(\text{mortal}) \\ &= 0.9*0.1 + 0.05*0.9 \\ &=0.135 \end{aligned} \]

\[\begin{aligned} \Pr(\text{vampire} | +) &= \Pr(+ | \text{vampire}) \Pr(\text{vampire}) / \Pr(+) \\ &= 0.9*0.1 / 0.135 \end{aligned} \]

Bayes’ theorem

A theorem about conditional probabilities.
$\Pr(B \mid A) = \displaystyle{\frac{ \Pr(A \mid B) \; \Pr(B)}{\Pr(A)}}$

Easy to mess up with letters. Might be easier to remember when written like this:

\[ \Pr(\text{hypothesis} \mid \text{data}) = \frac{ \Pr(\text{data} \mid \text{hypothesis}) \; \Pr(\text{hypothesis})}{\Pr(\text{data})} \]

The “hypothesis” is typically something unobserved or unknown. It’s what you want to learn about using the data.
For regression models, the “hypothesis” is a parameter (intercept, slopes or error terms).
Bayes theorem tells you the probability of the hypothesis given the data.

Frequentist versus Bayesian

Typical stats problems involve estimating parameter $\theta$ with available data.
The frequentist approach (maximum likelihood estimation – MLE) assumes that the parameters are fixed, but have unknown values to be estimated.
Classical estimates generally provide a point estimate of the parameter of interest.
The Bayesian approach assumes that the parameters are not fixed but have some fixed unknown distribution - a distribution for the parameter.

What is the Bayesian approach?

The approach is based upon the idea that the experimenter begins with some prior beliefs about the system.
And then updates these beliefs on the basis of observed data.
This updating procedure is based upon the Bayes’ Theorem:

\[\Pr(A \mid B) = \frac{\Pr(B \mid A) \; \Pr(A)}{\Pr(B)}\]

Schematically if $A = \theta$ and $B = \text{data}$, then
The Bayes’ theorem

\[\Pr(A \mid B) = \frac{\Pr(B \mid A) \; \Pr(A)}{\Pr(B)}\]

Translates into:

\[\Pr(\theta \mid \text{data}) = \frac{\Pr(\text{data} \mid \theta) \; \Pr(\theta)}{\Pr(\text{data})}\]

Bayes’ theorem

\[{\color{red}{\Pr(\theta \mid \text{data})}} = \frac{\color{blue}{\Pr(\text{data} \mid \theta)} \; \color{green}{\Pr(\theta)}}{\color{orange}{\Pr(\text{data})}}\]

$\color{red}{Posterior ~distribution}$: Represents what you know after having seen the data. The basis for inference, a distribution, possibly multivariate if more than one parameter ($\theta$).
$\color{blue}{Likelihood}$: We know that guy from before, same as in the MLE approach.
$\color{green}{Prior~ distribution}$: Represents what you know before seeing the data. The source of much discussion about the Bayesian approach.
$\color{orange}{\Pr(\text{data}) = \int L(\text{data} \mid \theta) \;\Pr(\theta) d\theta }$: Possibly high-dimensional integral, difficult if not impossible to calculate. This is one of the reasons why we need simulation (MCMC) methods - more soon.

\[ \color{dodgerblue}{\boxed{\underbrace{P(\theta \mid y)}_{\text{Posterior}} = \frac{\overbrace{P(y \mid \theta)}^{\text{Likelihood}} \cdot \overbrace{P(\theta)}^{\text{Prior}}}{\underbrace{P(y)}_{\text{Normalizing Constant}}}}} \]

Likelihood

Usually, when talking about probability distributions, we assume that we know the parameter values.
In the real world, it is usually the other way around.

A question of interest might be for example:

We have observed 3 births by a female during her 10 breeding attempts. What does this tell us about the true probability of getting a successful breeding attempt from this female? For the population?

We don’t know what the probability of a birth is.
But we can calculate the probability of getting our data for different values:

dbinom(x=3,size=10,prob=0.1)
#> [1] 0.05739563

We don’t know what the probability of a birth is.
But we can see what the probability of getting our data would be for different values:

dbinom(x=3,size=10,prob=0.9)
#> [1] 8.748e-06

We don’t know what the probability of a birth is.
But we can see what the probability of getting our data would be for different values:


dbinom(x=3,size=10,prob=0.25)
#> [1] 0.2502823

dbinom(x=3,size=10,prob=0.1)
#> [1] 0.05739563
dbinom(x=3,size=10,prob=0.9)
#> [1] 8.748e-06
dbinom(x=3,size=10,prob=0.25)
#> [1] 0.2502823

So we would be more likely to observe 3 births if the probability is 0.25 than 0.1 or 0.9.

The likelihood

This reasoning is so common in statistics that it has a special name:
The likelihood is the probability of observing the data under a certain model.
The data are known, we usually consider the likelihood as a function of the model parameters $\theta_1,\theta_2, \ldots, \theta_p$

\[L = P(\theta_1,\theta_2, \ldots, \theta_p \mid \text{data})\]

A Formula for the Likelihood Function

Let $f(x; \theta)$ denote the pdf of a random variable $X$ with associated parameter $\theta$. Suppose $X_1, X_2, \ldots , X_n$ are random samples from this distribution, and $x_1, x_2, \ldots , x_n$ are the corresponding observed values.

\[\color{dodgerblue}{\boxed{L(\theta \mid x_1, x_2, \ldots , x_n) = f(x_1; \theta) f(x_2; \theta) \ldots f(x_n; \theta) = \prod_{i=1}^n f(x_i; \theta).}}\]

Likelihood functions in R

We may create a function to calculate a likelihood:

lik.fun <- function(parameter){
  ll <- dbinom(x=3, size=10, prob=parameter)
  return(ll)
}

lik.fun(0.3)
#> [1] 0.2668279

lik.fun(0.6)
#> [1] 0.04246733

Maximize the likelihood (3 successes ot of 10 attempts)

The maximum of the likelihood is at value $0.3$

Maximum likelihood estimation

There is always a set of parameters that gives you the highest likelihood of observing the data, and this is the MLE.
These can be calculated using:
- Trial and error (not efficient!).
- Compute the maximum of a function by hand (rarely doable in practice).
- An iterative optimization algorithm: ?optim in R.

By hand:

compute MLE of $p$ from $Y \sim \text{Bin}(N=10,p)$ with $k=3$ successes

$P(Y=k) = {{k}\choose{N}} p^k (1-p)^{N-k} = L(p)$.
$\log(L(p)) = \text{cte} + k \log(p) + (N-k) \log(1-p)$.
We are searching for the maximum of $L$, or equivalently that of $\log(L)$.
Compute derivate w.r.t. $p$: $\displaystyle{{{d\log(L)}\over{dp}} = {{k}\over{p}} - {{(N-k)}\over{(1-p)}}}$.
Then solve $\displaystyle{{{d\log(L)}\over{dp}}=0}$; the MLE is $\displaystyle{\hat{p} = {{k}\over{N}}={{3}\over{10}}=0.3}$.
Here, the MLE is the proportion of observed successes.

Using a computer:

MLE of $p$ from $Y \sim \text{Bin}(N=10,p)$ with $k=3$ successes

lik.fun <- function(parameter) dbinom(x=3, size=10, prob=parameter)
# ?optimize
optimize(lik.fun,c(0,1),maximum=TRUE)
#> $maximum
#> [1] 0.3000157
#> 
#> $objective
#> [1] 0.2668279

Use optim when the number of parameters is $> 1$.

The Bayesian approach

The Bayesian starts off with a prior.
Now, the one thing we know about $\theta$ is that is a continuous random variable and that it lies between zero and one.
Thus, a suitable prior distribution might be the Beta defined on $[0,1]$.
What is the Beta distribution?

What is the Beta distribution?

\[ q(\theta \mid \alpha, \beta) = \frac{1}{\text{Beta}(\alpha, \beta)}{\theta^{\alpha-1}} {(1-\theta)^{\theta-1}} \]

with $\text{Beta}(\alpha, \beta) = \displaystyle{\frac{\Gamma(\alpha)\Gamma(\beta)}{\Gamma(\alpha+\beta)}}$ and $\Gamma(n) = (n-1)!$

Beta Beta Conjugate Prior

We assume a priori that (let $\theta =p$)

\[Prior: ~p\sim Beta(\alpha,\beta)\] so that $\Pr(p) = p^{\alpha-1} (1 - p)^{\beta-1}$

also :

\[Likelihood~:X|p \sim Binomial(n,p)\]

Then we have:

\[L(P) ={\color{blue}{\binom{n}{k}p^k(1-p)^{n-k}}}\]

\[P(p) = {\color{green}{\frac{p^{\alpha-1} (1 - p)^{\beta-1}}{Beta(\alpha,\beta)}}}\]

\[ \begin{aligned} {\color{red}{Pr(p \mid k)}} & \propto {\color{blue}{\binom{n}{k}p^k(1-p)^{n-k}}} \; {\color{green}{p^{a-1} (1 - p)^{\beta-1}}}\\ & \propto {p^{(\alpha+k)-1}} {(1-p)^{(\beta+n-k)-1}} \end{aligned} \]

That is:

\[ p \mid k \sim Beta(\alpha+k,\beta+n-k)\]

Take a Beta prior with a Binomial likelihood, you get a Beta posterior (conjugacy)

Elicitation

Often, one has a belief about the distribution of one’s data. You may think that your data come from a binomial distribution and in that case you typically know the $n$, the number of trials but you usually do not know $p$, the probability of success. Or you may think that your data come from a normal distribution. But you do not know the mean $\mu$ or the standard deviation $\sigma$ of the normal. Beside to knowing the distribution of one’s data, you may also have beliefs about the unknown $p$ in the binomial or the unknown mean $\mu$ in the normal.

Bayesians express their belief in terms of personal probabilities. These personal probabilities encapsulate everything a Bayesian knows or believes about the problem. But these beliefs must obey the laws of probability, and be consistent with everything else the Bayesian knows.

So a Bayesian will seek to express their belief about the value of $p$ through a probability distribution, and a very flexible family of distributions for this purpose is the beta family. A member of the beta family is specified by two parameters, $\alpha$ and $\beta$; we denote this as $p \sim \text{beta}(\alpha, \beta)$. The probability density function is

\[\begin{equation} f(p) = \frac{\Gamma(\alpha+\beta)}{\Gamma(\alpha)\Gamma(\beta)} p^{\alpha-1} (1-p)^{\beta-1} \end{equation}\] where $0 \leq p \leq 1, \alpha>0, \beta>0$, and $\Gamma$ is a factorial:

\[\Gamma(n) = (n-1)! = (n-1) \times (n-2) \times \cdots \times 1\] Here is a summary of the key ideas:

Bayesians express their uncertainty through probability distributions.
One can think about the situation and self-elicit a probability distribution that approximately reflects his/her personal probability.
One’s personal probability should change according Bayes’ rule, as new data are observed.
The beta family of distribution can describe a wide range of prior beliefs.

Conjugacy

Next, let’s introduce the concept of conjugacy in Bayesian statistics.

Suppose we have the prior beliefs about the data as below:

Binomial distribution $\text{Bin}(n,p)$ with $n$ known and $p$ unknown
Prior belief about $p$ is $\text{beta}(\alpha,\beta)$

Then we observe $x$ success in $n$ trials, and it turns out the Bayes’ rule implies that our new belief about the probability density of $p$ is also the beta distribution, but with different parameters. In mathematical terms,

\[\begin{equation} p|x \sim \text{beta}(\alpha+x, \beta+n-x). \end{equation}\]

This is an example of conjugacy. Conjugacy occurs when the posterior distribution is in the same family of probability density functions as the prior belief, but with new parameter values, which have been updated to reflect what we have learned from the data.

Why are the beta binomial families conjugate? Here is a mathematical explanation.

Recall the discrete form of the Bayes’ rule:

\[P(A_i|B) = \frac{P(B|A_i)P(A_i)}{\sum^n_{j=1}P(B|A_j)P(A_j)}\]

However, this formula does not apply to continuous random variables, such as the $p$ which follows a beta distribution, because the denominator sums over all possible values (must be finitely many) of the random variable.

But the good news is that the $p$ has a finite range – it can take any value only between 0 and 1. Hence we can perform integration, which is a generalization of the summation. The Bayes’ rule can also be written in continuous form as:

\[\pi^*(p|x) = \frac{P(x|p)\pi(p)}{\int^1_0 P(x|p)\pi(p) dp}.\]

This is analogus to the discrete form, since the integral in the denominator will also be equal to some constant, just like a summation. This constant ensures that the total area under the curve, i.e. the posterior density function, equals 1.

Note that in the numerator, the first term, $P(x|p)$, is the data likelihood – the probability of observing the data given a specific value of $p$. The second term, $\pi(p)$, is the probability density function that reflects the prior belief about $p$.

In the beta-binomial case, we have $P(x|p)=\text{Bin}(n,p)$ and $\pi(p)=\text{beta}(\alpha,\beta)$.

Plugging in these distributions, we get

$$ \[\begin{aligned} \pi^*(p|x) &= \frac{1}{\text{some number}} \times P(x|p)\pi(p) \\ &= \frac{1}{\text{some number}} [\left( \begin{array}{c} n \\ x \end{array} \right) p^x (1-p)^{n-x}] [\frac{\Gamma(\alpha+\beta)}{\Gamma(\alpha)\Gamma(\beta)} p^{\alpha-1} (1-p)^{\beta-1}] \\ &= \frac{\Gamma(\alpha+\beta+n)}{\Gamma(\alpha+x)\Gamma(\beta+n-x)} \times p^{\alpha+x-1} (1-p)^{\beta+n-x-1} \end{aligned}\]

Let $\alpha^* = \alpha + x$ and $\beta^* = \beta+n-x$, and we get

\[\pi^*(p|x) = \text{beta}(\alpha^*,\beta^*) = \text{beta}(\alpha+x, \beta+n-x),\]

same as the posterior formula in Equation.

We can recognize the posterior distribution from the numerator $p^{\alpha+x-1}$ and $(1-p)^{\beta+n-x-1}$. Everything else are just constants, and they must take the unique value, which is needed to ensure that the area under the curve between 0 and 1 equals 1. So they have to take the values of the beta, which has parameters $\alpha+x$ and $\beta+n-x$.

This is a cute trick. We can find the answer without doing the integral simply by looking at the form of the numerator.

Without conjugacy, one has to do the integral. Often, the integral is impossible to evaluate. That obstacle is the primary reason that most statistical theory in the 20th century was not Bayesian. The situation did not change until modern computing allowed researchers to compute integrals numerically.

Three Conjugate Families

In this section, the three conjugate families are beta-binomial, gamma-Poisson, and normal-normal pairs. Each of them has its own applications in everyday life.

The Gamma-Poisson Conjugate Families

The Poisson distribution has a single parameter $\lambda$, and it is denoted as $X \sim \text{Pois}(\lambda)$ with $\lambda>0$. The probability mass function is

\[P(X=k) = \frac{\lambda^k}{k!} e^{-\lambda} \text{ for } k=0,1,\cdots,\]

where $k! = k \times (k-1) \times \cdots \times 1$. This gives the probability of observing a random variable equal to $k$.

Note that $\lambda$ is both the mean and the variance of the Poisson random variable. It is obvious that $\lambda$ must be greater than zero, because it represents the mean number of counts, and the variance should be greater than zero (except for constants, which have zero variance).

The gamma family is flexible, and Figure below illustrates a wide range of gamma shapes.

The probability density function for the gamma is indexed by shape $k$ and scale $\theta$, denoted as $\text{Gamma}(k,\theta)$ with $k,\theta > 0$. The mathematical form of the density is

\[\begin{equation} f(x) = \dfrac{1}{\Gamma(k)\theta^k} x^{k-1} e^{-x/\theta} \end{equation}\] where

\[\begin{equation} \Gamma(z) = \int^{\infty}_0 x^{z-1} e^{-x} dx. \end{equation}\] $\Gamma(z)$, the gamma function, is simply a constant that ensures the area under curve between 0 and 1 sums to 1, just like in the beta probability distribution case of Equation @ref(eq:beta). A special case is that $\Gamma(n) = (n-1)!$ when $n$ is a positive integer.

However, some books parameterize the gamma distribution in a slightly different way with shape $\alpha = k$ and rate (inverse scale) $\beta=1/\theta$:

\[f(x) = \frac{\beta^{\alpha}}{\Gamma(\alpha)} x^{\alpha-1} e^{-\beta x}\]

For this example, we use the $k$-$\theta$ parameterization, but you should always check which parameterization is being used. For example, $\mathsf{R}$ uses the $\alpha$-$\beta$ parameterization by default.
In the the later material we find that using the rate parameterization is more convenient.

For the gamma Poisson conjugate family, suppose we observed data $x_1, x_2, \cdots, x_n$ that follow a Poisson distribution.Then similar to the previous section, we would recognize the kernel of the gamma when using the gamma-Poisson family. The posterior $\text{Gamma}(k^*, \theta^*)$ has parameters

\[k^* = k + \sum^n_{i=1} x_i \text{ and } \theta^* = \frac{\theta}{(n\theta+1)}.\]

The Normal-Normal Conjugate Families

There are other conjugate families, and one is the normal-normal pair. If your data come from a normal distribution with known variance $\sigma^2$ but unknown mean $\mu$, and if your prior on the mean $\mu$, has a normal distribution with self-elicited mean $\nu$ and self-elicited variance $\tau^2$, then your posterior density for the mean, after seeing a sample of size $n$ with sample mean $\bar{x}$, is also normal. In mathematical notation, we have

\[\begin{aligned} x|\mu &\sim N(\mu,\sigma^2) \\ \mu &\sim N(\nu, \tau^2) \end{aligned}\]

As a practical matter, one often does not know $\sigma^2$, the standard deviation of the normal from which the data. But there are cases in which it is reasonable to treat the $\sigma^2$ as known.

For the normal-normal conjugate families, assume the prior on the unknown mean follows a normal distribution, i.e. $\mu \sim N(\nu, \tau^2)$. We also assume that the data $x_1,x_2,\cdots,x_n$ are independent and come from a normal with variance $\sigma^2$.

Then the posterior distribution of $\mu$ is also normal, with mean as a weighted average of the prior mean and the sample mean. We have

\[\mu|x_1,x_2,\cdots,x_n \sim N(\nu^*, \tau^{*2}),\]

where

\[\nu^* = \frac{\nu\sigma^2 + n\bar{x}\tau^2}{\sigma^2 + n\tau^2} \text{ and } \tau^* = \sqrt{\frac{\sigma^2\tau^2}{\sigma^2 + n\tau^2}}.\]

Example: Numerical Example

Suppose we have a sample of size $n=10$ with sample $\overline{X}=5$ , known variance $\sigma^2=4$ , prior mean $\mu_0=3$ and prior variance $\tau^2=2$. The posterior mean and variance are:

Credible and Confidence Intervals

credible intervals, the Bayesian alternative to confidence intervals. Confidence intervals,are the frequentist way to express uncertainty about an estimate of a population mean, a population proportion or some other parameter.

A confidence interval has the form of an upper and lower bound.

\[L, U = \text{pe} \pm \text{se} \times \text{cv}\]

where $L$, $U$ are the lower bound and upper bound of the confidence interval respectively, $\text{pe}$ represents “point estimates”, $\text{se}$ is the standard error, and $\text{cv}$ is the critical value.

Most importantly, the interpretation of a 95% confidence interval on the mean is that “95% of similarly constructed intervals will contain the true mean”, not “the probability that true mean lies between $L$ and $U$ is 0.95”.

The reason for this frequentist wording is that a frequentist may not express his uncertainty as a probability. The true mean is either within the interval or not, so the probability is zero or one. The problem is that the frequentist does not know which is the case.

On the other hand, Bayesians have no such qualms. It is fine for us to say that “the probability that the true mean is contained within a given interval is 0.95”. To distinguish our intervals from confidence intervals, we call them credible intervals.

Definitin 1

A confidence interval of X% for a parameter is an interval $(a, b)$ generated by a repeated sampling procedure has probability X% of containing the true value of the parameter, for all possible values of the parameter.

Definition 2

Say you performed a statistical analysis to compare the efficacy of a public policy between two groups and you obtain a difference between the mean of these groups.You can express this difference as a confidence interval. Often we choose 95% confidence.This means that $\textbf{95 studies out of 100}$, that uses the $\textbf{same sample size and target population}$, performing the $\textbf{same statistical test}$,will expect to find a result of the mean difference between groups inside the confidence interval.