Probability and likelihood

Rules of probability, examples with a 10-side die:

• The probability of any single outcome: \(Pr(X = x) = \frac{1}{10}\) for \(x = 1, 2, \dots, 10\).
• The probability of an event (e.g., rolling an even number): \(Pr(X \text{ is even}) = \frac{5}{10} = \frac{1}{2}\).

Probability and likelihood

Rules of probability, examples with a 10-side die:

• The probability of any single outcome: \(Pr(X = x) = \frac{1}{10}\) for \(x = 1, 2, \dots, 10\).
• The probability of an event (e.g., rolling an even number): \(Pr(X \text{ is even}) = \frac{5}{10} = \frac{1}{2}\).

That even number could be rewritten as a bunch of “or” statements: \[Pr(X \text{ is even}) = Pr(X=2 \text{ or } X=4 \text{ or } X=6 \text{ or } X=8 \text{ or } X=10)\] \[\quad = \frac{1}{10} + \frac{1}{10} + \frac{1}{10} + \frac{1}{10} + \frac{1}{10}\] \[\quad = \frac{5}{10}\] \[\quad = \frac{1}{2}\]

Probability and likelihood

What about 2 and 4 and 6 and 8 and 10 as successive rolls independent rolls?

The probability of rolling 2 and 4 and 6 and 8 and 10 on successive independent rolls is: \(Pr(X_1=2 \text{ and } X_2=4 \text{ and } X_3=6 \text{ and } X_4=8 \text{ and } X_5=10) = \frac{1}{10} \cdot \frac{1}{10} \cdot \frac{1}{10} \cdot \frac{1}{10} \cdot \frac{1}{10} = \frac{1}{100000}\)

Probability and likelihood

All probabilities have models underlying them. The probability of rolling a 2 on a 10-sided die is \(Pr(X=2) = \frac{1}{10}\) because the die is fair and has 10 sides.

We could write this as: \[Pr(X=2 \mid \text{10-sided die}) = \frac{1}{10}\]

A probability predicts the outcome of a random event given a model. A likelihood evaluates the plausibility of a model given an observed outcome - written as \(L(\theta \mid X)\).

Probability and likelihood

Likelihood - your relative level of suprise of an observation given a model.

Let \(X = \{2,1,1,4\}\)

\[L(\text{10-sided die} \mid X) = Pr(X \mid \text{10-sided die}) = \frac{1}{10} \cdot \frac{1}{10} \cdot \frac{1}{10} \cdot \frac{1}{10} = \frac{1}{10000}\]

\[L(\text{10-sided die} \mid X) = f(X \mid \text{10-sided die}) = \prod_{i=1}^{n} f(x_i \mid \text{10-sided die}) = \frac{1}{10000}\]

\[L(\theta \mid X) = f(X \mid \theta) = \prod_{i=1}^{n} f(x_i \mid \theta) = \frac{1}{10000}\]

Probability and likelihood

Likelihood - your relative level of suprise of an observation given a model.

Let \(X = \{2,1,1,4\}\) and \(\theta\) be a four-sided die.

\(L(\theta \mid X) = f(X \mid \theta) = \prod_{i=1}^{n} f(x_i \mid \theta) = \frac{1}{256}\)

Changing the model to a four-sided die increases the likelihood of observing \(X\) from \(\frac{1}{10000}\) to \(\frac{1}{256}\). It is likely that a four-sided die was used to generate the data.

Likelihood provide a way to estimate model parameters

A likelihood treats the data as fixed and \(\theta\) as the variable:

\[\mathcal{L}(\theta \mid X = k) = P(X = k \mid \theta)\]

Same formula — different question being asked.

# Observed: k = 3 mutations on a branch. Model: X ~ Poisson(λ).
lambda <- seq(0.01, 10, length.out = 500)  # range of possible λ values
lik    <- dpois(3, lambda = lambda)          # likelihood at each λ
plot(lambda, lik, type = "l", lwd = 2, col = "#0A2342",
     xlab = expression(lambda),
     ylab = expression(L(lambda ~ "|" ~ k == 3)),
     main = "Likelihood of \u03bb given 3 observed mutations")
abline(v = 3, lty = 2, col = "#F15025", lwd = 2)  # mark the MLE
legend("topright", legend = expression(hat(lambda)[MLE] == 3),
       lty = 2, col = "#F15025", bty = "n")

The MLE \(\hat\lambda = k = 3\) — the observed count is the best estimate of the Poisson rate.

