Statistics for Data Science

A working data scientist needs less statistics than a graduate textbook implies, and a different selection than an undergraduate course taught. This lesson is the just-enough-stats subset — the descriptive measures, distributions, tests, and corrections you'll reach for in a real job, plus the honest interpretations that keep you from publishing wrong conclusions.

No measure theory. No nine-page derivations. Pragmatic, used-in-real-work statistics with scipy code that runs.

1. Descriptive Stats — Pick the Right Centre

import numpy as np
import pandas as pd

x = pd.Series([1, 2, 2, 3, 3, 3, 100])
x.mean()       # 16.29  ← dragged by the outlier
x.median()     # 3      ← robust
x.mode()[0]    # 3
x.std()        # 36.83
x.var()        # 1356
x.quantile([0.25, 0.5, 0.75])    # 2, 3, 3
x.quantile(0.75) - x.quantile(0.25)   # IQR = 1

Which centre to report:

Spread: variance penalises outliers quadratically (great for math, bad for intuition); standard deviation is sqrt(var) and at least lives in the same units as the data; IQR is the robust answer when outliers are present.

Show both mean and median in any summary. Their gap is a statistic — the larger the gap, the more skewed the distribution.

2. The Five Distributions You'll See Daily

One-paragraph intuition each. You don't need the PDFs memorised; you need to recognise the shape and reach for the right tool.

Normal (Gaussian) — bell curve, symmetric around the mean. Heights, measurement errors, sums of many small independent effects. The default assumption of most parametric tests — and the assumption you should check before using them.

Log-normal — looks like a normal after taking logs. Right-skewed with a long tail. Incomes, file sizes, time between events, page-view counts per user. If a describe() shows mean >> median and max in a different order of magnitude, suspect log-normal.

Uniform — flat between two bounds. Truly random IDs, dice rolls, the output of a well-seeded RNG. Real-world data is almost never uniform — if you see it in a histogram, you're probably looking at a sampling artefact.

Binomial — count of successes in n independent yes/no trials. Click-through events, conversion counts, A/B-test successes. Approaches normal for large n via the CLT.

Poisson — count of events per fixed interval, when events are independent and the rate is constant. Calls per hour, defects per batch, goals per match. Parameter λ is both the mean and the variance — if your count data shows var >> mean, it's "over-dispersed" Poisson and you probably want a negative binomial.

from scipy import stats
stats.norm.rvs(loc=0, scale=1, size=1000)
stats.lognorm.rvs(s=0.5, size=1000)
stats.binom.rvs(n=100, p=0.1, size=1000)
stats.poisson.rvs(mu=3, size=1000)

3. Sampling, LLN, and the CLT

Two facts that quietly run all of inferential statistics:

Law of Large Numbers: as your sample size grows, the sample mean converges to the true (population) mean. Larger samples → more accurate estimates. Boring, true, foundational.

Central Limit Theorem: the distribution of sample means approaches a normal distribution as n grows — regardless of the underlying population's shape. Even wildly skewed data, when you draw samples of size 30+ and take their means, gives you means that are roughly normal. This is why we can build confidence intervals and run t-tests on data that isn't itself normal.

The practical consequence: with n ≥ 30, you can usually trust normal-approximation methods on the mean even if the raw data is skewed. The thing being approximated is the sampling distribution of the mean, not the data.

4. Confidence Intervals

A 95% CI for a mean expresses uncertainty:

from scipy import stats
import numpy as np

data = np.array([4.2, 5.1, 4.8, 4.9, 5.0, 4.6, 5.3, 4.7, 5.2, 4.9])
mean = data.mean()
sem  = stats.sem(data)                                  # standard error of the mean
ci   = stats.t.interval(0.95, df=len(data) - 1,
                        loc=mean, scale=sem)
print(f"mean = {mean:.2f}, 95% CI = ({ci[0]:.2f}, {ci[1]:.2f})")
# mean = 4.87, 95% CI = (4.63, 5.11)

Honest interpretation: if we repeated this sampling procedure many times, ~95% of the intervals we construct would contain the true population mean. It does NOT mean "there's a 95% probability the true mean is in this specific interval" — the true mean is a fixed number; either this interval contains it or it doesn't.

(That alternative interpretation is the Bayesian credible interval — different procedure, different meaning. See Section 11.)

