Resampling
The bootstrap, permutation tests and the rank tests (Mann–Whitney, Kruskal–Wallis, Wilcoxon signed rank), computed without loops over replicates.
This module gives distribution-free inference: the bootstrap for standard errors, bias and percentile intervals; permutation tests for a difference in means; and the classical rank tests. None of them assumes a particular distribution for the data.
Randomness is supplied by the caller. A bootstrap takes an index matrix i of shape B×n: each of its B rows is one replicate and holds n whole numbers drawn uniformly from \(\{0, \dots, n-1\}\) with replacement. A permutation test takes a matrix p of the same shape whose rows are shuffles of \(0, \dots, n-1\). Both are passed as ordinary numbers (a double[,] from C#). Every replicate is computed at once, as one tensor of shape B×n (and an intermediate B×n×n comparison), so memory grows with \(B\,n^2\).
The rank tests use one-based midranks, the convention of printed tables: tied values share the average of the positions they occupy. Groups are given as a one-hot membership matrix g of shape n×G over the pooled sample. The statistics are returned as raw values (U, H, W); the module does not convert them to p-values, and the Kruskal–Wallis statistic has no correction for ties.
open Prelude.Resampling
Gathering
A resample is a gather: drawn(i, x) replaces every index in i with the value of x at that position. It is written as a comparison of each index against every position, followed by one contraction, so all replicates are gathered together.
drawn
drawn(i, x) · i indices matrix B×n · x data vector n → resampled values matrix B×n
Let \(x = (x_0, \dots, x_{n-1}) \in \mathbb{R}^n\) be the data, indexed by positions \(0, 1, \dots, n-1\), and let \(i \in \{0, \dots, n-1\}^{B \times m}\) be an index matrix (usually \(m = n\)). Then
\[\operatorname{drawn}(i, x)_{bk} = \sum_{j=0}^{n-1} [\, i_{bk} = j \,]\, x_j = x_{i_{bk}}, \qquad b = 1, \dots, B,\; k = 1, \dots, m .\]The data gathered at the given indices: row b is replicate b. The indices must be whole numbers in \([0, n)\); an index outside that range gathers 0.
Intuition The bootstrap treats your sample as a stand-in for the population: drawing n values from it with replacement imitates collecting a fresh sample. Each row of i is one such imaginary sample, and this function looks the values up for all of them in a single step.
Use it when you want the resampled data themselves, to apply a statistic this module does not provide (a median, a trimmed mean, a ratio) to every replicate, or to relabel a pooled sample with permutation rows. Draw the indices in C# from a seeded generator so the result is reproducible. The intermediate one-hot tensor is B×n×n, so a few thousand replicates of a few thousand points is already large.
simulationuncertainty-quantificationnon-normal-data
def drawn(i, x) ≔ einsum("bkj,j->bk", picked3(i, x), x)
Building blocks
drawnAt(i, x)— Let \(i \in \mathbb{R}^{B \times m}\) and \(x \in \mathbb{R}^n\). Then \(\operatorname{drawnAt}(i, x)_{bkj} = i_{bk}\) for \(j = 0, \dots, n-1\). The index matrix lifted to B×m×n, repeating each drawn index along the positions.placeAt(i, x)— Let \(i \in \mathbb{R}^{B \times m}\) and \(x \in \mathbb{R}^n\). Then \(\operatorname{placeAt}(i, x)_{bkj} = j\) for \(j = 0, \dots, n-1\). The positions lifted to the same B×m×n shape.picked3(i, x)— Let \(i \in \mathbb{R}^{B \times m}\) and \(x \in \mathbb{R}^n\). Then \(\operatorname{picked3}(i, x)_{bkj} = [\, i_{bk} = j \,]\). The one-hot B×m×n comparison of drawn index against position.
The bootstrap
Each row of the gathered matrix is a bootstrap sample. The replicate means estimate the sampling distribution of the mean: their standard deviation is the bootstrap standard error, their mean minus the sample mean is the bias, and their quantiles give the percentile interval
\[\left[\, Q_{a/2}(\bar{x}^*),\; Q_{1-a/2}(\bar{x}^*) \,\right]\]where a is the significance level (0.05 for a 95% interval). The quantile interpolates linearly between order statistics, as quantile does in Prelude.Statistics.
drawnMeans
drawnMeans(i, x) · i indices matrix B×n · x data vector n → replicate means vector B
Let \(x \in \mathbb{R}^n\) and let \(i \in \{0, \dots, n-1\}^{B \times m}\) be an index matrix, and write \(\tilde{x} = \operatorname{drawn}(i, x) \in \mathbb{R}^{B \times m}\). Then
\[\operatorname{drawnMeans}(i, x)_b = \bar{x}^*_b = \frac{1}{m} \sum_{k=1}^{m} \tilde{x}_{bk}, \qquad b = 1, \dots, B .\]The mean of each bootstrap replicate.
Intuition The mean you would have got in each imaginary repeat of the study. Their histogram approximates the sampling distribution of the mean, read straight off the data instead of assumed to be normal.
Use it when you want to see that distribution (plot it to check skew) or feed it to your own summaries. bootError, bootBias, bootLower and bootUpper are all summaries of this vector, so for those use the functions directly. Hundreds of replicates are enough for a standard error; percentile intervals need a few thousand for stable tails.
uncertainty-quantificationsimulation
def drawnMeans(i, x) ≔ einsum("bk->b", drawn(i, x)) / columns(i)
drawnSpreads
drawnSpreads(i, x) · i indices matrix B×n · x data vector n → replicate variances vector B
Let \(x \in \mathbb{R}^n\), let \(i \in \{0, \dots, n-1\}^{B \times m}\) be an index matrix with \(B = m\) (otherwise the shapes do not broadcast), and write \(\tilde{x} = \operatorname{drawn}(i, x)\) and \(\bar{x}^* = \operatorname{drawnMeans}(i, x)\). Then the code computes
\[\operatorname{drawnSpreads}(i, x)_b = \frac{1}{m - 1} \sum_{k=1}^{m} \bigl(\tilde{x}_{bk} - \bar{x}^*_k\bigr)^2 ,\]with the mean of replicate \(k\), not of replicate \(b\), subtracted from entry \(k\). The sample variance of replicate \(b\) would subtract \(\bar{x}^*_b\) throughout.
Intended as the sample variance of each bootstrap replicate. Caution: as written, the replicate means (length B) are subtracted along the last axis of the B×n matrix. The call therefore fails with a shape error when B differs from n, and returns incorrect values when B equals n. Do not rely on it; the other bootstrap functions are unaffected.
Intuition The intent is the spread within each imaginary resample, whose distribution would show how uncertain the sample variance is, just as drawnMeans does for the mean.
Use it when it has been fixed. As written, do not use it: per the caution above, it fails when B differs from n and is wrong when they are equal. If you need bootstrap variances, take drawn(i, x) and compute each row’s sample variance yourself. For uncertainty about the mean, bootError and the percentile interval work correctly.
uncertainty-quantificationsimulation
def drawnSpreads(i, x) ≔ einsum("bk,bk->b", drawn(i, x) − drawnMeans(i, x), drawn(i, x) − drawnMeans(i, x)) / (columns(i) − 1)
bootLower
bootLower(i, x, a) · i indices matrix B×n · x data vector n · a significance level scalar → lower limit scalar
Let \(x \in \mathbb{R}^n\), let \(i \in \{0, \dots, n-1\}^{B \times m}\) be an index matrix with replicate means \(\bar{x}^* = \operatorname{drawnMeans}(i, x) \in \mathbb{R}^B\), and let \(0 < a < 1\) be the significance level. Then
\[\operatorname{bootLower}(i, x, a) = Q_{a/2}(\bar{x}^*) ,\]where, for \(v \in \mathbb{R}^B\) with order statistics \(v_{(0)} \le \dots \le v_{(B-1)}\) (positions counted from zero), \(h = q\,(B - 1)\) and \(\phi = h - \lfloor h \rfloor\),
\[Q_q(v) = (1 - \phi)\, v_{(\lfloor h \rfloor)} + \phi\, v_{(\lceil h \rceil)} .\]The lower limit of the bootstrap percentile interval for the mean: the \(a/2\) quantile of the replicate means. Note that a is the significance level, not the coverage.
Intuition Sort the replicate means and cut off the lowest \(a/2\) fraction: the interval’s lower end is where the imaginary repeats start to become common. No normal shape is assumed, so a skewed sampling distribution gives a lopsided interval, as it should.
Use it when you want a confidence interval for a mean without trusting the t formula, typically for skewed data such as incomes, latencies or costs. Pass a = 0.05 for 95%, not 0.95. Use at least a few thousand replicates, since the tail quantile rests on few of them. With very small samples percentile intervals tend to be too narrow; the t interval lowerBound in Inference may be safer if the data are not strongly skewed.
confidence-intervalsnon-normal-data
def bootLower(i, x, a) ≔ quantile(drawnMeans(i, x), a / 2)
bootUpper
bootUpper(i, x, a) · i indices matrix B×n · x data vector n · a significance level scalar → upper limit scalar
Let \(x \in \mathbb{R}^n\), let \(i \in \{0, \dots, n-1\}^{B \times m}\) be an index matrix with replicate means \(\bar{x}^* = \operatorname{drawnMeans}(i, x)\), let \(0 < a < 1\) be the significance level, and let \(Q_q\) be the interpolated quantile of bootLower. Then
The upper limit of the percentile interval: the \(1 - a/2\) quantile of the replicate means.
Intuition The mirror of bootLower: the value above which only an \(a/2\) fraction of the replicate means fall. With right-skewed data it usually sits farther from the sample mean than the lower limit does.
Use it when you report the upper end of a bootstrap interval, or a one-sided upper bound for a mean (a call with a gives a bound at level \(1 - a/2\)). Always pass the same index matrix i and the same a as to bootLower, so both limits come from one set of replicates.
confidence-intervalsnon-normal-data
def bootUpper(i, x, a) ≔ quantile(drawnMeans(i, x), 1 − a / 2)
bootError
bootError(i, x) · i indices matrix B×n · x data vector n → standard error scalar
Let \(x \in \mathbb{R}^n\) and let \(i \in \{0, \dots, n-1\}^{B \times m}\) be an index matrix, \(B \ge 2\), with replicate means \(\bar{x}^*_b = \operatorname{drawnMeans}(i, x)_b\). Then
\[\operatorname{bootError}(i, x) = \sqrt{\frac{1}{B - 1} \sum_{b=1}^{B} \bigl(\bar{x}^*_b - \bar{\bar{x}}^*\bigr)^2}, \qquad \bar{\bar{x}}^* = \frac{1}{B} \sum_{b=1}^{B} \bar{x}^*_b .\]The bootstrap standard error of the mean: the sample standard deviation of the replicate means.
Intuition How much the mean moves between imaginary repeats of the study. For the mean it should come out close to the textbook \(s/\sqrt{n}\) (slightly smaller, by about \(\sqrt{(n-1)/n}\)), which makes it a good way to see and trust what the bootstrap does.
Use it when you want a standard error with no formula, or check the textbook one on awkward data. Its real value is for statistics with no simple formula, which you would compute from drawn yourself; this function covers the mean only. A few hundred replicates usually suffice. For an interval rather than a single spread, use bootLower and bootUpper.
uncertainty-quantificationnon-normal-data
def bootError(i, x) ≔ std(drawnMeans(i, x))
bootBias
bootBias(i, x) · i indices matrix B×n · x data vector n → bias scalar
Let \(x \in \mathbb{R}^n\) and let \(i \in \{0, \dots, n-1\}^{B \times m}\) be an index matrix with replicate means \(\bar{x}^*_b\). Then
\[\operatorname{bootBias}(i, x) = \frac{1}{B} \sum_{b=1}^{B} \bar{x}^*_b - \frac{1}{n} \sum_{j=0}^{n-1} x_j .\]The bootstrap estimate of bias: the mean of the replicate means minus the sample mean.
Intuition If the estimate computed on resamples systematically lands above or below the estimate on the original sample, the same offset probably separates the original estimate from the truth. The sample mean is unbiased, so here the answer should be near zero, and what remains is Monte Carlo noise.
Use it when you check that the replicates are well behaved: a value that is large relative to bootError means too few replicates or badly drawn indices. The method matters for biased statistics (ratios, correlations, variances), which you would bootstrap from drawn yourself; subtracting a bias estimate adds variance, so correct only when the bias is large relative to the standard error.
bias-correctionsanity-checking
def bootBias(i, x) ≔ mean(drawnMeans(i, x)) − mean(x)
Building blocks
replicates(i)— Let \(i \in \mathbb{R}^{B \times m}\). Then \(\operatorname{replicates}(i) = B\). The number of replicates, the row count of the index matrix.
Permutation
A permutation test uses the same gathering with shuffles instead of draws. Put the two samples one after the other in a single pooled vector x, the first sample occupying the first k positions. Each row of p relabels the pooled data; the difference between the mean of the first k values and the mean of the rest, over all rows, is the null distribution of the difference in means. The p-value counts the observed statistic among the shuffles:
so it can never be exactly zero. The test is two-sided.
permDifference
permDifference(p, x, k) · p permutations matrix B×n · x pooled data vector n · k size of the first sample scalar → differences vector B
Let \(x \in \mathbb{R}^n\) be the pooled data, indexed by positions \(0, \dots, n-1\), let \(p \in \{0, \dots, n-1\}^{B \times n}\) be a matrix whose rows are permutations, let \(0 < k < n\), and write \(\tilde{x} = \operatorname{drawn}(p, x)\). Then
\[\operatorname{permDifference}(p, x, k)_b = \delta_b = \frac{1}{k} \sum_{j=0}^{k-1} \tilde{x}_{bj} - \frac{1}{n - k} \sum_{j=k}^{n-1} \tilde{x}_{bj}, \qquad b = 1, \dots, B .\]For every shuffle, the mean of the first k relabelled values minus the mean of the remaining \(n - k\).
Intuition If the group labels made no difference, any reshuffling of them would be as plausible as the real one. Shuffling many times and recomputing the gap shows what differences pure chance produces, building the null distribution from the data themselves.
Use it when you want to see that null distribution, plot it next to the observed difference, or summarise it in your own way. For the p-value itself use permP. The two samples must be laid out first-k-then-rest in x, and the permutation rows must be genuine shuffles of \(0, \dots, n-1\).
comparing-groupshypothesis-testingnon-normal-data
def permDifference(p, x, k) ≔ headMean(p, x, k) − tailMean(p, x, k)
permP
permP(p, x, k, d) · p permutations matrix B×n · x pooled data vector n · k size of the first sample scalar · d observed difference scalar → p-value scalar
Let \(x \in \mathbb{R}^n\), \(p \in \{0, \dots, n-1\}^{B \times n}\) and \(0 < k < n\) be as in permDifference, with shuffled differences \(\delta_b\), and let \(d \in \mathbb{R}\) be the observed difference. Then
The two-sided permutation p-value for the observed difference d, which you compute from the unshuffled data (the mean of the first k values minus the mean of the rest).
Intuition The fraction of shuffles whose mean gap is at least as large as the real one. The observed labelling counts as one of the shuffles, which is why 1 is added on top and bottom: a p-value is never claimed to be exactly zero from finitely many shuffles.
Use it when you compare two groups’ means without assuming normality, particularly for small or oddly distributed samples, A/B tests with heavy tails, or when a t-test’s assumptions are doubtful. It assumes the observations are exchangeable under “no difference”, which means the same distribution in both groups, so very unequal spreads can distort it. The smallest possible p-value is \(1/(B+1)\), so use at least 1000 shuffles for a test at 0.05 and many more for small thresholds. For a rank-based alternative, mannWhitney; for the t version, pWelch in Inference.
comparing-groupshypothesis-testingnon-normal-datasmall-samples
def permP(p, x, k, d) ≔ (permCount(p, x, k, d) + 1) / (replicates(p) + 1)
Building blocks
firstOf(x, k)— Let \(x \in \mathbb{R}^n\) and \(k \in \mathbb{R}\). Then \(\operatorname{firstOf}(x, k)_j = [\, j < k \,]\) for positions \(j = 0, \dots, n-1\). A 0/1 mask over positions, 1 for the first k.headMean(p, x, k)— Let \(p \in \mathbb{R}^{B \times n}\), \(x \in \mathbb{R}^n\), \(k > 0\) and \(\tilde{x} = \operatorname{drawn}(p, x)\). Then \(\operatorname{headMean}(p, x, k)_b = \frac{1}{k} \sum_{j=0}^{n-1} \tilde{x}_{bj}\, [\, j < k \,]\). The mean of the first k relabelled values in each shuffle.tailMean(p, x, k)— Let \(p \in \mathbb{R}^{B \times n}\), \(x \in \mathbb{R}^n\), \(k < n\) and \(\tilde{x} = \operatorname{drawn}(p, x)\). Then \(\operatorname{tailMean}(p, x, k)_b = \frac{1}{n - k} \sum_{j=0}^{n-1} \tilde{x}_{bj}\, (1 - [\, j < k \,])\). The mean of the remaining relabelled values in each shuffle.permCount(p, x, k, d)— Let \(p, x, k\) be as inpermDifference, with differences \(\delta_b\), and let \(d \in \mathbb{R}\). Then \(\operatorname{permCount}(p, x, k, d) = \sum_{b=1}^{B} [\, \lvert \delta_b \rvert \ge \lvert d \rvert \,]\). How many shuffles give a difference at least as large in absolute value as d.
The rank tests
All three tests rank the pooled data once and add the ranks up. With \(R_g\) the rank sum and \(n_g\) the size of group \(g\), and \(N\) the pooled size:
\[U_g = R_g - \frac{n_g (n_g + 1)}{2}, \qquad H = \frac{12}{N(N+1)} \sum_g \frac{R_g^2}{n_g} - 3(N+1)\]For two groups the two U statistics add up to \(n_1 n_2\). The Wilcoxon signed-rank statistic \(W^+\) ranks the absolute paired differences and sums the ranks of the positive ones.
mannWhitney
mannWhitney(g, x) · g one-hot membership matrix n×2 · x pooled data vector n → U statistics vector 2
Let \(x \in \mathbb{R}^n\) be the pooled data and let \(g \in \{0, 1\}^{n \times 2}\) be a one-hot membership matrix. Let \(r_i = \tfrac{1}{2}\bigl(\#\{j : x_j < x_i\} + \#\{j : x_j \le x_i\} + 1\bigr)\) be the one-based midrank of \(x_i\), \(n_k = \sum_i g_{ik}\) and \(R_k = \sum_i g_{ik}\, r_i\). Then
\[\operatorname{mannWhitney}(g, x)_k = R_k - \frac{n_k (n_k + 1)}{2}, \qquad k = 1, 2 .\]The Mann–Whitney U statistic for both groups at once. Entry g is the number of pairs (one observation from group g, one from the other) in which group g’s observation is larger, with ties counting one half.
Intuition Pair every observation from one group with every observation from the other and count who wins. Only the ordering matters, not the size of the gap, so one enormous outlier counts as a single win, not as a huge pull on a mean.
Use it when you compare two independent groups whose data are skewed, have outliers, or are only ordinal (ratings, ranks, scores): “does group A tend to produce larger values than B?”. The module gives U, not a p-value. For a large sample, with no ties, U is approximately normal with mean \(n_1 n_2 / 2\) and variance \(n_1 n_2 (n + 1)/12\), so a z-score gives the p-value. Divide by mannWhitneyTotal for an effect size. For more than two groups use kruskal.
comparing-groupshypothesis-testingnon-normal-data
def mannWhitney(g, x) ≔ rankSums(g, x) − memberCounts(g) · (memberCounts(g) + 1) / 2
mannWhitneyTotal
mannWhitneyTotal(g) · g one-hot membership matrix n×2 → total scalar
Let \(g \in \{0, 1\}^{n \times G}\) be a one-hot membership matrix with group sizes \(n_k = \sum_i g_{ik}\). Then
\[\operatorname{mannWhitneyTotal}(g) = \frac{1}{2} \Bigl(\sum_{k=1}^{G} n_k\Bigr)^2 - \frac{1}{2} \sum_{k=1}^{G} n_k^2 = \sum_{k < l} n_k\, n_l ,\]which for two groups is \(n_1 n_2\).
The sum of the two U statistics, \(n_1 n_2\): the number of cross-group pairs. Dividing a U statistic by it gives the probability of superiority, which is the area under the ROC curve (see auc in Prelude.Metrics).
Intuition The number of head-to-head contests between the two groups. Every contest has one winner (or a tie split in half), so the two U statistics must share this total, and a U of half of it means neither group tends to win.
Use it when you turn a U statistic into an effect size: U / mannWhitneyTotal(g) is the probability that a random member of one group beats a random member of the other, 0.5 meaning no effect. It also serves as a check that the ranking was shared: the two entries of mannWhitney should add up to it. It depends only on group sizes, not on the data.
effect-sizesanity-checking
def mannWhitneyTotal(g) ≔ einsum("g->", memberCounts(g)) · einsum("g->", memberCounts(g)) / 2 − einsum("g,g->", memberCounts(g), memberCounts(g)) / 2
kruskal
kruskal(g, x) · g one-hot membership matrix n×G · x pooled data vector n → H statistic scalar
Let \(x \in \mathbb{R}^n\) be the pooled data and let \(g \in \{0, 1\}^{n \times G}\) be a one-hot membership matrix with every group non-empty, and let \(r_i\), \(n_k\) and \(R_k = \sum_i g_{ik}\, r_i\) be the one-based midranks, group sizes and rank sums of mannWhitney. Then
The Kruskal–Wallis H statistic, the rank analogue of one-way ANOVA for any number of groups. No tie correction is applied. For moderate group sizes H is approximately chi-squared on \(G - 1\) degrees of freedom.
Intuition Rank everything together; if the groups come from the same distribution, each group’s average rank should sit near the overall average \((n+1)/2\). H measures how far the group average ranks stray from it, weighted by group size, just as ANOVA does with means.
Use it when you compare three or more independent groups whose data are skewed, contain outliers or are ordinal, where anovaF in Inference would be unreliable. Compare H with a chi-square on \(G - 1\) degrees of freedom when every group has at least about five members. Without a tie correction, heavily tied data give an H that is too small, so the test is conservative. A significant H says some group differs; use pairwise mannWhitney with a correction to find which.
comparing-groupshypothesis-testingnon-normal-data
def kruskal(g, x) ≔ 12 / (count(x) · (count(x) + 1)) · einsum("g,g->", rankSums(g, x) / memberCounts(g), rankSums(g, x)) − 3 · (count(x) + 1)
signedRanks
signedRanks(d) · d paired differences vector n → signed ranks vector n
Let \(d \in \mathbb{R}^n\) be paired differences and let \(r_i = \tfrac{1}{2}\bigl(\#\{j : \lvert d_j \rvert < \lvert d_i \rvert\} + \#\{j : \lvert d_j \rvert \le \lvert d_i \rvert\} + 1\bigr)\) be the one-based midrank of \(\lvert d_i \rvert\). Then
\[\operatorname{signedRanks}(d)_i = r_i\, [\, d_i > 0 \,], \qquad i = 1, \dots, n .\]The rank of each \(\lvert d_i \rvert\) among all the absolute differences, kept where \(d_i > 0\) and 0 elsewhere. The ranking is done before masking, so negative differences still occupy their ranks. Zero differences are ranked too rather than dropped.
Intuition Rank how big each change was, ignoring its direction, then keep the credit only for the changes that went up. Large increases earn large ranks; if there is no real effect, ups and downs should hold ranks of similar total size.
Use it when you want the per-pair contributions to wilcoxon, for example to see which pairs drive the statistic. The ranks of the negative differences come from signedRanks(−d). Remove zero differences before calling, since they occupy ranks here.
before-after-comparisonnon-normal-data
def signedRanks(d) ≔ ranking(abs(d)) · (d > 0)
wilcoxon
wilcoxon(d) · d paired differences vector n → W⁺ scalar
Let \(d \in \mathbb{R}^n\) be paired differences with midranks \(r_i\) of \(\lvert d_i \rvert\), as in signedRanks. Then
The Wilcoxon signed-rank statistic: the sum of the ranks of the positive differences. Because zero differences are not removed first, remove them yourself before calling if your data contain any.
Intuition Under “no change”, each difference is as likely to be positive as negative, so the positive ones should collect about half of all ranks, \(n(n+1)/4\). A \(W^+\) far above that says the changes lean upward, and one far below says downward. Using ranks rather than values keeps one wild pair from dominating.
Use it when you have paired measurements (before and after, two methods on the same items) and the differences are skewed or have outliers, so pPaired in Inference is doubtful. It assumes the differences are roughly symmetric about their centre. The module returns \(W^+\), not a p-value; for about 20 or more pairs a normal approximation with variance \(n(n+1)(2n+1)/24\) (no ties) is standard.
before-after-comparisonhypothesis-testingnon-normal-data
def wilcoxon(d) ≔ einsum("i->", signedRanks(d))
Building blocks
ranking(x)— Let \(x \in \mathbb{R}^n\). Then \(\operatorname{ranking}(x)_i = \tfrac{1}{2}\bigl(\#\{j : x_j < x_i\} + \#\{j : x_j \le x_i\} + 1\bigr)\). One-based midranks: the zero-basedrankfromPrelude.Statisticsplus one.memberCounts(g)— Let \(g \in \{0, 1\}^{n \times G}\). Then \(\operatorname{memberCounts}(g)_k = n_k = \sum_{i=1}^{n} g_{ik}\). The size of each group.rankSums(g, x)— Let \(g \in \{0, 1\}^{n \times G}\) and \(x \in \mathbb{R}^n\) with one-based midranks \(r_i\). Then \(\operatorname{rankSums}(g, x)_k = R_k = \sum_{i=1}^{n} g_{ik}\, r_i\). The sum of the pooled ranks within each group.
Examples
Bootstrap a mean
Five replicates of a four-value sample. In practice B would be in the hundreds or thousands and the indices would come from a random generator; here they are fixed so the answer can be checked by hand.
open Prelude.Resampling
fn resampled(i, x) ↦ drawn(i, x)
fn replicateMeans(i, x) ↦ drawnMeans(i, x)
fn intervalLow(i, x, a) ↦ bootLower(i, x, a)
fn intervalHigh(i, x, a) ↦ bootUpper(i, x, a)
fn standardError(i, x) ↦ bootError(i, x)
fn biasOf(i, x) ↦ bootBias(i, x)
using Examples;
var x = new[] { 2d, 4, 6, 8 };
var i = new double[,] { { 0, 0, 1, 3 }, { 3, 2, 2, 1 }, { 1, 1, 1, 1 }, { 0, 1, 2, 3 }, { 3, 3, 2, 0 } };
var drawn = Resample.Resampled(i, x);
var means = Resample.ReplicateMeans(i, x);
var low = Resample.IntervalLow(i, x, 0.2);
var high = Resample.IntervalHigh(i, x, 0.2);
var se = Resample.StandardError(i, x);
var bias = Resample.BiasOf(i, x);
Output (verified):
drawn = shape [5, 4] values [2, 2, 4, 8, 8, 6, 6, 4, 4, 4, 4, 4, 2, 4, 6, 8, 8, 8, 6, 2]
means = [4, 6, 4, 5, 6]
low = 4
high = 6
se = 1
bias = 0
The first replicate draws positions 0, 0, 1, 3, that is 2, 2, 4, 8, with mean 4. The 80% interval is the 0.1 and 0.9 quantiles of the five replicate means.
A permutation test for a difference in means
The first three values (7, 9, 8) are one group and the last three (1, 2, 3) the other; the observed difference is \(8 - 2 = 6\). Four shuffles are given, the first of them the identity.
open Prelude.Resampling
fn shuffledGap(p, x, k) ↦ permDifference(p, x, k)
fn shuffleP(p, x, k, d) ↦ permP(p, x, k, d)
using Examples;
var pooled = new[] { 7d, 9, 8, 1, 2, 3 };
var p = new double[,] { { 0,1,2,3,4,5 }, { 3,4,5,0,1,2 }, { 0,3,1,4,2,5 }, { 5,4,3,2,1,0 } };
var gaps = Resample.ShuffledGap(p, pooled, 3.0);
var pValue = Resample.ShuffleP(p, pooled, 3.0, 6.0);
Output (verified, rounded to 6 decimal places):
gaps = [6, -6, 1.333333, -6]
permP = 0.8
Three of the four shuffles reach \(\lvert d \rvert \ge 6\), so \(p = (3 + 1)/(4 + 1)\).
Mann–Whitney U for two groups
open Prelude.Resampling
fn uStatistics(g, x) ↦ mannWhitney(g, x)
fn uTotal(g) ↦ mannWhitneyTotal(g)
using Examples;
var g = new double[,] { {1,0},{1,0},{1,0},{0,1},{0,1},{0,1},{0,1} };
var s = new[] { 1.2d, 3.4, 2.2, 5.0, 2.8, 6.1, 4.4 };
var u = Resample.UStatistics(g, s);
var total = Resample.UTotal(g);
Output (verified):
U = [1, 11]
n1n2 = 12
The first group’s values have ranks 1, 4 and 2, a rank sum of 7, so \(U_1 = 7 - 6 = 1\); the two statistics add up to \(3 \cdot 4 = 12\).
Kruskal–Wallis and Wilcoxon
open Prelude.Resampling
fn kruskalH(g, x) ↦ kruskal(g, x)
fn signedRankSum(d) ↦ wilcoxon(d)
using Examples;
var g3 = 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 h = Resample.KruskalH(g3, new[] { 1d, 2, 3, 4, 5, 6, 7, 8, 9 });
var w = Resample.SignedRankSum(new[] { 1.5d, -0.5, 2.0, -3.0, 4.0 });
Output (verified):
H = 7.2
W = 10
The rank sums are 6, 15 and 24, so \(H = \tfrac{12}{90}(12 + 75 + 192) - 30 = 7.2\). For the differences, the absolute values rank 2, 1, 3, 4, 5, and the positive ones (1.5, 2.0, 4.0) have ranks 2 + 3 + 5 = 10.