Random variables & expectations

What is a random variable?

Random Variables & Expectations

A random variable (RV) is a variable whose possible values are numerical outcomes of a random phenomenon. 1. discrete RVs take on a countable number of distinct values (e.g., 0, 1, 2, …). 2. continuous RVs take on an uncountable number of values, typically intervals of real numbers (e.g., \([0,1]\)).

The behavior of discrete RVs is described by a probability mass function (PMF), which gives the probability of each possible value.

\(Pr(X = x) = f(x)\)

The behavior of continuous RVs is described by a probability density function (PDF), which gives probabilities over intervals of values. Any single value has probability 0, but the area under the curve over an interval gives the probability of falling in that interval.

\(Pr( a \le X \le b) = \int_a^b f(x)\,dx\)

Both discrete and continuous RVs have cumulative distribution functions (CDFs), which give the probability of falling below a given value:

\(F(x) = Pr(X \le x)\)

Expectation and variance

It is often useful to summarize a random variable with a single number. The mean of a distribution is also the expectation or expected value of a random variable \(X\):

For discrete RVs:

\[\bar{x} = \sum_x x\,Pr(X = x)\]

and the variance (or expected squared deviation from the mean) is given by:

\[\mathrm{V}_{x} = \sum_x (x - \bar{x})^2 Pr(X = x)\]

Variance measures the scatter of a variable around its mean. The standard deviation is the square root of the variance and often used when showing data since it is in the same units as the variable itself.

Expectation and variance

The notation of expectations is sometimes encountered in the wild.

Mean:

\[\mathbb{E}[X] = \bar{x}\]

Variance can apply linearity to achieve a shortcut:

\[\mathbb{E}[(X - \bar{X})^2] = \mathbb{E}[X^2] - (\mathbb{E}[X])^2 = \mathrm{Var}(X)\]

The mean and variance are derived from a moment generating function as the first and second moments of a distribution. Moment generating functions are beyond our scope, but you might encounter the term “moments” in the literature.

Expectation and variance

The mean and variance of continuous RVs:

\[\mathbb{E}[X] = \int_{-\infty}^{\infty} x\,f(x)\,dx\]

Variance: \[\mathrm{Var}(X) = \int_{-\infty}^{\infty} (x - \mathbb{E}[X])^2 f(x)\,dx\]

In practice though, we are not evaluating integrals analytically, but rather using numerical methods to approximate the mean and variance of continuous distributions. Approximation means sum!

Expectation and variance

The mean and variance of continuous RVs:

\[\hat{\bar{x}} = \frac{1}{n} \sum_{i=1}^{n} x_i\]

\[\hat{\mathrm{V}}_{x} = \frac{1}{n-1} \sum_{i=1}^{n} (x_i - \hat{\bar{x}})^2\]

Sneaking in the \(1/(n-1)\) term is a correction for the fact that we are estimating the population variance from a sample. This is to ensure our estimator is unbiased. if the difference between \(n\) and \(n-1\) ever matters to you, then you are probably up to no good anyway.

Random Variables & Expectations Summary

• A random variable \(X\) maps possible outcomes from a model to real numbers.
• Discrete RVs: probability mass function (PMF) \(Pr(X = k) \ge 0\), \(\sum_k Pr(X=k) = 1\).
• Continuous RVs: probability density function (PDF) \(f(x) \ge 0\), \(\int_{-\infty}^{\infty} f(x)\,dx = 1\).
• Expectation: \(\mathbb{E}[X] = \sum_k k\,P(X=k)\) or \(\int x\,f(x)\,dx\).
• Variance: \(\mathrm{Var}(X) = \mathbb{E}[X^2] - \bigl(\mathbb{E}[X]\bigr)^2\).
• Linearity: \(\mathbb{E}[aX + bY] = a\,\mathbb{E}[X] + b\,\mathbb{E}[Y]\) always (even if \(X, Y\) are dependent).

Statistical independence: AND means multiply

If events \(A\) and \(B\) are independent, observing one gives no information about the other:

\[P(A \text{ and } B) = P(A) \times P(B)\]

Example — probability that two randomly drawn alleles are both \(A_1\) at frequency \(p\):

\[P(\text{draw } A_1) \times P(\text{draw } A_1) = p \times p = p^2\]

This is the Hardy–Weinberg homozygote frequency.

Statistical independence: OR means add

For mutually exclusive events (cannot both occur):