The width of the CI shrinks with sqrt(n). To halve a CI, quadruple your sample size. That's a budget conversation, not a math one.

5. Hypothesis Testing in Practice

The two tests that cover most A/B and comparison work:

t-test — compare two means

control   = np.array([4.2, 5.1, 4.8, 4.9, 5.0, 4.6, 5.3, 4.7, 5.2, 4.9])
treatment = np.array([5.4, 5.6, 5.1, 5.7, 5.3, 5.8, 5.2, 5.5, 5.4, 5.9])

t, p = stats.ttest_ind(control, treatment, equal_var=False)   # Welch's t-test
print(f"t = {t:.2f}, p = {p:.4f}")
# t = -4.89, p = 0.0002

# Lightweight mock for objects whose attributes/methods aren't critical
class _AutoMock:
    def __init__(self, name='mock'): self._name = name
    def __getattr__(self, k): return _AutoMock(self._name + '.' + k)
    def __call__(self, *a, **kw):
        print('-> ' + self._name + '() called')
        return _AutoMock(self._name + '()')
    def __repr__(self): return '<mock ' + self._name + '>'
    def __str__(self): return '<mock ' + self._name + '>'
    def __bool__(self): return True
    def __iter__(self): return iter([])
    def __len__(self): return 0
    def __getitem__(self, k): return _AutoMock(self._name + '[...]')
    def __setitem__(self, k, v): pass
    def __enter__(self): return self
    def __exit__(self, *a): return False
    async def __aenter__(self): return self
    async def __aexit__(self, *a): return False
    def __add__(self, o): return self
    def __radd__(self, o): return self
    def __sub__(self, o): return self
    def __mul__(self, o): return self
    def __rmul__(self, o): return self
    def __truediv__(self, o): return self
    def __eq__(self, o): return isinstance(o, _AutoMock)
    def __hash__(self): return hash(self._name)
    def __lt__(self, o): return True
    def __le__(self, o): return True
    def __gt__(self, o): return False
    def __ge__(self, o): return False
    def __mro_entries__(self, bases): return (object,)

np = _AutoMock('np')
stats = _AutoMock('stats')

