Bayesian Machine Learning

Philosophical Foundations and Bayes’ Theorem¶

Interpretations of probability¶

Frequentist definition: Probability of event $E$ is the long-run relative frequency of $E$ in infinitely many repetitions of the same experiment. Statements are about data, not parameters ($P(E)$ is defined via an idealized limit of repetitions).

Bayesian definition: Probability is a subjective degree of belief about events or parameters, given current knowledge ($P(E)$ quantifies how plausible you think $E$ is).

Key philosophical contrast¶

Frequentist: uncertainty is in the world (randomness of experiments).

Bayesian: uncertainty is in the head of the observer (state of information; different observers with different prior knowledge can have different probabilities).

Sampling vs. inference¶

Sampling (forward problem): Assume parameter is known, ask: “If $p$ is the probability of success, what is $P(Y = k \mid p)$?”. For example: Number of heads in 10 tosses with known $p$ ~ Binomial($n=10, p$).

Inference (inverse problem): Assume data is observed, but the parameters are unknown, ask: “Given I observed $k$ heads in 10 tosses, what do I believe about $p$?”

Frequentist: Estimate single “best” value (e.g. MLE $\hat p = k/n$)

Bayesian: derive a posterior distribution over $p$: $p(p \mid \text{data})$

Bayes’ theorem and roles of likelihood¶

Bayes’ theorem (discrete): $P(A \mid B) = \frac{P(B\mid A)P(A)}{P(B)}$

Belief-updater view: $\underbrace{P(A\mid d)}_{\text{posterior}}\propto\underbrace{P(d\mid A)}_{\text{likelihood}}\times\underbrace{P(A)}_{\text{prior}}$

Prior: beliefs before data.

Likelihood: how plausible data are under different parameter values.

Posterior: updated beliefs after data.

Likelihood in both schools:

In frequentist statistics, the likelihood function is used to perform Maximum Likelihood Estimation (MLE). This identifies the single parameter value that maximizes the probability of having produced the observed data.

In Bayesian statistics, the likelihood function serves as a belief updater. It is multiplied by the prior distribution to calculate the posterior, allowing you to estimate a full distribution of all plausible values rather than just one “best” number.

When Bayesian methods are preferable¶

Small sample sizes / rare events: frequentist asymptotics unreliable.

Strong prior information available (expert opinion, previous studies) and should be integrated.

Need full uncertainty quantification (posterior distributions, credible intervals) rather than point estimates or p-values.

Complex hierarchical or graphical models, where likelihood-based methods get complicated but Bayesian methods remain conceptually coherent.

The computational curse of Bayesian statistics¶

Core difficulty: the evidence (marginal likelihood) is an often intractable integral: $p(d) = \int p(d\mid \theta)p(\theta)\,d\theta$

In simple conjugate families this integral is analytic.

In realistic models (many parameters, non-Gaussian likelihoods), high‑dimensional integrals are impossible to compute in closed form.

This computational curse motivates Markov Chain Monte Carlo (MCMC) and Variational Inference.

Conjugate Families and Continuous Inference¶

Conjugate families are pedagogically useful to understand how priors and data trade off.

Beta-Binomial¶

Likelihood: $P(Y \mid \pi) \sim \text{Binomial}(n, \pi) \sim P(Y=k \mid \pi) = \binom{n}{k}\pi^k(1-\pi)^{n-k}$

$\mathbb{E}[Y] = n\pi \quad \operatorname{Var}[Y] = n\pi(1-\pi)$

$\operatorname{Mode}[Y] = \begin{cases} \lfloor (n+1)\pi \rfloor & \text{if } (n+1)\pi \notin \mathbb{Z} \\ (n+1)\pi, (n+1)\pi - 1 & \text{if } (n+1)\pi \in \mathbb{Z} \end{cases}$

$\binom{n}{k}$ is the binomial coefficient, representing the number of ways to choose $k$ successes from $n$ trials: $\binom{n}{k} = \frac{n!}{k!(n-k)!}$

$\pi^k$ is the probability of getting exactly $k$ successes

$(1-\pi)^{n-k}$ is the probability of getting exactly $n-k$ failures