\[P(A \text{ or } B) = P(A) + P(B)\]

For non-mutually exclusive events, subtract the overlap:

\[P(A \text{ or } B) = P(A) + P(B) - P(A \text{ and } B)\]

Example — probability that a diploid is a heterozygote (one copy of \(A_1\), one of \(A_2\)):

\[P(A_1 A_2) + P(A_2 A_1) = pq + qp = 2pq\]

The two allele-draw orderings are mutually exclusive events, so we add.

Joint and marginal probabilities

A joint probability \(P(A, B)\) assigns probability to combinations of two variables simultaneously.

Marginal probabilities are recovered by summing over the other variable — “collapsing” one dimension:

\[P(A) = \sum_b P(A,\, B = b) \qquad \text{(law of total probability)}\]

Conditional probability links joint and marginal:

\[P(A \mid B) = \frac{P(A,\, B)}{P(B)}\]

Rearranging gives Bayes’ theorem: \(P(A \mid B) \propto P(B \mid A)\,P(A)\), the foundation of every Bayesian inference method in the course.

A joint probability table

Pollen production level crossed with season:

• \(P(\text{High pollen}) = 0.30 + 0.10 = \mathbf{0.40}\) — sum across columns.
• \(P(\text{Wet season}) = 0.30 + 0.20 = \mathbf{0.50}\) — sum across rows.
• \(P(\text{High} \mid \text{Wet}) = 0.30\,/\,0.50 = \mathbf{0.60}\) — conditional on column.

Are pollen level and season independent?

\(P(\text{High}) \times P(\text{Wet}) = 0.40 \times 0.50 = 0.20 \ne 0.30\).

No — they are associated.

Binomial: definition

\(X \sim \text{Binomial}(n,\, p)\) — count of successes in \(n\) independent Bernoulli trials, each with success probability \(p\).

\[P(X = k) = \binom{n}{k} p^k (1-p)^{n-k}, \quad k = 0, 1, \ldots, n\]

• Parameters: \(n \in \mathbb{Z}^+\) (trials), \(p \in [0,1]\) (success probability).
• Support: \(\{0, 1, \ldots, n\}\).

Biological reading: \(n\) = gene copies sampled from a population, \(k\) = copies that are allele \(A_1\), \(p\) = allele frequency. One round of Wright–Fisher drift is one binomial draw.

Binomial: E[X] and Var(X)

Indicator variable trick. Let \(X_i = 1\) if trial \(i\) succeeds, 0 otherwise. Then:

\[X = X_1 + X_2 + \cdots + X_n, \qquad \mathbb{E}[X_i] = p\]

By linearity of expectation:

\[\boxed{\mathbb{E}[X] = \sum_{i=1}^{n} \mathbb{E}[X_i] = np}\]

Since trials are independent, \(\mathrm{Var}(X_i) = \mathbb{E}[X_i^2] - p^2 = p - p^2 = p(1-p)\), so:

\[\boxed{\mathrm{Var}(X) = np(1-p)}\]

Variance is maximised at \(p = 0.5\) and collapses to zero when \(p \in \{0, 1\}\).

Binomial in R

par(mfrow = c(1, 3), mar = c(4, 4, 3, 1))  # 3 panels side by side
ns <- c(10, 20, 50)  # number of trials to compare
p  <- 0.3            # success probability
for (n in ns) {
  k <- 0:n  # all possible count outcomes
  barplot(dbinom(k, size = n, prob = p), names.arg = k,
          main = paste0("Binom(n=", n, ", p=0.3)"),
          xlab = "k", ylab = "P(X = k)",
          col = "#2CA58D", border = NA, cex.names = 0.6)
}

Increasing \(n\) concentrates the distribution around \(\mathbb{E}[X] = np\); shape approaches Normal for large \(n\) (Central Limit Theorem).

Poisson: definition

\(X \sim \text{Poisson}(\lambda)\) — count of independent events in a fixed interval when events arrive at constant rate \(\lambda\).

\[P(X = k) = \frac{e^{-\lambda}\,\lambda^k}{k!}, \quad k = 0, 1, 2, \ldots\]

• Parameter: \(\lambda > 0\) (rate \(\times\) time = expected count).
• Support: \(\{0, 1, 2, \ldots\}\) — unbounded above.

Biological reading: \(\lambda = \mu t\) where \(\mu\) is the per-site substitution rate and \(t\) is branch length. The number of mutations on a branch \(\sim \text{Poisson}(\mu t)\).

