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

Module 1 started from known parameters ("defects occur at 2%", "calls arrive at 3 per hour") and computed probabilities. Real work runs the other way: you have data and need the parameter. What defect rate best explains today's inspections? What failure rate fits these component lifetimes? **Maximum likelihood estimation (MLE)** is the standard answer. It is how logistic regression, most machine-learning loss functions, survival models and actuarial frequency models are fitted.

## Estimators and their qualities

An **estimator** is a rule that turns data into a guess for a parameter, such as the sample mean `x̄` for the population mean `μ`. Two properties matter most:

- **Bias**: on average over many samples, does the estimator hit the true value? The sample variance with `n - 1` is unbiased; with `n` it is slightly biased low.
- **Variance**: how much does it bounce around from sample to sample? Its square root is the standard error.

A good estimator has both small. Their combined effect is the **mean squared error**: `MSE = bias² + variance`.

## Likelihood

The **likelihood** `L(θ)` is the probability (or density) of the data you actually observed, viewed as a function of the parameter `θ`. For independent observations it is a product:

```text
L(θ) = P(x1 | θ) × P(x2 | θ) × ... × P(xn | θ)
```

The **maximum-likelihood estimate** `θ̂` is the value of `θ` that makes the observed data most probable.

Products of many small numbers underflow on a computer, so in practice we maximise the **log-likelihood**, a sum:

```text
ℓ(θ) = log L(θ) = sum of log P(xi | θ)
```

The log is increasing, so the same `θ` maximises both.

## Worked example: a coin or a pass rate

You observe 7 successes in 10 trials. Under a binomial model,

```text
ℓ(p) = 7 log p + 3 log(1 - p) + constant
```

Try a few values: `ℓ(0.5) ≈ -6.93`, `ℓ(0.7) ≈ -6.11`, `ℓ(0.9) ≈ -7.64` (dropping the constant). The peak is at `p̂ = 0.7`. Setting the derivative to zero confirms it: `7/p - 3/(1-p) = 0` gives `p = 7/10`. The MLE is the sample proportion, which is reassuring.

## Worked example: Poisson defects per day

Daily defect counts over 8 days: `2, 3, 1, 0, 4, 2, 3, 1`. For a Poisson model with rate `λ`,

```text
log P(k | λ) = k log λ - λ - log(k!)
ℓ(λ) = (sum of k) log λ - n λ - sum of log(k!)
```

Setting the derivative to zero gives `λ̂ = (sum of k) / n = x̄ = 16 / 8 = 2.0`. Again the MLE is the sample mean.

You do not always get a closed form, so it is worth knowing how to find the peak numerically. A **grid search** is the simplest method:

```python
import math

data = [2, 3, 1, 0, 4, 2, 3, 1]

def loglik(lam):
    return sum(k * math.log(lam) - lam - math.lgamma(k + 1) for k in data)

grid = [i / 100 for i in range(1, 801)]     # 0.01 to 8.00
best = max(grid, key=loglik)
print(best)                                  # 2.0
```

`math.lgamma(k + 1)` is `log(k!)` computed without overflow.

## Common closed-form MLEs

| model | parameter | MLE |
|---|---|---|
| Bernoulli / binomial | `p` | successes / trials |
| Poisson | `λ` | sample mean |
| Exponential | rate `λ` | `1 / x̄` |
| Normal | `μ` | `x̄` |
| Normal | `σ²` | average squared deviation, dividing by `n` (slightly biased) |

The last row is a useful reminder that MLEs are not always unbiased. They are, however, **consistent** (they converge to the truth as `n` grows) and, for large samples, approximately normal with the smallest possible variance. That is why they are the default.

## In libraries

Real models have several parameters, and you maximise with an optimiser instead of a grid. A typical pattern is to minimise the **negative** log-likelihood:

```python
from scipy.optimize import minimize_scalar
res = minimize_scalar(lambda lam: -loglik(lam), bounds=(0.01, 10), method="bounded")
print(res.x)        # about 2.0
```