Prior: $P(\pi) \sim \text{Beta}(\pi \mid \alpha,\beta) = \frac{\Gamma(\alpha + \beta)}{\Gamma(\alpha)\Gamma(\beta)} \pi^{\alpha-1}(1-\pi)^{\beta-1}$

with: $\Gamma(n) = (n - 1)! \quad \forall n \in \mathbb{N}$

$\mathbb{E}[\pi] = \frac{\alpha}{\alpha + \beta} \quad \operatorname{Var}[\pi] = \frac{\alpha\beta}{(\alpha+\beta)^2(\alpha+\beta+1)}$

$\operatorname{Mode}[\pi] = \dfrac{\alpha - 1}{\alpha + \beta - 2} \quad \text{if } \alpha, \beta > 1$

Posterior: $P(\pi \mid Y=k) = \text{Beta}(\alpha + k, \beta + n - k)$

Key properties: Prior “pseudo-counts”: $\alpha-1$ successes, $\beta-1$ failures

With more data, likelihood dominates; with little data, prior dominates

Order of data does not matter; only sufficient statistics $(n,k)$

Poisson–Gamma: For Counts¶

Likelihood: $P(Y \mid \lambda) \sim \text{Poisson}(\lambda) \sim P(Y = y \mid \lambda) = \frac{e^{-\lambda}\lambda^y}{y!}$

$\mathbb{E}[Y] = \lambda \quad \text{Var}[Y] = \lambda$

$\operatorname{Mode}[\lambda] = \begin{cases} \lfloor \lambda \rfloor & \text{if } \lambda \text{ is not an integer} \\ \lambda, \lambda - 1 & \text{if } \lambda \text{ is an integer} \end{cases}$

• $\lambda^y$ is the contribution from $y$ events occurring

• $e^{-\lambda}$ is the probability of no other events occurring

• $\frac{1}{y!}$ normalizes for the ordering of events

Prior: $P(\lambda) \sim \text{Gamma}(s, r) \sim \frac{r^s}{\Gamma(s)}\lambda^{s-1}e^{-r\lambda}$

with $\Gamma(s) = (s-1)! \quad \forall n \in \mathbb{N}$

$\mathbb{E}[\lambda] = \frac{s}{r} \quad \text{Var}[\lambda] = \frac{s}{r^2} \quad \operatorname{Mode}[\lambda] = \frac{s - 1}{r} \quad \text{for } s \geq 1$

Posterior: $P(\lambda \mid Y = y) \sim \text{Gamma}\left(s + \sum_i y_i, r + n\right)$

Key properties: Prior “pseudo-counts”: $s$ represents prior event count, $r$ represents prior observation periods.

With more data, likelihood dominates; with little data, prior dominates.

Order of data does not matter; only sufficient statistic $\sum y_i$.

Negative Binomial For Overdispersed Counts¶

Unlike the Poisson distribution, the variance is not equal to the mean; it is larger (overdispersion). As the dispersion parameter $\alpha \to \infty$, the Negative Binomial distribution converges to the Poisson distribution. The Negative Binomial likelihood does not have a simple closed-form conjugate posterior - posterior inference must be done via MCMC simulation.

Likelihood: $P(Y \mid r, p) \sim \text{NegBin}(r, p)$

$\mathbb{E}[k \mid r,p] = \frac{r(1-p)}{p} \quad \text{and} \quad \mathbb{E}[k \mid r, \lambda] = \lambda$

$\text{Var}[k \mid r,p] = \frac{r(1-p)}{p^2} \quad \text{and} \quad \text{Var}[k \mid r, \lambda] = \lambda\left(1 + \frac{\lambda}{r}\right)$

• $p$ probability of success

• $r$ is the number of successes (-1 because last success is fixed as its the stopping criterion)

• $k$ is the number of failures

Prior: Negative Binomial is typically used in models where conjugacy does not hold. Therefore, we often use:

• in simple models $\lambda \sim \text{Gamma}(s, r)$ or in regression $\log(\lambda) = \beta_0 + \beta_1 X$ with $\beta \sim \text{Normal}$