Poisson: E[X] and Var(X)

\[\mathbb{E}[X] = \sum_{k=0}^{\infty} k \cdot \frac{e^{-\lambda}\lambda^k}{k!} = e^{-\lambda} \sum_{k=1}^{\infty} \frac{\lambda^k}{(k-1)!}\]

Let \(j = k - 1\):

\[= e^{-\lambda}\,\lambda \sum_{j=0}^{\infty} \frac{\lambda^j}{j!} = e^{-\lambda}\,\lambda\,e^{\lambda}\]

\[\boxed{\mathbb{E}[X] = \lambda}\]

Using \(\mathbb{E}[X(X-1)] = \lambda^2\) (same re-indexing, \(j = k-2\)) gives \(\mathbb{E}[X^2] = \lambda^2 + \lambda\), so:

\[\boxed{\mathrm{Var}(X) = \lambda}\]

Key property: mean \(=\) variance \(= \lambda\). Overdispersion (\(\mathrm{Var} > \mathrm{Mean}\)) signals a non-Poisson process.

Poisson in R

par(mfrow = c(1, 3), mar = c(4, 4, 3, 1))  # 3 panels side by side
lambdas <- c(1, 5, 15)  # three rate values to compare
for (lam in lambdas) {
  # cover the realistic range: mean ± 4 standard deviations
  k <- 0:ceiling(lam + 4 * sqrt(lam))
  barplot(dpois(k, lambda = lam), names.arg = k,
          main = bquote(Poisson(lambda == .(lam))),
          xlab = "k", ylab = "P(X = k)",
          col = "#D5A021", border = NA, cex.names = 0.6)
}

As \(\lambda\) grows the Poisson converges to a Normal with mean and variance both equal to \(\lambda\).

Exponential: definition

\(X \sim \text{Exponential}(\lambda)\) — continuous waiting time until the first event of a Poisson(\(\lambda\)) process.

\[f(x) = \lambda\,e^{-\lambda x}, \quad x \ge 0 \qquad F(x) = 1 - e^{-\lambda x}\]

• Parameter: \(\lambda > 0\) (rate; larger \(\lambda\) → shorter expected waits).
• Memoryless: \(P(X > s + t \mid X > s) = P(X > t)\) — the process has no memory of how long it has waited.

Biological reading: time until coalescence of two lineages in a population of size \(N_e\) is \(\text{Exponential}(1/N_e)\) (in units of generations). Waiting time to a mutation at rate \(\mu\) is \(\text{Exponential}(\mu)\).

Exponential: E[X] and Var(X)

\[\mathbb{E}[X] = \int_0^{\infty} x\,\lambda\,e^{-\lambda x}\,dx\]

Integration by parts with \(u = x\), \(dv = \lambda e^{-\lambda x}dx\):

\[= \Bigl[-x\,e^{-\lambda x}\Bigr]_0^{\infty} + \int_0^{\infty} e^{-\lambda x}\,dx = 0 + \frac{1}{\lambda}\]

\[\boxed{\mathbb{E}[X] = \frac{1}{\lambda}}\]

Applying the same method to \(\mathbb{E}[X^2]\) (with \(u = x^2\)) gives \(\mathbb{E}[X^2] = 2/\lambda^2\), so:

\[\boxed{\mathrm{Var}(X) = \frac{1}{\lambda^2}}\]

Exponential in R

par(mfrow = c(1, 3), mar = c(4, 4, 3, 1))  # 3 panels side by side
rates <- c(0.5, 1, 3)            # three rate values to compare
x     <- seq(0, 8, length.out = 400)  # x-axis values for the density curve
for (lam in rates) {
  plot(x, dexp(x, rate = lam), type = "l", lwd = 2, col = "#6F1A07",
       main = bquote(Exp(lambda == .(lam))),
       xlab = "x", ylab = "f(x)", ylim = c(0, 3.2))
  # dashed vertical line at E[X] = 1/lambda
  abline(v = 1/lam, lty = 2, col = "#0A2342", lwd = 1.5)
  # label the expected value in the legend
  legend("topright", legend = bquote(E[X] == .(round(1/lam, 2))),
         lty = 2, col = "#0A2342", bty = "n")
}

Dashed line marks \(\mathbb{E}[X] = 1/\lambda\). High rate → short expected wait; low rate → long tail.

Geometric: definition

\(X \sim \text{Geometric}(p)\) — number of trials until (and including) the first success, where each trial succeeds independently with probability \(p\).

