Statistics for Data Science

A working data scientist needs less statistics than a graduate textbook suggests, and different topics from a typical undergraduate course. This lesson covers the parts you will use in a real job: descriptive measures, distributions, tests and corrections, plus the correct way to interpret each result, so you don't publish wrong conclusions.

There is no measure theory and there are no long derivations. It is practical statistics, used in real work, with scipy code that runs.

1. Descriptive Stats: Pick the Right Centre

python

import numpy as np
import pandas as pd

x = pd.Series([1, 2, 2, 3, 3, 3, 100])
x.mean()       # 16.29  (pulled up by the outlier)
x.median()     # 3      (the outlier does not affect it)
x.mode()[0]    # 3
x.std()        # 36.92
x.var()        # 1363.24
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:

Statistic	When it's right
Mean	Symmetric distribution, no extreme outliers. Useful for budgets, totals (mean × n = sum).
Median	Skewed data, outliers present. Income, page-view counts, response times.
Mode	Categorical or discrete data where "most common" is the question.

Spread: variance squares each distance from the mean, so outliers have a very large effect (useful for maths, but hard to picture). Standard deviation is sqrt(var), so at least it is in the same units as the data. IQR is the robust choice when there are outliers.

Show both the mean and the median in any summary. The gap between them tells you something too: the larger the gap, the more skewed the distribution.

2. The Five Distributions You'll See Daily

Here is one short paragraph on each. You don't need to memorise the formulas. You need to recognise the shape and choose the right tool.

Normal (Gaussian): a bell curve, symmetric around the mean. Heights, measurement errors, sums of many small independent effects. Most parametric tests assume this distribution, so check that assumption before you use them.

Log-normal: looks like a normal distribution after you take logs. Right-skewed with a long tail. Incomes, file sizes, time between events, page-view counts per user. If describe() shows a mean much larger than the median, and a max many times bigger than both, suspect log-normal.

Uniform: flat between two bounds. Truly random IDs, dice rolls, the output of a well-seeded RNG. Real-world measurements are almost never uniform. If you see a flat histogram, it is probably caused by how the data was collected.

Binomial: the count of successes in n independent yes/no trials. Click-through events, conversion counts, A/B-test successes. For large n, it gets close to a normal distribution (because of the CLT).

Poisson: the count of events in a fixed interval, when events are independent and the rate is constant. Calls per hour, defects per batch, goals per match. The parameter λ is both the mean and the variance. If the variance of your count data is much larger than the mean, the data is "over-dispersed", and you probably want a negative binomial distribution instead.

python

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 all of inferential statistics is built on:

Law of Large Numbers: as your sample size grows, the sample mean converges to the true (population) mean. Larger samples give more accurate estimates. This sounds obvious, but everything else depends on it.

Central Limit Theorem: the distribution of sample means gets closer to a normal distribution as n grows, whatever the shape of the population. Even with very skewed data, if you draw samples of size 30 or more and take their means, the means 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:

python

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.64, 5.10)

Correct 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 that 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 second interpretation belongs to the Bayesian credible interval, which is a different procedure with a different meaning. See Section 11.)

The width of the CI shrinks with sqrt(n). To halve a CI, quadruple your sample size. Whether that is worth it is a question of budget, not of maths.

5. Hypothesis Testing in Practice

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

t-test: compare two means

python

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.76, p = 0.0002

+ setup added so this can run · re-runs 3 earlier blocks on this page first

# 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')

# Earlier blocks on this page that this one builds on, re-run quietly first.
try:
    import numpy as np
    import pandas as pd
    from scipy import stats
except Exception:
    pass