• $r \sim \text{Gamma}$ or $r \sim \text{Exponential}$ (must be positive)

Normal–Normal (variance needs to be known!)¶

Likelihood: $P(Y \mid \mu, \sigma) \sim \text{Normal}(\mu, \sigma^2) \sim \frac{1}{\sqrt{2\pi\sigma^2}}\exp\left(-\frac{(y-\mu)^2}{2\sigma^2}\right)$

$\mathbb{E}[Y] = \mu \quad \text{Var}[Y] = \sigma^2 \quad \operatorname{Mode}[Y] = \mu$

• $\mu$ is the mean (location parameter)

• $\sigma^2$ is the variance (assumed known)

• $(y-\mu)^2$ measures squared deviation from the mean

Prior: $P(\mu) \sim \text{Normal}(\mu_0, \tau_0^2)$

$\mathbb{E}[\mu] = \mu_0 \quad \text{Var}[\mu] = \tau_0^2$

• $\mu_0$ is the prior mean

• $\tau_0^2$ is the prior variance (reflects uncertainty about $\mu$)

Posterior: $P(\mu \mid Y = y) \sim \text{Normal}(\mu_n, \tau_n^2)$

Key properties: Posterior mean is a precision-weighted average of prior mean $\mu_0$ and data $y$.

Precision = $\frac{1}{\text{variance}}$: higher precision = more influence.

Posterior variance is always smaller than both prior and data variance (information accumulation).

With more data ($n$ observations), data precision $\frac{n}{\sigma^2}$ dominates.

Multivariate Problems¶

Covariance: Quantifying Relationship Between Variables¶

• $\text{Cov}(X, X) = \text{Var}[X] \geq 0$: A variable’s covariance with itself is its variance

• $\text{Cov}(X, Y) \approx 0$: $X$ and $Y$ are unrelated; on average there is no “co-variance”

For observed data (empirical), the sample covariance is:

• $\bar{x}$ and $\bar{y}$: the sample means

• $n$ is the number of observations (the denominator is $n-1$ to account for the loss of one degree of freedom when estimating the means)

Pearson Correlation = Standardized Covariance¶

Pearson correlation normalizes covariance by the standard deviations of both variables, producing a unitless measure bounded between 0 and 1. By standardizing covariance, correlation becomes interpretable on a common scale where $r = 1$ indicates perfect positive correlation, $r = 0$ indicates no correlation, and $r = -1$ indicates perfect negative correlation.

• $\sigma_X$ and $\sigma_Y$ are the standard deviations of $X$ and $Y$ respectively

For observed data (empirical), the Pearson Correlation is:

• $\hat{\sigma}_x$ and $\hat{\sigma}_y$ are the sample standard deviations

Multinomial–Dirichlet: For Categorical Counts (Generalized Beta-Binomial)¶

Use whenever we want to model experiments with > 2 mutually exclusive outcomes!

Likelihood: $P(\mathbf{K} \mid \boldsymbol{\pi}) \sim \text{Multinomial}(n, \boldsymbol{\pi})$

• $\prod_{i=1}^{r} \pi_i^{k_i}$ is the contribution from observed category counts

• $\frac{n!}{\prod_{i=1}^{r} k_i!}$ normalizes for the number of orderings

• Constraint: $\sum_{i=1}^{r} k_i = n$ and $\sum_{i=1}^{r} \pi_i = 1$

Prior: $P(\boldsymbol{\pi}) \sim \text{Dirichlet}(\boldsymbol{\alpha})$

with $\alpha_0 = \sum_{i=1}^{r} \alpha_i$

Posterior: $P(\boldsymbol{\pi} \mid \mathbf{K} = \mathbf{k}) \sim \text{Dirichlet}(\boldsymbol{\alpha} + \mathbf{k})$

Key properties: Prior “pseudo-counts”: $\alpha_i$ represents prior observations in category $i$. Setting $\alpha_i = 1$ (flat prior) means no prior preference.

With more data, likelihood dominates; with little data, prior dominates.

Order of data does not matter; only sufficient statistic $\mathbf{k} = (k_1, \ldots, k_r)$.

Conjugate family ensures posterior is also Dirichlet with simple parameter update.