\[P(X = k) = (1-p)^{k-1}\,p, \quad k = 1, 2, 3, \ldots\]

• Parameter: \(p \in (0,1)\) — success probability per trial.
• Discrete analogue of the Exponential; also memoryless.

Biological reading: in a haploid population of size \(N\), two lineages coalesce in any given generation with probability \(1/N\). The generation of coalescence \(\sim \text{Geometric}(1/N)\).

Geometric: E[X] and Var(X)

\[\mathbb{E}[X] = \sum_{k=1}^{\infty} k\,(1-p)^{k-1}\,p = p \cdot \frac{d}{dq}\!\left[\sum_{k=0}^{\infty} q^k\right]_{q=1-p}\]

Using the geometric series \(\sum_{k=0}^\infty q^k = \frac{1}{1-q}\) and differentiating:

\[\frac{d}{dq}\frac{1}{1-q} = \frac{1}{(1-q)^2} \implies \mathbb{E}[X] = p \cdot \frac{1}{p^2}\]

\[\boxed{\mathbb{E}[X] = \frac{1}{p}}\]

Using \(\mathbb{E}[X(X-1)]\) (same technique), \(\mathbb{E}[X^2] = (2-p)/p^2\), giving:

\[\boxed{\mathrm{Var}(X) = \frac{1-p}{p^2}}\]

Geometric in R

# R's dgeom(k, prob) counts *failures* before first success (support 0, 1, 2, …).
# Add 1 to shift to the "trial of first success" convention used above.
par(mfrow = c(1, 3), mar = c(4, 4, 3, 1))
probs <- c(0.1, 0.3, 0.7)
for (p in probs) {
  k <- 0:40
  barplot(dgeom(k, prob = p), names.arg = k + 1,
          main = bquote(Geom(p == .(p))),
          xlab = "k (trial of first success)", ylab = "P(X = k)",
          col = "#F15025", border = NA, cex.names = 0.55)
}

Normal: definition

\(X \sim \mathcal{N}(\mu,\, \sigma^2)\) — the bell curve; the limiting distribution of sums of independent random variables (Central Limit Theorem).

\[f(x) = \frac{1}{\sigma\sqrt{2\pi}}\exp\!\left(-\frac{(x-\mu)^2}{2\sigma^2}\right), \quad x \in \mathbb{R}\]

• Parameters: \(\mu \in \mathbb{R}\) (mean/location), \(\sigma^2 > 0\) (variance/spread).
• \(\mathbb{E}[X] = \mu\), \(\quad\mathrm{Var}(X) = \sigma^2\) — parameters are the moments.
• Symmetric around \(\mu\); 68 / 95 / 99.7% of mass within \(1/2/3\,\sigma\).

Biological reading: trait values under polygenic additive models; log-transformed allele frequencies near fixation; large-sample approximation to Binomial and Poisson.

Normal in R

par(mfrow = c(1, 3), mar = c(4, 4, 3, 1))  # 3 panels side by side
mus    <- c(0, 0, 3)   # mean values
sigmas <- c(1, 2, 1)   # standard deviations
labels <- c("N(0,1)", "N(0,4)", "N(3,1)")
x <- seq(-6, 9, length.out = 400)  # x-axis values
for (i in seq_along(labels)) {
  mu <- mus[i]     # mean
  s  <- sigmas[i]  # standard deviation
  plot(x, dnorm(x, mean = mu, sd = s), type = "l", lwd = 2, col = "#2CA58D",
       main = labels[i], xlab = "x", ylab = "f(x)", ylim = c(0, 0.45))
  # dashed line at the mean mu
  abline(v = mu, lty = 2, col = "#6F1A07", lwd = 1.5)
}

Dashed line marks \(\mu\). Changing \(\sigma\) scales the spread without shifting the centre.

Normal: polygenic traits

Many loci each contributing a small additive effect → trait values converge to Normal (CLT).

set.seed(1)
par(mfrow = c(1, 3), mar = c(4, 4, 3, 1))
for (L in c(5, 20, 100)) {
  effects <- matrix(sample(c(-1, 1), 5000 * L, replace = TRUE), nrow = 5000)
  trait   <- rowSums(effects)
  hist(trait, breaks = 40, col = "#2CA58D", border = NA, freq = FALSE,
       main = paste0(L, " additive loci"),
       xlab = "Trait value", ylab = "Density")
  curve(dnorm(x, 0, sqrt(L)), add = TRUE, col = "#6F1A07", lwd = 2)
}