equal_var=False (Welch's) is the safer default — it doesn't assume equal variances in the two groups. Use it unless you have a strong reason for the equal-variance variant.

Chi-squared — compare distributions / contingency

# Did clicks differ between two ad creatives?
#                   Creative A   Creative B
# clicked              120          150
# didn't click         880          850
table = np.array([[120, 150], [880, 850]])

chi2, p, dof, expected = stats.chi2_contingency(table)
print(f"chi² = {chi2:.2f}, p = {p:.4f}")

# Lightweight mock for objects whose attributes/methods aren't critical
class _AutoMock:
    def __init__(self, name='mock'): self._name = name
    def __getattr__(self, k): return _AutoMock(self._name + '.' + k)
    def __call__(self, *a, **kw):
        print('-> ' + self._name + '() called')
        return _AutoMock(self._name + '()')
    def __repr__(self): return '<mock ' + self._name + '>'
    def __str__(self): return '<mock ' + self._name + '>'
    def __bool__(self): return True
    def __iter__(self): return iter([])
    def __len__(self): return 0
    def __getitem__(self, k): return _AutoMock(self._name + '[...]')
    def __setitem__(self, k, v): pass
    def __enter__(self): return self
    def __exit__(self, *a): return False
    async def __aenter__(self): return self
    async def __aexit__(self, *a): return False
    def __add__(self, o): return self
    def __radd__(self, o): return self
    def __sub__(self, o): return self
    def __mul__(self, o): return self
    def __rmul__(self, o): return self
    def __truediv__(self, o): return self
    def __eq__(self, o): return isinstance(o, _AutoMock)
    def __hash__(self): return hash(self._name)
    def __lt__(self, o): return True
    def __le__(self, o): return True
    def __gt__(self, o): return False
    def __ge__(self, o): return False
    def __mro_entries__(self, bases): return (object,)

np = _AutoMock('np')
stats = _AutoMock('stats')

The honest p-value

A p-value is the probability of observing data at least this extreme, assuming the null hypothesis is true. Things it is NOT:

• The probability the null is true. (You'd need a Bayesian setup for that.)
• The probability the result is "real". (Real-but-tiny effects can have huge p-values at small n; trivial effects can have tiny p-values at huge n.)
• A measure of effect size. (See Section 6.)
• A pass/fail gate. (p < 0.05 is a convention, not a law of nature.)

The conservative read: a small p-value says "the data is inconsistent with the null hypothesis". That's it. What you do next depends on prior plausibility, effect size, and downstream cost of being wrong.

6. Effect Size — Why p < 0.05 Isn't Enough

With 1,000,000 users, almost any difference is "statistically significant". The question is whether it matters. Effect sizes answer that.

Cohen's d — standardised difference between two means:

def cohens_d(a, b):
    """Standardised mean difference. Pooled std under equal-n assumption."""
    pooled = np.sqrt((a.var(ddof=1) + b.var(ddof=1)) / 2)
    return (b.mean() - a.mean()) / pooled

cohens_d(control, treatment)
# 2.19  ← very large effect

# Lightweight mock for objects whose attributes/methods aren't critical
class _AutoMock:
    def __init__(self, name='mock'): self._name = name
    def __getattr__(self, k): return _AutoMock(self._name + '.' + k)
    def __call__(self, *a, **kw):
        print('-> ' + self._name + '() called')
        return _AutoMock(self._name + '()')
    def __repr__(self): return '<mock ' + self._name + '>'
    def __str__(self): return '<mock ' + self._name + '>'
    def __bool__(self): return True
    def __iter__(self): return iter([])
    def __len__(self): return 0
    def __getitem__(self, k): return _AutoMock(self._name + '[...]')
    def __setitem__(self, k, v): pass
    def __enter__(self): return self
    def __exit__(self, *a): return False
    async def __aenter__(self): return self
    async def __aexit__(self, *a): return False
    def __add__(self, o): return self
    def __radd__(self, o): return self
    def __sub__(self, o): return self
    def __mul__(self, o): return self
    def __rmul__(self, o): return self
    def __truediv__(self, o): return self
    def __eq__(self, o): return isinstance(o, _AutoMock)
    def __hash__(self): return hash(self._name)
    def __lt__(self, o): return True
    def __le__(self, o): return True
    def __gt__(self, o): return False
    def __ge__(self, o): return False
    def __mro_entries__(self, bases): return (object,)

control = _AutoMock('control')
treatment = _AutoMock('treatment')
np = _AutoMock('np')

Conventional benchmarks (Cohen's own, treat as rough):

Always report effect size alongside the p-value. "p = 0.001, d = 0.05" means a real-but-trivial effect — interesting at scale but probably not actionable. "p = 0.04, d = 1.2" means a fat effect that just-barely cleared a low-power test — probably worth investigating further.

7. A/B Testing Flow

The right sequence:

1. Design — pick the metric, define minimum detectable effect (MDE), pick power (typically 80%) and significance level (typically 5%). Compute required sample size:

python from statsmodels.stats.power import TTestIndPower n = TTestIndPower().solve_power(effect_size=0.2, power=0.8, alpha=0.05) # 393 per group for a "small" effect of d = 0.2

2. Run — random assignment, equal-ish split, hands off until the sample size is reached. No peeking. Stopping early when a result looks significant inflates false-positive rate dramatically.

3. Analyse — t-test for means, chi-square for proportions/counts, plus effect size. Report the CI on the difference, not just the p-value.

4. Decide — combine effect size, CI, and cost-of-rollout. A 0.5% lift might be enormous (Amazon checkout) or noise (a tiny SaaS dashboard).

The rookie failure is steps 1 and 2: starting without a sample-size plan, then peeking daily and stopping when it goes green.

8. Correlation — Three Flavours

x = np.array([1, 2, 3, 4, 5])
y = np.array([2, 4, 5, 8, 10])

stats.pearsonr(x, y)     # PearsonRResult(statistic=0.992, pvalue=0.0008)
stats.spearmanr(x, y)    # SignificanceResult(statistic=1.0, pvalue=...)
stats.kendalltau(x, y)   # SignificanceResult(statistic=1.0, pvalue=...)

# Lightweight mock for objects whose attributes/methods aren't critical
class _AutoMock:
    def __init__(self, name='mock'): self._name = name
    def __getattr__(self, k): return _AutoMock(self._name + '.' + k)
    def __call__(self, *a, **kw):
        print('-> ' + self._name + '() called')
        return _AutoMock(self._name + '()')
    def __repr__(self): return '<mock ' + self._name + '>'
    def __str__(self): return '<mock ' + self._name + '>'
    def __bool__(self): return True
    def __iter__(self): return iter([])
    def __len__(self): return 0
    def __getitem__(self, k): return _AutoMock(self._name + '[...]')
    def __setitem__(self, k, v): pass
    def __enter__(self): return self
    def __exit__(self, *a): return False
    async def __aenter__(self): return self
    async def __aexit__(self, *a): return False
    def __add__(self, o): return self
    def __radd__(self, o): return self
    def __sub__(self, o): return self
    def __mul__(self, o): return self
    def __rmul__(self, o): return self
    def __truediv__(self, o): return self
    def __eq__(self, o): return isinstance(o, _AutoMock)
    def __hash__(self): return hash(self._name)
    def __lt__(self, o): return True
    def __le__(self, o): return True
    def __gt__(self, o): return False
    def __ge__(self, o): return False
    def __mro_entries__(self, bases): return (object,)

np = _AutoMock('np')
stats = _AutoMock('stats')

Pearson is the default everyone reaches for; Spearman is the smart default when you don't know the relationship's shape; Kendall is more conservative and slower but plays better with small samples.

Correlation ≠ Causation

The cautionary classics:

• US spending on science correlates almost perfectly with suicides by hanging (both rose monotonically 1999-2009). One does not cause the other.
• Ice cream sales correlate with drowning rates. Confounder: summer.
• Nicolas Cage films per year correlates with swimming-pool drownings. Coincidence at large scales is the default, not the exception.

Two correlated variables can share a cause, can have a coincidental match, or one can cause the other. Correlation alone can't tell you which. Causal inference (instrumental variables, regression discontinuity, randomised experiments) is how you go further — outside the scope of this lesson but well worth a follow-up.

9. The Multiple-Testing Problem

Run 20 independent tests at α = 0.05. Expected false positives even if nothing is real: 1. Run 100 tests, expect 5. This is why "we tested 47 features and found one significant" is not a discovery — it's the math working as designed.

Two corrections:

Bonferroni — divide α by the number of tests. Simple, very conservative.

alpha_corrected = 0.05 / n_tests

n_tests = 1

Benjamini-Hochberg (FDR) — controls expected proportion of false positives among rejections. Less conservative, more power:

from statsmodels.stats.multitest import multipletests
pvals = [0.001, 0.04, 0.03, 0.21, 0.5]
reject, p_adj, _, _ = multipletests(pvals, alpha=0.05, method="fdr_bh")
print(reject)        # [ True  True  True False False]
print(p_adj)         # adjusted p-values

Default to BH for exploratory work where you're scanning many features. Reserve Bonferroni for cases where any false positive is catastrophic (drug safety, hard-science publication).

10. Bootstrap — Get a CI for Anything

The bootstrap principle: if you don't know the sampling distribution of your statistic, resample your data with replacement many times, compute the statistic each time, and use the empirical spread.

import numpy as np

def bootstrap_median_ci(data, n_iter=10000, ci=95):
    rng = np.random.default_rng(0)
    boot_stats = [
        np.median(rng.choice(data, size=len(data), replace=True))
        for _ in range(n_iter)
    ]
    lo = np.percentile(boot_stats, (100 - ci) / 2)
    hi = np.percentile(boot_stats, 100 - (100 - ci) / 2)
    return lo, hi

data = np.array([2, 3, 5, 7, 11, 13, 17, 19, 23, 29])
bootstrap_median_ci(data)
# (4.0, 19.0) approximately

This works for any statistic — median, 95th percentile, ratio of means, AUC, custom KPI. No closed-form formula needed. The price is compute (n_iter resamples) — but with n_iter = 10,000 and modern hardware, it's seconds.

Limits: bootstrap needs the original sample to be a reasonable representation of the population. With n = 5, it's a fancy way of measuring noise. With n = 500, it's usually reliable.

11. Bayesian Intro — One Example

The Bayesian setup: combine prior belief with likelihood of observed data to get a posterior belief.

Example — a possibly-biased coin. We flip it 10 times and see 7 heads. Is it biased?

import numpy as np
from scipy import stats

# Prior: Beta(2, 2) — vague preference for fair coin, but open to evidence
# Likelihood: 7 heads / 10 flips → Beta(7, 3) update
# Posterior: Beta(2 + 7, 2 + 3) = Beta(9, 5) — conjugate update

posterior = stats.beta(9, 5)
print(f"Posterior mean: {posterior.mean():.2f}")
# 0.64

# 95% credible interval — actual probability the parameter is in this range
ci = posterior.interval(0.95)
print(f"95% credible interval: ({ci[0]:.2f}, {ci[1]:.2f})")
# (0.39, 0.86)

The Bayesian credible interval can be interpreted as "95% probability the true bias is between 0.39 and 0.86" — because in the Bayesian framework, the parameter is a random variable with a probability distribution. That's the interpretation people want from a frequentist CI and don't get.

The price: you have to specify a prior, and conjugate updates only work for nice prior/likelihood pairs. For non-trivial models, reach for pymc or numpyro.

12. The Practical Stats Checklist

Before any analysis:

• [ ] Look at the distribution — histogram, KDE, or at least describe(). Skewed? Multimodal? Outliers?
• [ ] Check sample size — n < 30 means be careful with normality-assuming tests; n < 5 means the test is mostly hope.
• [ ] Note base rates — a 10% conversion rate jumping to 11% (relative 10% lift) is very different from 10% jumping to 50%.
• [ ] Plan tests in advance — including how many you'll run, so you can correct correctly.
• [ ] Report effect size with every p-value — p < 0.05 is necessary, not sufficient.

After:

• [ ] CI on the difference, not just on each group separately.
• [ ] Sanity-check direction — did the metric move the way the theory predicts?
• [ ] Sensitivity analysis — does the conclusion change if you remove a handful of extreme rows?

13. Common Mistakes

1. Mean on heavy-tailed data
Reporting "mean response time = 4.2s" when 99% of responses are <1s and a long tail of timeouts skews the mean. Report median + 95th percentile.

2. p-hacking
Run 20 tests, report the one with p = 0.04 as if you'd planned to test it. By design, ~1 of 20 independent tests "significant" under the null. Pre-register hypotheses, or correct for multiple testing.

3. Ignoring effect size
With n = 10,000,000 a 0.01% difference is "highly significant" (p < 0.0001) and almost always meaningless. Report Cohen's d, relative lift, or the absolute difference in the metric you care about.

4. Peeking during an A/B test
Each look at the data is an implicit test. After 20 looks, your effective alpha is much higher than 0.05. Either commit to one analysis at the planned sample size, or use sequential testing methods (mSPRT, group-sequential, Bayesian).

5. Correlation as causation
"X correlates with Y" is the start of an investigation, not its end. Plausible confounders to consider every time: time/seasonality, selection bias, reverse causation, lurking third variables.

6. Bonferroni when you have 1000 tests
With many tests, Bonferroni is so conservative you'll miss every real effect. Use Benjamini-Hochberg for FDR control.

🎯 Your Turn — A General Bootstrap CI

Implement bootstrap_ci(data, statistic_fn, n_iter=1000, ci=95) that works for any statistic — not just the mean or median. The function:

• Resamples data with replacement, n_iter times.
• Computes statistic_fn(resample) each time.
• Returns (low, high) percentile-based CI bounds.

data = np.array([2, 3, 5, 7, 11, 13, 17, 19, 23, 29])

bootstrap_ci(data, np.median)                         # CI for median
bootstrap_ci(data, np.mean)                           # CI for mean
bootstrap_ci(data, lambda x: np.percentile(x, 90))    # CI for the 90th percentile
bootstrap_ci(data, lambda x: x.max() - x.min())       # CI for range — yes, it works

# Lightweight mock for objects whose attributes/methods aren't critical
class _AutoMock:
    def __init__(self, name='mock'): self._name = name
    def __getattr__(self, k): return _AutoMock(self._name + '.' + k)
    def __call__(self, *a, **kw):
        print('-> ' + self._name + '() called')
        return _AutoMock(self._name + '()')
    def __repr__(self): return '<mock ' + self._name + '>'
    def __str__(self): return '<mock ' + self._name + '>'
    def __bool__(self): return True
    def __iter__(self): return iter([])
    def __len__(self): return 0
    def __getitem__(self, k): return _AutoMock(self._name + '[...]')
    def __setitem__(self, k, v): pass
    def __enter__(self): return self
    def __exit__(self, *a): return False
    async def __aenter__(self): return self
    async def __aexit__(self, *a): return False
    def __add__(self, o): return self
    def __radd__(self, o): return self
    def __sub__(self, o): return self
    def __mul__(self, o): return self
    def __rmul__(self, o): return self
    def __truediv__(self, o): return self
    def __eq__(self, o): return isinstance(o, _AutoMock)
    def __hash__(self): return hash(self._name)
    def __lt__(self, o): return True
    def __le__(self, o): return True
    def __gt__(self, o): return False
    def __ge__(self, o): return False
    def __mro_entries__(self, bases): return (object,)

def bootstrap_ci(*_a, **_kw):
    print('-> bootstrap_ci() called')
    return _AutoMock('bootstrap_ci()')
np = _AutoMock('np')

Skeleton:

import numpy as np

def bootstrap_ci(data, statistic_fn, n_iter=1000, ci=95, random_state=0):
    rng = np.random.default_rng(random_state)
    data = np.asarray(data)

    # TODO 1: For n_iter iterations:
    #   - draw a resample of len(data) from data WITH replacement
    #   - compute statistic_fn on the resample
    #   - append to a list

    # TODO 2: Compute the (100-ci)/2 and 100-(100-ci)/2 percentiles
    #         of the bootstrap distribution
    ...

import numpy as np

def bootstrap_ci(data, statistic_fn, n_iter=1000, ci=95, random_state=0):
    """Percentile bootstrap CI for any statistic computed by `statistic_fn`.

    Parameters
    ----------
    data : array-like
        1D sample.
    statistic_fn : callable
        Function that takes an array and returns a scalar (e.g. np.median).
    n_iter : int
        Number of bootstrap resamples. 1000 is enough for most uses;
        10000 if you need stable estimates of extreme percentiles.
    ci : float
        Confidence level as a percentage (default 95).
    random_state : int
        Seed for reproducibility.

    Returns
    -------
    (low, high) : tuple of floats
        Percentile bootstrap interval.
    """
    rng = np.random.default_rng(random_state)
    data = np.asarray(data)
    n = len(data)

    boot_stats = np.empty(n_iter)
    for i in range(n_iter):
        resample = rng.choice(data, size=n, replace=True)
        boot_stats[i] = statistic_fn(resample)

    half = (100 - ci) / 2
    low  = np.percentile(boot_stats, half)
    high = np.percentile(boot_stats, 100 - half)
    return low, high


data = np.array([2, 3, 5, 7, 11, 13, 17, 19, 23, 29])

print("median:     ", bootstrap_ci(data, np.median))
print("mean:       ", bootstrap_ci(data, np.mean))
print("90th pct:   ", bootstrap_ci(data, lambda x: np.percentile(x, 90)))
print("range:      ", bootstrap_ci(data, lambda x: x.max() - x.min()))
print("trimmed mean:",
      bootstrap_ci(data,
                   lambda x: np.mean(np.sort(x)[1:-1])))   # drop min + max
# median:      (3.0, 19.0)
# mean:        (7.6, 18.1)
# 90th pct:    (16.4, 29.0)
# range:       (10.0, 27.0)
# trimmed mean: (6.625, 17.0)

What You Learned

• Descriptive stats: mean vs median vs mode — pick by data shape, report both mean and median.
• Five distributions: normal, log-normal, uniform, binomial, Poisson — recognise the shape, reach for the right tool.
• LLN + CLT: more data → more accurate; sample means become normal even when data isn't.
• Confidence intervals: stats.t.interval for the mean; the interpretation is about the procedure, not the specific interval.
• Hypothesis tests: Welch's t-test for two means, chi-squared for contingency tables. A p-value is "probability of data this extreme under the null".
• Effect size: Cohen's d. Report alongside every p-value. Big sample + tiny effect is significant but useless.
• A/B testing: design → run → analyse → decide. Sample size in advance, no peeking.
• Correlations: Pearson (linear), Spearman (monotonic), Kendall (robust). Correlation is not causation.
• Multiple testing: Bonferroni (conservative) or Benjamini-Hochberg (FDR, default for exploration).
• Bootstrap: percentile resampling gives a CI for any statistic without a closed-form formula.
• Bayesian basics: prior + likelihood = posterior; credible intervals support the intuitive interpretation.

Next: with cleaning, features, and stats in hand, you're ready for the modelling lessons — regression, trees, and clustering. Or jump straight into the ML track: ml-overview.