Multivariate Normal¶

• $x\in\mathbb R^d$, $x\sim\mathcal N(\mu,\Sigma)$ with mean $\mu \in \mathbb{R}^d$ and covariance matrix $\Sigma \in \mathbb{R}^{d\times d}$

• Closed under linear transformations, addition, marginalization, conditioning

• Covariance controls scale and correlation; must be symmetric and positive definite

• LKJ prior for correlation matrices: flexible prior on correlation structures; together with priors on standard deviations gives LKJCov prior for $\Sigma$

Multinomial–Dirichlet¶

• Likelihood: counts over $K$ categories, $Y\sim\text{Multinomial}(n,\boldsymbol\pi)$.

• Prior: Dirichlet$(\boldsymbol\alpha)$, $\boldsymbol\alpha=(\alpha_1,\dots,\alpha_K)$.

• Posterior: Dirichlet$(\boldsymbol\alpha + \mathbf k)$, where $k_j$ counts in category $j$.

Used for compositional data (election results, topic proportions).

Markov Chain Monte Carlo and the Metropolis Algorithm¶

Goal: obtain samples from $p(\theta\mid d)$ when no closed form is available. MCMC constructs a Markov chain whose stationary distribution is the posterior.

Metropolis algorithm (random-walk variant)¶

For target density $\pi(\theta) = p(\theta\mid d)$ and symmetric proposal distribution $q(\theta' \mid \theta)$, the Metropolis steps are:

1. Initialize $\theta^{(0)}$ (e.g. some reasonable value).

2. For iterations $t = 0,1,\dots$:

3. Propose new state $\theta' \sim q(\theta' \mid \theta^{(t)})$. E.g. random walk: $q = \mathcal N(\theta^{(t)}, \sigma_{\text{prop}}^2)$.