Red curve: \(\mathcal{N}(0, L)\). With 100 loci the histogram is already nearly indistinguishable from the Normal — the basis of the infinitesimal model of quantitative genetics.

Log-Normal: definition

\(X \sim \text{LogNormal}(\mu,\, \sigma^2)\) — \(X > 0\) and \(\ln X \sim \mathcal{N}(\mu, \sigma^2)\).

\[f(x) = \frac{1}{x\,\sigma\sqrt{2\pi}}\exp\!\left(-\frac{(\ln x - \mu)^2}{2\sigma^2}\right), \quad x > 0\]

• Parameters: \(\mu, \sigma^2\) are the mean and variance of the log-scale.
• \(\mathbb{E}[X] = e^{\mu + \sigma^2/2}\), \(\quad\mathrm{Var}(X) = (e^{\sigma^2}-1)\,e^{2\mu+\sigma^2}\).
• Right-skewed; naturally bounded at zero; multiplicative processes generate it.

Biological reading: body sizes, gene expression levels, species abundances, substitution rates across sites — all tend to be log-normally distributed because they arise from multiplicative growth processes.

Log-Normal in R

par(mfrow = c(1, 3), mar = c(4, 4, 3, 1))  # 3 panels side by side
mus    <- c(0,    0,   1)    # log-scale means
var_s  <- c(0.25, 1,   0.5)  # log-scale variances
labels <- c("LN(0, 0.25)", "LN(0, 1)", "LN(1, 0.5)")
x <- seq(0.001, 10, length.out = 400)  # x-axis values (must be > 0)
for (i in seq_along(labels)) {
  mu <- mus[i]    # log-scale mean
  s2 <- var_s[i]  # log-scale variance
  plot(x, dlnorm(x, meanlog = mu, sdlog = sqrt(s2)), type = "l", lwd = 2,
       col = "#D5A021", main = labels[i], xlab = "x", ylab = "f(x)")
  # dashed line at E[X] = exp(mu + sigma^2/2)
  abline(v = exp(mu + s2/2), lty = 2, col = "#0A2342", lwd = 1.5)
}

Dashed line marks \(\mathbb{E}[X] = e^{\mu + \sigma^2/2}\). Increasing \(\sigma^2\) pulls the mean right and increases right-skew.

Log-Normal: Ks distributions after WGD

Synonymous divergence (Ks) between duplicated gene pairs is log-normally distributed — log-transforming reveals the underlying Normal.

set.seed(42)
mu_log <- -2   # log-scale mean
sd_log <- 0.8  # log-scale standard deviation
ks_sim <- rlnorm(800, meanlog = mu_log, sdlog = sd_log)  # simulate 800 Ks values
par(mfrow = c(1, 3), mar = c(4, 4, 3, 1))
hist(ks_sim, breaks = 50, col = "#D5A021", border = NA, freq = FALSE,
     main = "Ks (raw)", xlab = "Ks", ylab = "Density")
hist(log(ks_sim), breaks = 40, col = "#0A2342", border = NA, freq = FALSE,
     main = "log(Ks) — underlying Normal", xlab = "log(Ks)", ylab = "Density")
curve(dnorm(x, mu_log, sd_log), add = TRUE, col = "#F15025", lwd = 2)
plot(sort(ks_sim), seq_along(ks_sim) / length(ks_sim), type = "l", lwd = 2, col = "#D5A021",
     main = "CDF vs. Log-Normal theory", xlab = "Ks", ylab = "P(Ks ≤ x)")
curve(plnorm(x, mu_log, sd_log), add = TRUE, lty = 2, col = "#6F1A07", lwd = 2)
legend("bottomright", c("empirical", "theoretical"), lty = c(1,2),
       col = c("#D5A021", "#6F1A07"), bty = "n", cex = 0.8)

Left: raw Ks is right-skewed. Middle: log(Ks) is approximately Normal. Right: empirical CDF closely tracks Log-Normal theory. → Macro 03 (genome architecture & polyploidy).

Gamma: definition

\(X \sim \text{Gamma}(\alpha,\, \beta)\) — sum of \(\alpha\) independent \(\text{Exponential}(\beta)\) waiting times.

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