`scipy.stats` distributions also have a `.fit(data)` method that returns MLEs directly, and machine-learning "log loss" or "cross-entropy" is exactly a negative log-likelihood.

## Common mistakes

- **Multiplying probabilities instead of adding logs**, which underflows to 0 for large data.
- **Searching outside the valid range.** `log(0)` fails; keep `p` in (0, 1) and `λ > 0`.
- **A grid that is too coarse or too narrow.** If the best value sits at the edge of the grid, widen it.
- **Treating the MLE as the truth.** It is an estimate with a standard error; report one.

## Recap

- Likelihood: probability of the observed data as a function of the parameter.
- Maximise the log-likelihood (a sum) instead of the likelihood (a product).
- Many MLEs are familiar: the sample proportion, the sample mean, `1 / x̄`.
- When no formula exists, search numerically (grid search, then an optimiser).

## Try it

In the **Do** tab, write a Poisson log-likelihood and find its maximum by grid search.


Write a log-likelihood and find the maximum-likelihood estimate of a parameter, numerically and in closed form.

Estimation and MLE

Maximum Likelihood, clearly explained!!!

Maximum Likelihood Estimation (MLE) with Examples

## Why this matters

Statistics Fundamentals introduced the logic of a hypothesis test: assume no effect (the null hypothesis), measure how surprising the data is, and report a p-value. In practice the hard part is choosing **which** test. Did a new catalyst change yield? Did a redesign change the sign-up rate? Is defect type related to shift? Each question has a standard test, and picking the wrong one is one of the most common errors reviewers flag.

## Step 1: what kind of data?

- **A numeric outcome** (yield, time, concentration) compared across one or two groups: a **t-test**.
- **Counts in categories** (pass/fail by shift, defect type by machine): a **chi-square test**.

## The t-test family

All t-tests share one shape:

```text
t = (observed difference - difference under H0) / standard error
```

and compare `t` to a t distribution with some degrees of freedom (df). The versions differ in what the difference and standard error are.

| situation | test | in SciPy |
|---|---|---|
| One group vs a target value (is mean fill 500 ml?) | one-sample t | `stats.ttest_1samp(x, 500)` |
| Same units measured twice (before/after on each patient or machine) | paired t | `stats.ttest_rel(before, after)` |
| Two independent groups (catalyst A vs catalyst B batches) | Welch two-sample t | `stats.ttest_ind(a, b, equal_var=False)` |

A **paired** test is just a one-sample t-test on the differences. Pairing removes unit-to-unit variation, so when the design is paired, use the paired test: treating paired data as independent throws away power.

For two independent groups, prefer **Welch's** test, which does not assume equal variances. The older "Student" version (`equal_var=True`) is only safe when the variances really are similar, and Welch loses almost nothing when they are.

### Welch's formulas

```text
SE = sqrt(s1²/n1 + s2²/n2)
t  = (x̄1 - x̄2) / SE
df = (s1²/n1 + s2²/n2)² / [ (s1²/n1)²/(n1-1) + (s2²/n2)²/(n2-1) ]
```

### Worked example: two catalysts

Yields (%) from 8 batches with catalyst A and 7 with catalyst B:

```text
A: 12.1 11.4 13.0 12.7 11.9 12.5 13.3 12.2    mean 12.39
B: 11.2 11.8 10.9 11.5 12.0 11.1 11.6         mean 11.44
```

Welch's test gives `t ≈ 3.58` with `df ≈ 12.0`, and a two-sided p-value of about 0.004. The data is hard to explain if the catalysts were equivalent. Always report the effect alongside the p-value: A's mean yield is about 0.95 percentage points higher.

Assumptions: independent observations, and either roughly normal data in each group or samples large enough for the CLT (module 2). With small, heavily skewed samples, consider a rank test (`stats.mannwhitneyu`) or a bootstrap.

## The chi-square test of independence

For two categorical variables, arrange counts in a **contingency table**. Under the null hypothesis of independence, the expected count in each cell is

```text
expected = (row total × column total) / grand total
```

and the statistic adds up how far the observed counts are from expected:

```text
χ² = sum over cells of (observed - expected)² / expected
df = (rows - 1) × (columns - 1)
```

### Worked example: pass rate by shift

|  | pass | fail | total |
|---|---|---|---|
| Day shift | 30 | 20 | 50 |
| Night shift | 15 | 35 | 50 |
| total | 45 | 55 | 100 |

Expected "day, pass" = 50 × 45 / 100 = 22.5, and similarly 27.5, 22.5, 27.5 for the other cells. Then

```text
χ² = (30-22.5)²/22.5 + (20-27.5)²/27.5 + (15-22.5)²/22.5 + (35-27.5)²/27.5 ≈ 9.09
df = (2-1)(2-1) = 1        p ≈ 0.0026
```

Pass rate depends on shift (60% vs 30%). In SciPy: `stats.chi2_contingency(table)`. Note that for 2×2 tables SciPy applies a continuity correction by default, so its statistic (about 7.92 here) is slightly smaller; pass `correction=False` to match the hand calculation.

Assumption: expected counts should be at least about 5 in each cell. For small tables with small counts, use Fisher's exact test (`stats.fisher_exact`).

A related test, **chi-square goodness of fit** (`stats.chisquare`), compares one categorical variable to expected proportions, such as checking that a die is fair.

## A quick decision guide

1. Categories and counts? Chi-square (or Fisher for small counts).
2. Numbers, one group vs a target? One-sample t.
3. Numbers, same units twice? Paired t.
4. Numbers, two separate groups? Welch t.
5. More than two groups? ANOVA, covered in Regression & Experimental Design.

## Common mistakes

- **Unpaired test on paired data.** You lose the benefit of pairing.
- **Running a t-test on proportions or counts in categories.** Use chi-square.
- **Using percentages in a chi-square table.** It needs raw counts.
- **Reporting only "p < 0.05".** Report the effect size, a confidence interval and the exact p-value.

## Recap

- t-tests compare means: one-sample, paired, or Welch for two independent groups.
- Chi-square tests compare observed and expected counts; `df = (r-1)(c-1)`.
- Check assumptions: independence, normality or large n for t; expected counts ≥ 5 for chi-square.

## Try it

In the **Do** tab, compute Welch's t and the chi-square statistic yourself and check them against the worked examples.


Choose between one-sample, paired, Welch and chi-square tests, and compute their test statistics.

t-tests and chi-square

StatQuickie: Which t test to use

Chi-Square Tests: Crash Course Statistics #29

## Why this matters

Two errors ruin studies. The first is running a study too small to detect the effect you care about, then concluding "no effect". The second is running so many tests that some come out "significant" by luck. A/B testing teams, lab scientists screening thousands of genes, and analysts slicing a dashboard by every segment meet both problems every week. This lesson gives you the tools to plan against the first and correct for the second.

## Errors and power

| | H0 true (no effect) | H0 false (real effect) |
|---|---|---|
| Reject H0 | Type I error (false positive), rate `α` | Correct: probability = **power** |
| Do not reject | Correct | Type II error (false negative), rate `β` |

**Power = 1 - β** is the probability that your study detects an effect of a given size if it is really there. A common target is 80%.

Power goes up when:

- the true effect `δ` is larger,
- the noise `σ` is smaller (better instruments, paired designs),
- the sample size `n` is larger,
- `α` is larger (a looser threshold, which costs more false positives).

## Computing power for two groups

For comparing two means with `n` per group, the normal approximation is

```text
SE    = σ × sqrt(2 / n)
power ≈ Φ( |δ| / SE - z(1 - α/2) )
```

where `Φ` is the standard normal CDF and `z(1 - α/2)` is 1.96 for `α = 0.05`. (The tiny chance of rejecting in the wrong direction is ignored.)

**Worked example.** A new process might raise mean strength by `δ = 5` MPa; the SD between specimens is `σ = 10` MPa.

```text
n = 20 per group:  SE = 10 × sqrt(2/20) = 3.16   power ≈ Φ(5/3.16 - 1.96) = Φ(-0.38) ≈ 0.35
n = 63 per group:  SE = 1.78                       power ≈ Φ(2.81 - 1.96) = Φ(0.85) ≈ 0.80
```

