Inference
t-tests, one-way ANOVA, confidence intervals, the correlation test and multiple-comparison corrections, each returning a statistic or a p-value.
This module turns data into the numbers a hypothesis test reports. Every procedure follows the same pattern: compute one statistic from the data (a t ratio, an F ratio), then read it against a distribution from Prelude.Distributions. A p-value is a cumulative distribution function read backwards, \(2\,(1 - F(\lvert t \rvert))\) for a two-sided test and \(1 - F(t)\) for a one-sided one; a confidence interval uses the quantile function to find the critical value.
Conventions. Samples x and y are vectors. Variances are sample variances (divided by \(n - 1\)), exactly as in Prelude.Statistics. A coverage c is a probability such as 0.95; a significance level a is a probability such as 0.05. For analysis of variance the groups are given as a one-hot membership matrix g of shape n×G: row i has a 1 in the column of the group observation i belongs to, and zeros elsewhere. A collection of p-values for a correction is a vector.
Limitations. The p-values and critical values are only as accurate as the Student t, F and normal functions in Prelude.Distributions. The proportion interval is the Wald interval, which is known to cover poorly when the proportion is near 0 or 1 or the sample is small. The Welch test uses the Welch–Satterthwaite degrees of freedom, which are fractional; the Student t distribution accepts a real number of degrees of freedom, so no rounding takes place.
open Prelude.Inference
One sample and paired samples
The one-sample t statistic compares a sample mean with a hypothesised mean \(\mu\), in units of the standard error:
\[t = \frac{\bar{x} - \mu}{s / \sqrt{n}}, \qquad \nu = n - 1\]One degree of freedom is spent on the mean, for the same reason the sample variance divides by \(n - 1\). A paired test is a one-sample test of the differences \(x_i - y_i\) against zero.
twoSided
twoSided(v, t) · v degrees of freedom scalar · t t statistic scalar → p-value scalar
Let \(\nu > 0\) (the argument v) be a number of degrees of freedom and let \(t \in \mathbb{R}\). Let \(F_{\nu}(u) = \operatorname{studentCdf}(\nu, u)\) be the Student t distribution function on \(\nu\) degrees of freedom. Then
The two-sided p-value of a t statistic on v degrees of freedom: the probability of a value at least as far from zero as t, in either direction.
Intuition A p-value asks: if there were really no effect, how often would chance alone produce a statistic this extreme? “Extreme” here counts both tails, so the area beyond \(\lvert t \rvert\) on one side is doubled. Fewer degrees of freedom mean heavier tails, so the same t gives a larger p-value in a small sample.
Use it when you already have a t statistic and its degrees of freedom from anywhere, for example a regression coefficient’s statistics with \(n-p\), and you have no advance reason to expect one direction. It is the default choice. If the hypothesis fixed the direction before seeing the data, use oneSided.
hypothesis-testingsmall-samples
def twoSided(v, t) ≔ 2 · (1 − studentCdf(v, abs(t)))
oneSided
oneSided(v, t) · v degrees of freedom scalar · t t statistic scalar → p-value scalar
Let \(\nu > 0\) (the argument v) be a number of degrees of freedom, let \(t \in \mathbb{R}\), and let \(F_{\nu}\) be the Student t distribution function, as in twoSided. Then
The upper-tail p-value: the probability of a value at least as large as t. For the lower tail, pass -t.
Intuition Only one tail counts as evidence: the probability of a value at least as large as t. When the effect goes in the predicted direction the p-value is half the two-sided one, and when it goes the other way it is close to 1.
Use it when the question is directional and fixed in advance: “is the new version faster?”, “does the dose raise the level?”. Choosing the side after looking at the data halves the p-value dishonestly, so if in doubt use twoSided. The pOneSample, pPooled, pWelch and pCorrelation functions are all two-sided; for a one-sided version feed their t statistic and degrees of freedom to this function.
hypothesis-testingdirectional-hypotheses
def oneSided(v, t) ≔ 1 − studentCdf(v, t)
tOneSample
tOneSample(x, m) · x sample vector n · m hypothesised mean scalar → t statistic scalar
Let \(x = (x_1, \dots, x_n) \in \mathbb{R}^n\) be a sample with \(n \ge 2\), and let \(\mu \in \mathbb{R}\) (the argument m) be the hypothesised mean. Then
where \(\bar{x} = \frac{1}{n} \sum_{i=1}^{n} x_i\) and \(s^2 = \frac{1}{n-1} \sum_{i=1}^{n} (x_i - \bar{x})^2\).
The one-sample t statistic: the distance of the sample mean from m, divided by the standard error of the mean.
Intuition The gap between what you measured and what you expected, measured in units of how much the sample mean itself would wobble from sample to sample. The standard error shrinks as \(1/\sqrt{n}\), so the same gap becomes more convincing with more data.
Use it when you compare one sample’s mean with a fixed reference value: a target, a specification, a published figure. It assumes independent observations and a roughly normal mean, which holds for moderate n unless the data are very skewed. For the p-value use pOneSample. For skewed small samples a rank or resampling method (wilcoxon, bootLower in Resampling) is safer.
hypothesis-testingcomparing-to-a-referencesmall-samples
def tOneSample(x, m) ≔ (mean(x) − m) / error(x)
pOneSample
pOneSample(x, m) · x sample vector n · m hypothesised mean scalar → p-value scalar
Let \(x \in \mathbb{R}^n\) with \(n \ge 2\), let \(\mu \in \mathbb{R}\) (the argument m), let \(t = \operatorname{tOneSample}(x, \mu)\), and let \(F_{\nu}\) be the Student t distribution function on \(\nu\) degrees of freedom. Then
The two-sided p-value of the one-sample t-test on \(n - 1\) degrees of freedom.
Intuition Turns tOneSample into a probability: if the true mean really were m, how often would a sample of this size land at least this far from it? A small value says the data are hard to reconcile with m.
Use it when you need a yes/no style decision against a reference mean, typically rejecting at 0.05. A p-value is not the size of the effect: with large n a trivial difference is “significant”, so report the interval from lowerBound and upperBound alongside it. For a directional test use oneSided(n − 1, tOneSample(x, m)).
hypothesis-testingcomparing-to-a-reference
def pOneSample(x, m) ≔ twoSided(freedomOf(x), tOneSample(x, m))
tPaired
tPaired(x, y) · x first measurements vector n · y second measurements vector n → t statistic scalar
Let \(x, y \in \mathbb{R}^n\), \(n \ge 2\), be paired measurements, and let \(d = x - y\). Then
\[\operatorname{tPaired}(x, y) = \frac{\bar{d}}{s_d / \sqrt{n}} ,\]where \(\bar{d} = \frac{1}{n} \sum_{i=1}^{n} d_i\) and \(s_d^2 = \frac{1}{n-1} \sum_{i=1}^{n} (d_i - \bar{d})^2\).
The paired t statistic: the one-sample statistic of the differences x − y against zero. The two vectors must be the same length and in matching order.
Intuition When the same subject is measured twice, most of the variation is between subjects, not between conditions. Subtracting within each pair cancels that shared part, so only the change is left, and the test is far more sensitive than comparing the two columns as unrelated groups.
Use it when each value in x has a natural partner in y: before and after treatment, two methods on the same specimens, left and right, twins. Do not use it for unrelated groups; use tWelch. If the differences are strongly skewed or have outliers, the rank-based wilcoxon in Resampling is more robust.
before-after-comparisoncomparing-groupshypothesis-testing
def tPaired(x, y) ≔ tOneSample(x − y, 0)
pPaired
pPaired(x, y) · x first measurements vector n · y second measurements vector n → p-value scalar
Let \(x, y \in \mathbb{R}^n\), \(n \ge 2\), be paired measurements and let \(t = \operatorname{tPaired}(x, y)\). Then
\[\operatorname{pPaired}(x, y) = 2\,\bigl(1 - F_{n-1}(\lvert t \rvert)\bigr) .\]The two-sided p-value of the paired t-test, on \(n - 1\) degrees of freedom.
Intuition The probability of seeing an average change at least this large, in either direction, if the two conditions really made no difference. Degrees of freedom come from the number of pairs, not the number of values.
Use it when you decide whether a before/after or matched-pairs change is real. To report how large the change is, compute an interval on the differences with lowerBound(x − y, c) and upperBound(x − y, c). For unrelated samples use pWelch; for a rank-based alternative on the differences, wilcoxon in Resampling.
before-after-comparisonhypothesis-testing
def pPaired(x, y) ≔ pOneSample(x − y, 0)
Building blocks
freedomOf(x)— Let \(x \in \mathbb{R}^n\). Then \(\operatorname{freedomOf}(x) = n - 1\). The degrees of freedom of one sample; used bypOneSample, the Welch test and the intervals.
Two independent samples
Two tests compare the means of independent samples. The pooled (Student) test assumes the two populations share one variance and estimates it from both samples, on \(n_1 + n_2 - 2\) degrees of freedom:
\[s_p^2 = \frac{(n_1 - 1)\,s_1^2 + (n_2 - 1)\,s_2^2}{n_1 + n_2 - 2}, \qquad t = \frac{\bar{x} - \bar{y}}{\sqrt{s_p^2\,(1/n_1 + 1/n_2)}}\]The Welch test drops that assumption. It divides the same difference by a different standard error, and its degrees of freedom come from the Welch–Satterthwaite approximation:
\[t = \frac{\bar{x} - \bar{y}}{\sqrt{s_1^2/n_1 + s_2^2/n_2}}, \qquad \nu = \frac{\left(s_1^2/n_1 + s_2^2/n_2\right)^2}{\dfrac{(s_1^2/n_1)^2}{n_1 - 1} + \dfrac{(s_2^2/n_2)^2}{n_2 - 1}}\]When in doubt about equal variances, prefer Welch.
tPooled
tPooled(x, y) · x first sample vector n₁ · y second sample vector n₂ → t statistic scalar
Let \(x \in \mathbb{R}^{n_1}\) and \(y \in \mathbb{R}^{n_2}\) be independent samples with \(n_1, n_2 \ge 2\), means \(\bar{x}, \bar{y}\) and sample variances \(s_1^2 = \frac{1}{n_1 - 1} \sum_{i} (x_i - \bar{x})^2\) and \(s_2^2 = \frac{1}{n_2 - 1} \sum_{i} (y_i - \bar{y})^2\). Then
\[\operatorname{tPooled}(x, y) = \frac{\bar{x} - \bar{y}}{\sqrt{s_p^2 \left( \frac{1}{n_1} + \frac{1}{n_2} \right)}} , \qquad s_p^2 = \frac{(n_1 - 1)\, s_1^2 + (n_2 - 1)\, s_2^2}{n_1 + n_2 - 2} .\]The pooled two-sample t statistic. The samples may have different lengths.
Intuition The difference between two group means, in units of how much that difference would vary by chance. If both groups share one true variance, the best estimate of it averages the two sample variances weighted by their degrees of freedom, and every observation helps estimate the noise.
Use it when you compare two independent groups and have good reason to believe their variances are equal, for example the same measurement process on randomised groups. It is the classic Student test. If the spreads or group sizes differ noticeably, it can give badly wrong p-values; tWelch costs almost nothing when variances are equal and is safe when they are not.
comparing-groupshypothesis-testing
def tPooled(x, y) ≔ (mean(x) − mean(y)) / sqrt(pooledVariance(x, y) · (1 / count(x) + 1 / count(y)))
pPooled
pPooled(x, y) · x first sample vector n₁ · y second sample vector n₂ → p-value scalar
Let \(x \in \mathbb{R}^{n_1}\) and \(y \in \mathbb{R}^{n_2}\) be independent samples with \(n_1, n_2 \ge 2\), and let \(t = \operatorname{tPooled}(x, y)\). Then
\[\operatorname{pPooled}(x, y) = 2\,\bigl(1 - F_{n_1 + n_2 - 2}(\lvert t \rvert)\bigr) .\]The two-sided p-value of the pooled t-test, on \(n_1 + n_2 - 2\) degrees of freedom.
Intuition How often two groups drawn from the same population, with the same variance, would show a mean gap this large. Pooling both samples gives the most degrees of freedom any two-sample t-test can have, which is where its extra power comes from.
Use it when the equal-variance assumption is justified and you want the textbook Student test, for instance to match a published analysis. Otherwise prefer pWelch. It is also the two-group special case of one-way ANOVA: with two groups pAnova gives the same p-value. For non-normal data with outliers, mannWhitney or permP in Resampling.
comparing-groupshypothesis-testing
def pPooled(x, y) ≔ twoSided(freedomPooled(x, y), tPooled(x, y))
tWelch
tWelch(x, y) · x first sample vector n₁ · y second sample vector n₂ → t statistic scalar
Let \(x \in \mathbb{R}^{n_1}\) and \(y \in \mathbb{R}^{n_2}\) be independent samples with \(n_1, n_2 \ge 2\), means \(\bar{x}, \bar{y}\) and sample variances \(s_1^2, s_2^2\) (divided by \(n_1 - 1\) and \(n_2 - 1\)). Then
\[\operatorname{tWelch}(x, y) = \frac{\bar{x} - \bar{y}}{\sqrt{s_1^2 / n_1 + s_2^2 / n_2}} .\]The Welch t statistic, which does not assume equal variances.
Intuition Each mean carries its own uncertainty, \(s^2/n\), and the uncertainty of their difference is just the two added together. No shared variance is assumed, so a small noisy group cannot borrow false precision from a large quiet one.
Use it when you compare the means of two independent groups: A/B tests, treatment versus control, two machines or two batches. It is the sensible default for two samples. The observations must be independent within and between groups; for matched measurements use tPaired. It pairs with freedomWelch for its degrees of freedom, and pWelch gives the p-value directly.
comparing-groupshypothesis-testingunequal-variance
def tWelch(x, y) ≔ (mean(x) − mean(y)) / sqrt(spreadOf(x) + spreadOf(y))
freedomWelch
freedomWelch(x, y) · x first sample vector n₁ · y second sample vector n₂ → degrees of freedom scalar
Let \(x \in \mathbb{R}^{n_1}\) and \(y \in \mathbb{R}^{n_2}\) be samples with \(n_1, n_2 \ge 2\) and sample variances \(s_1^2, s_2^2\), and let \(a = s_1^2 / n_1\) and \(b = s_2^2 / n_2\). Then
\[\operatorname{freedomWelch}(x, y) = \frac{(a + b)^2}{\dfrac{a^2}{n_1 - 1} + \dfrac{b^2}{n_2 - 1}} .\]The Welch–Satterthwaite degrees of freedom. The result is generally not a whole number.
Intuition How much information the combined standard error really rests on. When one group’s variance term dominates, the error is effectively estimated from that group alone, and the value falls toward its \(n - 1\); when both contribute equally it rises toward \(n_1 + n_2 - 2\).
Use it when you need the reference distribution for tWelch, for example to call oneSided or to build a Welch confidence interval with critical(freedomWelch(x, y), c) times \(\sqrt{s_1^2/n_1 + s_2^2/n_2}\). A value much smaller than \(n_1 + n_2 - 2\) is itself a sign that the pooled test would be unreliable.
unequal-varianceconfidence-intervalssmall-samples
def freedomWelch(x, y) ≔ (spreadOf(x) + spreadOf(y)) · (spreadOf(x) + spreadOf(y)) / (spreadOf(x) · spreadOf(x) / freedomOf(x) + spreadOf(y) · spreadOf(y) / freedomOf(y))
pWelch
pWelch(x, y) · x first sample vector n₁ · y second sample vector n₂ → p-value scalar
Let \(x \in \mathbb{R}^{n_1}\) and \(y \in \mathbb{R}^{n_2}\) be independent samples with \(n_1, n_2 \ge 2\), let \(t = \operatorname{tWelch}(x, y)\) and \(\nu = \operatorname{freedomWelch}(x, y)\). Then
\[\operatorname{pWelch}(x, y) = 2\,\bigl(1 - F_{\nu}(\lvert t \rvert)\bigr) .\]The two-sided p-value of Welch’s t-test.
Intuition The chance that two independent groups with the same true mean, but possibly different spreads, would show a mean gap at least this large. It keeps its stated error rate whether or not the variances match.
Use it when you want one p-value for “do these two groups differ on average?”, which is the usual two-sample question. When there are more than two groups, testing every pair inflates false positives: run pAnova first, or correct the pairwise p-values with bonferroni or discoveries. For heavy-tailed data, mannWhitney or permP in Resampling.
comparing-groupshypothesis-testingunequal-variance
def pWelch(x, y) ≔ twoSided(freedomWelch(x, y), tWelch(x, y))
Building blocks
pooledVariance(x, y)— Let \(x \in \mathbb{R}^{n_1}\) and \(y \in \mathbb{R}^{n_2}\) with sample variances \(s_1^2, s_2^2\). Then \(\operatorname{pooledVariance}(x, y) = \frac{(n_1 - 1)\, s_1^2 + (n_2 - 1)\, s_2^2}{n_1 + n_2 - 2}\). The variance estimated from both samples together; used bytPooled.freedomPooled(x, y)— Let \(x \in \mathbb{R}^{n_1}\) and \(y \in \mathbb{R}^{n_2}\). Then \(\operatorname{freedomPooled}(x, y) = n_1 + n_2 - 2\). The pooled degrees of freedom; used bypPooled.spreadOf(x)— Let \(x \in \mathbb{R}^n\) with sample variance \(s^2\). Then \(\operatorname{spreadOf}(x) = s^2 / n\). The squared standard error of one mean; used by the Welch test.
Confidence intervals
A t interval for a mean is the sample mean plus or minus a margin, where the margin is the critical value times the standard error:
\[\bar{x} \pm t_{\nu,\,1 - (1-c)/2} \cdot \frac{s}{\sqrt{n}}\]The critical value is the Student t quantile. For a proportion of \(k\) successes in \(n\) trials the module gives the Wald margin on the normal distribution, \(z_{1-(1-c)/2}\sqrt{\hat{p}(1-\hat{p})/n}\).
critical
critical(v, c) · v degrees of freedom scalar · c coverage, for example 0.95 scalar → critical value scalar
Let \(\nu > 0\) (the argument v) be a number of degrees of freedom and let \(0 < c < 1\) be the coverage. Let \(Q_{\nu}(q) = \operatorname{studentQuantile}(\nu, q)\) be the Student t quantile function on \(\nu\) degrees of freedom. Then
The two-sided Student t critical value for coverage c.
Intuition How many standard errors wide an interval must be to catch the true value with probability c. It is the t version of the familiar 1.96: larger for small samples, because the standard error is itself estimated and may be too small, and approaching 1.96 (at c = 0.95) as degrees of freedom grow.
Use it when you build any t-based interval yourself: a regression coefficient plus or minus critical(n − p, c) times its errors, a Welch interval, or a rejection threshold for a t statistic. For the mean of one sample, margin already combines it with the standard error. It comes from a Newton-inverted quantile, so it is only as accurate as studentQuantile.
confidence-intervalssmall-samples
def critical(v, c) ≔ studentQuantile(v, 1 − (1 − c) / 2)
margin
margin(x, c) · x sample vector n · c coverage scalar → half-width scalar
Let \(x \in \mathbb{R}^n\), \(n \ge 2\), be a sample with sample standard deviation \(s\) (divided by \(n - 1\)), and let \(0 < c < 1\). Then
\[\operatorname{margin}(x, c) = Q_{n-1}\!\left(1 - \frac{1 - c}{2}\right) \frac{s}{\sqrt{n}} ,\]where \(Q_{\nu}\) is the Student t quantile function, as in critical.
The half-width of the t interval for the mean of x: the critical value on \(n - 1\) degrees of freedom times the standard error.
Intuition The “plus or minus” in “5.47 ± 0.44”: how far the true mean could plausibly be from the sample mean. It shrinks like \(1/\sqrt{n}\), so halving it takes about four times the data.
Use it when you report a mean with its precision, compare precision across experiments, or plan a study (run a pilot, then pick n so the margin meets a target). It assumes independent observations and a roughly normal sample mean. For the interval’s endpoints use lowerBound and upperBound; for proportions use proportionMargin; for statistics other than the mean, or strongly skewed data, a bootstrap interval (bootLower in Resampling).
confidence-intervalsuncertainty-quantificationsample-size-planning
def margin(x, c) ≔ critical(freedomOf(x), c) · error(x)
lowerBound
lowerBound(x, c) · x sample vector n · c coverage scalar → lower limit scalar
Let \(x \in \mathbb{R}^n\), \(n \ge 2\), be a sample with mean \(\bar{x}\) and sample standard deviation \(s\), and let \(0 < c < 1\). Then
\[\operatorname{lowerBound}(x, c) = \bar{x} - Q_{n-1}\!\left(1 - \frac{1 - c}{2}\right) \frac{s}{\sqrt{n}} .\]The lower limit of the t confidence interval for the mean.
Intuition The smallest mean the data comfortably allow. If you repeated the experiment many times, a fraction c of intervals built this way would contain the true mean; the interval is the random thing, not the mean.
Use it when you report the interval’s lower end, or check a one-sided requirement such as “the average strength is at least 40”: with coverage \(c = 1 - 2\alpha\) this limit is a one-sided lower bound at level \(1 - \alpha\) (0.90 gives a 95% lower bound). An interval that excludes a value m agrees with pOneSample(x, m) falling below \(1 - c\). Pair it with upperBound.
confidence-intervalsuncertainty-quantification
def lowerBound(x, c) ≔ mean(x) − margin(x, c)
upperBound
upperBound(x, c) · x sample vector n · c coverage scalar → upper limit scalar
Let \(x \in \mathbb{R}^n\), \(n \ge 2\), be a sample with mean \(\bar{x}\) and sample standard deviation \(s\), and let \(0 < c < 1\). Then
\[\operatorname{upperBound}(x, c) = \bar{x} + Q_{n-1}\!\left(1 - \frac{1 - c}{2}\right) \frac{s}{\sqrt{n}} .\]The upper limit of the t confidence interval for the mean.
Intuition The largest mean the data comfortably allow, the mirror image of lowerBound about the sample mean. The two together bracket the values of m that a two-sided test at level \(1 - c\) would not reject.
Use it when you report the upper end of the interval, or check a ceiling such as “the average defect rate stays below the limit”: with coverage \(c = 1 - 2\alpha\) it is a one-sided upper bound at level \(1 - \alpha\). The interval is symmetric, so for a quantity bounded below (times, concentrations) with skewed data, a bootstrap interval (bootUpper in Resampling) may fit better.
confidence-intervalsuncertainty-quantification
def upperBound(x, c) ≔ mean(x) + margin(x, c)
proportionMargin
proportionMargin(k, n, c) · k number of successes scalar · n number of trials scalar · c coverage scalar → half-width scalar
Let \(k \ge 0\) be a number of successes in \(n > 0\) trials, let \(\hat{p} = k / n\), and let \(0 < c < 1\). Let \(\Phi^{-1} = \operatorname{standardQuantile}\) be the standard normal quantile function. Then
\[\operatorname{proportionMargin}(k, n, c) = \Phi^{-1}\!\left(1 - \frac{1 - c}{2}\right) \sqrt{\frac{\hat{p}\,(1 - \hat{p})}{n}} .\]The half-width of the Wald interval for a proportion. The interval itself is \(k/n\) plus or minus this margin; the result is not clamped to \([0, 1]\).
Intuition Each trial is a coin flip with success chance \(p\), whose variance is \(p(1-p)\); averaging n flips divides that by n. The margin is that standard error, estimated by plugging in \(\hat p\), times a normal critical value. It is widest at \(\hat p = 0.5\) and shrinks like \(1/\sqrt{n}\).
Use it when you report a rate from counts: a poll, a conversion rate, a defect fraction, a model’s accuracy on a test set. It is reliable only when both successes and failures are reasonably numerous (a common rule is at least about 10 of each). Near 0 or 1, or for small n, it undercovers and can extend below 0 or above 1; do not use it for rare events.
confidence-intervalsrate-estimationsample-size-planning
def proportionMargin(k, n, c) ≔ standardQuantile(1 − (1 − c) / 2) · sqrt(proportion(k, n) · (1 − proportion(k, n)) / n)
Building blocks
proportion(k, n)— Let \(k \in \mathbb{R}\) and \(n \neq 0\). Then \(\operatorname{proportion}(k, n) = k / n\). The observed proportion; used byproportionMargin.
The correlation test
A Pearson correlation \(r\) from \(n\) pairs becomes a t statistic on \(n - 2\) degrees of freedom, which tests whether the population correlation is zero:
\[t = r\,\sqrt{\frac{n - 2}{1 - r^2}}\]tCorrelation
tCorrelation(x, y) · x first variable vector n · y second variable vector n → t statistic scalar
Let \(x, y \in \mathbb{R}^n\) with \(n \ge 3\), and let
\[r = \frac{\sum_{i=1}^{n} (x_i - \bar{x})(y_i - \bar{y})}{(n - 1)\, s_x\, s_y}\]be their Pearson correlation, where \(s_x, s_y\) are the sample standard deviations; let \(\lvert r \rvert < 1\). Then
\[\operatorname{tCorrelation}(x, y) = r \sqrt{\frac{n - 2}{1 - r^2}} .\]The t statistic of the Pearson correlation between x and y.
Intuition A sample correlation is never exactly zero, even for unrelated variables. This rescales \(r\) so it can be read against a t distribution: the same \(r\) is weak evidence from 10 pairs and strong evidence from 1000, and the factor \(\sqrt{n-2}\) is what captures that.
Use it when you need the statistic itself, for a one-sided test with oneSided(n − 2, t) or to rank many variable pairs by evidence. It detects only linear association and assumes roughly normal variables; a single outlier can create or hide a correlation. It divides by \(1 - r^2\), so \(\lvert r \rvert = 1\) fails. For the two-sided p-value use pCorrelation.
correlation-analysishypothesis-testing
def tCorrelation(x, y) ≔ correlation(x, y) · sqrt((count(x) − 2) / (1 − correlation(x, y) · correlation(x, y)))
pCorrelation
pCorrelation(x, y) · x first variable vector n · y second variable vector n → p-value scalar
Let \(x, y \in \mathbb{R}^n\) with \(n \ge 3\), and let \(t = \operatorname{tCorrelation}(x, y)\). Then
\[\operatorname{pCorrelation}(x, y) = 2\,\bigl(1 - F_{n-2}(\lvert t \rvert)\bigr) .\]The two-sided p-value for the null hypothesis of zero correlation.
Intuition If x and y were truly unrelated, how often would a sample of this size show a correlation at least this strong, positive or negative? It answers “is there any linear relationship at all?”, not “is it strong?”.
Use it when you screen whether two measured variables move together, or confirm a relationship before modelling it. With large n, tiny correlations become significant, so look at \(r\) itself for strength. Screening many pairs needs bonferroni or discoveries. Correlation is not causation, and a strong non-linear relationship can still give a large p-value. It is equivalent to the slope test of a simple regression.
correlation-analysishypothesis-testing
def pCorrelation(x, y) ≔ twoSided(count(x) − 2, tCorrelation(x, y))
One-way analysis of variance
One-way ANOVA asks whether \(G\) group means differ. Every per-group quantity is a contraction against the membership matrix g, so there is no loop over groups. The total sum of squares splits exactly into a between-group part and a within-group part, and the F statistic compares their mean squares:
on \(G - 1\) and \(n - G\) degrees of freedom.
anovaF
anovaF(g, y) · g one-hot group membership matrix n×G · y observations vector n → F statistic scalar
Let \(y \in \mathbb{R}^n\) be the observations and let \(g \in \{0, 1\}^{n \times G}\) be a one-hot membership matrix, \(\sum_{k=1}^{G} g_{ik} = 1\) for every \(i\), with every group non-empty and \(1 < G < n\). Write \(n_k = \sum_{i} g_{ik}\), \(\bar{y}_k = \frac{1}{n_k} \sum_{i} g_{ik}\, y_i\) and \(\bar{y} = \frac{1}{n} \sum_{i} y_i\). Then
\[\operatorname{anovaF}(g, y) = \frac{\dfrac{1}{G - 1} \displaystyle\sum_{k=1}^{G} n_k\, (\bar{y}_k - \bar{y})^2}{\dfrac{1}{n - G} \displaystyle\sum_{i=1}^{n} \Bigl(y_i - \sum_{k=1}^{G} g_{ik}\, \bar{y}_k\Bigr)^2} .\]The one-way ANOVA F statistic. Every group should have at least one member; the number of groups is the number of columns of g.
Intuition Compare two estimates of the noise. One comes from how far the group means sit from the overall mean; the other from how far points sit from their own group mean. If the groups really share one mean, both estimate the same variance and \(F\) is near 1; real differences inflate only the first, pushing \(F\) up.
Use it when you compare the means of three or more independent groups in one test, instead of running every pairwise t-test and inflating false positives. It assumes roughly normal data with similar variance in every group. It says only that some means differ, not which; follow up with pairwise tests corrected by bonferroni. For the p-value use pAnova; for skewed data, kruskal in Resampling.
comparing-groupshypothesis-testing
def anovaF(g, y) ≔ (betweenSquares(g, y) / (groups(g) − 1)) / (withinSquares(g, y) / (count(y) − groups(g)))
pAnova
pAnova(g, y) · g one-hot group membership matrix n×G · y observations vector n → p-value scalar
Let \(g \in \{0, 1\}^{n \times G}\) and \(y \in \mathbb{R}^n\) be as in anovaF, and let \(f = \operatorname{anovaF}(g, y)\). Then
where \(\operatorname{fCdf}(d_1, d_2, \cdot)\) is the distribution function of the F distribution on \(d_1\) and \(d_2\) degrees of freedom.
The p-value of the F statistic: the upper tail of the F distribution on \(G - 1\) and \(n - G\) degrees of freedom.
Intuition How often groups that all share the same true mean would produce an F ratio this large. Only the upper tail counts, because group differences can only make the between-group estimate bigger.
Use it when you want a single decision on “do any of these groups differ?”, for example several treatments, suppliers or configurations. A small p-value licenses pairwise follow-up comparisons (with pWelch and a correction such as bonferroni or discoveries). With exactly two groups it matches pPooled. For data far from normal or with outliers, use kruskal in Resampling or a permutation test.
comparing-groupshypothesis-testingmultiple-comparisons
def pAnova(g, y) ≔ 1 − fCdf(groups(g) − 1, count(y) − groups(g), anovaF(g, y))
Building blocks
groupCounts(g)— Let \(g \in \{0, 1\}^{n \times G}\). Then \(\operatorname{groupCounts}(g)_k = n_k = \sum_{i=1}^{n} g_{ik}\) for \(k = 1, \dots, G\). The number of observations in each group, a vector of length G.groupSums(g, y)— Let \(g \in \{0, 1\}^{n \times G}\) and \(y \in \mathbb{R}^n\). Then \(\operatorname{groupSums}(g, y)_k = \sum_{i=1}^{n} g_{ik}\, y_i\) for \(k = 1, \dots, G\). The sum of the observations in each group.groupAverages(g, y)— Let \(g \in \{0, 1\}^{n \times G}\) with every \(n_k = \sum_i g_{ik} > 0\), and let \(y \in \mathbb{R}^n\). Then \(\operatorname{groupAverages}(g, y)_k = \bar{y}_k = \frac{1}{n_k} \sum_{i=1}^{n} g_{ik}\, y_i\). The mean of each group.groupFitted(g, y)— Let \(g \in \{0, 1\}^{n \times G}\) and \(y \in \mathbb{R}^n\), with group means \(\bar{y}_k\). Then \(\operatorname{groupFitted}(g, y)_i = \sum_{k=1}^{G} g_{ik}\, \bar{y}_k\) for \(i = 1, \dots, n\). Each observation’s own group mean, a vector of length n.betweenSquares(g, y)— Let \(g \in \{0, 1\}^{n \times G}\) and \(y \in \mathbb{R}^n\), with group sizes \(n_k\), group means \(\bar{y}_k\) and overall mean \(\bar{y} = \frac{1}{n} \sum_i y_i\). Then \(\operatorname{betweenSquares}(g, y) = \sum_{k=1}^{G} n_k\, (\bar{y}_k - \bar{y})^2\). The between-group sum of squares.withinSquares(g, y)— Let \(g \in \{0, 1\}^{n \times G}\) and \(y \in \mathbb{R}^n\), with group means \(\bar{y}_k\). Then \(\operatorname{withinSquares}(g, y) = \sum_{i=1}^{n} \bigl(y_i - \sum_{k} g_{ik}\, \bar{y}_k\bigr)^2\). The within-group (residual) sum of squares.groups(g)— Let \(g \in \mathbb{R}^{n \times G}\). Then \(\operatorname{groups}(g) = G\). The number of groups, the column count ofg.
Many tests at once
When \(m\) tests are run together, some p-values will look small by chance. The Bonferroni correction multiplies every p-value by \(m\) and caps it at 1, which controls the family-wise error rate. The Benjamini–Hochberg procedure controls the false discovery rate instead: sort the p-values, compare the \(k\)-th smallest (counting from 1) against \(k\,a/m\), find the largest one that passes, and declare it and every smaller p-value a discovery.
bonferroni
bonferroni(p) · p p-values vector m → adjusted p-values vector m
Let \(p = (p_1, \dots, p_m) \in [0, 1]^m\) be p-values. Then
\[\operatorname{bonferroni}(p)_i = \min(m\, p_i,\; 1), \qquad i = 1, \dots, m .\]Bonferroni-adjusted p-values, in the same order as the input. Compare them with the usual significance level.
Intuition Run 20 tests at 0.05 on pure noise and you expect one “significant” result. Multiplying each p-value by the number of tests m is the same as testing each at \(a/m\), which keeps the chance of even one false positive across the whole family at or below \(a\), whatever the dependence between the tests.
Use it when a single false positive is costly and the number of tests is modest: a handful of pairwise comparisons after pAnova, several endpoints in a trial, safety checks. It becomes very conservative as m grows and misses real effects; with many tests (screening genes, features, metrics) use discoveries.
multiple-comparisonshypothesis-testing
def bonferroni(p) ≔ min(p · count(p), 1)
discoveries
discoveries(p, a) · p p-values vector m · a false discovery rate, for example 0.05 scalar → 0/1 mask vector m
Let \(p = (p_1, \dots, p_m) \in [0, 1]^m\) be p-values and let \(0 < a < 1\) be the false discovery rate. Let \(p_{(1)} \le p_{(2)} \le \dots \le p_{(m)}\) be the same values in ascending order. Then
\[\operatorname{discoveries}(p, a)_i = [\,p_i \le c\,]\,[\,c > 0\,], \qquad i = 1, \dots, m,\]where
\[c = \max_{1 \le j \le m} \; p_{(j)} \, \Bigl[\, p_{(j)} \le \frac{j\, a}{m} \,\Bigr]\]is the largest sorted p-value that passes its own threshold.
The Benjamini–Hochberg decision for each test, in the input order: 1 where the null hypothesis is rejected, 0 where it is not. The p-values need not be sorted. When no p-value passes, every entry is 0. Because the cutoff must be strictly positive, a p-value of exactly 0 does not count on its own: if 0 is the only value that passes, every entry is 0, where the textbook procedure would reject that one test.
Intuition Instead of guarding against any false positive, it controls the expected fraction of false ones among everything you declare a finding. The threshold grows with rank, \(j\,a/m\), so once several p-values are small the bar relaxes for the next ones: many real effects make each other easier to detect.
Use it when you screen many hypotheses and can tolerate a known share of false leads: feature or gene screening, many A/B metrics, anomaly flags. With a = 0.05, about 5% of the discoveries are expected to be false. It is valid for independent or positively correlated tests. When any single false positive is unacceptable, use bonferroni. Clamp p-values of exactly 0 to a tiny positive number first (see above).
multiple-comparisonsfeature-selectionhypothesis-testing
def discoveries(p, a) ≔ (p ≤ cutoff(p, a)) · (cutoff(p, a) > 0)
Building blocks
steps(p)— Let \(p \in \mathbb{R}^m\). Then \(\operatorname{steps}(p)_j = j / m\) for \(j = 1, \dots, m\). The Benjamini–Hochberg fractions, built from the positions \(0, 1, \dots, m - 1\) plus one.passing(p, a)— Let \(p \in \mathbb{R}^m\) with ascending order \(p_{(1)} \le \dots \le p_{(m)}\), and let \(a \in \mathbb{R}\). Then \(\operatorname{passing}(p, a)_j = [\,p_{(j)} \le j\, a / m\,]\) for \(j = 1, \dots, m\). A 0/1 mask over the sorted p-values: which fall at or below their own threshold.cutoff(p, a)— Let \(p\) and \(a\) be as inpassing. Then \(\operatorname{cutoff}(p, a) = \max_{1 \le j \le m} p_{(j)}\, [\,p_{(j)} \le j\, a / m\,]\). The largest sorted p-value that passes, or 0 when none does.
Examples
Test a mean and give its 95% interval
Six measurements; is their mean different from 5.0?
open Prelude.Inference
fn studentT(x, m) ↦ tOneSample(x, m)
fn studentP(x, m) ↦ pOneSample(x, m)
fn meanLow(x, c) ↦ lowerBound(x, c)
fn meanHigh(x, c) ↦ upperBound(x, c)
using Examples;
var x = new[] { 5.1d, 4.9, 5.6, 5.8, 6.0, 5.4 };
var t = Tests.StudentT(x, 5.0);
var p = Tests.StudentP(x, 5.0);
var low = Tests.MeanLow(x, 0.95);
var high = Tests.MeanHigh(x, 0.95);
Output (verified, rounded to 6 decimal places):
t = 2.735126
p = 0.041029
low = 5.028075
high = 5.905259
The mean is 5.4667 and the standard error 0.1706, so \(t = 0.4667 / 0.1706\); the interval uses the critical value 2.571 on 5 degrees of freedom.
Compare two independent samples
open Prelude.Inference
fn welchP(x, y) ↦ pWelch(x, y)
fn welchDegrees(x, y) ↦ freedomWelch(x, y)
fn pooledP(x, y) ↦ pPooled(x, y)
using Examples;
var a = new[] { 1d, 2, 3, 4, 5 };
var b = new[] { 3d, 4, 5, 6, 7, 8, 9 };
var welch = Tests.WelchP(a, b);
var degrees = Tests.WelchDegrees(a, b);
var pooled = Tests.PooledP(a, b);
Output (verified, rounded to 6 decimal places):
welch p = 0.019594
welch v = 9.966102
pooled p = 0.025236
The variances are 2.5 and 4.667, so the Welch degrees of freedom fall between 4 and 10 rather than being the pooled 10.
One-way analysis of variance on three groups
Three groups of three observations, with group means 2, 5 and 8.
open Prelude.Inference
fn anovaStatistic(g, y) ↦ anovaF(g, y)
fn anovaP(g, y) ↦ pAnova(g, y)
using Examples;
var g = new double[,] { {1,0,0},{1,0,0},{1,0,0},{0,1,0},{0,1,0},{0,1,0},{0,0,1},{0,0,1},{0,0,1} };
var y = new[] { 1d, 2, 3, 4, 5, 6, 7, 8, 9 };
var f = Tests.AnovaStatistic(g, y);
var p = Tests.AnovaP(g, y);
Output (verified):
F = 27
pF = 0.001
The between-group sum of squares is 54 on 2 degrees of freedom and the within-group sum is 6 on 6, so \(F = 27 / 1\). For two numerator degrees of freedom the tail has the closed form \((1 + 2F/6)^{-3} = 10^{-3}\).
Correct four p-values for multiple testing
open Prelude.Inference
fn adjusted(p) ↦ bonferroni(p)
fn findings(p, a) ↦ discoveries(p, a)
using Examples;
var p = new[] { 0.01d, 0.04, 0.03, 0.20 };
var bonf = Tests.Adjusted(p);
var bh = Tests.Findings(p, 0.05);
Output (verified):
bonf = [0.04, 0.16, 0.12, 0.8]
bh = [1, 0, 0, 0]
Sorted, the p-values are 0.01, 0.03, 0.04, 0.20 against thresholds 0.0125, 0.025, 0.0375, 0.05. Only 0.01 passes, so only the first test is a discovery.
A margin of error for a poll
40 of 100 respondents said yes.
open Prelude.Inference
fn waldMargin(k, n, c) ↦ proportionMargin(k, n, c)
using Examples;
var margin = Tests.WaldMargin(40.0, 100.0, 0.95);
Output (verified, rounded to 6 decimal places):
margin = 0.096018
That is \(1.96 \sqrt{0.4 \cdot 0.6 / 100}\), so the interval is roughly 0.304 to 0.496.