• Parameters: \(\alpha > 0\) (shape), \(\beta > 0\) (rate).
• \(\mathbb{E}[X] = \alpha/\beta\), \(\quad\mathrm{Var}(X) = \alpha/\beta^2\).
• Special cases: \(\text{Gamma}(1, \beta) = \text{Exponential}(\beta)\); \(\text{Gamma}(k/2, 1/2) = \chi^2_k\).

Biological reading: the Gamma model of rate variation across sites is central to phylogenetics — site rates are drawn from a Gamma distribution with shape parameter \(\alpha\), which controls how much rate heterogeneity there is. Small \(\alpha\) → high variance; \(\alpha \to \infty\) → constant rates.

Gamma in R

par(mfrow = c(1, 3), mar = c(4, 4, 3, 1))  # 3 panels side by side
alphas <- c(0.5, 2, 5)  # shape values to compare
betas  <- c(1,   1, 1)  # rate values (all 1 here)
labels <- c("Gamma(0.5, 1)", "Gamma(2, 1)", "Gamma(5, 1)")
x <- seq(0.001, 15, length.out = 500)  # x-axis values (must be > 0)
for (i in seq_along(labels)) {
  a <- alphas[i]  # shape alpha
  b <- betas[i]   # rate beta
  plot(x, dgamma(x, shape = a, rate = b), type = "l", lwd = 2, col = "#F15025",
       main = labels[i], xlab = "x", ylab = "f(x)", ylim = c(0, 1.5))
  # dashed line at E[X] = alpha/beta
  abline(v = a/b, lty = 2, col = "#0A2342", lwd = 1.5)
}

Dashed line marks \(\mathbb{E}[X] = \alpha/\beta\). The shape parameter \(\alpha\) controls skewness: small \(\alpha\) → highly right-skewed (many slow sites, few fast); large \(\alpha\) → near-Normal.

Gamma: among-site rate variation

The \(+\Gamma\) model draws a relative substitution rate for each site from Gamma\((\alpha, \alpha)\) (mean \(= 1\)). Small \(\alpha\) creates extreme rate heterogeneity; ignoring it biases branch-length estimates.

set.seed(7)
par(mfrow = c(1, 3), mar = c(4, 4, 3, 1))  # 3 panels side by side
alphas <- c(0.3, 1, 5)  # three shape values to compare
for (a in alphas) {
  # setting rate = alpha normalises the mean to 1 (mean = shape/rate)
  rates <- rgamma(2000, shape = a, rate = a)
  hist(rates, breaks = 60, col = "#F15025", border = NA, freq = FALSE,
       xlim = c(0, 5), main = bquote(alpha == .(a) ~ "(mean = 1)"),
       xlab = "Relative site rate", ylab = "Density")
  # overlay the theoretical Gamma density
  curve(dgamma(x, shape = a, rate = a), add = TRUE, col = "#0A2342", lwd = 2)
  # dashed line at the mean (= 1 by construction)
  abline(v = 1, lty = 2, col = "#6F1A07")
}

\(\alpha = 0.3\): most sites evolve slowly, a few extremely fast (hard to model without \(+\Gamma\)). \(\alpha = 5\): rates nearly uniform. → Speciation 05 (phylogenetics theory).

Beta: definition

\(X \sim \text{Beta}(\alpha,\, \beta)\) — a flexible distribution for probabilities and proportions bounded on \([0, 1]\).

\[f(x) = \frac{x^{\alpha-1}(1-x)^{\beta-1}}{B(\alpha,\beta)}, \quad x \in [0,1]\]

where \(B(\alpha,\beta) = \Gamma(\alpha)\Gamma(\beta)/\Gamma(\alpha+\beta)\) is the beta function.

• \(\mathbb{E}[X] = \frac{\alpha}{\alpha+\beta}\), \(\quad\mathrm{Var}(X) = \frac{\alpha\beta}{(\alpha+\beta)^2(\alpha+\beta+1)}\).
• \(\alpha = \beta = 1\) → Uniform\((0,1)\); \(\alpha = \beta > 1\) → symmetric bell; \(\alpha \ne \beta\) → skewed.
• Conjugate prior for the Binomial likelihood in Bayesian inference.

Biological reading: allele frequencies, heterozygosity, admixture proportions — any quantity naturally constrained to \([0,1]\).

Beta in R