_pm_earlier = [
    'import numpy as np\nimport pandas as pd\n\nx = pd.Series([1, 2, 2, 3, 3, 3, 100])\nx.mean()       # 16.29  (pulled up by the outlier)\nx.median()     # 3      (the outlier does not affect it)\nx.mode()[0]    # 3\nx.std()        # 36.92\nx.var()        # 1363.24\nx.quantile([0.25, 0.5, 0.75])    # 2, 3, 3\nx.quantile(0.75) - x.quantile(0.25)   # IQR = 1',
    'from scipy import stats\nstats.norm.rvs(loc=0, scale=1, size=1000)\nstats.lognorm.rvs(s=0.5, size=1000)\nstats.binom.rvs(n=100, p=0.1, size=1000)\nstats.poisson.rvs(mu=3, size=1000)',
    'from scipy import stats\nimport numpy as np\n\ndata = np.array([4.2, 5.1, 4.8, 4.9, 5.0, 4.6, 5.3, 4.7, 5.2, 4.9])\nmean = data.mean()\nsem  = stats.sem(data)                                  # standard error of the mean\nci   = stats.t.interval(0.95, df=len(data) - 1,\n                        loc=mean, scale=sem)\nprint(f"mean = {mean:.2f}, 95% CI = ({ci[0]:.2f}, {ci[1]:.2f})")\n# mean = 4.87, 95% CI = (4.64, 5.10)',
]
def _pm_run_earlier(_srcs, _g):
    import sys as _s, io as _io, warnings as _w
    for _src in _srcs:
        _o, _e = _s.stdout, _s.stderr
        _s.stdout = _s.stderr = _io.StringIO()
        try:
            with _w.catch_warnings():
                _w.simplefilter('ignore')
                exec(compile(_src, '<earlier block>', 'exec'), _g)
        except (Exception, SystemExit):
            pass
        finally:
            _s.stdout, _s.stderr = _o, _e
    if 'matplotlib.pyplot' in _s.modules:
        _s.modules['matplotlib.pyplot'].close('all')   # their plots belong to their own blocks
_pm_run_earlier(_pm_earlier, globals())
del _pm_earlier, _pm_run_earlier

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

Chi-squared: compare distributions / contingency tables

python

# 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}")

+ setup added so this can run · re-runs 4 earlier blocks on this page first

# 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')

# Earlier blocks on this page that this one builds on, re-run quietly first.
try:
    import numpy as np
    import pandas as pd
    from scipy import stats
except Exception:
    pass
_pm_earlier = [
    'import numpy as np\nimport pandas as pd\n\nx = pd.Series([1, 2, 2, 3, 3, 3, 100])\nx.mean()       # 16.29  (pulled up by the outlier)\nx.median()     # 3      (the outlier does not affect it)\nx.mode()[0]    # 3\nx.std()        # 36.92\nx.var()        # 1363.24\nx.quantile([0.25, 0.5, 0.75])    # 2, 3, 3\nx.quantile(0.75) - x.quantile(0.25)   # IQR = 1',
    'from scipy import stats\nstats.norm.rvs(loc=0, scale=1, size=1000)\nstats.lognorm.rvs(s=0.5, size=1000)\nstats.binom.rvs(n=100, p=0.1, size=1000)\nstats.poisson.rvs(mu=3, size=1000)',
    'from scipy import stats\nimport numpy as np\n\ndata = np.array([4.2, 5.1, 4.8, 4.9, 5.0, 4.6, 5.3, 4.7, 5.2, 4.9])\nmean = data.mean()\nsem  = stats.sem(data)                                  # standard error of the mean\nci   = stats.t.interval(0.95, df=len(data) - 1,\n                        loc=mean, scale=sem)\nprint(f"mean = {mean:.2f}, 95% CI = ({ci[0]:.2f}, {ci[1]:.2f})")\n# mean = 4.87, 95% CI = (4.64, 5.10)',
    'control   = np.array([4.2, 5.1, 4.8, 4.9, 5.0, 4.6, 5.3, 4.7, 5.2, 4.9])\ntreatment = np.array([5.4, 5.6, 5.1, 5.7, 5.3, 5.8, 5.2, 5.5, 5.4, 5.9])\n\nt, p = stats.ttest_ind(control, treatment, equal_var=False)   # Welch\'s t-test\nprint(f"t = {t:.2f}, p = {p:.4f}")\n# t = -4.76, p = 0.0002',
]
def _pm_run_earlier(_srcs, _g):
    import sys as _s, io as _io, warnings as _w
    for _src in _srcs:
        _o, _e = _s.stdout, _s.stderr
        _s.stdout = _s.stderr = _io.StringIO()
        try:
            with _w.catch_warnings():
                _w.simplefilter('ignore')
                exec(compile(_src, '<earlier block>', 'exec'), _g)
        except (Exception, SystemExit):
            pass
        finally:
            _s.stdout, _s.stderr = _o, _e
    if 'matplotlib.pyplot' in _s.modules:
        _s.modules['matplotlib.pyplot'].close('all')   # their plots belong to their own blocks