4. Compute acceptance probability $\alpha = \min\left\{1,\;\frac{\pi(\theta')}{\pi(\theta^{(t)})}\right\} = \min\left\{1,\;\frac{p(d\mid \theta')p(\theta')}{p(d\mid \theta^{(t)})p(\theta^{(t)})}\right\}$

5. Accept or reject: Draw $u \sim \text{Uniform}(0,1)$, If $u < \alpha$, set $\theta^{(t+1)}=\theta'$ (move), Else set $\theta^{(t+1)}=\theta^{(t)}$ (stay).

6. After burn‑in and thinning, use remaining samples as approximate draws from the posterior.

Proposal distribution: the distribution from which proposed new parameters are drawn ($q(\theta' \mid \theta)$). Choice of its scale strongly affects mixing and acceptance rate.

MCMC with PyMC¶

Diagnosing convergence and mixing¶

Key visual and numeric diagnostics:

Trace plots: Plot each chain’s sampled values over iterations

For Metropolis specifically, problematic trace plots often indicate: Proposal variance too small (tiny steps, high acceptance, poor exploration). Proposal variance too large (proposals rejected, chain rarely moves). Tunable by adjusting $\sigma_{\text{prop}}$.

Chains explore the same region (“fuzzy caterpillars”)

No systematic upward or downward drift.

No chains stuck in a narrow band.

Different chains overlap well.

Run multiple chains: Posterior histograms or kernel density estimates for each chain separately

Distributions from all chains almost perfectly overlap.

No unique behavior in one chain.

Autocorrelation plots: Autocorrelation function (ACF) of each chain as a function of lag

Autocorrelation decays quickly with lag.

Very slow decay or persistent high correlation means poor mixing and low effective sample size.

Rank plots / rank bars: For each parameter, pool all samples from all chains, compute their ranks, then see how ranks are distributed per chain (as in ArviZ’s rank plots)

Ranks within each chain are roughly uniformly distributed.

All chains show similar, flat-ish rank patterns.

Strong deviations (U‑shape, spikes, or clearly different patterns per chain) indicate non‑convergence or poor mixing.

$\hat R$ (R‑hat, Gelman–Rubin statistic): Compares total chain data variance (concatenated) vs. the mean variance within individual chains. Summarizes the overlaid density plots of multiple chains and provides a numerical measure for their diagreement.

$\hat R \approx 1.00$ (e.g. < 1.01–1.05) means chains have likely converged to same target distribution.

$\hat R > 1.05$ means potential non‑convergence; increase iterations, needs better initialization, or reparameterization.

ESS (effective sample size): Estimates how many independent draws your autocorrelated chain is equivalent to

ESS large relative to the number of draws, giving stable posterior summaries. Rule of thumb: ESS / N > 10%

Very small ESS (especially compared to total draws) means strong autocorrelation, poor mixing, unreliable summaries.

ess_bulk (Bulk ESS): Use this to assess reliability of posterior means, medians, and central credible intervals. Estimates the effective sample size for the central part of the distribution (around the median/mean). Interpretation: Higher is better. If ess_bulk is very small relative to the total number of draws, the chain is autocorrelated or has poor mixing in the center of the distribution.

ess_tail (Tail ESS): Use this to assess reliability of HDI/credible interval bounds and tail probabilities. Estimates the effective sample size specifically for the tails of the distribution (extreme quantiles, typically 5th and 95th percentiles). Interpretation: Higher is better. Can be substantially lower than ess_bulk if the tails are poorly explored.

Posterior Summaries, Prediction, and Uncertainty¶

Point Summaries: Mean, Variance, Maximum A Posteriori (MAP / Posterior Mode)¶

Given posterior samples $\theta_1,\dots,\theta_N$:

Posterior mean: $\hat\theta_{\text{mean}} = \frac{1}{N}\sum_{i=1}^N \theta_i$

Posterior variance: $\widehat{\text{Var}}(\theta) = \frac{1}{N-1}\sum_{i=1}^N(\theta_i - \hat\theta_{\text{mean}})^2$

Maximum A Posteriori (MAP) or Posterior Mode: $\hat{\theta}_{\text{MAP}}=\arg\max_{\theta} P(\theta \mid D) = \arg\max_{\theta} P(D\mid\theta) P(\theta) = \log P(D\mid\theta) + \log P(\theta)$; To compute: Use closed form Posterior Mode if available; or Differentiate with respect to $\theta$, set the derivative to zero, and solve for $\theta$; or use numerical optimization when there is no closed form.

Credible Intervals vs. Confidence Intervals¶

Credible interval (Bayesian): For posterior $p(\theta\mid d)$, an 80% credible interval $[a,b]$ satisfies: $P(a \le \theta \le b \mid d) = 0.8$. It’s a direct statement about the parameter: “Given data and model, we believe with 80% plausibility that $\theta$ lies in $[a,b]$.”

HDI (Highest Density Interval): Among all 80% intervals, choose the one where posterior density is everywhere at least as high as outside. It minimizes width for a given mass; It contains the most probable values by construction.

Confidence interval (frequentist): 80% CI: in 80% of repeated experiments, the interval constructed by a given procedure will contain the fixed true $\theta$. It does not allow the statement “$\theta$ is in this specific interval with probability 0.8”.

Aleatoric vs. Epistemic Uncertainty¶

Law of total variance:

Aleatoric Uncertainty: irreducible data noise (e.g. binomial or Gaussian observation noise).

Epistemic Uncertainty: uncertainty about model parameters; reducible with more/better data.

If Epistemic is bigger than Aleatoric, it means collecting more data can meaningfully reduce predictive uncertainty.

If Aleatoric is bigger than Epistemic, it means that even infinite data cannot substantially reduce predictive spread (inherent randomness).

Hypothesis Testing¶

Many branches of science depend on hypothesis testing to build knowledge. Hypothesis tests are often used for feature selection (e.g. t-test p-value for features in frequentist linear regression). In Bayesian statistics, one method for model selection is based on hypothesis testing (a model is a hypothesis).

Frequentist vs. Bayesian Hypothesis Tests¶

Frequentist hypothesis testing can only compute the probability of data under a given hypothesis. Set up a null hypothesis and discard it if the probability $P(d \mid H_0)$ is too low. This approach is dangerous when we measure unlikely data or are looking at unlikely hypotheses.

Bayesian hypothesis testing can compute $P(H \mid d)$, i.e. our posterior belief in $H$ given data $d$! Does not need a null hypothesis - can compare any two or more competing hypotheses! Bayesian and frequentist hypothesis tests return different results whenever there is heavy prior belief in $H_0$ and we have measured unlikely data that are in favor of $H_1$: The Bayesian approach is more robust in this case because it supports the philosophy of “one random accident should not overthrow an established theory”.

One-Sided Hypothesis Tests¶

Frequentist Hypothesis Test: Null Hypothesis Significance Testing (NHST)

1. Define $H_0$ and $H_1$

2. Define an estimator $\hat\pi$ and significance level $\alpha$ (false positive rate, usually $\alpha = 0.05$)

3. Define a test statistic (i.e. z-score in case of binominal distribution, type of test statistic depends on estimator!): $z = \frac{\hat\pi - \pi_0}{\sigma} = \frac{\hat\pi - \pi_0}{\sqrt{\frac{\hat\pi_0(1-\hat\pi_0)}{n}}}$

4. Determine p-value: $p = P(Z > z) = 1 - \Phi(z)$

5. If $p < \alpha$, reject $H_0$$

Bayesian Hypothesis Test: Compare posterior probability of hypotheses

1. Define $H_0$ and $H_1$:

   • $H_0:$ 80% of people or less in the local population are able to roll their tongue

   • $H_1:$ More than 80% of people in the local population are able to roll their tongue

2. Define a prior distribution for $\pi$ (e.g. Beta(1,1) for flat prior)

3. Define a likelihood function (e.g. Binomial likelihood)

4. Sample posterior $\pi_i$

5. Approximate posterior distribution of $\pi$ using $H_0$ and $H_1$

pH1d = np.mean(trace.posterior.pi.values > 0.8)
pH1d, 1.0 - pH1d

Bayes Factor and Odds¶

Prior odds for $H_1$ vs $H_0$ (how much more likely is $H_1$ than $H_0$): $\frac{P(H_1)}{P(H_0)}$, flip the ratio if you are interested in how much more likely is $H_0$ than $H_1$

Posterior odds after data: $\frac{P(H_1\mid d)}{P(H_0\mid d)}$

Bayes factor $BF = \frac{P(d\mid H_1)}{P(d\mid H_0)} = \frac{\text{posterior odds}}{\text{prior odds}}$, describes how much my believe changed or how much I learned from the data. Rules of Thumb: 1–3: barely worth mentioning, 3–10: substantial, 10–100: strong, bigger than 100: decisive evidence in favour of $H_1$.

Two-Sided Hypothesis Tests¶

Region of Practical Equivalence (ROPE)¶

In Bayesian inference, testing whether a parameter exactly equals a specific value is statistically problematic because the probability of a continuous parameter equaling any single point is essentially zero. Instead, ROPE defines a range of values considered “practically equivalent” to the null hypothesis. For testing whether a parameter is practically equal to some value $\pi_0$, define a ROPE, e.g. $[\pi_0-\delta, \pi_0+\delta]$.

• $H_0$: Exactly 93% of people in the local population are able to roll their tongue $\implies \pi_0 \in [\pi_0-\delta,\pi_0+\delta]$ for $\pi_0 = 0.93$ and some $\delta > 0$.

• $H_1$: More or less than 93% of people in the local population are able to roll their tongue $\implies \pi_0$ outside that interval.

pi0 = 0.93; delta = 0.03
pi_values = trace.posterior.pi.values
pH0d = np.mean((pi_values > pi0 - delta) &
    (pi_values < pi0 + delta))
pH1d = 1 - pH0d
pH0d, pH1d, (pH0d / pH1d)

If most posterior mass lies inside the ROPE, you accept practical equivalence; otherwise, you conclude a meaningful difference. If the posterior distribution is fairly diffuse, with substantial mass both inside and outside this range a narrow ROPE {small $\delta$) can create an inconclusive result. A larger ROPE (bigger $\delta$) reflects what differences you actually care about practically: If you only need to know whether the true proportion is “around 93%” rather than “exactly 93%”, a wider ROPE makes sense.

Posterior Predictive Checks¶

1. Fit model and obtain posterior samples.

2. For each posterior draw, simulate an entire replicated dataset from the model (posterior predictive).

3. Compare distribution of replicated data to observed data (histograms, kernel density plots, summary statistics).

If simulated data systematically differ (e.g. lower variance, wrong tail behaviour), the model is likely misspecified (e.g. wrong likelihood family).

Model Checking, Selection, and Multivariate Distributions¶

Model selection criteria¶

The Fully Bayesian way: Use Bayes Factor (based on marginal likelihood). This is conceptually clean, but marginal likelihood can be hard to compute.

Use Predictive Criteria instead: These criteria trade off fit vs complexity; Bayesian evidence automatically penalizes overly complex models (“built-in Occam’s razor”).

• RMSE, MAE (point prediction error using posterior predictive means)

• ELPD / LOO (Expected Log Predictive Density with leave-one-out cross-validation, Computed with PSIS-LOO, higher is better)

• Bayesian $R^2$ for regression

Bayesian Linear Regression¶

What is the posterior distribution of the linear regression parameters $\beta_0$, $\beta_1$, and $\sigma$? Use Bayes’ theorem (and condition on $\mathbf{x}$):

• Likelihood: $p(\mathbf{y} \mid \beta_0, \beta_1, \sigma, \mathbf{x})$ — probability of observing the data given parameters

• Prior: $p(\beta_0) p(\beta_1) p(\sigma)$ — treat parameters as independent, no direct dependence on $\mathbf{x}$

• Evidence: $p(\mathbf{y} \mid \mathbf{x})$ — marginal likelihood computed via MCMC

Then define priors and likelihood distributions:

Probabilistic Graphical Model Notation¶

Compact graphical way to describe a model without the distributions of the individual variables. Every unobserved variable that has no arrows going into it does need a prior. Bayesian inference will get us probability distributions for the unobserved variables.

Implementation and interpretation (PyMC/Bambi)¶

Robust Linear Regression¶

Motivation: Outliers can have a big impact on the resulting regression line.

Reason: The Gaussian likelihood does not like points far away from its center.

Solution: Can be mitigated using a Student’s t-distribution with a heavier tail. Bayesians often prefer to change the likelihood instead of the data because its more reproducible.

Gaussian Processes (GPs)¶

Intuition and definition¶

Gaussian Processes are prior distributions over functions: $f(x) \sim \mathcal{GP}(m(x), k(x,x'))$. A Gaussian Process is fully characterized by: Mean $m(x) = E[f(x)]$ and a Covariance function $k(x,x') = cov[f(x), f(x')]$

Covariance function $k(x,x')$ must be: Symmetric and Positive semi-definite: for any finite set of inputs, the matrix $K$ must produce non-negative quadratic forms $z^\top K z \ge 0$. If using a “kernel” that is not positive semi-definite (e.g. $K(x,x') = 2 + |x-x'|$), you can obtain negative predictive variances in the formula $v - k^\top C^{-1}k$. Such a function is not a valid covariance function, because a covariance matrix must be positive semi-definite, implying predictive variance must be $\ge 0$.

Key Intuition about Gaussian Processes: If things are close to each other, they should behave similar (i.e. be stronger correlated). That implies that the covariance between outputs is given in terms of the inputs.

Conditioning Rule: (“whats the prediction given the training data”) $p(f_* \mid f) = \mathcal{N}\left(f_* \mid K_{f_* f} K_{ff}^{-1} f, K_{f_* f_*} - K_{f_* f} K_{ff}^{-1} K_{f_* f}^T\right)$

Posterior Distribution for $f_*$ given $y$, including observation noise:

Inverting the covariance matrix of the noisy observations $C = (K_{ff} + \sigma^2_{obs}I) \in \mathbb{R}^{n \times n}$ is the most computationally expensive step, resulting in a complexity of $O(N^3)$. Therefore, large datasets are challenging for standard GPs.

GP regression: predictive mean and variance¶

Predictive distribution for response $y_*$:

• Training outputs $y \in \mathbb R^N$.

• Covariance matrix for training outputs (including noise): $C = K_{ff} + \sigma_{\text{obs}}^2 I = k(X, X) + \sigma_{\text{obs}}^2 I$

• Covariances between test point and training points: $k = k(x_*, X) \in \mathbb R^N$

• Prior variance at test point (for $y_*$): $v = k(x_*,x_*) + \sigma_{\text{obs}}^2$

Connection to kernels and PSD requirement¶

Common kernels and hyperparameters:

• Squared Exponential / RBF: $k(x,x') = \alpha \exp\left(-\frac{\|x-x'\|^2}{2\ell^2}\right)$, ($\alpha$: output variance, $\ell$: length scale (smoothness))

• Matérn (controls differentiability via $\nu$):

• Rational Quadratic (mixture over length scales).

Bayesian Networks and Causal Inference¶

Structural Causal Models (SCMs) and Bayesian Networks¶

SCM formalism (from BN lecture notes):

• Endogenous (ovservered) variables $V$: determined by structural equations.

• Exogenous (unobserved) variables $U$: external noise.

• Structural equations: each variable $V_i$ is a function of its parents (endogenous causes) and noise (exogenous causes):

• Joint distribution defined via $P(U)$ and structural equations.

A Bayesian Network (BN) is a DAG over variables $V$ such that joint distribution factorizes as:

• Graph encodes conditional independencies.

• Each node has a Conditional Probability Table (CPT) or parameterized distribution.

Basic graph motifs and conditional independence¶

Three canonical patterns:

1. Chain (serial): $X \to Z \to Y$

   • $X$ and $Y$ correlated marginally.

   • Conditional independence given the middle node: $X \perp Y \mid Z$.

2. Fork (common cause): $X \leftarrow Z \to Y$

   • $X$ and $Y$ correlated marginally (through common cause $Z$).

   • Conditional independence given cause: $X \perp Y \mid Z$.

3. Collider (v‑structure): $X \to Z \leftarrow Y$

   • $X$ and $Y$ independent marginally.

   • Conditioning on collider (or its descendants) induces dependence (explaining away).

   • Explaining away: In collider $X \to Z \leftarrow Y$, observing the effect $Z$ makes causes $X$ and $Y$ compete.

     • Example: “Smoker” and “Covid” both cause “Hospitalization”; given hospitalization, learning that patient is a smoker reduces belief that Covid is cause, and vice versa.

Computing Probabilities from Joint Tables¶

Given full joint distribution table $P(H,S,C)$:

Marginal: $P(C=1) = \sum_{s,h} P(C=1, H=h, S=s) = P(C=1, H=0, S=0) + P(C=1, H=0, S=1) + \dots$

Conditional: $P(C=1 \mid H=1) = \dfrac{\sum_s P(C=1,H=1,S=s)}{\sum_{c,s}P(C=c,H=1,S=s)}$

Example: $P(C = 1 \mid H = 1, S = 0) = \dfrac{P(C = 1, H = 1, S = 0)}{P(H = 1, S = 0)} = \dfrac{P(C = 1, H = 1, S = 0)}{\sum_c P(C = c, H = 1, S = 0)}$

Interventions and the do‑operator¶

Conditioning: $P(Y\mid X=x)$ describes correlations in observed data; $X$ may itself be caused by other variables (confounded).

Intervention: $P(Y\mid \text{do}(X=x))$ describes what happens if we force $X$ to $x$ (surgery on the graph):

• All incoming edges into $X$ are cut.

• $X$ is set externally.

• The rest of the graph and conditional distributions remain unchanged.

For a treatment $X$, outcome $Y$, and set of observed confounders $Z$ that block all backdoor paths (backdoor criterion), adjustment formula:

Expectation Maximization (EM)¶

The EM algorithm iteratively maximizes the log-likelihood of the observed data. $H$ is hidden, $E = e$ is observed.

1. Initialize $\vartheta$ randomly

2. Repeat until convergence:

E-step:

• Compute $q(h) = p(H = h \mid E = e; \vartheta)$ for each $h$ (probabilistic inference)

• Create fully-observed weighted examples: $(h, e)$ with weight $q(h)$

M-step:

• Perform maximum likelihood (count and normalize) on weighted examples to update $\vartheta$