Model uncertainty with random variables and distributions, then estimate, test and update beliefs from data, with Python simulations you run in the browser.

Probability & Statistical Inference

## Why this matters

A lab test is 90% sensitive and 95% specific, and a patient tests positive. What is the chance they have the condition? Most people, including many professionals, say about 90%. For a rare condition the right answer can be closer to 15%. **Bayes' rule** is the tool that gets this right. It drives diagnostic testing, spam filters, fraud alerts, quality inspection and the whole Bayesian approach to statistics that the rest of this module builds.

## The rule

For a hypothesis `H` and observed evidence `E`:

```text
P(H | E) = P(E | H) × P(H) / P(E)
P(E)     = P(E | H) P(H) + P(E | not H) P(not H)
```

- `P(H)` is the **prior**: how plausible `H` was before the evidence.
- `P(E | H)` is the **likelihood**: how expected the evidence is if `H` is true.
- `P(H | E)` is the **posterior**: the updated belief.

The denominator just makes sure the posterior probabilities of `H` and `not H` add to 1.

## Worked example: a screening test

- Prevalence (prior): 1% of people tested have the condition.
- **Sensitivity** `P(+ | condition) = 0.90`.
- **Specificity** `P(- | no condition) = 0.95`, so the false-positive rate is 0.05.

Picture 10,000 people:

```text
               condition (100)    no condition (9,900)
test +              90                 495
test -              10               9,405
```

Of the 585 positives, only 90 have the condition:

```text
P(condition | +) = 90 / 585 ≈ 0.154
```

With the formula: `0.90 × 0.01 / (0.90 × 0.01 + 0.05 × 0.99) = 0.009 / 0.0585 ≈ 0.154`.

The test is good, but the condition is rare, so false positives from the large healthy group outnumber true positives. Ignoring the prior like this is called **base-rate neglect**.

## Updating again

Today's posterior is tomorrow's prior. If the same patient takes a second, independent test and is positive again, start from 0.154:

```text
P(condition | +, +) = 0.90 × 0.154 / (0.90 × 0.154 + 0.05 × 0.846) ≈ 0.766
```

A negative result updates too, using `P(- | condition) = 0.10` and `P(- | no condition) = 0.95`. After one positive and then one negative, the probability drops to about 0.019.

This is why screening programmes confirm positives with a second test, and why the order of independent evidence does not matter: you reach the same posterior either way.

## The odds form

Bayes' rule is easier to do in your head with **odds**:

```text
posterior odds = prior odds × likelihood ratio
likelihood ratio for a positive (LR+) = sensitivity / (1 - specificity) = 0.90 / 0.05 = 18
```

Prior odds of 1 : 99, times 18, gives 18 : 99, which is a probability of 18 / 117 ≈ 0.154. Each independent piece of evidence multiplies the odds by its own likelihood ratio.

## In code

```python
def update(prior, p_e_given_h, p_e_given_not_h):
    numerator = p_e_given_h * prior
    return numerator / (numerator + p_e_given_not_h * (1 - prior))

p = update(0.01, 0.90, 0.05)     # positive test
print(round(p, 3))               # 0.154
p = update(p, 0.90, 0.05)        # second positive
print(round(p, 3))               # 0.766
```

## Beyond diagnostics

The same calculation appears everywhere:

- **Spam filters** multiply likelihood ratios for words like "invoice" or "winner".
- **Quality inspection**: a vision system flags a part; how likely is it really defective, given a 0.5% defect rate?
- **Fraud and security alerts**: a 99%-accurate rule on 1-in-10,000 events produces mostly false alarms.

## Bayesian vs frequentist, briefly

Earlier modules used **frequentist** methods: parameters are fixed, and probability describes the data you would get over repeated samples. The **Bayesian** view treats an unknown parameter as uncertain and describes that uncertainty with a probability distribution, updated with Bayes' rule. Both are mainstream; you will use each where it fits. The next two lessons apply the Bayesian view to rates.

## Common mistakes

- **Confusing P(+ | condition) with P(condition | +).** Sensitivity is not the chance a positive is real.
- **Ignoring the base rate.** The prior matters most when it is small.
- **Treating dependent evidence as independent.** Repeating the *same* test on the *same* sample does not give two independent likelihood ratios.
- **Using a negative result's likelihoods for a positive one.** For a negative, use `1 - sensitivity` and `specificity`.

## Recap

- Posterior ∝ likelihood × prior; normalise so the probabilities add to 1.
- Rare conditions mean many positives are false, even with a good test.
- Update sequentially: each posterior becomes the next prior.
- Odds form: posterior odds = prior odds × likelihood ratio.

## Try it

In the **Do** tab, write an update function and apply it to a sequence of test results.


Update the probability of a hypothesis after a test result, one piece of evidence at a time.

Bayes' rule in practice

Bayes theorem, the geometry of changing beliefs

You Know I'm All About that Bayes: Crash Course Statistics #24

## Why this matters