_pm_run_earlier(_pm_earlier, globals())
del _pm_earlier, _pm_run_earlier

What a p-value means

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". (A real but small effect can have a large p-value when n is small. A trivial effect can have a tiny p-value when n is huge.)
A measure of effect size. (See Section 6.)
A pass/fail gate. (p < 0.05 is a convention, not a law of nature.)

The careful reading: a small p-value says "the data doesn't fit the null hypothesis well". Nothing more. What you do next depends on how plausible the effect was to begin with, the effect size, and the 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 is the standardised difference between two means:

python

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.13  (a very large effect)

+ setup added so this can run · re-runs 5 earlier blocks on this page first

# 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')
stats = _AutoMock('stats')

# Earlier blocks on this page that this one builds on, re-run quietly first.
try:
    import numpy as np
    import pandas as pd
    from scipy import stats
except Exception:
    pass
_pm_earlier = [
    'import numpy as np\nimport pandas as pd\n\nx = pd.Series([1, 2, 2, 3, 3, 3, 100])\nx.mean()       # 16.29  (pulled up by the outlier)\nx.median()     # 3      (the outlier does not affect it)\nx.mode()[0]    # 3\nx.std()        # 36.92\nx.var()        # 1363.24\nx.quantile([0.25, 0.5, 0.75])    # 2, 3, 3\nx.quantile(0.75) - x.quantile(0.25)   # IQR = 1',
    'from scipy import stats\nstats.norm.rvs(loc=0, scale=1, size=1000)\nstats.lognorm.rvs(s=0.5, size=1000)\nstats.binom.rvs(n=100, p=0.1, size=1000)\nstats.poisson.rvs(mu=3, size=1000)',
    'from scipy import stats\nimport numpy as np\n\ndata = np.array([4.2, 5.1, 4.8, 4.9, 5.0, 4.6, 5.3, 4.7, 5.2, 4.9])\nmean = data.mean()\nsem  = stats.sem(data)                                  # standard error of the mean\nci   = stats.t.interval(0.95, df=len(data) - 1,\n                        loc=mean, scale=sem)\nprint(f"mean = {mean:.2f}, 95% CI = ({ci[0]:.2f}, {ci[1]:.2f})")\n# mean = 4.87, 95% CI = (4.64, 5.10)',
    'control   = np.array([4.2, 5.1, 4.8, 4.9, 5.0, 4.6, 5.3, 4.7, 5.2, 4.9])\ntreatment = np.array([5.4, 5.6, 5.1, 5.7, 5.3, 5.8, 5.2, 5.5, 5.4, 5.9])\n\nt, p = stats.ttest_ind(control, treatment, equal_var=False)   # Welch\'s t-test\nprint(f"t = {t:.2f}, p = {p:.4f}")\n# t = -4.76, p = 0.0002',
    '# Did clicks differ between two ad creatives?\n#                   Creative A   Creative B\n# clicked              120          150\n# didn\'t click         880          850\ntable = np.array([[120, 150], [880, 850]])\n\nchi2, p, dof, expected = stats.chi2_contingency(table)\nprint(f"chi² = {chi2:.2f}, p = {p:.4f}")',
]
def _pm_run_earlier(_srcs, _g):
    import sys as _s, io as _io, warnings as _w
    for _src in _srcs:
        _o, _e = _s.stdout, _s.stderr
        _s.stdout = _s.stderr = _io.StringIO()
        try:
            with _w.catch_warnings():
                _w.simplefilter('ignore')
                exec(compile(_src, '<earlier block>', 'exec'), _g)
        except (Exception, SystemExit):
            pass
        finally:
            _s.stdout, _s.stderr = _o, _e
    if 'matplotlib.pyplot' in _s.modules:
        _s.modules['matplotlib.pyplot'].close('all')   # their plots belong to their own blocks