par(mfrow = c(1, 3), mar = c(4, 4, 3, 1))  # 3 panels side by side
alphas <- c(0.5, 2, 5)  # alpha (shape1)
betas  <- c(0.5, 5, 2)  # beta  (shape2)
labels <- c("Beta(0.5, 0.5)", "Beta(2, 5)", "Beta(5, 2)")
x <- seq(0.001, 0.999, length.out = 400)  # x in (0,1)
for (i in seq_along(labels)) {
  a <- alphas[i]  # alpha (shape1)
  b <- betas[i]   # beta  (shape2)
  plot(x, dbeta(x, shape1 = a, shape2 = b), type = "l", lwd = 2, col = "#6F1A07",
       main = labels[i], xlab = "x", ylab = "f(x)")
  # dashed line at E[X] = alpha/(alpha+beta)
  abline(v = a/(a+b), lty = 2, col = "#0A2342", lwd = 1.5)
}

Dashed line marks \(\mathbb{E}[X] = \alpha/(\alpha+\beta)\). Beta(2,5) skews toward low allele frequencies; Beta(5,2) toward high.

Beta: allele frequency drift

Under neutral drift, the stationary distribution of allele frequencies is Beta\((\theta/2,\,\theta/2)\) where \(\theta = 4N_e\mu\). Low \(\theta\) → most sites near fixation or loss; high \(\theta\) → intermediate frequencies.

set.seed(3)
par(mfrow = c(1, 3), mar = c(4, 4, 3, 1))  # 3 panels side by side
for (theta in c(0.5, 2, 10)) {
  # stationary Beta distribution: both shape parameters = theta/2
  freqs <- rbeta(5000, theta/2, theta/2)
  hist(freqs, breaks = 40, col = "#6F1A07", border = NA, freq = FALSE,
       main = bquote(theta == .(theta)),
       xlab = "Allele frequency p", ylab = "Density",
       xlim = c(0, 1))
  # overlay the theoretical Beta density
  curve(dbeta(x, theta/2, theta/2), add = TRUE, col = "#0A2342", lwd = 2)
}

Left (\(\theta=0.5\)): most SNPs are near fixation or loss — typical of small populations. Right (\(\theta=10\)): many polymorphic sites at intermediate frequency. → Popgen 01 (genetic diversity) and Popgen 02 (population structure).

Uniform: definition

\(X \sim \text{Uniform}(a,\, b)\) — all values in \([a, b]\) are equally likely.

\[f(x) = \frac{1}{b-a}, \quad x \in [a, b]\]

• \(\mathbb{E}[X] = \frac{a+b}{2}\), \(\quad\mathrm{Var}(X) = \frac{(b-a)^2}{12}\).
• Maximum entropy distribution on a bounded interval — the least informative prior when only bounds are known.
• Discrete analogue: \(\text{Uniform}\{1, \ldots, n\}\) with \(\mathbb{E}[X] = (n+1)/2\).

Biological reading: non-informative (flat) Bayesian priors on bounded parameters, e.g. \(\pi \sim \text{Uniform}(0,1)\) for base frequency priors in phylogenetic models; random starting values in optimisation.

Uniform in R

par(mfrow = c(1, 2), mar = c(4, 4, 3, 1))  # 2 panels side by side

# Panel 1: continuous Uniform density
x <- seq(-0.1, 1.1, length.out = 400)
plot(x, dunif(x, min = 0, max = 1), type = "l", lwd = 2, col = "#2CA58D",
     main = "Uniform(0, 1) = Beta(1,1)", xlab = "x", ylab = "f(x)", ylim = c(0, 1.5))

# Panel 2: discrete Uniform by simulation
set.seed(42)
draws <- sample(1:6, 10000, replace = TRUE)  # roll a fair 6-sided die 10,000 times
barplot(table(draws)/10000, col = "#D5A021", border = NA,
        main = "Discrete Uniform{1,...,6}", xlab = "x", ylab = "Proportion")

Uniform\((0,1)\) is a special case of Beta\((1,1)\) — a flat prior that says all values in \([0,1]\) are equally plausible.

Evolution is a stochastic process

• A stochastic process = a random variable indexed by time.
• Poisson process — events at constant rate \(\lambda\); interarrival times \(\sim \text{Exponential}(\lambda)\); counts \(\sim \text{Poisson}(\lambda t)\).
• Markov chains — memoryless state change; the substitution CTMC (phylo) and drift.

Review Questions

• What is an RV?
• What is the difference between discrete and continuous RVs?
• What are the expectations of RVs?
• Can you apply basic AND and OR rules of probability?
• How is a likelihood different from a probability?
• Can you recall the binomial, Poisson, geometric, exponential distributions and their biological uses?