Concept · intermediate
Hypothesis tests
Ask whether a difference between two groups is bigger than chance alone would make it.
You grow 8 tomato plants with compost and 8 without, and the compost plants yield more. But plants vary, so some gap would show up even if compost did nothing. A hypothesis test asks how surprising your gap would be if chance were the only thing at work.
The dataset
data.csv holds the yield of 16 plants in whole grams, with a group
column of control or compost. The data is synthetic, and we know the
truth: make_data.py gives compost plants an extra 50 g on average. The
control plants average 403.875 g and the compost plants 453.625 g, a gap
of 49.75 g.
The question, turned around
The test starts by assuming the dull explanation, called the null hypothesis: compost makes no difference. If that’s true, the labels are arbitrary. Each plant would have yielded the same whatever group it was in, so any 8 of the 16 could just as well have been the compost group.
That gives us a way to see what chance alone produces. Choose 8 plants, call them compost, and compute the gap. There are exactly
ways to choose them, few enough to try every one. That’s what makes this an exact permutation test: no random shuffling and no table of a probability distribution, just counting.
The method
- Observe the real gap: compost total minus control total.
- Relabel every possible way, 12,870 in all, and compute each gap the same way.
- Count the relabelings whose gap is at least as large as the real one, in either direction. Counting both directions makes the test two-sided: we’d be just as surprised if compost cut the yield.
- Divide by 12,870. That fraction is the p-value.
Both groups have 8 plants, so comparing totals is the same as comparing means. Totals are whole numbers, though, so every language compares them exactly. That matters: 6 relabelings give a gap of exactly 49.75 g, and they must all count as “at least as large”.
Reading the result
398 of the 12,870 relabelings give a gap of 49.75 g or more, so
If compost did nothing, a gap this large would turn up about 3% of the time. By the common convention of rejecting the null hypothesis when is below 0.05, we’d conclude compost helps.
Two cautions. First, is not the chance that compost does nothing. It’s the chance of a gap this large if compost does nothing.
Second, small samples swing. We drew 20 datasets from the same recipe, each with the same true 50 g effect, and chose this one as a clear example. Their p-values ran from 0.0002 to 0.77, and only 8 of the 20 fell below 0.05. A large p-value means the data can’t rule chance out, not that the effect is absent. With 8 plants a group, a real 50 g effect is easy to miss.
What to watch for across languages
None of these implementations calls a ready-made test. SciPy’s
permutation_test, R’s coin package, and Apache Commons Math have them,
but they aren’t installed. The interest is in how each language lists the
12,870 splits.
Some hand you the list. Python’s itertools.combinations, Ruby’s
combination, R’s combn, and Octave’s nchoosek give every choice of 8
plants in one call. Haskell’s subsequences and Prolog’s subseq come
close, but they list all 65,536 subsets, and you still have to keep the
ones of size 8.
Others count in binary. Each number from 0 to 65,535 is a 16-bit mask of
which plants are chosen, and you keep the masks with exactly eight 1s:
Java’s Integer.bitCount, Julia’s count_ones, and Common Lisp’s
logcount count them for you. C++ goes the other way, with
std::next_permutation stepping through the arrangements of eight 0s and
eight 1s. APL builds the whole 16 by 65,536 table of masks at once.
The rest write the recursion out, taking or skipping one plant at a time,
and that’s where the traps are. In awk, a function’s only local variables
are extra parameters. Leave the loop counter out of the list, and every
call shares one counter: the program reports 9 splits instead of 12,870.
In Go, a closure declared with := can’t call itself, so it’s declared
with var first. BASIC, COBOL, and Fortran skip recursion and step from
one split to the next in dictionary order, and bwBASIC adds a trap of its
own: it runs a FOR loop’s body once even when the loop should run zero
times.
Finally, watch the division that prints the means. Lua’s and Elixir’s /
always give a float. Tcl’s rounds two integers down, giving 403, so it
converts with double() first, and Common Lisp’s gives the exact fraction
3231/8, so it divides by a double instead.
Data preview
| group | grams |
|---|---|
| control | 459 |
| control | 406 |
| control | 416 |
| control | 447 |
| control | 357 |
| control | 391 |
| control | 387 |
| control | 368 |
First 8 of 16 rows · Download data.csv
Expected results
Every implementation must print these, to within a relative tolerance of 10⁻⁴.
control_count | 8 |
compost_count | 8 |
control_mean | 403.875000 |
compost_mean | 453.625000 |
observed_difference | 49.750000 |
splits | 12870 |
as_extreme | 398 |
p_value | 0.030925 |
Implementations
"""Exact two-sided permutation test for a difference in means, using only
the standard library.
If compost made no difference, the "compost" and "control" labels are
arbitrary, and any 8 of the 16 plants could have carried the compost label.
We try every one of those splits and count how many give a difference at
least as large as the one we observed. That fraction is the p-value.
"""
import csv
import sys
from itertools import combinations
def read(path):
groups = {"control": [], "compost": []}
with open(path, newline="") as f:
for row in csv.DictReader(f):
groups[row["group"]].append(int(row["grams"]))
return groups["control"], groups["compost"]
def main():
control, compost = read(sys.argv[1])
n = len(compost)
plants = control + compost
total = sum(plants)
# Both groups have n plants, so comparing totals is the same as comparing
# means, and totals stay whole numbers: no floating-point near-ties.
observed = sum(compost) - sum(control)
splits = 0
as_extreme = 0
for chosen in combinations(range(len(plants)), n):
relabeled = sum(plants[i] for i in chosen)
difference = relabeled - (total - relabeled)
splits += 1
if abs(difference) >= abs(observed): # two-sided: either direction
as_extreme += 1
print(f"control_count {len(control)}")
print(f"compost_count {n}")
print(f"control_mean {sum(control) / len(control):.6f}")
print(f"compost_mean {sum(compost) / n:.6f}")
print(f"observed_difference {observed / n:.6f}")
print(f"splits {splits}")
print(f"as_extreme {as_extreme}")
print(f"p_value {as_extreme / splits:.6f}")
if __name__ == "__main__":
main()We never shuffle at random. With 8 plants in each group, there are exactly
12,870 ways to choose which 8 carry the compost label, so we try them all
and the p-value is exact. itertools.combinations lists each choice of 8
plant positions once, in a fixed order.
Because both groups have 8 plants, comparing totals is the same as
comparing means, and the totals stay whole numbers, so >= never meets a
near-tie. Python’s standard library has no permutation test; SciPy’s
permutation_test does, but SciPy isn’t installed.
-- Exact two-sided permutation test for a difference in means, using only
-- Ada's standard library.
--
-- If compost made no difference, the "compost" and "control" labels are
-- arbitrary, and any 8 of the 16 plants could have carried the compost
-- label. We try every one of those splits and count how many give a
-- difference at least as large as the one we observed. That fraction is
-- the p-value.
with Ada.Text_IO; use Ada.Text_IO;
with Ada.Long_Float_Text_IO; use Ada.Long_Float_Text_IO;
with Ada.Command_Line;
with Ada.Strings.Fixed; use Ada.Strings.Fixed;
with Ada.Containers.Vectors;
procedure Main is
package Integer_Vectors is new Ada.Containers.Vectors (Positive, Integer);
use Integer_Vectors;
-- The text of field number Col (counting from 1) in a CSV line.
function Field (Line : String; Col : Positive) return String is
Start : Positive := Line'First;
begin
for I in 1 .. Col - 1 loop
Start := Index (Line (Start .. Line'Last), ",") + 1;
end loop;
declare
Stop : constant Natural := Index (Line (Start .. Line'Last), ",");
begin
return (if Stop = 0 then Line (Start .. Line'Last)
else Line (Start .. Stop - 1));
end;
end Field;
-- The number of the header field called Name.
function Column (Header, Name : String) return Positive is
begin
for C in 1 .. Ada.Strings.Fixed.Count (Header, ",") + 1 loop
if Field (Header, C) = Name then
return C;
end if;
end loop;
raise Constraint_Error with "column not found: " & Name;
end Column;
function Total (V : Vector) return Integer is
Sum : Integer := 0;
begin
for X of V loop
Sum := Sum + X;
end loop;
return Sum;
end Total;
procedure Put_Value (Name : String; X : Long_Float) is
begin
Put (Name & " ");
Put (X, Fore => 1, Aft => 6, Exp => 0);
New_Line;
end Put_Value;
File : File_Type;
Control, Compost : Vector;
begin
Open (File, In_File, Ada.Command_Line.Argument (1));
declare
Header : constant String := Get_Line (File);
Group_C : constant Positive := Column (Header, "group");
Grams_C : constant Positive := Column (Header, "grams");
begin
while not End_Of_File (File) loop
declare
Line : constant String := Get_Line (File);
Grams : constant Integer := Integer'Value (Field (Line, Grams_C));
begin
if Field (Line, Group_C) = "control" then
Control.Append (Grams);
else
Compost.Append (Grams);
end if;
end;
end loop;
end;
Close (File);
declare
N : constant Positive := Natural (Compost.Length);
Plants : constant Vector := Control & Compost; -- & joins vectors
All_Grams : constant Integer := Total (Plants);
-- Both groups have N plants, so comparing totals is the same as
-- comparing means, and Integer totals stay whole numbers: no
-- floating-point near-ties.
Observed : constant Integer := Total (Compost) - Total (Control);
Splits : Natural := 0;
As_Extreme : Natural := 0;
-- Visit every way to choose Left more plants from position From
-- onward. Each plant is either chosen (add it, and choose one fewer
-- from the rest) or not (choose them all from the rest). A nested
-- procedure can update Splits and As_Extreme in the block around it.
procedure Choose (From : Positive; Left : Natural; Chosen : Integer) is
begin
if Left = 0 then
declare
Difference : constant Integer :=
Chosen - (All_Grams - Chosen);
begin
Splits := Splits + 1;
-- Two-sided: a difference counts in either direction.
if abs Difference >= abs Observed then
As_Extreme := As_Extreme + 1;
end if;
end;
elsif Plants.Last_Index - From + 1 >= Left then
Choose (From + 1, Left - 1, Chosen + Plants (From));
Choose (From + 1, Left, Chosen);
end if;
end Choose;
begin
Choose (1, N, 0);
Put_Line ("control_count" & Control.Length'Image);
Put_Line ("compost_count" & N'Image);
Put_Value ("control_mean",
Long_Float (Total (Control)) / Long_Float (Control.Length));
Put_Value ("compost_mean", Long_Float (Total (Compost)) / Long_Float (N));
Put_Value ("observed_difference", Long_Float (Observed) / Long_Float (N));
Put_Line ("splits" & Splits'Image);
Put_Line ("as_extreme" & As_Extreme'Image);
Put_Value ("p_value", Long_Float (As_Extreme) / Long_Float (Splits));
end;
end Main;We never shuffle at random. A recursive procedure, Choose, walks the
plants in order and either chooses each one or skips it, so it reaches
all 12,870 ways to pick 8 of the 16, and the p-value is exact. Choose is
nested inside the block that declares Splits and As_Extreme, so it
can update those counters directly, with no globals or extra parameters.
The & operator from Ada.Containers.Vectors joins the two groups into
one vector. Totals are Integer, so >= compares whole numbers, and the
6 exact ties all count.
⍝ Exact two-sided permutation test for a difference in means, from scratch
⍝ in GNU APL. If compost made no difference, the "compost" and "control"
⍝ labels are arbitrary, and any 8 of the 16 plants could have carried the
⍝ compost label. We try every one of those splits and count how many give a
⍝ difference at least as large as the one we observed.
⍝ APL reads right to left, with no operator precedence: 2×3+4 is 2×(3+4).
⍝ ⍝ starts a comment; ← assigns; ⍵ is a function's right argument.
∇NAME SAY N
⍝ Print "name value" for a whole number. ⍕ formats a number as text.
⎕←NAME,' ',⍕N
∇
∇NAME SAY6 X;T
⍝ Print "name value" with six decimal places.
T←0 6⍕|X ⍝ |X is the size of X; 0 6⍕ gives 6 decimals
T←(∨\T≠' ')/T ⍝ drop the leading blank ⍕ adds
T←((T[1]='.')/'0'),T ⍝ ⍕ writes .03 for 0.03, so put the 0 back
⎕←NAME,' ',((X<0)/'-'),T ⍝ sign as -, since APL's own minus is ¯
∇
⍝ Read the file. ⎕ARG is the command line; ↑¯1↑ takes its last item, the
⍝ data path. ⎕FIO[26] reads the whole file as one string.
TEXT←⎕FIO[26]↑¯1↑⎕ARG
LINES←(TEXT≠⎕UCS 10)⊂TEXT ⍝ ⊂ splits at the 0s of the mask: one item per line
⍝ Split each (¨) line at commas; ⊃ stacks them into a table.
ROWS←⊃{(⍵≠',')⊂⍵}¨LINES
⍝ A column by header name (⍳), minus the header (1↓).
COLUMN←{1↓ROWS[;ROWS[1;]⍳⊂⍵]}
GRAMS←⍎¨COLUMN 'grams' ⍝ ⍎ turns text into a number
COMPOST←(COLUMN 'group')≡¨⊂'compost' ⍝ 1 for each compost plant
N←+/COMPOST ⍝ the compost group's size
CONTROLSUM←+/(~COMPOST)/GRAMS ⍝ ~ is not; / keeps where the mask is 1
COMPOSTSUM←+/COMPOST/GRAMS
TOTAL←+/GRAMS
⍝ Both groups have the same size, so comparing totals is the same as
⍝ comparing means, and totals stay whole numbers: no floating-point near-ties.
OBSERVED←COMPOSTSUM-CONTROLSUM
⍝ Every split at once. (P⍴2)⊤ writes each whole number from 0 to 2*P - 1 in
⍝ base 2, one column per number and one row per plant: a 16 by 65536 table
⍝ of 0s and 1s, where 1 means "labeled compost".
P←≢GRAMS
MASKS←(P⍴2)⊤¯1+⍳2*P
⍝ +⌿ sums each column; keep (/) the columns with exactly N ones.
MASKS←(N=+⌿MASKS)/MASKS
⍝ +.× is the matrix product: each split's chosen grams in one step.
CHOSEN←GRAMS+.×MASKS
DIFFERENCE←CHOSEN-TOTAL-CHOSEN ⍝ right to left: CHOSEN-(TOTAL-CHOSEN)
SPLITS←≢CHOSEN
ASEXTREME←+/(|DIFFERENCE)≥|OBSERVED ⍝ two-sided: | is absolute value
'control_count' SAY +/~COMPOST
'compost_count' SAY N
'control_mean' SAY6 CONTROLSUM÷N
'compost_mean' SAY6 COMPOSTSUM÷N
'observed_difference' SAY6 OBSERVED÷N
'splits' SAY SPLITS
'as_extreme' SAY ASEXTREME
'p_value' SAY6 ASEXTREME÷SPLITS
)OFFAPL tests every split at once, without a loop. (P⍴2)⊤¯1+⍳2*P writes each
number from 0 to 65,535 in base 2, giving a 16 by 65,536 table of 0s and
1s, one column per labeling. (N=+⌿MASKS)/MASKS keeps the 12,870 columns
with exactly 8 ones.
The inner product GRAMS+.×MASKS then multiplies and sums down every
column in one step, giving each split’s chosen total. The grams are whole
numbers, so ≥ compares exact integers, and the 6 splits that tie the
observed difference all count.
# Exact two-sided permutation test for a difference in means, in POSIX awk.
#
# If compost made no difference, the "compost" and "control" labels are
# arbitrary, and any 8 of the 16 plants could have carried the compost label.
# We try every one of those splits and count how many give a difference at
# least as large as the one we observed. That fraction is the p-value. Rows
# are stored as they are read; the work happens in END.
BEGIN { FS = "," }
# The header line: find each column by name.
NR == 1 {
for (i = 1; i <= NF; i++) {
if ($i == "group") c_group = i
if ($i == "grams") c_grams = i
}
if (!c_group || !c_grams) {
print "missing a column (group, grams)" > "/dev/stderr"
failed = 1
exit 1
}
next
}
# Every plant goes into one array; the group totals are kept as we go.
{
n_plants++
plant[n_plants] = $c_grams + 0
if ($c_group == "control") { n_control++; control_total += $c_grams }
if ($c_group == "compost") { n_compost++; compost_total += $c_grams }
}
function abs(x) { return x < 0 ? -x : x }
# Walk every way to pick `left` more plants from position start onward,
# carrying the running total of the plants picked so far. awk functions can
# call themselves, and i is local (declared as an extra parameter), so each
# call keeps its own loop counter.
function choose(start, left, chosen_total, i, difference) {
if (left == 0) {
difference = chosen_total - (total - chosen_total)
splits++
if (abs(difference) >= abs(observed)) as_extreme++ # two-sided
return
}
# Stop early enough that `left` plants still remain to pick.
for (i = start; i <= n_plants - left + 1; i++)
choose(i + 1, left - 1, chosen_total + plant[i])
}
END {
if (failed) exit 1
total = control_total + compost_total
# Both groups have the same size, so comparing totals is the same as
# comparing means. The totals are whole numbers, and awk's doubles hold
# whole numbers this size exactly, so >= never meets a near-tie.
observed = compost_total - control_total
choose(1, n_compost, 0)
printf "control_count %d\n", n_control
printf "compost_count %d\n", n_compost
printf "control_mean %.6f\n", control_total / n_control
printf "compost_mean %.6f\n", compost_total / n_compost
printf "observed_difference %.6f\n", observed / n_compost
printf "splits %d\n", splits
printf "as_extreme %d\n", as_extreme
printf "p_value %.6f\n", as_extreme / splits
}We never shuffle at random. A recursive function, choose, walks every way
to pick 8 of the 16 plants, carrying a running total, and counts all 12,870
splits. awk has no local variables except extra parameters, so the loop
counter i sits after the real ones. Leave it out and every call shares one
i: on this data the run then reports 9 splits instead of 12,870.
Every awk number is a double, but whole numbers this size are exact, so the
totals compare without near-ties. mawk has no bit operations, since
and(5, 3) is an undefined function, so the splits come from recursion.
10 REM Exact two-sided permutation test for a difference in means, in
20 REM classic BASIC.
30 REM
40 REM If compost made no difference, the group labels are arbitrary, and
50 REM any K of the N plants could have carried the compost label. We try
60 REM every one of those splits and count how many give a difference at
70 REM least as large as the one we observed. That fraction is the p-value.
80 REM
90 REM run.sh passes the data path in the DATA environment variable,
100 REM because bwbasic can't read command-line arguments.
110 DIM G(100), IX(100), FL$(20) : REM grams per plant, chosen positions
120 GOSUB 1000 : REM read the plants into G(1..N), group sums SA and SB
130 K = NB : REM both groups have K plants
140 T = SA + SB
150 REM Both groups have K plants, so comparing totals is the same as
160 REM comparing means, and totals stay whole numbers: no near-ties.
170 OB = ABS(SB - SA)
180 REM IX(1..K) holds the chosen positions in increasing order. It starts
190 REM at 1, 2, ..., K and steps through every choice in dictionary order.
200 FOR J = 1 TO K
210 IX(J) = J
220 NEXT J
230 SP = 0 : AE = 0
240 R = 0
250 FOR J = 1 TO K
260 R = R + G(IX(J))
270 NEXT J
280 SP = SP + 1
290 REM Two-sided: ABS counts a difference in either direction.
300 IF ABS(R - (T - R)) >= OB THEN AE = AE + 1
310 REM Find the rightmost position that can still move up. Position J can
320 REM go as high as N - K + J, leaving room for the ones after it.
330 J = K
340 IF J < 1 THEN 400
350 IF IX(J) < N - K + J THEN 360
355 J = J - 1 : GOTO 340
360 IX(J) = IX(J) + 1
365 REM Refill the positions after J in a run. bwbasic runs a FOR body at
370 REM least once, even FOR M = 4 TO 3, so skip the loop when J = K.
375 IF J = K THEN 240
380 FOR M = J + 1 TO K
385 IX(M) = IX(M - 1) + 1
390 NEXT M
395 GOTO 240
400 PRINT "control_count"; NA
410 PRINT "compost_count"; NB
420 L$ = "control_mean" : Y = SA / NA : GOSUB 3000
430 L$ = "compost_mean" : Y = SB / NB : GOSUB 3000
440 L$ = "observed_difference" : Y = (SB - SA) / K : GOSUB 3000
450 PRINT "splits"; SP
460 PRINT "as_extreme"; AE
470 L$ = "p_value" : Y = AE / SP : GOSUB 3000
480 SYSTEM : REM without SYSTEM, bwbasic waits at its prompt
1000 REM Read every plant into G, and total each group as we go.
1010 F$ = ENVIRON$("DATA")
1020 OPEN "I", #1, F$
1030 LINE INPUT #1, L$
1040 GOSUB 4000 : REM find which field of the header holds each column
1050 FOR J = 1 TO NF
1060 IF FL$(J) = "group" THEN CG = J
1070 IF FL$(J) = "grams" THEN CM = J
1080 NEXT J
1090 N = 0 : NA = 0 : NB = 0 : SA = 0 : SB = 0
1100 IF EOF(1) THEN 1190
1110 LINE INPUT #1, L$
1120 GOSUB 4000
1130 N = N + 1 : G(N) = VAL(FL$(CM))
1140 IF FL$(CG) = "control" THEN NA = NA + 1 : SA = SA + G(N)
1150 IF FL$(CG) = "compost" THEN NB = NB + 1 : SB = SB + G(N)
1160 GOTO 1100
1190 CLOSE #1
1200 RETURN
3000 REM Print L$ and Y with six decimals.
3010 PRINT L$; " ";
3020 PRINT USING "#.######"; Y
3030 RETURN
4000 REM Split L$ at its commas into FL$(1) to FL$(NF).
4010 NF = 0 : P = 1
4020 Q = INSTR(P, L$, ",")
4030 NF = NF + 1
4040 IF Q = 0 THEN FL$(NF) = MID$(L$, P) : RETURN
4050 FL$(NF) = MID$(L$, P, Q - P) : P = Q + 1
4060 GOTO 4020We step through the splits in place, with no recursion. IX(1) to IX(K)
hold the chosen positions in increasing order, starting at 1 to 8. Each
step finds the rightmost position that can still move up, adds one to it,
and refills the positions after it in a run. That visits all 12,870 choices
in dictionary order.
One quirk shaped the code: bwbasic runs a FOR body at least once, even
for FOR M = 4 TO 3, so line 375 skips the refill loop by hand. Totals are
whole numbers, so the >= test has no rounding to trip on.
/* Exact two-sided permutation test for a difference in means, using only
* the C standard library.
*
* If compost made no difference, the "compost" and "control" labels are
* arbitrary, and any 8 of the 16 plants could have carried the compost
* label. We try every one of those splits and count how many give a
* difference at least as large as the one we observed. That fraction is
* the p-value. */
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#define MAX_PLANTS 64
#define MAX_LINE 256
/* Every plant's grams, controls first, then compost, in file order. */
static long plants[MAX_PLANTS];
static int plant_count;
/* Running tallies for the enumeration below. */
static long total, threshold;
static long splits, as_extreme;
/* Cut a line at its commas, in place, and return the number of fields.
* strtok replaces each comma with '\0', so fields[i] points into line. */
static int split(char *line, char *fields[], int max) {
int n = 0;
for (char *f = strtok(line, ",\r\n"); f && n < max; f = strtok(NULL, ",\r\n"))
fields[n++] = f;
return n;
}
/* Choose `left` more plants from positions start onward, with `chosen`
* grams already chosen. Each call either takes plant `start` or skips it,
* so every set of positions is reached exactly once. */
static void choose(int start, int left, long chosen) {
if (left == 0) {
long difference = chosen - (total - chosen);
splits++;
/* labs is absolute value for a long; whole numbers compare exactly. */
if (labs(difference) >= threshold) as_extreme++; /* two-sided */
return;
}
if (plant_count - start < left) return; /* too few plants remain */
choose(start + 1, left - 1, chosen + plants[start]); /* take it */
choose(start + 1, left, chosen); /* skip it */
}
int main(int argc, char **argv) {
if (argc < 2) { fprintf(stderr, "usage: %s data.csv\n", argv[0]); return 1; }
FILE *file = fopen(argv[1], "r");
if (!file) { perror(argv[1]); return 1; }
char line[MAX_LINE], *fields[4];
if (!fgets(line, sizeof line, file)) { fprintf(stderr, "empty file\n"); return 1; }
long control_sum = 0, compost_sum = 0;
int control_count = 0, compost_count = 0;
while (fgets(line, sizeof line, file) && plant_count < MAX_PLANTS) {
if (split(line, fields, 4) != 2) continue;
long grams = strtol(fields[1], NULL, 10);
if (strcmp(fields[0], "compost") == 0) {
compost_sum += grams;
compost_count++;
} else {
control_sum += grams;
control_count++;
}
plants[plant_count++] = grams;
}
fclose(file);
if (compost_count == 0 || compost_count != control_count) {
fprintf(stderr, "need two groups of equal size\n");
return 1;
}
/* Both groups have the same size, so comparing totals is the same as
* comparing means, and totals stay whole numbers: no floating-point
* near-ties. */
total = control_sum + compost_sum;
long observed = compost_sum - control_sum;
threshold = labs(observed);
choose(0, compost_count, 0);
printf("control_count %d\n", control_count);
printf("compost_count %d\n", compost_count);
printf("control_mean %.6f\n", (double)control_sum / control_count);
printf("compost_mean %.6f\n", (double)compost_sum / compost_count);
printf("observed_difference %.6f\n", (double)observed / compost_count);
printf("splits %ld\n", splits);
printf("as_extreme %ld\n", as_extreme);
printf("p_value %.6f\n", (double)as_extreme / splits);
return 0;
}C’s standard library has no combinations helper, so a recursive function,
choose, builds the splits. At each plant it either takes the plant, adding
its grams to the running total, or skips it, and it stops when 8 are taken,
so each of the 12,870 ways to pick 8 of the 16 positions is reached exactly
once. It also gives up early when too few plants remain to finish a split.
Every total is a long, so the test compares whole numbers with labs and
>=, and the 6 splits that tie the observed difference count exactly. We
convert to double only to print the means and the p-value.
// Exact two-sided permutation test for a difference in means, using only the
// C++ standard library.
//
// If compost made no difference, the "compost" and "control" labels are
// arbitrary, and any 8 of the 16 plants could have carried the compost label.
// We try every one of those splits and count how many give a difference at
// least as large as the one we observed. That fraction is the p-value.
#include <algorithm>
#include <cstdlib>
#include <fstream>
#include <numeric>
#include <print>
#include <stdexcept>
#include <string>
#include <vector>
int main(int argc, char* argv[]) {
if (argc != 2) {
std::println(stderr, "usage: program data.csv");
return 1;
}
// Read group,grams lines. Controls come first in the file, then compost.
std::ifstream file(argv[1]);
std::string line;
std::getline(file, line); // skip the header
std::vector<long> control, compost;
while (std::getline(file, line)) {
auto comma = line.find(',');
if (comma == std::string::npos) {
continue;
}
long grams = std::stol(line.substr(comma + 1));
(line.substr(0, comma) == "compost" ? compost : control).push_back(grams);
}
if (compost.empty() || compost.size() != control.size()) {
throw std::runtime_error("need two groups of equal size");
}
std::vector<long> plants = control;
plants.insert(plants.end(), compost.begin(), compost.end());
long control_sum = std::reduce(control.begin(), control.end(), 0L);
long compost_sum = std::reduce(compost.begin(), compost.end(), 0L);
long total = control_sum + compost_sum;
// Both groups have the same size, so comparing totals is the same as
// comparing means, and totals stay whole numbers: no floating-point
// near-ties.
long observed = compost_sum - control_sum;
// A mask of 0s and 1s, one per plant: 1 means "labeled compost". It
// starts sorted, all 0s then all 1s, and std::next_permutation steps it
// through every distinct arrangement once, returning false after the last.
std::vector<int> mask(plants.size(), 0);
std::fill(mask.begin() + control.size(), mask.end(), 1);
long splits = 0, as_extreme = 0;
do {
// transform_reduce multiplies each plant by its 0 or 1 and adds them up.
long chosen = std::transform_reduce(plants.begin(), plants.end(), mask.begin(), 0L);
long difference = chosen - (total - chosen);
splits++;
if (std::labs(difference) >= std::labs(observed)) { // two-sided
as_extreme++;
}
} while (std::next_permutation(mask.begin(), mask.end()));
double n = static_cast<double>(compost.size());
std::println("control_count {}", control.size());
std::println("compost_count {}", compost.size());
std::println("control_mean {:.6f}", control_sum / n);
std::println("compost_mean {:.6f}", compost_sum / n);
std::println("observed_difference {:.6f}", observed / n);
std::println("splits {}", splits);
std::println("as_extreme {}", as_extreme);
std::println("p_value {:.6f}", static_cast<double>(as_extreme) / splits);
}We enumerate the splits with a mask: a vector of eight 0s then eight 1s,
where 1 means “labeled compost”. std::next_permutation steps the mask to
its next arrangement in sorted order and returns false after the last one.
Because equal values are interchangeable, it visits each distinct
arrangement once: 12,870 of them, not 16 factorial.
std::transform_reduce multiplies each plant’s grams by its 0 or 1 and adds
the products, giving the chosen total. Totals stay long integers, so >=
counts the 6 exact ties correctly.
*> Exact two-sided permutation test for a difference in means,
*> in COBOL.
*>
*> If compost made no difference, the "compost" and "control"
*> labels are arbitrary, and any 8 of the 16 plants could have
*> carried the compost label. We try every one of those splits
*> and count how many give a difference at least as large as the
*> one we observed. That fraction is the p-value.
IDENTIFICATION DIVISION.
PROGRAM-ID. HYPOTHESIS-TESTS.
ENVIRONMENT DIVISION.
INPUT-OUTPUT SECTION.
FILE-CONTROL.
SELECT DATA-FILE ASSIGN TO DYNAMIC WS-FILE-NAME
ORGANIZATION IS LINE SEQUENTIAL.
DATA DIVISION.
FILE SECTION.
FD DATA-FILE.
01 DATA-LINE PIC X(200).
WORKING-STORAGE SECTION.
01 WS-FILE-NAME PIC X(256).
01 WS-EOF PIC X VALUE "N".
88 END-OF-FILE VALUE "Y".
01 WS-GROUP PIC X(40).
01 WS-GRAMS-TEXT PIC X(40).
*> Every plant's grams, in file order. Grams are whole numbers,
*> so every total below is an exact integer.
01 WS-PLANT-COUNT PIC 99 VALUE 0.
01 WS-PLANTS.
05 WS-GRAMS PIC 9(6) OCCURS 64 TIMES.
01 WS-CONTROL-COUNT PIC 99 VALUE 0.
01 WS-COMPOST-COUNT PIC 99 VALUE 0.
01 WS-CONTROL-SUM PIC 9(9) VALUE 0.
01 WS-COMPOST-SUM PIC 9(9) VALUE 0.
01 WS-TOTAL PIC 9(9) VALUE 0.
01 WS-OBSERVED PIC S9(9).
01 WS-THRESHOLD PIC 9(9).
*> The current split: the positions of the plants labeled
*> compost, in increasing order, such as 1 2 3 ... 8.
01 WS-PICKS.
05 WS-PICK PIC 99 OCCURS 64 TIMES.
01 WS-CHOSEN PIC 9(9).
01 WS-DIFFERENCE PIC S9(9).
01 WS-SPLITS PIC 9(9) VALUE 0.
01 WS-AS-EXTREME PIC 9(9) VALUE 0.
01 WS-DONE PIC X VALUE "N".
88 ALL-SPLITS-DONE VALUE "Y".
01 I PIC 99.
01 J PIC 99.
01 WS-RESULT PIC S9(9)V9(12).
01 WS-RESULT-OUT PIC -(9)9.9(6).
01 WS-COUNT-OUT PIC Z(8)9.
PROCEDURE DIVISION.
MAIN.
ACCEPT WS-FILE-NAME FROM ARGUMENT-VALUE
OPEN INPUT DATA-FILE
*> The first line is the header, group,grams; skip it.
READ DATA-FILE
AT END SET END-OF-FILE TO TRUE
END-READ
PERFORM UNTIL END-OF-FILE
READ DATA-FILE
AT END SET END-OF-FILE TO TRUE
NOT AT END PERFORM ADD-PLANT
END-READ
END-PERFORM
CLOSE DATA-FILE
IF WS-COMPOST-COUNT = 0
OR WS-COMPOST-COUNT NOT = WS-CONTROL-COUNT
DISPLAY "need two groups of equal size" UPON SYSERR
MOVE 1 TO RETURN-CODE
STOP RUN
END-IF
*> Both groups have the same size, so comparing totals is the
*> same as comparing means.
COMPUTE WS-TOTAL = WS-CONTROL-SUM + WS-COMPOST-SUM
COMPUTE WS-OBSERVED = WS-COMPOST-SUM - WS-CONTROL-SUM
COMPUTE WS-THRESHOLD = FUNCTION ABS(WS-OBSERVED)
*> Start from the first split, positions 1 to N, and step
*> through every split in increasing order.
PERFORM VARYING I FROM 1 BY 1 UNTIL I > WS-COMPOST-COUNT
MOVE I TO WS-PICK(I)
END-PERFORM
PERFORM UNTIL ALL-SPLITS-DONE
PERFORM SCORE-SPLIT
PERFORM NEXT-SPLIT
END-PERFORM
MOVE WS-CONTROL-COUNT TO WS-COUNT-OUT
DISPLAY "control_count " FUNCTION TRIM(WS-COUNT-OUT)
MOVE WS-COMPOST-COUNT TO WS-COUNT-OUT
DISPLAY "compost_count " FUNCTION TRIM(WS-COUNT-OUT)
COMPUTE WS-RESULT ROUNDED =
WS-CONTROL-SUM / WS-CONTROL-COUNT
COMPUTE WS-RESULT-OUT ROUNDED = WS-RESULT
DISPLAY "control_mean " FUNCTION TRIM(WS-RESULT-OUT)
COMPUTE WS-RESULT ROUNDED =
WS-COMPOST-SUM / WS-COMPOST-COUNT
COMPUTE WS-RESULT-OUT ROUNDED = WS-RESULT
DISPLAY "compost_mean " FUNCTION TRIM(WS-RESULT-OUT)
COMPUTE WS-RESULT ROUNDED = WS-OBSERVED / WS-COMPOST-COUNT
COMPUTE WS-RESULT-OUT ROUNDED = WS-RESULT
DISPLAY "observed_difference " FUNCTION TRIM(WS-RESULT-OUT)
MOVE WS-SPLITS TO WS-COUNT-OUT
DISPLAY "splits " FUNCTION TRIM(WS-COUNT-OUT)
MOVE WS-AS-EXTREME TO WS-COUNT-OUT
DISPLAY "as_extreme " FUNCTION TRIM(WS-COUNT-OUT)
COMPUTE WS-RESULT ROUNDED = WS-AS-EXTREME / WS-SPLITS
COMPUTE WS-RESULT-OUT ROUNDED = WS-RESULT
DISPLAY "p_value " FUNCTION TRIM(WS-RESULT-OUT)
STOP RUN.
*> Store one group,grams line.
ADD-PLANT.
MOVE SPACES TO WS-GROUP WS-GRAMS-TEXT
UNSTRING DATA-LINE DELIMITED BY ","
INTO WS-GROUP WS-GRAMS-TEXT
END-UNSTRING
IF WS-GROUP = SPACES
EXIT PARAGRAPH
END-IF
ADD 1 TO WS-PLANT-COUNT
COMPUTE WS-GRAMS(WS-PLANT-COUNT) =
FUNCTION NUMVAL(WS-GRAMS-TEXT)
IF WS-GROUP = "compost"
ADD 1 TO WS-COMPOST-COUNT
ADD WS-GRAMS(WS-PLANT-COUNT) TO WS-COMPOST-SUM
ELSE
ADD 1 TO WS-CONTROL-COUNT
ADD WS-GRAMS(WS-PLANT-COUNT) TO WS-CONTROL-SUM
END-IF.
*> Count the current split, and whether it is as extreme as the
*> observed one, in either direction.
SCORE-SPLIT.
MOVE 0 TO WS-CHOSEN
PERFORM VARYING I FROM 1 BY 1 UNTIL I > WS-COMPOST-COUNT
ADD WS-GRAMS(WS-PICK(I)) TO WS-CHOSEN
END-PERFORM
COMPUTE WS-DIFFERENCE = WS-CHOSEN - (WS-TOTAL - WS-CHOSEN)
ADD 1 TO WS-SPLITS
IF FUNCTION ABS(WS-DIFFERENCE) >= WS-THRESHOLD
ADD 1 TO WS-AS-EXTREME
END-IF.
*> Step to the next split in increasing order: find the last pick
*> that can still move right, move it one place, and put the
*> picks after it just behind it. Pick I can go no further than
*> position PLANT-COUNT - N + I, so when even pick 1 is there,
*> that was the last split.
NEXT-SPLIT.
MOVE WS-COMPOST-COUNT TO I
PERFORM UNTIL I = 0
OR WS-PICK(I) <
WS-PLANT-COUNT - WS-COMPOST-COUNT + I
SUBTRACT 1 FROM I
END-PERFORM
IF I = 0
SET ALL-SPLITS-DONE TO TRUE
ELSE
ADD 1 TO WS-PICK(I)
PERFORM VARYING J FROM I BY 1
UNTIL J >= WS-COMPOST-COUNT
COMPUTE WS-PICK(J + 1) = WS-PICK(J) + 1
END-PERFORM
END-IF.We step through the splits in increasing order with a table of 8 plant
positions, starting at 1 through 8. NEXT-SPLIT finds the last position
that can still move right, moves it one place, and packs the positions
after it just behind it. When even the first position can’t move, all
12,870 splits are done.
The grams are whole-number PIC 9(6) fields, so every total is an exact
integer, and FUNCTION ABS with >= counts the 6 ties exactly. The means
are the group totals divided by 8, not FUNCTION MEAN.
;;;; Exact two-sided permutation test for a difference in means, using only
;;;; standard Common Lisp.
;;;;
;;;; If compost made no difference, the "compost" and "control" labels are
;;;; arbitrary, and any 8 of the 16 plants could have carried the compost
;;;; label. We try every one of those splits and count how many give a
;;;; difference at least as large as the one we observed. That fraction is
;;;; the p-value.
(defun split-on-comma (line)
"The fields of LINE, split at every comma."
(loop for start = 0 then (1+ end)
for end = (position #\, line :start start)
collect (subseq line start end)
while end))
(defun read-plants (path)
"The rows of the CSV file at PATH, as (group . grams) pairs in file order."
(with-open-file (in path)
(let* ((header (split-on-comma (read-line in)))
(group-at (or (position "group" header :test #'string=)
(error "no column group")))
(grams-at (or (position "grams" header :test #'string=)
(error "no column grams"))))
(loop for line = (read-line in nil)
while line
unless (string= line "")
collect (let ((fields (split-on-comma line)))
;; parse-integer reads digits only; it never runs code.
(cons (nth group-at fields)
(parse-integer (nth grams-at fields))))))))
(defun grams-of (group rows)
"The grams of every row in GROUP."
(loop for (g . grams) in rows
when (string= g group) collect grams))
(let* ((rows (read-plants (second sb-ext:*posix-argv*)))
(control (grams-of "control" rows))
(compost (grams-of "compost" rows))
(n (length compost))
(plants (coerce (append control compost) 'vector))
(total (reduce #'+ plants))
;; Both groups have N plants, so comparing totals is the same as
;; comparing means, and integer totals stay whole numbers: no
;; floating-point near-ties.
(observed (- (reduce #'+ compost) (reduce #'+ control)))
(splits 0)
(as-extreme 0))
;; Each whole number from 0 below 2^16 is a set of plants: bit I is on
;; when plant I is chosen. (ash 1 16) is 2^16. logcount counts the 1 bits,
;; so masks with exactly N of them are the splits, and logbitp tests one bit.
(dotimes (mask (ash 1 (length plants)))
(when (= (logcount mask) n)
(let* ((chosen (loop for i below (length plants)
when (logbitp i mask) sum (aref plants i)))
(difference (- chosen (- total chosen))))
(incf splits)
;; Two-sided: a difference counts in either direction.
(when (>= (abs difference) (abs observed))
(incf as-extreme)))))
;; ~d prints an integer, ~,6f a number with six decimals, ~% a newline.
;; (/ 3231 8) would be the exact fraction 3231/8; dividing by a double
;; gives a double.
(format t "control_count ~d~%" (length control))
(format t "compost_count ~d~%" n)
(format t "control_mean ~,6f~%" (/ (reduce #'+ control) (float (length control) 1d0)))
(format t "compost_mean ~,6f~%" (/ (reduce #'+ compost) (float n 1d0)))
(format t "observed_difference ~,6f~%" (/ observed (float n 1d0)))
(format t "splits ~d~%" splits)
(format t "as_extreme ~d~%" as-extreme)
(format t "p_value ~,6f~%" (/ as-extreme (float splits 1d0))))We never shuffle at random. Each whole number from 0 to 65,535 stands for
a set of the 16 plants, one bit per plant, and logcount counts its 1
bits. The 12,870 numbers with exactly 8 are the splits, and logbitp
tells us which plants each one chose, so the p-value is exact.
The grams are read with parse-integer, and the totals stay integers, so
>= compares whole numbers and the 6 exact ties all count. The means
divide by a double, since (/ 3231 8) would give the fraction 3231/8.
# Exact two-sided permutation test for a difference in means, using only
# Elixir's standard library.
#
# If compost made no difference, the "compost" and "control" labels are
# arbitrary, and any 8 of the 16 plants could have carried the compost
# label. We try every one of those splits and count how many give a
# difference at least as large as the one we observed. That fraction is the
# p-value.
defmodule HypothesisTest do
# Every row as a map, keyed by the header's column names.
def read_rows(path) do
[header | rows] =
path
|> File.read!()
|> String.split("\n", trim: true)
|> Enum.map(&String.split(&1, ","))
Enum.map(rows, fn fields -> header |> Enum.zip(fields) |> Map.new() end)
end
# The grams of every row in one group, as integers.
def grams(rows, group) do
for %{"group" => ^group, "grams" => text} <- rows, do: String.to_integer(text)
end
# The total of every way to choose k of the plants. The first plant is
# either chosen (add it, and choose k - 1 from the rest) or not (choose
# all k from the rest). The list has one total per split.
def chosen_totals(0, _plants), do: [0]
def chosen_totals(_k, []), do: []
def chosen_totals(k, [plant | rest]) do
Enum.map(chosen_totals(k - 1, rest), &(&1 + plant)) ++ chosen_totals(k, rest)
end
end
[path] = System.argv()
rows = HypothesisTest.read_rows(path)
control = HypothesisTest.grams(rows, "control")
compost = HypothesisTest.grams(rows, "compost")
n = length(compost)
total = Enum.sum(control) + Enum.sum(compost)
# Both groups have n plants, so comparing totals is the same as comparing
# means, and integer totals stay whole numbers: no floating-point near-ties.
observed = Enum.sum(compost) - Enum.sum(control)
differences =
for chosen <- HypothesisTest.chosen_totals(n, control ++ compost) do
chosen - (total - chosen)
end
splits = length(differences)
# Two-sided: a difference counts in either direction.
as_extreme = Enum.count(differences, &(abs(&1) >= abs(observed)))
# / always returns a float in Elixir, even for 3232 / 8.
format = &:erlang.float_to_binary(&1, decimals: 6)
IO.puts("control_count #{length(control)}")
IO.puts("compost_count #{n}")
IO.puts("control_mean #{format.(Enum.sum(control) / length(control))}")
IO.puts("compost_mean #{format.(Enum.sum(compost) / n)}")
IO.puts("observed_difference #{format.(observed / n)}")
IO.puts("splits #{splits}")
IO.puts("as_extreme #{as_extreme}")
IO.puts("p_value #{format.(as_extreme / splits)}")We never shuffle at random. Enum.shuffle/1 would give one random
relabeling, but with 8 plants in each group we can try all 12,870 splits,
so the p-value is exact. Enum has no combinations function, so the
recursive chosen_totals/2 either adds the first plant and chooses one
fewer from the rest, or skips it.
The pattern %{"group" => ^group} in a for keeps only one group’s rows;
the pin ^ matches the existing value. Totals stay integers, so >=
compares whole numbers, and the 6 exact ties all count.
! Exact two-sided permutation test for a difference in means, using only
! standard Fortran.
!
! If compost made no difference, the "compost" and "control" labels are
! arbitrary, and any 8 of the 16 plants could have carried the compost label.
! We try every one of those splits and count how many give a difference at
! least as large as the one we observed. That fraction is the p-value.
program hypothesis_test
implicit none
integer, parameter :: dp = kind(1.0d0) ! double precision
integer, allocatable :: control(:), compost(:), plants(:), chosen(:)
integer :: n, m, total, observed, relabeled, difference, i, j
integer :: splits, as_extreme
call read_plants(control, compost)
n = size(compost)
plants = [control, compost] ! the array constructor joins the two groups
m = size(plants)
total = sum(plants)
! Both groups have n plants, so comparing totals is the same as comparing
! means, and integer totals stay whole numbers: no floating-point near-ties.
observed = sum(compost) - sum(control)
! chosen holds the positions of the n relabeled plants, in increasing
! order, starting from 1, 2, ..., n. Each pass moves to the next set in
! dictionary order, so every set comes up exactly once.
chosen = [(i, i = 1, n)] ! an implied-do loop builds [1, 2, ..., n]
splits = 0
as_extreme = 0
do
! plants(chosen) is a vector subscript: the elements at those positions.
relabeled = sum(plants(chosen))
difference = relabeled - (total - relabeled)
splits = splits + 1
if (abs(difference) >= abs(observed)) as_extreme = as_extreme + 1 ! two-sided
! Find the rightmost position that can still move up. Position i can go
! no higher than m - n + i, or the positions after it would run out.
i = n
do while (i >= 1)
if (chosen(i) < m - n + i) exit
i = i - 1
end do
if (i == 0) exit ! the last set, [m - n + 1, ..., m], is done
! Move it up by one, and set the positions after it to the next numbers
! in a row: if chosen(i) becomes 5, the ones after are 6, 7, and so on.
chosen(i:) = chosen(i) + [(j, j = 1, n - i + 1)]
end do
print "(a, 1x, i0)", "control_count", size(control)
print "(a, 1x, i0)", "compost_count", n
call say("control_mean", real(sum(control), dp) / size(control))
call say("compost_mean", real(sum(compost), dp) / n)
call say("observed_difference", real(observed, dp) / n)
print "(a, 1x, i0)", "splits", splits
print "(a, 1x, i0)", "as_extreme", as_extreme
call say("p_value", real(as_extreme, dp) / splits)
contains
subroutine read_plants(control, compost)
integer, allocatable, intent(out) :: control(:), compost(:)
character(len=256) :: path, line, text
integer :: unit, status, group_col, grams_col, grams
call get_command_argument(1, path)
open(newunit=unit, file=path, status="old", action="read")
read(unit, "(a)") line
group_col = field_index(line, "group")
grams_col = field_index(line, "grams")
allocate(control(0), compost(0))
do
read(unit, "(a)", iostat=status) line
if (status /= 0) exit
if (len_trim(line) == 0) cycle
! An internal read parses a number. It must read from a variable,
! not straight from the result of field(), so we copy into text.
text = field(line, grams_col)
read(text, *) grams
if (field(line, group_col) == "control") then
control = [control, grams] ! grows the array by one
else
compost = [compost, grams]
end if
end do
close(unit)
end subroutine read_plants
! The position of a name in a comma-separated header line.
integer function field_index(header, name)
character(len=*), intent(in) :: header, name
integer :: i, fields
fields = 1
do i = 1, len_trim(header)
if (header(i:i) == ",") fields = fields + 1
end do
do field_index = 1, fields
if (field(header, field_index) == name) return
end do
error stop "column not found: " // name
end function field_index
! The k-th comma-separated field of a line.
function field(line, k) result(text)
character(len=*), intent(in) :: line
integer, intent(in) :: k
character(len=:), allocatable :: text
integer :: start, i, comma
start = 1
do i = 1, k - 1
start = start + index(line(start:), ",")
end do
comma = index(line(start:), ",")
if (comma == 0) then
text = trim(line(start:))
else
text = line(start:start + comma - 2)
end if
end function field
! Print "name value" with six decimals. The f0.6 format drops the zero
! before the point (0.5 prints as .500000), so we put it back.
subroutine say(name, x)
character(len=*), intent(in) :: name
real(dp), intent(in) :: x
character(len=40) :: text
write(text, "(f0.6)") x
if (text(1:1) == ".") text = "0" // text
if (text(1:2) == "-.") text = "-0" // text(2:)
print "(a, 1x, a)", name, trim(text)
end subroutine say
end program hypothesis_testWe never shuffle at random. The array chosen holds the positions of the
8 relabeled plants, starting at 1 to 8. Each pass finds the rightmost
position that can still move up, adds one to it, and resets the positions
after it, so the loop visits all 12,870 sets once, in dictionary order,
with no recursion.
plants(chosen) is a vector subscript, the elements at those positions,
so sum(plants(chosen)) is the relabeled total. Totals are integers, so
the 6 exact ties all count.
// Exact two-sided permutation test for a difference in means, using only the
// standard library.
//
// If compost made no difference, the "compost" and "control" labels are
// arbitrary, and any 8 of the 16 plants could have carried the compost label.
// We try every one of those splits and count how many give a difference at
// least as large as the one we observed. That fraction is the p-value.
package main
import (
"encoding/csv"
"fmt"
"os"
"strconv"
)
func read(path string) (control, compost []int, err error) {
f, err := os.Open(path)
if err != nil {
return nil, nil, err
}
defer f.Close()
records, err := csv.NewReader(f).ReadAll()
if err != nil {
return nil, nil, err
}
for _, record := range records[1:] { // skip the header
grams, err := strconv.Atoi(record[1])
if err != nil {
return nil, nil, err
}
switch record[0] {
case "control":
control = append(control, grams)
case "compost":
compost = append(compost, grams)
}
}
return control, compost, nil
}
func sum(xs []int) int {
total := 0
for _, x := range xs {
total += x
}
return total
}
// abs for ints; math.Abs takes and returns float64.
func abs(x int) int {
if x < 0 {
return -x
}
return x
}
func main() {
control, compost, err := read(os.Args[1])
if err != nil {
fmt.Fprintln(os.Stderr, err)
os.Exit(1)
}
n := len(compost)
plants := append(append([]int{}, control...), compost...)
total := sum(plants)
// Both groups have n plants, so comparing totals is the same as comparing
// means, and totals stay whole numbers: no floating-point near-ties.
observed := sum(compost) - sum(control)
splits, asExtreme := 0, 0
// choose walks every way to pick `left` more plants from positions start
// onward, carrying the running total of the plants picked so far. A
// recursive closure must be declared before it is assigned.
var choose func(start, left, chosenTotal int)
choose = func(start, left, chosenTotal int) {
if left == 0 {
difference := chosenTotal - (total - chosenTotal)
splits++
if abs(difference) >= abs(observed) { // two-sided
asExtreme++
}
return
}
// Stop early enough that `left` plants still remain to pick.
for i := start; i <= len(plants)-left; i++ {
choose(i+1, left-1, chosenTotal+plants[i])
}
}
choose(0, n, 0)
fmt.Printf("control_count %d\n", len(control))
fmt.Printf("compost_count %d\n", n)
fmt.Printf("control_mean %.6f\n", float64(sum(control))/float64(len(control)))
fmt.Printf("compost_mean %.6f\n", float64(sum(compost))/float64(n))
fmt.Printf("observed_difference %.6f\n", float64(observed)/float64(n))
fmt.Printf("splits %d\n", splits)
fmt.Printf("as_extreme %d\n", asExtreme)
fmt.Printf("p_value %.6f\n", float64(asExtreme)/float64(splits))
}We never shuffle at random. A recursive closure, choose, walks every way
to pick 8 of the 16 plants for the compost label, carrying a running total,
and counts all 12,870 splits. A closure can’t call itself if it’s declared
with :=, because the name doesn’t exist yet on the right-hand side, so we
write var choose func(...) first and assign it after.
The totals are int, so >= never meets a near-tie. math.Abs takes and
returns float64, so a three-line abs for ints keeps the comparison in
whole numbers.
-- Exact two-sided permutation test for a difference in means, using only
-- Haskell's base library.
--
-- If compost made no difference, the "compost" and "control" labels are
-- arbitrary, and any 8 of the 16 plants could have carried the compost
-- label. We try every one of those splits and count how many give a
-- difference at least as large as the one we observed. That fraction is the
-- p-value.
import Data.List (elemIndex)
import System.Environment (getArgs)
import Text.Printf (printf)
-- Split a line at every comma. break returns the text before the first comma
-- and the rest, which still starts with the comma.
splitOnComma :: String -> [String]
splitOnComma line = case break (== ',') line of
(field, ',' : rest) -> field : splitOnComma rest
(field, _) -> [field]
-- The position of a column in the header, or a stop with a message.
columnIndex :: [String] -> String -> Int
columnIndex header name = case elemIndex name header of
Just index -> index
Nothing -> error ("no column " ++ name)
-- Every row as (group, grams), in file order.
readPlants :: String -> [(String, Int)]
readPlants text = case lines text of
[] -> error "empty file"
headerLine : rows ->
[ (field "group", read (field "grams"))
| row <- rows
, not (null row)
, let fields = splitOnComma row
field name = fields !! columnIndex header name
]
where
header = splitOnComma headerLine
-- The total of every way to choose k of the plants. Each plant is either in
-- the chosen group (add it and choose k - 1 from the rest) or not (choose
-- all k from the rest). The list has one total per split.
chosenTotals :: Int -> [Int] -> [Int]
chosenTotals 0 _ = [0]
chosenTotals _ [] = []
chosenTotals k (x : xs) = map (+ x) (chosenTotals (k - 1) xs) ++ chosenTotals k xs
main :: IO ()
main = do
[path] <- getArgs
text <- readFile path
let plants = readPlants text
control = [grams | ("control", grams) <- plants]
compost = [grams | ("compost", grams) <- plants]
n = length compost
total = sum control + sum compost
-- Both groups have n plants, so comparing totals is the same as
-- comparing means, and Int totals stay whole numbers: no
-- floating-point near-ties.
observed = sum compost - sum control
differences = [chosen - (total - chosen) | chosen <- chosenTotals n (control ++ compost)]
splits = length differences
-- Two-sided: a difference counts in either direction.
asExtreme = length (filter (\d -> abs d >= abs observed) differences)
-- fromIntegral turns an Int into a Double for division.
mean xs = fromIntegral (sum xs) / fromIntegral (length xs) :: Double
printf "control_count %d\n" (length control)
printf "compost_count %d\n" n
printf "control_mean %.6f\n" (mean control)
printf "compost_mean %.6f\n" (mean compost)
printf "observed_difference %.6f\n" (fromIntegral observed / fromIntegral n :: Double)
printf "splits %d\n" splits
printf "as_extreme %d\n" asExtreme
printf "p_value %.6f\n" (fromIntegral asExtreme / fromIntegral splits :: Double)We never shuffle at random. A four-line recursive function, chosenTotals,
lists the total of every way to choose 8 of the 16 plants: each plant is
either in the chosen group or not. That gives exactly 12,870 totals, so the
p-value is exact.
The totals are Int, so abs d >= abs observed compares whole numbers
and the 6 exact ties all count. Data.List has subsequences, but it
lists all 65,536 subsets of the 16 plants, so we’d still have to filter
by size.
// Exact two-sided permutation test for a difference in means, using only the
// Java standard library.
//
// If compost made no difference, the "compost" and "control" labels are
// arbitrary, and any 8 of the 16 plants could have carried the compost label.
// We try every one of those splits and count how many give a difference at
// least as large as the one we observed. That fraction is the p-value.
import java.io.IOException;
import java.nio.file.Files;
import java.nio.file.Path;
import java.util.ArrayList;
import java.util.List;
import java.util.Locale;
class Main {
// Reading a file can fail with an IOException, and Java makes us say so
// with "throws".
public static void main(String[] args) throws IOException {
List<String> lines = Files.readAllLines(Path.of(args[0]));
List<Integer> control = new ArrayList<>(), compost = new ArrayList<>();
for (String line : lines.subList(1, lines.size())) { // skip the header
if (line.isBlank()) {
continue;
}
String[] fields = line.split(",");
// Strings are compared with equals; == would compare objects.
(fields[0].equals("compost") ? compost : control).add(Integer.parseInt(fields[1]));
}
int n = compost.size();
if (n == 0 || control.size() != n) {
throw new IllegalArgumentException("need two groups of equal size");
}
List<Integer> plants = new ArrayList<>(control);
plants.addAll(compost);
long controlSum = 0, compostSum = 0;
for (int grams : control) controlSum += grams;
for (int grams : compost) compostSum += grams;
long total = controlSum + compostSum;
// Both groups have the same size, so comparing totals is the same as
// comparing means, and totals stay whole numbers: no floating-point
// near-ties.
long observed = compostSum - controlSum;
// Each whole number from 0 to 2^16 - 1 is a mask: bit i set means
// plant i is labeled compost. Integer.bitCount counts the set bits,
// so the masks with exactly n of them are the splits.
int splits = 0, asExtreme = 0;
for (int mask = 0; mask < 1 << plants.size(); mask++) {
if (Integer.bitCount(mask) != n) {
continue;
}
long chosen = 0;
for (int i = 0; i < plants.size(); i++) {
if ((mask & (1 << i)) != 0) {
chosen += plants.get(i);
}
}
long difference = chosen - (total - chosen);
splits++;
if (Math.abs(difference) >= Math.abs(observed)) { // two-sided
asExtreme++;
}
}
Locale root = Locale.ROOT;
System.out.printf(root, "control_count %d%n", control.size());
System.out.printf(root, "compost_count %d%n", n);
System.out.printf(root, "control_mean %.6f%n", (double) controlSum / n);
System.out.printf(root, "compost_mean %.6f%n", (double) compostSum / n);
System.out.printf(root, "observed_difference %.6f%n", (double) observed / n);
System.out.printf(root, "splits %d%n", splits);
System.out.printf(root, "as_extreme %d%n", asExtreme);
System.out.printf(root, "p_value %.6f%n", (double) asExtreme / splits);
}
}We treat every whole number from 0 to 2^16 − 1 as a mask: bit i set means
plant i is labeled compost. Integer.bitCount counts the set bits, so the
12,870 masks with exactly 8 of them are the splits, and the other 52,666
are skipped. mask & (1 << i) tests one bit at a time to add up the chosen
grams.
The sums are long integers, so Math.abs and >= compare whole numbers
and the 6 splits that tie the observed difference all count. We divide only
when printing. Apache
Commons Math has statistical tests, but it isn’t installed.
// Exact two-sided permutation test for a difference in means, using only
// Node's built-in modules.
//
// If compost made no difference, the "compost" and "control" labels are
// arbitrary, and any 8 of the 16 plants could have carried the compost label.
// We try every one of those splits and count how many give a difference at
// least as large as the one we observed. That fraction is the p-value.
import { readFileSync } from "node:fs";
// Node has no CSV parser. No field here is quoted, so we split on commas.
function read(path) {
const groups = { control: [], compost: [] };
const [, ...lines] = readFileSync(path, "utf8").trim().split("\n");
for (const line of lines) {
const [group, grams] = line.split(",");
groups[group].push(Number(grams));
}
return [groups.control, groups.compost];
}
const sum = (xs) => xs.reduce((acc, x) => acc + x, 0);
const [control, compost] = read(process.argv[2]);
const n = compost.length;
const plants = [...control, ...compost];
const total = sum(plants);
// Both groups have n plants, so comparing totals is the same as comparing
// means, and totals stay whole numbers: no floating-point near-ties.
const observed = sum(compost) - sum(control);
let splits = 0;
let asExtreme = 0;
// Walk every way to pick `left` more plants from positions start onward,
// carrying the running total of the plants picked so far.
function choose(start, left, chosenTotal) {
if (left === 0) {
const difference = chosenTotal - (total - chosenTotal);
splits++;
if (Math.abs(difference) >= Math.abs(observed)) asExtreme++; // two-sided
return;
}
// Stop early enough that `left` plants still remain to pick.
for (let i = start; i <= plants.length - left; i++) {
choose(i + 1, left - 1, chosenTotal + plants[i]);
}
}
choose(0, n, 0);
console.log(`control_count ${control.length}`);
console.log(`compost_count ${n}`);
console.log(`control_mean ${(sum(control) / control.length).toFixed(6)}`);
console.log(`compost_mean ${(sum(compost) / n).toFixed(6)}`);
console.log(`observed_difference ${(observed / n).toFixed(6)}`);
console.log(`splits ${splits}`);
console.log(`as_extreme ${asExtreme}`);
console.log(`p_value ${(asExtreme / splits).toFixed(6)}`);We never shuffle at random. A recursive function, choose, walks every way
to pick 8 of the 16 plants for the compost label, carrying a running total
as it goes, so we add up each split once and visit all 12,870 of them. Its
loop stops while enough plants remain to fill the group, so no branch is a
dead end.
Because both groups have 8 plants, comparing totals is the same as comparing
means, and the totals stay whole numbers, so >= never meets a near-tie.
# Exact two-sided permutation test for a difference in means, using only
# Base Julia.
#
# If compost made no difference, the "compost" and "control" labels are
# arbitrary, and any 8 of the 16 plants could have carried the compost label.
# We try every one of those splits and count how many give a difference at
# least as large as the one we observed. That fraction is the p-value.
using Printf
function read_groups(path)
groups = Dict("control" => Int[], "compost" => Int[])
for line in Iterators.drop(eachline(path), 1) # skip the header
group, grams = split(line, ",")
push!(groups[group], parse(Int, grams))
end
return groups["control"], groups["compost"]
end
function main(path)
control, compost = read_groups(path)
n = length(compost)
plants = vcat(control, compost)
total = sum(plants)
# Both groups have n plants, so comparing totals is the same as comparing
# means, and totals stay whole numbers: no floating-point near-ties.
observed = sum(compost) - sum(control)
# Each split is a number whose binary digits mark the chosen plants: bit
# i - 1 set means plant i carries the compost label. Of the 2^16 numbers,
# the ones with exactly n bits set are the splits.
splits = 0
as_extreme = 0
for mask in 0:(1 << length(plants)) - 1
count_ones(mask) == n || continue
relabeled = sum(plants[i] for i in eachindex(plants) if (mask >> (i - 1)) & 1 == 1)
difference = relabeled - (total - relabeled)
splits += 1
if abs(difference) >= abs(observed) # two-sided
as_extreme += 1
end
end
println("control_count $(length(control))")
println("compost_count $n")
@printf("control_mean %.6f\n", sum(control) / length(control))
@printf("compost_mean %.6f\n", sum(compost) / n)
@printf("observed_difference %.6f\n", observed / n)
println("splits $splits")
println("as_extreme $as_extreme")
@printf("p_value %.6f\n", as_extreme / splits)
end
main(ARGS[1])We never shuffle at random. Each split is a number from 0 to 65,535 whose
binary digits mark the chosen plants, and count_ones(mask) == n keeps the
12,870 numbers with exactly 8 bits set. Julia arrays start at 1, so plant
i is bit i - 1, read with (mask >> (i - 1)) & 1.
Because both groups have 8 plants, comparing totals is the same as comparing
means, and the totals stay Int, so >= never meets a near-tie. Base Julia
has no combinations function or permutation test. The Combinatorics and
HypothesisTests packages have them, but neither is installed.
-- Exact two-sided permutation test for a difference in means, using only
-- Lua's standard library.
--
-- If compost made no difference, the "compost" and "control" labels are
-- arbitrary, and any n of the plants could have carried the compost label.
-- We try every one of those splits and count how many give a difference at
-- least as large as the one we observed. That fraction is the p-value.
-- Returns a table from group name to the list of its grams.
local function read_groups(path)
local f = assert(io.open(path, "r"))
local col, i = {}, 0
for field in f:read("l"):gmatch("[^,]+") do
i = i + 1
col[field] = i -- column name -> position
end
local groups = {control = {}, compost = {}}
for line in f:lines() do
local fields = {}
for field in line:gmatch("[^,]+") do
fields[#fields + 1] = field
end
local list = groups[fields[col.group]]
-- "459" becomes the integer 459, not the float 459.0.
list[#list + 1] = tonumber(fields[col.grams])
end
f:close()
return groups.control, groups.compost
end
local function sum(values)
local total = 0
for _, v in ipairs(values) do total = total + v end
return total
end
local control, compost = read_groups(arg[1])
local n = #compost
local plants = {}
for _, g in ipairs(control) do plants[#plants + 1] = g end
for _, g in ipairs(compost) do plants[#plants + 1] = g end
local total = sum(plants)
-- Both groups have n plants, so comparing totals is the same as comparing
-- means, and totals stay whole numbers: no floating-point near-ties.
local observed = sum(compost) - sum(control)
local splits, as_extreme = 0, 0
-- Walk the plants from position next on, choosing left more of them to
-- carry the compost label; chosen is the total of those picked so far. Each
-- call either takes plant next or skips it, so every choice is reached once.
-- "local function" lets choose call itself and update the counters above.
local function choose(next, left, chosen)
if left == 0 then
splits = splits + 1
-- Two-sided: math.abs counts a difference in either direction.
if math.abs(chosen - (total - chosen)) >= math.abs(observed) then
as_extreme = as_extreme + 1
end
return
end
if #plants - next + 1 < left then return end -- too few plants left
choose(next + 1, left - 1, chosen + plants[next])
choose(next + 1, left, chosen)
end
choose(1, n, 0) -- Lua arrays start at 1
-- / always divides as a decimal, even between two integers.
print("control_count " .. #control)
print("compost_count " .. n)
print(string.format("control_mean %.6f", sum(control) / #control))
print(string.format("compost_mean %.6f", sum(compost) / n))
print(string.format("observed_difference %.6f", observed / n))
print("splits " .. splits)
print("as_extreme " .. as_extreme)
print(string.format("p_value %.6f", as_extreme / splits))A recursive function, choose, walks the plants in order. At each plant it
calls itself twice: once taking the plant into the compost group and once
skipping it. When it has taken 8, it counts the split. It’s declared with
local function, which lets it call itself; local choose = function
would fail, because the name isn’t in scope yet inside the body.
tonumber("459") gives an integer, so the totals and the >= test stay
exact. Lua’s / always returns a float, even for 8 / 4, which is what the
means want. LuaRocks isn’t installed.
(* Exact two-sided permutation test for a difference in means, using only
OCaml's standard library.
If compost made no difference, the "compost" and "control" labels are
arbitrary, and any 8 of the 16 plants could have carried the compost
label. We try every one of those splits and count how many give a
difference at least as large as the one we observed. That fraction is
the p-value. *)
(* Every row as (group, grams), with each column found by its header name. *)
let read_plants path =
In_channel.with_open_text path (fun ic ->
let lines = In_channel.input_lines ic in
let header = String.split_on_char ',' (List.hd lines) in
let index name =
match List.find_index (String.equal name) header with
| Some i -> i
| None -> failwith ("column not found: " ^ name)
in
List.filter_map
(fun line ->
if line = "" then None
else
let fields = String.split_on_char ',' line in
let field name = List.nth fields (index name) in
Some (field "group", int_of_string (field "grams")))
(List.tl lines))
let sum = List.fold_left ( + ) 0
(* Call visit with the total of every way to choose k of the plants from
position i onward. Each plant is either chosen (add it, and choose k - 1
from the rest) or not (choose all k from the rest). *)
let rec each_chosen_total plants i k total visit =
if k = 0 then visit total
else if Array.length plants - i >= k then begin
each_chosen_total plants (i + 1) (k - 1) (total + plants.(i)) visit;
each_chosen_total plants (i + 1) k total visit
end
let () =
let rows = read_plants Sys.argv.(1) in
let group name = List.filter_map (fun (g, x) -> if g = name then Some x else None) rows in
let control = group "control" and compost = group "compost" in
let n = List.length compost in
let plants = Array.of_list (control @ compost) in
let total = sum control + sum compost in
(* Both groups have n plants, so comparing totals is the same as comparing
means, and int totals stay whole numbers: no floating-point near-ties. *)
let observed = sum compost - sum control in
(* ref makes a mutable cell; ! reads it and := sets it. *)
let splits = ref 0 and as_extreme = ref 0 in
each_chosen_total plants 0 n 0 (fun chosen ->
let difference = chosen - (total - chosen) in
incr splits;
(* Two-sided: a difference counts in either direction. *)
if abs difference >= abs observed then incr as_extreme);
let mean xs = float_of_int (sum xs) /. float_of_int (List.length xs) in
Printf.printf "control_count %d\n" (List.length control);
Printf.printf "compost_count %d\n" n;
Printf.printf "control_mean %.6f\n" (mean control);
Printf.printf "compost_mean %.6f\n" (mean compost);
Printf.printf "observed_difference %.6f\n"
(float_of_int observed /. float_of_int n);
Printf.printf "splits %d\n" !splits;
Printf.printf "as_extreme %d\n" !as_extreme;
Printf.printf "p_value %.6f\n"
(float_of_int !as_extreme /. float_of_int !splits)We never shuffle at random. OCaml 5.3’s Array.shuffle would give a
random relabeling, but a sample only estimates the p-value, and with 8
plants in each group we can try all 12,870 splits. The recursive
each_chosen_total walks the plants in order, either choosing each one or
skipping it, and calls a function with every finished total.
Two ref cells count the splits and the extreme ones. The totals are
int, so abs difference >= abs observed compares whole numbers, and the
6 exact ties all count.
% Exact two-sided permutation test for a difference in means, using only
% core GNU Octave.
%
% If compost made no difference, the "compost" and "control" labels are
% arbitrary, and any n of the plants could have carried the compost label.
% We try every one of those splits and count how many give a difference at
% least as large as the one we observed. That fraction is the p-value.
% The named column of a CSV file, as a column cell array of strings. Wrap it
% in str2double to get numbers.
function values = read_column(path, name)
lines = strsplit(strtrim(fileread(path)), "\n");
header = strsplit(strtrim(lines{1}), ",");
index = find(strcmp(header, name));
if isempty(index)
error("no column %s", name);
end
values = cell(numel(lines) - 1, 1);
for row = 2:numel(lines)
fields = strsplit(strtrim(lines{row}), ",");
values{row - 1} = fields{index};
end
end
args = argv();
group = read_column(args{1}, "group");
grams = str2double(read_column(args{1}, "grams"));
% strcmp compares a whole cell array with one string, giving a logical mask.
control = grams(strcmp(group, "control"));
compost = grams(strcmp(group, "compost"));
n = numel(compost);
plants = [control; compost];
total = sum(plants);
% Both groups have n plants, so comparing totals is the same as comparing
% means, and totals stay whole numbers: no floating-point near-ties.
observed = sum(compost) - sum(control);
% nchoosek with a vector returns every choice of n positions, one per row.
% Indexing plants with that matrix gives a matrix of grams the same shape,
% and sum(..., 2) totals each row: one relabeled compost total per split.
chosen = nchoosek(1:numel(plants), n);
relabeled = sum(plants(chosen), 2);
difference = relabeled - (total - relabeled);
splits = rows(chosen);
as_extreme = sum(abs(difference) >= abs(observed)); % two-sided
printf("control_count %d\n", numel(control));
printf("compost_count %d\n", n);
printf("control_mean %.6f\n", mean(control));
printf("compost_mean %.6f\n", mean(compost));
printf("observed_difference %.6f\n", observed / n);
printf("splits %d\n", splits);
printf("as_extreme %d\n", as_extreme);
printf("p_value %.6f\n", as_extreme / splits);There’s no loop over the splits. nchoosek(1:16, 8) returns all 12,870
choices of 8 plant positions as rows of a 12,870-by-8 matrix, in dictionary
order. Indexing the grams column with that matrix gives a matrix of grams
the same shape, and sum(..., 2) totals each row, so one line yields every
relabeled compost total.
We compare whole-number totals rather than means, so >= never meets a
floating-point near-tie. nchoosek ships with core Octave and only lists
the choices; the counting is ours. Octave Forge packages, such as
statistics, aren’t installed.
{ Exact two-sided permutation test for a difference in means, using only
Free Pascal's run-time library.
If compost made no difference, the "compost" and "control" labels are
arbitrary, and any N of the plants could have carried the compost label.
We try every one of those splits and count how many give a difference at
least as large as the one we observed. That fraction is the p-value. }
program HypothesisTest;
{$mode objfpc}{$H+} // 32-bit Integer, and strings longer than 255 characters
uses
SysUtils;
type
TIntArray = array of Integer;
var
Control, Compost, Plants: TIntArray;
N, Total, Observed, Splits, AsExtreme: Integer;
// Read the grams of each group into Control and Compost.
procedure ReadPlants(const Path: String);
var
F: TextFile;
Line: String;
Fields: TStringArray;
GroupCol, GramsCol, I: Integer;
begin
AssignFile(F, Path);
Reset(F);
ReadLn(F, Line);
Fields := Line.Split(',');
GroupCol := -1;
GramsCol := -1;
for I := 0 to High(Fields) do
if Fields[I] = 'group' then GroupCol := I
else if Fields[I] = 'grams' then GramsCol := I;
if (GroupCol < 0) or (GramsCol < 0) then
raise Exception.Create('column not found');
while not Eof(F) do
begin
ReadLn(F, Line);
Fields := Line.Split(',');
// Insert(Value, Array, Length(Array)) appends.
if Fields[GroupCol] = 'control' then
Insert(StrToInt(Fields[GramsCol]), Control, Length(Control))
else
Insert(StrToInt(Fields[GramsCol]), Compost, Length(Compost));
end;
CloseFile(F);
end;
function Sum(const Values: TIntArray): Integer;
var
V: Integer;
begin
Result := 0;
for V in Values do
Result := Result + V;
end;
// Walk the plants from position Next on, choosing Left more of them to carry
// the compost label. Chosen is the total of those picked so far. Each call
// either takes plant Next or skips it, so every choice is reached once.
procedure Choose(Next, Left, Chosen: Integer);
begin
if Left = 0 then
begin
Inc(Splits);
// Two-sided: Abs counts a difference in either direction.
if Abs(Chosen - (Total - Chosen)) >= Abs(Observed) then
Inc(AsExtreme);
Exit;
end;
if Length(Plants) - Next < Left then
Exit; // too few plants left to finish this choice
Choose(Next + 1, Left - 1, Chosen + Plants[Next]);
Choose(Next + 1, Left, Chosen);
end;
begin
Control := nil;
Compost := nil;
ReadPlants(ParamStr(1));
N := Length(Compost);
Plants := Concat(Control, Compost);
Total := Sum(Plants);
// Both groups have N plants, so comparing totals is the same as comparing
// means, and totals stay whole numbers: no floating-point near-ties.
Observed := Sum(Compost) - Sum(Control);
Splits := 0;
AsExtreme := 0;
Choose(0, N, 0);
WriteLn('control_count ', Length(Control));
WriteLn('compost_count ', N);
WriteLn('control_mean ', Sum(Control) / Length(Control):0:6);
WriteLn('compost_mean ', Sum(Compost) / N:0:6);
WriteLn('observed_difference ', Observed / N:0:6);
WriteLn('splits ', Splits);
WriteLn('as_extreme ', AsExtreme);
WriteLn('p_value ', AsExtreme / Splits:0:6);
end.A recursive procedure, Choose, walks the plants in order. At each plant
it calls itself twice: once taking the plant into the compost group and
once skipping it. When it has taken 8, it counts the split. A branch with
too few plants left to reach 8 returns early, so the walk visits exactly
the 12,870 splits.
The totals are Integer, and so is the comparison, so >= never meets a
floating-point near-tie. Concat joins the two groups’ dynamic arrays
into one list of plants, and Abs makes the test two-sided.
# Exact two-sided permutation test for a difference in means, using only
# Perl and its core modules.
#
# If compost made no difference, the "compost" and "control" labels are
# arbitrary, and any n of the plants could have carried the compost label.
# We try every one of those splits and count how many give a difference at
# least as large as the one we observed. That fraction is the p-value.
use strict;
use warnings;
# Returns a hash reference from group name to an array reference of grams.
sub read_groups {
my ($path) = @_;
open my $fh, '<', $path or die "can't open $path: $!";
chomp(my $header = <$fh>);
my @fields = split /,/, $header;
my %col = map { $fields[$_] => $_ } 0 .. $#fields; # name -> position
my %groups = (control => [], compost => []);
while (my $line = <$fh>) {
chomp $line;
my @row = split /,/, $line;
push @{ $groups{ $row[ $col{group} ] } }, $row[ $col{grams} ];
}
close $fh;
return \%groups;
}
sub sum {
my $total = 0;
$total += $_ for @_;
return $total;
}
my $groups = read_groups($ARGV[0]);
my @control = @{ $groups->{control} };
my @compost = @{ $groups->{compost} };
my $n = @compost; # an array in scalar context is its length
my @plants = (@control, @compost);
my $total = sum(@plants);
# Both groups have n plants, so comparing totals is the same as comparing
# means, and totals stay whole numbers: no floating-point near-ties.
my $observed = sum(@compost) - sum(@control);
# With N plants, every subset is a whole number from 0 to 2**N - 1 whose
# bits say which plants it holds: bit i set means plant i is chosen. We keep
# the numbers with exactly n bits set. sprintf '%b' writes the bits as a
# string of 0s and 1s, and tr/1// counts the 1s without changing anything.
my ($splits, $as_extreme) = (0, 0);
for my $mask (0 .. 2**@plants - 1) {
next unless sprintf('%b', $mask) =~ tr/1// == $n;
my $relabeled = sum(map { $plants[$_] } grep { ($mask >> $_) & 1 } 0 .. $#plants);
$splits++;
# Two-sided: abs counts a difference in either direction.
$as_extreme++ if abs($relabeled - ($total - $relabeled)) >= abs($observed);
}
print 'control_count ', scalar @control, "\n";
print "compost_count $n\n";
printf "control_mean %.6f\n", sum(@control) / @control;
printf "compost_mean %.6f\n", sum(@compost) / $n;
printf "observed_difference %.6f\n", $observed / $n;
print "splits $splits\n";
print "as_extreme $as_extreme\n";
printf "p_value %.6f\n", $as_extreme / $splits;We number the subsets instead of recursing. With 16 plants, each whole
number from 0 to 65,535 is a subset: bit i set means plant i is chosen.
sprintf('%b', $mask) writes the bits as text, and tr/1// counts the
1s, so we keep the 12,870 masks with exactly 8 bits set and total those
plants’ grams.
That checks all 65,536 masks to use a fifth of them, which is still
instant here. Totals are whole numbers, so >= is exact. CPAN modules,
such as Algorithm::Combinatorics, aren’t installed.
<?php
// Exact two-sided permutation test for a difference in means, using only
// PHP and the extensions built into it.
//
// If compost made no difference, the "compost" and "control" labels are
// arbitrary, and any n of the plants could have carried the compost label.
// We try every one of those splits and count how many give a difference at
// least as large as the one we observed. That fraction is the p-value.
declare(strict_types=1);
// Returns [control grams, compost grams] as two lists of integers.
function read_groups(string $path): array
{
$lines = file($path, FILE_IGNORE_NEW_LINES | FILE_SKIP_EMPTY_LINES);
$col = array_flip(explode(',', array_shift($lines))); // name => position
$groups = ['control' => [], 'compost' => []];
foreach ($lines as $line) {
$row = explode(',', $line);
$groups[$row[$col['group']]][] = (int) $row[$col['grams']];
}
return [$groups['control'], $groups['compost']];
}
// Yields the compost total of every way to choose $left more plants from
// position $next on; $chosen is the total of those picked so far. Each step
// either takes plant $next or skips it, so every choice is reached once. A
// function with yield is a generator: foreach pulls one value at a time,
// and yield from passes along everything a nested call yields.
function chosen_totals(array $plants, int $next, int $left, int $chosen): Generator
{
if ($left === 0) {
yield $chosen;
return;
}
if (count($plants) - $next < $left) {
return; // too few plants left to finish this choice
}
yield from chosen_totals($plants, $next + 1, $left - 1, $chosen + $plants[$next]);
yield from chosen_totals($plants, $next + 1, $left, $chosen);
}
[$control, $compost] = read_groups($argv[1]);
$n = count($compost);
$plants = array_merge($control, $compost);
$total = array_sum($plants);
// Both groups have n plants, so comparing totals is the same as comparing
// means, and totals stay whole numbers: no floating-point near-ties.
$observed = array_sum($compost) - array_sum($control);
$splits = 0;
$as_extreme = 0;
foreach (chosen_totals($plants, 0, $n, 0) as $relabeled) {
$splits++;
// Two-sided: abs counts a difference in either direction.
if (abs($relabeled - ($total - $relabeled)) >= abs($observed)) {
$as_extreme++;
}
}
echo 'control_count ', count($control), "\n";
echo "compost_count $n\n";
printf("control_mean %.6f\n", array_sum($control) / count($control));
printf("compost_mean %.6f\n", array_sum($compost) / $n);
printf("observed_difference %.6f\n", $observed / $n);
echo "splits $splits\n";
echo "as_extreme $as_extreme\n";
printf("p_value %.6f\n", $as_extreme / $splits);A recursive generator, chosen_totals, walks the plants in order and
yields the compost total of each split. At each plant it either takes the
plant or skips it, and yield from passes along everything the nested
call yields. The main loop uses foreach to pull the 12,870 totals one at
a time, so no list of splits is ever built.
The grams are cast with (int), so totals and the >= test stay exact.
PHP’s / returns an integer when the division is exact, as in 8 / 4, and
a float otherwise. Composer packages aren’t installed.
-- Exact two-sided permutation test for a difference in means, in PostgreSQL.
--
-- If compost made no difference, the "compost" and "control" labels are
-- arbitrary, and any 8 of the 16 plants could have carried the compost label.
-- We try every one of those splits and count how many give a difference at
-- least as large as the one we observed. That fraction is the p-value.
WITH RECURSIVE plants AS (
-- run.sh loads grams as double precision; they are whole numbers, so we
-- turn them back into integers and every total below stays exact.
-- "group" is a reserved word in SQL, so its column name needs quotes.
SELECT row_id AS id, "group" AS label, grams::integer AS grams FROM data
),
summary AS (
SELECT count(*) FILTER (WHERE label = 'control') AS control_count,
count(*) FILTER (WHERE label = 'compost') AS compost_count,
sum(grams) FILTER (WHERE label = 'control') AS control_sum,
sum(grams) FILTER (WHERE label = 'compost') AS compost_sum,
sum(grams) AS total
FROM plants
),
picks AS (
-- A recursive CTE grows its rows step by step. Each row is a partial
-- split: how many plants it has chosen, their grams, and the last id
-- chosen. A step adds one plant with a larger id, so every set of ids is
-- built once, in increasing order, and never twice in another order.
SELECT id AS last_id, 1 AS size, grams AS chosen FROM plants
UNION ALL
SELECT plants.id, picks.size + 1, picks.chosen + plants.grams
FROM picks
JOIN plants ON plants.id > picks.last_id
CROSS JOIN summary
WHERE picks.size < summary.compost_count
),
splits AS (
-- The finished splits have as many plants as the compost group. Both
-- groups have the same size, so comparing totals is the same as
-- comparing means, and integers have no floating-point near-ties.
SELECT picks.chosen - (summary.total - picks.chosen) AS difference
FROM picks CROSS JOIN summary
WHERE picks.size = summary.compost_count
),
test AS (
SELECT count(*) AS splits,
count(*) FILTER (WHERE abs(difference) >= -- two-sided
(SELECT abs(compost_sum - control_sum) FROM summary))
AS as_extreme
FROM splits
)
SELECT name, value
FROM (
-- Dividing numeric values and rounding to 6 places prints 6 decimals.
SELECT 1 AS position, 'control_count' AS name, control_count::text AS value FROM summary
UNION ALL SELECT 2, 'compost_count', compost_count::text FROM summary
UNION ALL SELECT 3, 'control_mean', round(control_sum::numeric / control_count, 6)::text FROM summary
UNION ALL SELECT 4, 'compost_mean', round(compost_sum::numeric / compost_count, 6)::text FROM summary
UNION ALL SELECT 5, 'observed_difference',
round((compost_sum - control_sum)::numeric / compost_count, 6)::text FROM summary
UNION ALL SELECT 6, 'splits', splits::text FROM test
UNION ALL SELECT 7, 'as_extreme', as_extreme::text FROM test
UNION ALL SELECT 8, 'p_value', round(as_extreme::numeric / splits, 6)::text FROM test
) AS results
ORDER BY position;SQL has no loops, so a recursive CTE builds the splits. Each row of picks
is a partial split: how many plants it holds, their total grams, and the
last plant id. Each step joins on plants.id > picks.last_id, so a set of
plants is built in one order only. The CTE makes 39,202 partial rows, and
the 12,870 with 8 plants are the splits.
run.sh loads grams as double precision, so we cast them to integer and
every total is exact. group is a reserved word; unquoted, it’s a syntax
error.
% Exact two-sided permutation test for a difference in means, using
% SWI-Prolog and the libraries that come with it.
%
% If compost made no difference, the "compost" and "control" labels are
% arbitrary, and any 8 of the 16 plants could have carried the compost
% label. We try every one of those splits and count how many give a
% difference at least as large as the one we observed. That fraction is
% the p-value.
:- initialization(main, main).
% read_plants(+Path, -Control, -Compost): the grams of each group, in file
% order, with each column found by its header name.
read_plants(Path, Control, Compost) :-
read_file_to_string(Path, Text, []),
split_string(Text, "\n", "", [HeaderLine|Lines]),
split_string(HeaderLine, ",", "", Header),
nth0(G, Header, "group"),
nth0(X, Header, "grams"),
findall(Group-Grams,
( member(Line, Lines),
Line \== "",
split_string(Line, ",", "", Fields),
nth0(G, Fields, Group),
nth0(X, Fields, GramsText),
number_string(Grams, GramsText) ),
Rows),
findall(Grams, member("control"-Grams, Rows), Control),
findall(Grams, member("compost"-Grams, Rows), Compost).
% chosen_total(+K, +Plants, -Total): Total is the grams of K of the Plants.
% On backtracking it gives every way to choose them, once each: the first
% plant is either chosen (add it, and choose K - 1 from the rest) or not
% (choose all K from the rest).
chosen_total(0, _, 0).
chosen_total(K, [Plant|Rest], Total) :-
K > 0,
K1 is K - 1,
chosen_total(K1, Rest, Total0),
Total is Total0 + Plant.
chosen_total(K, [_|Rest], Total) :-
K > 0,
chosen_total(K, Rest, Total).
% difference(+N, +Plants, +All, -Difference): one split's chosen total
% minus the rest, for each split in turn.
difference(N, Plants, All, Difference) :-
chosen_total(N, Plants, Chosen),
Difference is Chosen - (All - Chosen).
main :-
current_prolog_flag(argv, [Path|_]),
read_plants(Path, Control, Compost),
length(Control, ControlCount),
length(Compost, N),
append(Control, Compost, Plants),
sum_list(Control, ControlTotal),
sum_list(Compost, CompostTotal),
All is ControlTotal + CompostTotal,
% Both groups have N plants, so comparing totals is the same as
% comparing means, and integer totals stay whole numbers: no
% floating-point near-ties.
Observed is CompostTotal - ControlTotal,
% aggregate_all(count, Goal, Count) counts the ways Goal succeeds.
aggregate_all(count, difference(N, Plants, All, _), Splits),
% Two-sided: a difference counts in either direction.
aggregate_all(count,
( difference(N, Plants, All, D), abs(D) >= abs(Observed) ),
AsExtreme),
% / of two integers gives a float, or an integer when it divides
% exactly (3232 / 8 is 404); ~6f prints either with six decimals.
ControlMean is ControlTotal / ControlCount,
CompostMean is CompostTotal / N,
ObservedMean is Observed / N,
PValue is AsExtreme / Splits,
format("control_count ~d~n", [ControlCount]),
format("compost_count ~d~n", [N]),
format("control_mean ~6f~n", [ControlMean]),
format("compost_mean ~6f~n", [CompostMean]),
format("observed_difference ~6f~n", [ObservedMean]),
format("splits ~d~n", [Splits]),
format("as_extreme ~d~n", [AsExtreme]),
format("p_value ~6f~n", [PValue]).We never shuffle at random. chosen_total/3 has two clauses for each
plant, one that chooses it and one that skips it, so backtracking walks
every way to choose 8 of the 16 plants, 12,870 in all, and the p-value is
exact. aggregate_all(count, ...) counts the splits, and counts again the
ones with abs(D) >= abs(Observed).
The totals are integers, so the 6 exact ties all count. library(lists)
has subseq/3, but it walks all 65,536 subsets, so we’d still filter by
size.
# Exact two-sided permutation test for a difference in means, using only
# base R.
#
# If compost made no difference, the "compost" and "control" labels are
# arbitrary, and any 8 of the 16 plants could have carried the compost label.
# We try every one of those splits and count how many give a difference at
# least as large as the one we observed. That fraction is the p-value.
args <- commandArgs(trailingOnly = TRUE)
data <- read.csv(args[1])
control <- data$grams[data$group == "control"]
compost <- data$grams[data$group == "compost"]
n <- length(compost)
plants <- c(control, compost)
total <- sum(plants)
# Both groups have n plants, so comparing totals is the same as comparing
# means, and totals stay whole numbers: no floating-point near-ties.
observed <- sum(compost) - sum(control)
# combn() lists every choice of n plants, one per column. With FUN = sum it
# returns each choice's total instead, so we get all splits in one vector.
relabeled <- combn(plants, n, FUN = sum)
difference <- relabeled - (total - relabeled)
splits <- length(difference)
as_extreme <- sum(abs(difference) >= abs(observed)) # two-sided
cat(sprintf("control_count %d\n", length(control)))
cat(sprintf("compost_count %d\n", n))
cat(sprintf("control_mean %.6f\n", mean(control)))
cat(sprintf("compost_mean %.6f\n", mean(compost)))
cat(sprintf("observed_difference %.6f\n", observed / n))
cat(sprintf("splits %d\n", splits))
cat(sprintf("as_extreme %d\n", as_extreme))
cat(sprintf("p_value %.6f\n", as_extreme / splits))We never shuffle at random. combn(plants, n) lists every choice of 8 of
the 16 plants, one per column, in an 8 by 12,870 matrix. Passing
FUN = sum gives each choice’s total instead, so one call returns all
12,870 totals as a vector, and the count of extreme splits is a single
vectorized sum() of a comparison. combn() lives in the utils package,
which ships with R and loads by default.
read.csv() reads the grams as integers, so the totals stay whole numbers
and >= never meets a near-tie. The coin package has permutation tests,
but it isn’t installed.
# Exact two-sided permutation test for a difference in means, using only the
# standard library.
#
# If compost made no difference, the "compost" and "control" labels are
# arbitrary, and any 8 of the 16 plants could have carried the compost label.
# We try every one of those splits and count how many give a difference at
# least as large as the one we observed. That fraction is the p-value.
# The csv library is a bundled gem since Ruby 3.4, so split lines by hand
# (no quoted fields here).
def read(path)
groups = { "control" => [], "compost" => [] }
File.readlines(path, chomp: true).drop(1).each do |line|
group, grams = line.split(",")
groups.fetch(group) << Integer(grams)
end
groups.values_at("control", "compost")
end
control, compost = read(ARGV[0])
n = compost.length
plants = control + compost
total = plants.sum
# Both groups have n plants, so comparing totals is the same as comparing
# means, and totals stay whole numbers: no floating-point near-ties.
observed = compost.sum - control.sum
splits = 0
as_extreme = 0
# combination(n) yields each choice of n plants once, in a fixed order.
plants.combination(n) do |chosen|
relabeled = chosen.sum
difference = relabeled - (total - relabeled)
splits += 1
as_extreme += 1 if difference.abs >= observed.abs # two-sided
end
puts "control_count #{control.length}"
puts "compost_count #{n}"
puts format("control_mean %.6f", control.sum.fdiv(control.length))
puts format("compost_mean %.6f", compost.sum.fdiv(n))
puts format("observed_difference %.6f", observed.fdiv(n))
puts "splits #{splits}"
puts "as_extreme #{as_extreme}"
puts format("p_value %.6f", as_extreme.fdiv(splits))We never shuffle at random. plants.combination(8) yields each choice of 8
plants once, in a fixed order, so the loop sees all 12,870 splits and the
p-value is exact. Called without a block, combination returns an
Enumerator instead, so nothing builds a list of all the splits at once.
Because both groups have 8 plants, comparing totals is the same as comparing
means, and the totals stay whole Integers, so >= never meets a near-tie.
fdiv turns them into floats only for printing.
//! Exact two-sided permutation test for a difference in means, using only the
//! standard library.
//!
//! If compost made no difference, the "compost" and "control" labels are
//! arbitrary, and any 8 of the 16 plants could have carried the compost label.
//! We try every one of those splits and count how many give a difference at
//! least as large as the one we observed. That fraction is the p-value.
use std::env;
use std::fs;
/// There is no CSV parser in the standard library; the file has no quoted
/// fields, so commas suffice.
fn read(path: &str) -> (Vec<i64>, Vec<i64>) {
let text = fs::read_to_string(path).expect("cannot read data file");
let (mut control, mut compost) = (Vec::new(), Vec::new());
for line in text.lines().skip(1) {
let (group, grams) = line.split_once(',').expect("expected two fields");
let grams: i64 = grams.parse().expect("grams must be a whole number");
match group {
"control" => control.push(grams),
"compost" => compost.push(grams),
other => panic!("unknown group {other}"),
}
}
(control, compost)
}
fn main() {
let path = env::args().nth(1).expect("usage: main <data.csv>");
let (control, compost) = read(&path);
let n = compost.len();
let plants: Vec<i64> = control.iter().chain(&compost).copied().collect();
let total: i64 = plants.iter().sum();
// Both groups have n plants, so comparing totals is the same as comparing
// means, and totals stay whole numbers: no floating-point near-ties.
let observed = compost.iter().sum::<i64>() - control.iter().sum::<i64>();
// Each split is a number whose binary digits mark the chosen plants: bit i
// set means plant i carries the compost label. Of the 2^16 numbers, the
// ones with exactly n bits set are the splits.
let (mut splits, mut as_extreme) = (0u64, 0u64);
for mask in 0u32..(1 << plants.len()) {
if mask.count_ones() as usize != n {
continue;
}
let relabeled: i64 = (0..plants.len())
.filter(|&i| mask >> i & 1 == 1)
.map(|i| plants[i])
.sum();
let difference = relabeled - (total - relabeled);
splits += 1;
if difference.abs() >= observed.abs() {
as_extreme += 1; // two-sided: either direction
}
}
let mean = |xs: &[i64]| xs.iter().sum::<i64>() as f64 / xs.len() as f64;
println!("control_count {}", control.len());
println!("compost_count {n}");
println!("control_mean {:.6}", mean(&control));
println!("compost_mean {:.6}", mean(&compost));
println!("observed_difference {:.6}", observed as f64 / n as f64);
println!("splits {splits}");
println!("as_extreme {as_extreme}");
println!("p_value {:.6}", as_extreme as f64 / splits as f64);
}We never shuffle at random. Each split is a u32 from 0 to 65,535 whose
binary digits mark the chosen plants, and count_ones() keeps the 12,870
numbers with exactly 8 bits set. For each one, an iterator chain filters the
plant positions whose bit is set and sums their grams.
The totals are i64, so >= never meets a near-tie, and floats appear only
in the printed means and p-value. The itertools crate has a combinations adapter, but no
crates are installed in this image, and the run uses rustc directly.
// Exact two-sided permutation test for a difference in means, using only
// Swift's standard library and Foundation, which ships with it.
//
// If compost made no difference, the "compost" and "control" labels are
// arbitrary, and any 8 of the 16 plants could have carried the compost label.
// We try every one of those splits and count how many give a difference at
// least as large as the one we observed. That fraction is the p-value.
import Foundation
// Read group,grams lines. split skips empty pieces, so a blank last line
// doesn't become a row. Int(_:) gives an optional, which ! unwraps.
let text = try String(contentsOfFile: CommandLine.arguments[1], encoding: .utf8)
var control: [Int] = []
var compost: [Int] = []
for line in text.split(separator: "\n").dropFirst() { // dropFirst skips the header
let fields = line.split(separator: ",")
let grams = Int(fields[1])!
if fields[0] == "compost" {
compost.append(grams)
} else {
control.append(grams)
}
}
let n = compost.count
precondition(n > 0 && control.count == n, "need two groups of equal size")
let plants = control + compost
let total = plants.reduce(0, +)
// Both groups have the same size, so comparing totals is the same as
// comparing means, and totals stay whole numbers: no floating-point near-ties.
let observed = compost.reduce(0, +) - control.reduce(0, +)
// Each whole number from 0 to 2^16 - 1 is a mask: bit i set means plant i is
// labeled compost. nonzeroBitCount counts the set bits, so the masks with
// exactly n of them are the splits.
var splits = 0
var asExtreme = 0
for mask in 0..<(1 << plants.count) where mask.nonzeroBitCount == n {
var chosen = 0
for i in plants.indices where mask & (1 << i) != 0 {
chosen += plants[i]
}
let difference = chosen - (total - chosen)
splits += 1
if abs(difference) >= abs(observed) { // two-sided
asExtreme += 1
}
}
func six(_ x: Double) -> String { String(format: "%.6f", x) }
print("control_count", control.count)
print("compost_count", n)
print("control_mean", six(Double(control.reduce(0, +)) / Double(n)))
print("compost_mean", six(Double(compost.reduce(0, +)) / Double(n)))
print("observed_difference", six(Double(observed) / Double(n)))
print("splits", splits)
print("as_extreme", asExtreme)
print("p_value", six(Double(asExtreme) / Double(splits)))Every whole number from 0 to 2^16 − 1 is a mask: bit i set means plant i is
labeled compost. A for loop with a where clause keeps only the masks
whose nonzeroBitCount equals the group size, 12,870 of the 65,536, and a
second where loop adds up the grams of the plants whose bits are set.
Int totals stay whole numbers, so abs and >= count the 6 ties with
the observed difference exactly; Double appears only in the printed
output. Arrays have no combinations(ofCount:) here: that method comes
from the Swift Algorithms package, which isn’t installed.
# Exact two-sided permutation test for a difference in means, using only
# Tcl's built-in commands.
#
# If compost made no difference, the "compost" and "control" labels are
# arbitrary, and any n of the plants could have carried the compost label.
# We try every one of those splits and count how many give a difference at
# least as large as the one we observed. That fraction is the p-value.
# Returns a dict from group name to the list of its grams.
proc read_groups {path} {
set f [open $path r]
set header [split [gets $f] ,]
set group_col [lsearch -exact $header group]
set grams_col [lsearch -exact $header grams]
if {$group_col < 0 || $grams_col < 0} { error "column not found" }
set groups [dict create control {} compost {}]
while {[gets $f line] >= 0} {
set fields [split $line ,]
dict lappend groups [lindex $fields $group_col] [lindex $fields $grams_col]
}
close $f
return $groups
}
proc sum {values} {
set total 0
foreach v $values { incr total $v }
return $total
}
# Walk the plants from index next on, choosing left more of them to carry
# the compost label; chosen is the total of those picked so far. Each call
# either takes plant next or skips it, so every choice is reached once.
proc choose {next left chosen} {
global plants total observed splits as_extreme
if {$left == 0} {
incr splits
# Two-sided: abs counts a difference in either direction.
if {abs($chosen - ($total - $chosen)) >= abs($observed)} {
incr as_extreme
}
return
}
# Too few plants left to finish this choice.
if {[llength $plants] - $next < $left} { return }
choose [expr {$next + 1}] [expr {$left - 1}] \
[expr {$chosen + [lindex $plants $next]}]
choose [expr {$next + 1}] $left $chosen
}
set groups [read_groups [lindex $argv 0]]
set control [dict get $groups control]
set compost [dict get $groups compost]
set n [llength $compost]
set plants [concat $control $compost]
set total [sum $plants]
# Both groups have n plants, so comparing totals is the same as comparing
# means, and totals stay whole numbers: no floating-point near-ties.
set observed [expr {[sum $compost] - [sum $control]}]
set splits 0
set as_extreme 0
choose 0 $n 0
# double() makes / divide as a decimal; between two integers, / rounds down.
puts "control_count [llength $control]"
puts "compost_count $n"
puts "control_mean [format %.6f [expr {double([sum $control]) / [llength $control]}]]"
puts "compost_mean [format %.6f [expr {double([sum $compost]) / $n}]]"
puts "observed_difference [format %.6f [expr {double($observed) / $n}]]"
puts "splits $splits"
puts "as_extreme $as_extreme"
puts "p_value [format %.6f [expr {double($as_extreme) / $splits}]]"A recursive proc, choose, walks the plants in order. At each plant it
calls itself twice: once taking the plant into the compost group and once
skipping it. When it has taken 8, it counts the split, so every one of the
12,870 splits is reached once. The counters are global, since a Tcl proc
sees only its own variables otherwise.
The totals are whole numbers, so >= never meets a floating-point
near-tie. In expr, / between two integers rounds down, so the means
use double() first. Tcllib, which has a statistics package, isn’t
installed.