_pm_run_earlier(_pm_earlier, globals())
del _pm_earlier, _pm_run_earlier

The usual benchmarks (from Cohen himself, so treat them as rough guides):

d	Magnitude
0.2	small
0.5	medium
0.8	large

Always report effect size alongside the p-value. "p = 0.001, d = 0.05" means a real but trivial effect: interesting at large scale, but probably not worth acting on. "p = 0.04, d = 1.2" means a large effect that only just passed a test with low power. It is 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: assign users randomly, split them roughly equally, and don't touch the test until the sample size is reached. Don't check the results early. Stopping early when a result looks significant greatly increases the false-positive rate.

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 the effect size, the CI and the cost of rolling out the change. A 0.5% lift might be worth a lot (on a very large checkout page) or be just noise (on a small SaaS dashboard).

Beginners usually fail at steps 1 and 2: they start without a sample-size plan, then check the results every day and stop as soon as the result looks significant.

8. Correlation: Three Kinds

python

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

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

+ setup added so this can run · re-runs 6 earlier blocks on this page first

# 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')
control = _AutoMock('control')
treatment = _AutoMock('treatment')

# Earlier blocks on this page that this one builds on, re-run quietly first.
try:
    import numpy as np
    import pandas as pd
    from scipy import stats
except Exception:
    pass
_pm_earlier = [
    'import numpy as np\nimport pandas as pd\n\nx = pd.Series([1, 2, 2, 3, 3, 3, 100])\nx.mean()       # 16.29  (pulled up by the outlier)\nx.median()     # 3      (the outlier does not affect it)\nx.mode()[0]    # 3\nx.std()        # 36.92\nx.var()        # 1363.24\nx.quantile([0.25, 0.5, 0.75])    # 2, 3, 3\nx.quantile(0.75) - x.quantile(0.25)   # IQR = 1',
    'from scipy import stats\nstats.norm.rvs(loc=0, scale=1, size=1000)\nstats.lognorm.rvs(s=0.5, size=1000)\nstats.binom.rvs(n=100, p=0.1, size=1000)\nstats.poisson.rvs(mu=3, size=1000)',
    'from scipy import stats\nimport numpy as np\n\ndata = np.array([4.2, 5.1, 4.8, 4.9, 5.0, 4.6, 5.3, 4.7, 5.2, 4.9])\nmean = data.mean()\nsem  = stats.sem(data)                                  # standard error of the mean\nci   = stats.t.interval(0.95, df=len(data) - 1,\n                        loc=mean, scale=sem)\nprint(f"mean = {mean:.2f}, 95% CI = ({ci[0]:.2f}, {ci[1]:.2f})")\n# mean = 4.87, 95% CI = (4.64, 5.10)',
    'control   = np.array([4.2, 5.1, 4.8, 4.9, 5.0, 4.6, 5.3, 4.7, 5.2, 4.9])\ntreatment = np.array([5.4, 5.6, 5.1, 5.7, 5.3, 5.8, 5.2, 5.5, 5.4, 5.9])\n\nt, p = stats.ttest_ind(control, treatment, equal_var=False)   # Welch\'s t-test\nprint(f"t = {t:.2f}, p = {p:.4f}")\n# t = -4.76, p = 0.0002',
    '# Did clicks differ between two ad creatives?\n#                   Creative A   Creative B\n# clicked              120          150\n# didn\'t click         880          850\ntable = np.array([[120, 150], [880, 850]])\n\nchi2, p, dof, expected = stats.chi2_contingency(table)\nprint(f"chi² = {chi2:.2f}, p = {p:.4f}")',
    'def cohens_d(a, b):\n    """Standardised mean difference. Pooled std under equal-n assumption."""\n    pooled = np.sqrt((a.var(ddof=1) + b.var(ddof=1)) / 2)\n    return (b.mean() - a.mean()) / pooled\n\ncohens_d(control, treatment)\n# 2.13  (a very large effect)',
]
def _pm_run_earlier(_srcs, _g):
    import sys as _s, io as _io, warnings as _w
    for _src in _srcs:
        _o, _e = _s.stdout, _s.stderr
        _s.stdout = _s.stderr = _io.StringIO()
        try:
            with _w.catch_warnings():
                _w.simplefilter('ignore')
                exec(compile(_src, '<earlier block>', 'exec'), _g)
        except (Exception, SystemExit):
            pass
        finally:
            _s.stdout, _s.stderr = _o, _e
    if 'matplotlib.pyplot' in _s.modules:
        _s.modules['matplotlib.pyplot'].close('all')   # their plots belong to their own blocks