The previous lesson updated a single yes/no hypothesis. Most real questions are about an unknown **rate**: a conversion rate, a defect rate, a cure rate, a claim frequency. Bayesian inference gives you a full probability distribution for that rate, so you can answer the questions decision-makers actually ask: "What is the chance version B is better than A?" or "How likely is the defect rate above 2%?" The beta-binomial model is the cleanest example, and it is used directly in A/B testing tools, reliability work and actuarial credibility theory.

## The beta distribution

The **beta distribution** lives on the interval 0 to 1, which makes it a natural model for a rate `p`. It has two positive parameters, `a` and `b`:

```text
mean = a / (a + b)
var  = a b / [ (a + b)² (a + b + 1) ]
```

A helpful reading is "`a` prior successes and `b` prior failures":

- `Beta(1, 1)`: flat (uniform). Every rate from 0 to 1 is equally plausible.
- `Beta(2, 8)`: centred near 0.2, fairly vague (like 10 observations' worth of information).
- `Beta(20, 80)`: also centred at 0.2, much more confident (100 observations' worth).

## The conjugate update

With a `Beta(a, b)` prior and `s` successes in `n` trials (binomial likelihood), the posterior is again a beta:

```text
posterior = Beta(a + s, b + n - s)
```

That is the whole calculation: add successes to `a` and failures to `b`. A prior and likelihood that give a posterior of the same family are called **conjugate**.

The posterior mean is a weighted compromise between the prior mean and the observed rate:

```text
posterior mean = (a + s) / (a + b + n)
```

With little data the prior dominates; with lots of data the observed rate dominates. Actuaries call this the **credibility** weighting.

## Worked example: a sign-up rate

A new landing page converted 18 of 100 visitors.

**Flat prior** `Beta(1, 1)`: posterior `Beta(19, 83)`, mean `19 / 102 ≈ 0.186`. A 95% **credible interval** (the middle 95% of the posterior) is about 0.117 to 0.267.

**Informed prior** from past pages, `Beta(20, 80)` (mean 0.20): posterior `Beta(38, 162)`, mean `0.19`, and a narrower interval because the prior adds information.

A credible interval has the reading people often (wrongly) give to confidence intervals: "given the model and prior, there is a 95% probability the rate lies in this range".

### Summaries in code

```python
import numpy as np
rng = np.random.default_rng(0)
draws = rng.beta(19, 83, size=100_000)
print(draws.mean())                         # about 0.186
print(np.percentile(draws, [2.5, 97.5]))    # about [0.117, 0.267]
print((draws > 0.25).mean())                # P(rate > 25%), about 0.05
```

With SciPy you can get exact values instead: `scipy.stats.beta(19, 83).ppf([0.025, 0.975])`.

## Worked example: an A/B test

Version A converted 18 of 100 visitors; version B converted 30 of 100. With flat priors:

```text
A: Beta(19, 83)        B: Beta(31, 71)
```

What is `P(p_B > p_A)`? There is no simple formula, but sampling makes it easy: draw from both posteriors and count how often B's draw is larger.

```python
a = rng.beta(19, 83, size=100_000)
b = rng.beta(31, 71, size=100_000)
print((b > a).mean())          # about 0.976
```

"There is about a 98% probability B has the higher conversion rate" is a statement product managers can act on. You can also look at the distribution of `b - a` to see **how much** better B is likely to be, which matters more than whether it is better at all.

## Choosing a prior

- Use a **weakly informative** prior (such as `Beta(1, 1)` or `Beta(2, 2)`) when you know little.
- Use an **informed** prior from historical data when you have it, and say where it came from.
- Run a **sensitivity check**: refit with a different reasonable prior. If the conclusion flips, the data is not yet strong enough to decide.

## Common mistakes

- **Mixing up `a` and `b`.** `a` counts successes, `b` failures.
- **Adding trials instead of failures to `b`.** The update is `b + (n - s)`.
- **A prior that is too strong.** `Beta(2000, 8000)` will barely move with 100 new visitors. Make sure the prior's "worth" in observations is honest.
- **Reading `P(B > A) = 0.98` as "B is 98% better".** It is the probability B is better at all, not the size of the improvement.

## Recap

- Beta prior + binomial data -> beta posterior: `Beta(a + s, b + n - s)`.
- Posterior mean = `(a + s) / (a + b + n)`, a compromise between prior and data.
- Credible intervals and probabilities like `P(B > A)` come straight from posterior draws.
- Check how sensitive conclusions are to the prior.

## Try it

In the **Do** tab, write the beta-binomial update and estimate `P(B > A)` by sampling with NumPy.


Use a beta prior with binomial data to get a beta posterior, summarise it, and compare two rates in an A/B test.

Priors and posteriors

Bayesian Updates and Conjugate Priors

Binomial distributions | Probabilities of probabilities, part 1

## Why this matters

The beta-binomial model has a neat closed-form posterior. Most real models do not: a logistic regression with a dozen predictors, a hierarchical model of many hospitals or machines, a pharmacokinetic curve. Their posteriors have no formula, and the normalising constant `P(data)` is an intractable integral. **Markov chain Monte Carlo (MCMC)** gets around this by drawing samples from the posterior without ever computing that constant. Tools such as PyMC and Stan, used across data science, biostatistics and actuarial modelling, are built on it.

## The key insight

Bayes' rule says

```text
posterior(θ) = likelihood(θ) × prior(θ) / P(data)
```

`P(data)` does not depend on `θ`. So the **ratio** of the posterior at two points needs only the numerator:

```text
posterior(θ_new) / posterior(θ_old) = [likelihood × prior](θ_new) / [likelihood × prior](θ_old)
```

If we can compare points, we can build a random walk that spends time in each region in proportion to its posterior probability. The positions it visits are then (dependent) samples from the posterior, and everything from module 2 applies: means, percentiles and probabilities come from counting samples.

## The Metropolis algorithm

1. Start at some value `θ`.
2. **Propose** a nearby value: `θ' = θ + noise`, with noise drawn from `Normal(0, step)`.
3. Compute the acceptance ratio `r = posterior(θ') / posterior(θ)`.
4. **Accept** the move with probability `min(1, r)`: always move uphill; sometimes move downhill.
5. Record the current position (the new one if accepted, otherwise the old one again).
6. Repeat thousands of times.

Work with **log** posteriors to avoid underflow: accept when `random < exp(min(0, log_post(θ') - log_post(θ)))`. A proposal outside the allowed range (such as `p < 0` for a rate) has log posterior `-∞`, so `exp(-∞) = 0` and it is always rejected.

## Worked example: a rate with a known answer

A coin (or a pass/fail process) gives 7 successes in 10 trials, with a `Beta(2, 2)` prior. From the previous lesson we know the exact posterior is `Beta(9, 5)` with mean `9/14 ≈ 0.643` and SD about `0.124`. That lets us check MCMC against the truth.

```python
import math, random

def log_post(p):
    if not 0 < p < 1:
        return -math.inf
    # log prior Beta(2, 2) + log likelihood (7 of 10), constants dropped
    return (1 + 7) * math.log(p) + (1 + 3) * math.log(1 - p)

rng = random.Random(42)
current, current_lp = 0.5, log_post(0.5)
samples, accepted = [], 0
for _ in range(30_000):
    proposal = current + rng.gauss(0, 0.1)
    proposal_lp = log_post(proposal)
    if rng.random() < math.exp(min(0, proposal_lp - current_lp)):
        current, current_lp = proposal, proposal_lp
        accepted += 1
    samples.append(current)

kept = samples[2_000:]                 # discard burn-in
print(sum(kept) / len(kept))           # about 0.64
print(accepted / 30_000)               # acceptance rate, about 0.77
```

The sample mean lands close to 0.643. Replace `log_post` with any model's log-likelihood plus log-prior and the same loop works.

## Tuning and diagnostics

- **Burn-in.** Early samples depend on the starting point. Discard the first part of the chain (here 2,000 steps).
- **Step size.** Too small and nearly every move is accepted but the chain crawls; too large and almost everything is rejected. For one parameter, an acceptance rate of very roughly 40 to 80% is workable; tools tune this automatically.
- **Autocorrelation.** Consecutive samples are dependent, so 28,000 kept samples carry less information than 28,000 independent ones. Libraries report an **effective sample size (ESS)**.
- **Convergence.** Run several chains from different starts. If they disagree (the `R-hat` diagnostic is well above 1.01), the chain has not converged and the results are not trustworthy.
- **Trace plots.** Plot samples against step number: a healthy chain looks like a fuzzy horizontal band, not a slow drift.

## In practice

You will rarely hand-write Metropolis for real work. PyMC (Python) and Stan use **Hamiltonian Monte Carlo (NUTS)**, a smarter proposal that handles many parameters efficiently:

```python
import pymc as pm
with pm.Model():
    p = pm.Beta("p", 2, 2)
    pm.Binomial("y", n=10, p=p, observed=7)
    trace = pm.sample(2000, chains=4, random_seed=1)
```

Writing the simple version once means you can read their diagnostics and know what can go wrong.

## Common mistakes

- **Forgetting to repeat the current value on rejection.** Recording only accepted moves biases the samples.
- **Using raw probabilities instead of logs**, which underflow for real datasets.
- **No burn-in and one short chain.** Always check convergence.
- **Treating autocorrelated samples as independent** when reporting precision.

## Recap

- MCMC samples a posterior using only `likelihood × prior`, never `P(data)`.
- Metropolis: propose nearby, accept with probability `min(1, ratio)`, record the current value every step.
- Discard burn-in, tune the step size, and check convergence with several chains.
- Production tools (PyMC, Stan) use the same idea with smarter proposals.

## Try it

In the **Do** tab, write the Metropolis acceptance step and check your samples against the exact Beta(9, 5) posterior.


Approximate a posterior you cannot compute in closed form by sampling it with the Metropolis algorithm.

Simulation-based Bayes

Metropolis - Hastings : Data Science Concepts

Introduction to Bayesian statistics, part 2: MCMC and the Metropolis–Hastings algorithm

Updating beliefs with Bayes' rule, conjugate beta-binomial posteriors for rates and A/B tests, and approximating posteriors by sampling with the Metropolis algorithm.