With 20 specimens per group you would miss a real 5 MPa improvement about two times in three. A non-significant result from that study says almost nothing.

### Solving for n

Rearranging gives the sample size per group for a target power:

```text
n = 2 × (z(1 - α/2) + z(power))² × σ² / δ²
  = 2 × (1.96 + 0.84)² × 100 / 25 ≈ 62.8  -> 63 per group
```

```python
from statistics import NormalDist
Z = NormalDist()
n = 2 * (Z.inv_cdf(0.975) + Z.inv_cdf(0.80)) ** 2 * 10**2 / 5**2
print(n)     # 62.79
```

Libraries such as `statsmodels.stats.power.TTestIndPower` do this with the exact t distribution and give a slightly larger `n`. Always round **up**, and plan for dropouts or failed runs.

The hard part is choosing `δ`: use the **smallest effect that would matter** to the business or science, not the effect you hope for.

## Multiple testing

Run one test at `α = 0.05` with no real effect and you have a 5% chance of a false positive. Run 20 independent tests and

```text
P(at least one false positive) = 1 - 0.95^20 ≈ 0.64
```

Testing 20 dashboard segments, 20 outcome measures or 20,000 genes guarantees "discoveries" unless you correct.

### Bonferroni: control the family-wise error rate

Test each of `m` hypotheses at `α / m`. With 10 tests at overall `α = 0.05`, each needs `p ≤ 0.005`. This keeps the chance of **any** false positive at or below 5%. It is simple and safe, but very strict when `m` is large.

### Benjamini-Hochberg: control the false discovery rate

When you screen many hypotheses and can tolerate a few false leads, control the **false discovery rate (FDR)**: the expected share of your discoveries that are false.

1. Sort the `m` p-values from smallest to largest: `p(1) ≤ p(2) ≤ ... ≤ p(m)`.
2. Find the **largest** rank `k` with `p(k) ≤ (k / m) × α`.
3. Reject the hypotheses with the `k` smallest p-values (all of them, even if an earlier one missed its own threshold).

**Worked example**, `m = 10`, `α = 0.05`:

```text
sorted p:   0.001  0.008  0.012  0.018  0.040  0.20 ...
threshold:  0.005  0.010  0.015  0.020  0.025  0.03 ...
passes?     yes    yes    yes    yes    no     no
```

The largest passing rank is 4, so BH rejects 4 hypotheses. Bonferroni (threshold 0.005 for all) rejects only 1. In Python: `statsmodels.stats.multitest.multipletests(p, method="fdr_bh")`, or SciPy's `scipy.stats.false_discovery_control`.

## Good practice

- **Pre-register** the primary outcome and the analysis before collecting data.
- **Don't peek and stop.** Checking an A/B test daily and stopping at the first p < 0.05 inflates false positives; use a fixed sample size or a sequential method designed for it.
- **Report every test you ran**, not just the significant ones.

## Common mistakes

- **"Not significant" read as "no effect"** from an underpowered study.
- **Power calculated after the fact** from the observed effect: it adds nothing beyond the p-value.
- **Correcting only the tests that looked interesting.** The correction must cover every test in the family.

## Recap

- Power = chance of detecting a real effect; plan for about 80% at the smallest effect that matters.
- For two means: `power ≈ Φ(|δ| / (σ sqrt(2/n)) - 1.96)`; `n ≈ 2(1.96 + 0.84)² σ²/δ²`.
- Many tests inflate false positives: Bonferroni controls any error, BH controls the false discovery rate.

## Try it

In the **Do** tab, write a power function and the Benjamini-Hochberg procedure.


Compute the power and sample size of a two-group study and control false positives across many tests with Bonferroni and Benjamini-Hochberg.

Power and multiple testing

Statistical Power, Clearly Explained!!!

False Discovery Rates, FDR, clearly explained

Estimating parameters by maximum likelihood, choosing and running t-tests and chi-square tests, and planning studies with power while controlling errors across many tests.