_pm_run_earlier(_pm_earlier, globals())
del _pm_earlier, _pm_run_earlier

Coefficient	Measures	Use when
Pearson r	Linear association	Both variables roughly normal, linear relationship expected.
Spearman ρ	Rank/monotonic association	Non-linear-but-monotonic, ordinal data, outliers present.
Kendall τ	Pairwise concordance	Small samples, lots of ties, you want a robust alternative to Spearman.

Pearson is the one most people use by default. Spearman is a better default when you don't know the shape of the relationship. Kendall is more conservative and slower, but works better with small samples.

Correlation ≠ Causation

Well-known warning examples:

US spending on science correlates almost perfectly with suicides by hanging (both rose from 1999 to 2009). One does not cause the other.
Ice cream sales correlate with drowning rates. Confounder: summer.
The number of Nicolas Cage films per year correlates with swimming-pool drownings. When you compare thousands of unrelated data series, some of them will match by chance.

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. It is outside the scope of this lesson, but worth learning next.

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 is exactly what chance alone would produce.

Two corrections:

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

python

alpha_corrected = 0.05 / n_tests

+ setup added so this can run · defines n_tests

n_tests = 1

Benjamini-Hochberg (FDR): controls the expected proportion of false positives among the results you call significant. It is less conservative and has more power:

python

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 False False False False]
print(p_adj)         # adjusted p-values

Default to BH for exploratory work where you're scanning many features. Use Bonferroni when even one false positive would be very costly (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 spread of those results.

python

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, a custom KPI. You don't need a formula. The cost is computing time (n_iter resamples), but with n_iter = 10,000 on a modern computer, it takes seconds.

Limits: bootstrap needs the original sample to be a reasonable representation of the population. With n = 5, it mostly measures 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 coin that might be biased. We flip it 10 times and see 7 heads. Is it biased?

python

import numpy as np
from scipy import stats

# Prior: Beta(2, 2), a weak belief that the coin is fair, which data can change
# Data: 7 heads in 10 flips, so we add 7 heads and 3 tails
# Posterior: Beta(2 + 7, 2 + 3) = Beta(9, 5), using the conjugate update rule

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

# 95% credible interval: a 95% probability that the true value 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 read as "95% probability that 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. People often read a frequentist CI this way, but that reading is wrong.

The cost: you have to choose a prior, and conjugate updates only work for certain prior/likelihood pairs. For more complex models, use 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: with n < 30, be careful with tests that assume normality. With n < 5, the test tells you very little.
[ ] 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 that you can apply the right correction.
[ ] Report effect size with every p-value. p < 0.05 alone is not enough.

After:

[ ] CI on the difference, not just on each group separately.
[ ] Check the 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. When nothing is real, about 1 in 20 independent tests will still look "significant". Pre-register hypotheses, or correct for multiple testing.

3. Ignoring effect size
With n = 10,000,000, a very small difference can be "highly significant" and still be meaningless in practice. 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 a hidden extra 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 the end. Possible explanations to check 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 only 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.

python

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 the range (this works too)

+ setup added so this can run · defines bootstrap_ci, np

# 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:

python

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
    ...

Hint 1: Resampling

rng.choice(data, size=len(data), replace=True) draws a bootstrap sample. replace=True is essential. Without it, every "resample" would only be the original data in a different order.

Hint 2: Vectorise if you want

The naïve loop works fine. For a faster version, pre-build a 2D array of indices with rng.integers(0, len(data), size=(n_iter, len(data))), index data with it, and call statistic_fn along axis=1, but only if statistic_fn is vectorised. The loop version works for any statistic. The vectorised version is faster for built-in NumPy statistics.

Show full solution

python

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 (for example, 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 round(float(low), 2), round(float(high), 2)


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 the min and the max
# median:      (5.0, 20.0)
# mean:        (8.0, 18.5)
# 90th pct:    (13.6, 29.0)
# range:       (14.98, 27.0)
# trimmed mean: (6.87, 19.0)

What the implementation gets right:

Works with any statistic: any callable that returns a scalar works. Median, mean, percentile, range, AUC, the gap between mean and median, anything.
Reproducible: random_state defaults to a fixed seed. Random methods that give different results on every run are very hard to debug.
Pre-allocated array instead of list.append, which is a little faster.
Percentile method is the simplest kind of bootstrap CI. Two refinements you might add later: bias-corrected (BC) and bias-corrected-accelerated (BCa). scipy.stats.bootstrap implements both and handles edge cases (degenerate resamples, small n) more carefully.

For production, scipy.stats.bootstrap(data, statistic_fn, confidence_level=0.95, method="BCa") is the standard choice. The function above is a teaching version: the same idea in less code, so you can see what happens.

A quick check: notice that the range statistic has a CI that doesn't include 0 (you can't have negative range). And the trimmed mean CI comes out slightly wider than the mean's here: trimming doesn't always tighten an interval, especially on 10 points. The bootstrap works with whatever statistic you give it, and you never have to derive a formula for its sampling distribution.

What You Learned

Descriptive stats: mean vs median vs mode. Choose by the shape of the data, and report both mean and median.
Five distributions: normal, log-normal, uniform, binomial, Poisson. Recognise the shape and choose the right tool.
LLN + CLT: more data gives more accurate estimates, and 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. A big sample with a tiny effect can be significant but useless.
A/B testing: design → run → analyse → decide. Fix the sample size in advance, and don't check the results early.
Correlations: Pearson (linear), Spearman (monotonic), Kendall (robust). Correlation is not causation.
Multiple testing: Bonferroni (conservative) or Benjamini-Hochberg (FDR, default for exploration).
Bootstrap: resampling with the percentile method gives a CI for any statistic without a closed-form formula.
Bayesian basics: prior + likelihood = posterior, and credible intervals can be read in the intuitive way.

Next: now that you know cleaning, features and statistics, you're ready for the modelling lessons: regression, trees, and clustering. Or jump straight into the ML track: ml-overview.

1 · The lesson

1. Descriptive Stats: Pick the Right Centre

2. The Five Distributions You'll See Daily

3. Sampling, LLN, and the CLT

4. Confidence Intervals

5. Hypothesis Testing in Practice

t-test: compare two means

Chi-squared: compare distributions / contingency tables

What a p-value means

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

7. A/B Testing Flow

8. Correlation: Three Kinds

Correlation ≠ Causation

9. The Multiple-Testing Problem

10. Bootstrap: Get a CI for Anything

11. Bayesian Intro: One Example

12. The Practical Stats Checklist

13. Common Mistakes

🎯 Your Turn: A General Bootstrap CI

What You Learned

Practice this