Descriptive statistics

"""Descriptive statistics on one CSV column, using only the standard library.

The statistics module would do most of this in one call each; it is avoided
on purpose so the arithmetic stays visible.
"""

import csv
import math
import sys

COLUMN = "minutes"


def read_column(path, column):
    with open(path, newline="") as f:
        return [float(row[column]) for row in csv.DictReader(f)]


def median(sorted_values):
    n = len(sorted_values)
    middle = n // 2
    if n % 2 == 1:
        return sorted_values[middle]
    return (sorted_values[middle - 1] + sorted_values[middle]) / 2


def main():
    values = read_column(sys.argv[1], COLUMN)
    values.sort()

    n = len(values)
    total = sum(values)
    mean = total / n
    # Sample variance divides by n - 1, not n (Bessel's correction).
    variance = sum((x - mean) ** 2 for x in values) / (n - 1)

    results = [
        ("sum", total),
        ("mean", mean),
        ("median", median(values)),
        ("min", values[0]),
        ("max", values[-1]),
        ("range", values[-1] - values[0]),
        ("variance", variance),
        ("std_dev", math.sqrt(variance)),
    ]

    print(f"count {n}")
    for name, value in results:
        print(f"{name} {value:.6f}")


if __name__ == "__main__":
    main()

"""Descriptive statistics on one CSV column, with the statistics module.

The built-in version of ../python/main.py: statistics.mean, median, variance,
and stdev do the arithmetic. variance and stdev are the sample versions,
dividing by n - 1.
"""

import csv
import statistics
import sys

COLUMN = "minutes"


def main():
    with open(sys.argv[1], newline="") as f:
        values = [float(row[COLUMN]) for row in csv.DictReader(f)]

    results = [
        ("sum", sum(values)),
        ("mean", statistics.mean(values)),
        ("median", statistics.median(values)),
        ("min", min(values)),
        ("max", max(values)),
        ("range", max(values) - min(values)),
        ("variance", statistics.variance(values)),
        ("std_dev", statistics.stdev(values)),
    ]

    print(f"count {len(values)}")
    for name, value in results:
        print(f"{name} {value:.6f}")


if __name__ == "__main__":
    main()

--  Descriptive statistics on one CSV column, using only Ada's standard
--  library.
--
--  Ada has no mean, median, or variance, so each is written out. Types are
--  strict: an Integer count must be converted with Long_Float (N) before it
--  can divide a Long_Float total.

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.Numerics.Long_Elementary_Functions;
use Ada.Numerics.Long_Elementary_Functions;
with Ada.Containers.Vectors;
with Ada.Containers.Generic_Array_Sort;

procedure Main is
   --  A growable list of numbers, indexed from 1.
   package Float_Vectors is new Ada.Containers.Vectors (Positive, Long_Float);

   type Float_Array is array (Positive range <>) of Long_Float;
   procedure Sort is new Ada.Containers.Generic_Array_Sort
     (Positive, Long_Float, Float_Array);

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

   function Read_Column (Path, Name : String) return Float_Array is
      File   : File_Type;
      Values : Float_Vectors.Vector;
      Col    : Natural := 0;
   begin
      Open (File, In_File, Path);
      declare
         Header : constant String := Get_Line (File);
      begin
         for C in 1 .. Ada.Strings.Fixed.Count (Header, ",") + 1 loop
            if Field (Header, C) = Name then
               Col := C;
            end if;
         end loop;
      end;
      if Col = 0 then
         raise Constraint_Error with "column not found: " & Name;
      end if;
      while not End_Of_File (File) loop
         Values.Append (Long_Float'Value (Field (Get_Line (File), Col)));
      end loop;
      Close (File);
      --  Copy the vector into a fixed-size array, which Sort works on.
      return Result : Float_Array (1 .. Natural (Values.Length)) do
         for I in Result'Range loop
            Result (I) := Values (I);
         end loop;
      end return;
   end Read_Column;

   procedure Say (Name : String; Value : Long_Float) is
   begin
      Put (Name & " ");
      Put (Value, Fore => 1, Aft => 6, Exp => 0);   --  six decimals, no E
      New_Line;
   end Say;

   Values : Float_Array := Read_Column (Ada.Command_Line.Argument (1), "minutes");
   N      : constant Positive := Values'Length;
   Total, Sum_Sq, Mean, Median, Variance : Long_Float := 0.0;
begin
   Sort (Values);

   for V of Values loop
      Total := Total + V;
   end loop;
   Mean := Total / Long_Float (N);

   --  N / 2 on two Integers rounds toward zero.
   if N mod 2 = 1 then
      Median := Values (N / 2 + 1);
   else
      Median := (Values (N / 2) + Values (N / 2 + 1)) / 2.0;
   end if;

   --  Sample variance divides by n - 1, not n (Bessel's correction).
   for V of Values loop
      Sum_Sq := Sum_Sq + (V - Mean) ** 2;
   end loop;
   Variance := Sum_Sq / Long_Float (N - 1);

   Put_Line ("count" & Integer'Image (N));
   Say ("sum", Total);
   Say ("mean", Mean);
   Say ("median", Median);
   Say ("min", Values (1));
   Say ("max", Values (N));
   Say ("range", Values (N) - Values (1));
   Say ("variance", Variance);
   Say ("std_dev", Sqrt (Variance));
end Main;

⍝ Descriptive statistics on one CSV column, from scratch in GNU APL.
⍝ APL reads right to left, with no operator precedence: 2×3+4 is 2×(3+4).
⍝ ⍝ starts a comment, and ← assigns.

∇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 .94 for 0.94, 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
COL←ROWS[1;]⍳⊂'minutes'     ⍝ ⍳ finds the column by its header name
⍝ Drop the header (1↓) and turn each text into a number (⍎¨).
X←⍎¨1↓ROWS[;COL]

⍝ ⍋ gives the indices that sort X; indexing with them sorts.
S←X[⍋X]
N←≢X                        ⍝ ≢ counts the items
MEAN←(+/X)÷N                ⍝ +/ puts + between the items: the sum
MID←(N+1)÷2                 ⍝ 60.5 for 120 values
MEDIAN←(+/S[(⌊MID),⌈MID])÷2 ⍝ average the values at the floor and ceiling of MID
⍝ Sample variance; right to left, so ÷N-1 divides by N-1.
VAR←(+/(X-MEAN)*2)÷N-1

'count' SAY N
'sum' SAY6 +/X
'mean' SAY6 MEAN
'median' SAY6 MEDIAN
'min' SAY6 S[1]             ⍝ APL counts from 1 (⎕IO is 1)
'max' SAY6 S[N]
'range' SAY6 S[N]-S[1]
'variance' SAY6 VAR
'std_dev' SAY6 VAR*0.5      ⍝ * is power, so *0.5 is the square root
)OFF

# Descriptive statistics on one CSV column, in POSIX awk.
#
# awk reads the file one line at a time, running each pattern { action } block
# on every line it matches. POSIX awk has no sort function, so values are kept
# in order as they arrive (insertion sort).

BEGIN { FS = ","; column = "minutes" }

# The header line: find which field holds the column.
NR == 1 {
    for (i = 1; i <= NF; i++)
        if ($i == column) col = i
    if (!col) {
        print "no column " column > "/dev/stderr"
        failed = 1
        exit 1
    }
    next
}

# Every other line: insert the value into the sorted array v[1..n].
{
    x = $col + 0
    for (i = n; i > 0 && v[i] > x; i--)
        v[i + 1] = v[i]
    v[i + 1] = x
    n++
    total += x
}

END {
    if (failed) exit 1

    mean = total / n
    middle = int(n / 2)
    median = (n % 2 == 1) ? v[middle + 1] : (v[middle] + v[middle + 1]) / 2

    # Sample variance divides by n - 1, not n (Bessel's correction).
    for (i = 1; i <= n; i++)
        squares += (v[i] - mean) ^ 2
    variance = squares / (n - 1)

    printf "count %d\n", n
    printf "sum %.6f\n", total
    printf "mean %.6f\n", mean
    printf "median %.6f\n", median
    printf "min %.6f\n", v[1]
    printf "max %.6f\n", v[n]
    printf "range %.6f\n", v[n] - v[1]
    printf "variance %.6f\n", variance
    printf "std_dev %.6f\n", sqrt(variance)
}

10 REM Descriptive statistics on one CSV column, in classic BASIC.
20 REM
30 REM There are no functions for sum, mean, median, variance, or sorting,
40 REM so each is written out. Subroutines are blocks of numbered lines
50 REM reached with GOSUB, and every variable is global.
60 REM
70 REM run.sh passes the data path in the DATA environment variable,
80 REM because bwbasic can't read command-line arguments.
100 COL$ = "minutes"
110 DIM X(1000), FL$(20) : REM arrays have a fixed size, set by DIM
120 GOSUB 1000 : REM read the column into X(1) to X(N)
130 GOSUB 2000 : REM sort it
140 T = 0
150 FOR I = 1 TO N
160 T = T + X(I)
170 NEXT I
180 M = T / N
190 REM Sample variance divides by N - 1, not N (Bessel's correction).
200 S = 0
210 FOR I = 1 TO N
220 S = S + (X(I) - M) ^ 2
230 NEXT I
240 V = S / (N - 1)
250 IF N MOD 2 = 1 THEN MED = X((N + 1) / 2) ELSE MED = (X(N / 2) + X(N / 2 + 1)) / 2
260 PRINT "count"; N
270 L$ = "sum" : Y = T : GOSUB 3000
280 L$ = "mean" : Y = M : GOSUB 3000
290 L$ = "median" : Y = MED : GOSUB 3000
300 L$ = "min" : Y = X(1) : GOSUB 3000
310 L$ = "max" : Y = X(N) : GOSUB 3000
320 L$ = "range" : Y = X(N) - X(1) : GOSUB 3000
330 L$ = "variance" : Y = V : GOSUB 3000
340 L$ = "std_dev" : Y = SQR(V) : GOSUB 3000
350 SYSTEM : REM without SYSTEM, bwbasic waits at its prompt
1000 REM Read the column named COL$ into X(1) to X(N).
1010 F$ = ENVIRON$("DATA")
1020 OPEN "I", #1, F$
1030 LINE INPUT #1, L$
1040 GOSUB 4000 : REM find which field of the header is COL$
1050 FOR K = 1 TO NF
1060 IF FL$(K) = COL$ THEN C = K
1070 NEXT K
1080 N = 0
1090 IF EOF(1) THEN 1140
1100 LINE INPUT #1, L$
1110 GOSUB 4000
1120 N = N + 1 : X(N) = VAL(FL$(C))
1130 GOTO 1090
1140 CLOSE #1
1150 RETURN
2000 REM Insertion sort: slide each value left past the larger ones.
2010 FOR I = 2 TO N
2020 Y = X(I) : J = I - 1
2030 IF J < 1 THEN 2060
2040 IF X(J) <= Y THEN 2060
2050 X(J + 1) = X(J) : J = J - 1 : GOTO 2030
2060 X(J + 1) = Y
2070 NEXT I
2080 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 4020

/* Descriptive statistics on one CSV column, using only the C standard library.
 *
 * C has no statistics functions to avoid, and no growable arrays or CSV
 * reader either: the program grows its own array with realloc and finds the
 * column by splitting each line on commas. */

#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>

#define COLUMN "minutes"
#define MAX_LINE 1024

/* Return the 0-based position of name in a header line, or -1. strtok
 * replaces each comma with '\0' as it goes, so it works on a copy. */
static int column_index(const char *header, const char *name) {
    char copy[MAX_LINE];
    strcpy(copy, header);
    int index = 0;
    for (char *field = strtok(copy, ",\r\n"); field; field = strtok(NULL, ",\r\n")) {
        if (strcmp(field, name) == 0) return index;
        index++;
    }
    return -1;
}

/* Read one column as doubles. The array doubles in size whenever it fills,
 * and the count comes back through *count. */
static double *read_column(const char *path, const char *name, size_t *count) {
    FILE *file = fopen(path, "r");
    if (!file) { perror(path); exit(1); }

    char line[MAX_LINE];
    if (!fgets(line, sizeof line, file)) { fprintf(stderr, "empty file\n"); exit(1); }
    int index = column_index(line, name);
    if (index < 0) { fprintf(stderr, "no column %s\n", name); exit(1); }

    size_t n = 0, capacity = 64;
    double *values = malloc(capacity * sizeof *values);
    while (fgets(line, sizeof line, file)) {
        char *field = strtok(line, ",");
        for (int i = 0; i < index && field; i++) field = strtok(NULL, ",");
        if (!field) continue;
        if (n == capacity) {
            capacity *= 2;
            values = realloc(values, capacity * sizeof *values);
        }
        values[n++] = strtod(field, NULL);
    }
    fclose(file);
    *count = n;
    return values;
}

/* qsort calls this to compare two elements, passed as untyped pointers. */
static int compare_doubles(const void *a, const void *b) {
    double x = *(const double *)a, y = *(const double *)b;
    return (x > y) - (x < y);
}

int main(int argc, char **argv) {
    if (argc < 2) { fprintf(stderr, "usage: %s data.csv\n", argv[0]); return 1; }

    size_t n;
    double *values = read_column(argv[1], COLUMN, &n);
    qsort(values, n, sizeof *values, compare_doubles);

    double total = 0;
    for (size_t i = 0; i < n; i++) total += values[i];
    double mean = total / n;

    /* Arrays count from 0: for n = 120 the middle pair is values[59], values[60]. */
    double median = n % 2 == 1 ? values[n / 2] : (values[n / 2 - 1] + values[n / 2]) / 2;

    /* Sample variance divides by n - 1, not n (Bessel's correction). */
    double squares = 0;
    for (size_t i = 0; i < n; i++) squares += (values[i] - mean) * (values[i] - mean);
    double variance = squares / (n - 1);

    printf("count %zu\n", n);
    printf("sum %.6f\n", total);
    printf("mean %.6f\n", mean);
    printf("median %.6f\n", median);
    printf("min %.6f\n", values[0]);
    printf("max %.6f\n", values[n - 1]);
    printf("range %.6f\n", values[n - 1] - values[0]);
    printf("variance %.6f\n", variance);
    printf("std_dev %.6f\n", sqrt(variance));

    free(values);
    return 0;
}

// Descriptive statistics on one CSV column, using only the C++ standard
// library.
//
// C++ has no mean, median, or variance. It does have std::accumulate, which
// adds up a range, and std::sort, so those do the general work here.

#include <algorithm>
#include <cmath>
#include <fstream>
#include <numeric>
#include <print>
#include <sstream>
#include <stdexcept>
#include <string>
#include <vector>

// Split one line at every comma.
std::vector<std::string> split(const std::string& line) {
    std::vector<std::string> fields;
    std::stringstream stream(line);
    std::string field;
    while (std::getline(stream, field, ',')) {
        fields.push_back(field);
    }
    return fields;
}

// The named column of a CSV file, as numbers. A vector grows as we
// push_back, so there's no need to know the row count first.
std::vector<double> read_column(const std::string& path, const std::string& name) {
    std::ifstream file(path);
    std::string line;
    std::getline(file, line);
    std::vector<std::string> header = split(line);
    auto found = std::find(header.begin(), header.end(), name);
    if (found == header.end()) {
        throw std::runtime_error("no column " + name);
    }
    std::size_t index = found - header.begin();

    std::vector<double> values;
    while (std::getline(file, line)) {
        if (!line.empty()) {
            values.push_back(std::stod(split(line).at(index)));
        }
    }
    return values;
}

int main(int argc, char* argv[]) {
    if (argc != 2) {
        std::println(stderr, "usage: program data.csv");
        return 1;
    }
    std::vector<double> values = read_column(argv[1], "minutes");
    std::sort(values.begin(), values.end());

    // size() is an unsigned std::size_t; n keeps it, and every division
    // below has a double on one side, so none of them rounds to a whole number.
    std::size_t n = values.size();
    // The 0.0 makes accumulate add doubles. A plain 0 would make it add
    // ints, cutting off the fraction at every step, and no warning says so.
    double sum = std::accumulate(values.begin(), values.end(), 0.0);
    double mean = sum / n;

    // Vectors count from 0: for n = 120 the middle pair is at 59 and 60.
    double median = n % 2 == 1
        ? values[n / 2]
        : (values[n / 2 - 1] + values[n / 2]) / 2;

    // Sample variance divides by n - 1, not n (Bessel's correction).
    double squares = 0.0;
    for (double value : values) {
        squares += (value - mean) * (value - mean);
    }
    double variance = squares / (n - 1);

    std::println("count {}", n);
    std::println("sum {:.6f}", sum);
    std::println("mean {:.6f}", mean);
    std::println("median {:.6f}", median);
    std::println("min {:.6f}", values.front());
    std::println("max {:.6f}", values.back());
    std::println("range {:.6f}", values.back() - values.front());
    std::println("variance {:.6f}", variance);
    std::println("std_dev {:.6f}", std::sqrt(variance));
}

      *> Descriptive statistics on one CSV column, in COBOL.
      *>
      *> COBOL has had MEAN, MEDIAN, VARIANCE, and STANDARD-DEVIATION
      *> intrinsic functions since 1989; they are avoided here so the
      *> arithmetic stays visible.
       IDENTIFICATION DIVISION.
       PROGRAM-ID. DESCRIPTIVE-STATISTICS.

       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-COLUMN-NAME           PIC X(40) VALUE "minutes".
       01  WS-EOF                   PIC X VALUE "N".
           88  END-OF-FILE          VALUE "Y".

      *> One CSV line split into fields, and where the column is.
       01  WS-FIELDS.
           05  WS-FIELD             PIC X(40) OCCURS 20 TIMES.
       01  WS-FIELD-COUNT           PIC 99.
       01  WS-POINTER               PIC 999.
       01  WS-LINE-LENGTH           PIC 999.
       01  WS-COLUMN                PIC 99 VALUE 0.

      *> The values: fixed-point decimals, 9 digits before the point
      *> and 6 after.
       01  WS-N                     PIC 9(4) VALUE 0.
       01  WS-VALUES.
           05  WS-VALUE             PIC S9(9)V9(6)
                                    OCCURS 1 TO 5000 TIMES
                                    DEPENDING ON WS-N.

       01  I                        PIC 9(4).
       01  WS-MIDDLE                PIC 9(4).
       01  WS-REMAINDER             PIC 9.
       01  WS-SUM                   PIC S9(12)V9(6) VALUE 0.
       01  WS-MEAN                  PIC S9(9)V9(12).
       01  WS-MEDIAN                PIC S9(9)V9(6).
       01  WS-SQUARES               PIC S9(12)V9(12) VALUE 0.
       01  WS-VARIANCE              PIC S9(9)V9(12).
       01  WS-STD-DEV               PIC S9(9)V9(12).

      *> Output: a name, then the value edited to six decimal places.
       01  WS-NAME                  PIC X(20).
       01  WS-RESULT                PIC S9(12)V9(12).
       01  WS-RESULT-OUT            PIC -(12)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

           READ DATA-FILE
               AT END SET END-OF-FILE TO TRUE
           END-READ
           PERFORM SPLIT-LINE
           PERFORM VARYING I FROM 1 BY 1 UNTIL I > WS-FIELD-COUNT
               IF WS-FIELD(I) = WS-COLUMN-NAME
                   MOVE I TO WS-COLUMN
               END-IF
           END-PERFORM
           IF WS-COLUMN = 0
               DISPLAY "no column " FUNCTION TRIM(WS-COLUMN-NAME)
                   UPON SYSERR
               MOVE 1 TO RETURN-CODE
               STOP RUN
           END-IF

           PERFORM UNTIL END-OF-FILE
               READ DATA-FILE
                   AT END SET END-OF-FILE TO TRUE
                   NOT AT END PERFORM ADD-VALUE
               END-READ
           END-PERFORM
           CLOSE DATA-FILE

      *>   SORT works on tables in working storage, not just files.
           SORT WS-VALUE ASCENDING

           COMPUTE WS-MEAN ROUNDED = WS-SUM / WS-N

           DIVIDE WS-N BY 2 GIVING WS-MIDDLE REMAINDER WS-REMAINDER
           IF WS-REMAINDER = 1
               MOVE WS-VALUE(WS-MIDDLE + 1) TO WS-MEDIAN
           ELSE
               COMPUTE WS-MEDIAN ROUNDED =
                   (WS-VALUE(WS-MIDDLE) + WS-VALUE(WS-MIDDLE + 1)) / 2
           END-IF

      *>   Sample variance divides by n - 1, not n (Bessel's correction).
           PERFORM VARYING I FROM 1 BY 1 UNTIL I > WS-N
               COMPUTE WS-SQUARES =
                   WS-SQUARES + (WS-VALUE(I) - WS-MEAN) ** 2
           END-PERFORM
           COMPUTE WS-VARIANCE ROUNDED = WS-SQUARES / (WS-N - 1)
           COMPUTE WS-STD-DEV ROUNDED = FUNCTION SQRT(WS-VARIANCE)

           MOVE WS-N TO WS-COUNT-OUT
           DISPLAY "count " FUNCTION TRIM(WS-COUNT-OUT)
           MOVE "sum" TO WS-NAME
           MOVE WS-SUM TO WS-RESULT
           PERFORM PRINT-RESULT
           MOVE "mean" TO WS-NAME
           MOVE WS-MEAN TO WS-RESULT
           PERFORM PRINT-RESULT
           MOVE "median" TO WS-NAME
           MOVE WS-MEDIAN TO WS-RESULT
           PERFORM PRINT-RESULT
           MOVE "min" TO WS-NAME
           MOVE WS-VALUE(1) TO WS-RESULT
           PERFORM PRINT-RESULT
           MOVE "max" TO WS-NAME
           MOVE WS-VALUE(WS-N) TO WS-RESULT
           PERFORM PRINT-RESULT
           MOVE "range" TO WS-NAME
           COMPUTE WS-RESULT = WS-VALUE(WS-N) - WS-VALUE(1)
           PERFORM PRINT-RESULT
           MOVE "variance" TO WS-NAME
           MOVE WS-VARIANCE TO WS-RESULT
           PERFORM PRINT-RESULT
           MOVE "std_dev" TO WS-NAME
           MOVE WS-STD-DEV TO WS-RESULT
           PERFORM PRINT-RESULT

           STOP RUN.

      *> Split DATA-LINE on commas into WS-FIELD(1), WS-FIELD(2), ...
       SPLIT-LINE.
           MOVE SPACES TO WS-FIELDS
           MOVE 0 TO WS-FIELD-COUNT
           MOVE 1 TO WS-POINTER
           MOVE FUNCTION LENGTH(FUNCTION TRIM(DATA-LINE TRAILING))
               TO WS-LINE-LENGTH
           PERFORM UNTIL WS-POINTER > WS-LINE-LENGTH
               ADD 1 TO WS-FIELD-COUNT
               UNSTRING DATA-LINE DELIMITED BY ","
                   INTO WS-FIELD(WS-FIELD-COUNT)
                   WITH POINTER WS-POINTER
               END-UNSTRING
           END-PERFORM.

      *> Store one data line's value and add it to the running sum.
       ADD-VALUE.
           PERFORM SPLIT-LINE
           ADD 1 TO WS-N
           COMPUTE WS-VALUE(WS-N) = FUNCTION NUMVAL(WS-FIELD(WS-COLUMN))
           ADD WS-VALUE(WS-N) TO WS-SUM.

      *> Print "name value", rounded to six decimal places.
       PRINT-RESULT.
           COMPUTE WS-RESULT-OUT ROUNDED = WS-RESULT
           DISPLAY FUNCTION TRIM(WS-NAME) " "
               FUNCTION TRIM(WS-RESULT-OUT).

;;;; Descriptive statistics on one CSV column, using only standard Common
;;;; Lisp.
;;;;
;;;; Every expression is a list with the operator first: (+ a b), (/ sum n).
;;;; Common Lisp has no mean, median, or variance, and no string split or
;;;; CSV reader, so we write those; sort, reduce, and sqrt are built in.

(defparameter *column* "minutes")

;; By default the reader turns 28.6 into a single float, which keeps only
;; about 7 digits; as a double that's 28.600000381469727. Read doubles.
;; eval-when makes the setting apply while the compiler reads this file, too,
;; not only once the program runs.
(eval-when (:compile-toplevel :load-toplevel :execute)
  (setf *read-default-float-format* 'double-float))

(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 parse-number (text)
  "The number written in TEXT, read by the Lisp reader. Binding *read-eval*
to nil stops the reader from running code written as #.(...)."
  (let ((*read-eval* nil))
    (coerce (read-from-string text) 'double-float)))

(defun read-column (path name)
  "The column called NAME in the CSV file at PATH, as a list of numbers."
  (with-open-file (in path)
    (let* ((header (split-on-comma (read-line in)))
           (index (or (position name header :test #'string=)
                      (error "no column ~a" name))))
      (loop for line = (read-line in nil)
            while line
            unless (string= line "")
              collect (parse-number (nth index (split-on-comma line)))))))

;; sort may take the list it's given apart, so we keep only what it returns.
;; With --script, the data path is the second item of *posix-argv*.
(let* ((minutes (sort (read-column (second sb-ext:*posix-argv*) *column*) #'<))
       (n (length minutes))
       (sum (reduce #'+ minutes))
       (mean (/ sum n))
       ;; Lists count from 0: for n = 120 the middle pair is at 59 and 60.
       ;; floor divides and rounds down, as integer division does elsewhere.
       (middle (floor n 2))
       (median (if (oddp n)
                   (nth middle minutes)
                   (/ (+ (nth (1- middle) minutes) (nth middle minutes)) 2)))
       ;; Sample variance divides by n - 1, not n (Bessel's correction).
       (variance (/ (reduce #'+ (mapcar (lambda (x) (expt (- x mean) 2)) minutes))
                    (1- n)))
       (low (first minutes))
       (high (car (last minutes))))
  ;; ~d prints an integer, ~,6f a number with six decimals, ~% a newline.
  (format t "count ~d~%" n)
  (format t "sum ~,6f~%" sum)
  (format t "mean ~,6f~%" mean)
  (format t "median ~,6f~%" median)
  (format t "min ~,6f~%" low)
  (format t "max ~,6f~%" high)
  (format t "range ~,6f~%" (- high low))
  (format t "variance ~,6f~%" variance)
  (format t "std_dev ~,6f~%" (sqrt variance)))

# Descriptive statistics on one CSV column, using only Elixir's standard
# library.
#
# Elixir has Enum.sum but no mean, median, or variance, so those are written
# out. Data never changes in place: each step takes a list and returns a new
# value, and |> passes a value into the next function as its first argument.

defmodule Stats do
  def read_column(path, name) do
    [header | rows] = path |> File.read!() |> String.split("\n", trim: true)
    col = header |> String.split(",") |> Enum.find_index(&(&1 == name))
    if col == nil, do: raise("column not found: #{name}")

    Enum.map(rows, fn row ->
      # Float.parse returns {number, rest}; "12" parses as 12.0 too.
      {value, ""} = row |> String.split(",") |> Enum.at(col) |> Float.parse()
      value
    end)
  end

  def median(sorted) do
    n = length(sorted)
    mid = div(n, 2)                    # div is integer division
    if rem(n, 2) == 1 do
      Enum.at(sorted, mid)
    else
      (Enum.at(sorted, mid - 1) + Enum.at(sorted, mid)) / 2
    end
  end

  def say(name, value) do
    IO.puts("#{name} #{:erlang.float_to_binary(value, decimals: 6)}")
  end
end

[path] = System.argv()
values = path |> Stats.read_column("minutes") |> Enum.sort()

n = length(values)
total = Enum.sum(values)
mean = total / n                       # / always gives a float

# Sample variance divides by n - 1, not n (Bessel's correction).
sumsq = values |> Enum.map(fn v -> (v - mean) ** 2 end) |> Enum.sum()
variance = sumsq / (n - 1)

min = List.first(values)
max = List.last(values)

IO.puts("count #{n}")
Stats.say("sum", total)
Stats.say("mean", mean)
Stats.say("median", Stats.median(values))
Stats.say("min", min)
Stats.say("max", max)
Stats.say("range", max - min)
Stats.say("variance", variance)
Stats.say("std_dev", :math.sqrt(variance))

! Descriptive statistics on one CSV column, using only standard Fortran.
!
! Fortran has intrinsic sum, minval, and maxval for arrays, but no mean,
! median, variance, or sort, so those are written out.

program descriptive_statistics
  implicit none
  integer, parameter :: dp = kind(1.0d0)   ! double precision
  character(len=*), parameter :: column = "minutes"
  real(dp), allocatable :: values(:)
  real(dp) :: total, mean, variance
  integer :: n

  values = read_column(column)
  call insertion_sort(values)

  n = size(values)
  total = sum(values)
  mean = total / n
  ! Sample variance divides by n - 1, not n (Bessel's correction).
  ! Whole-array arithmetic: values - mean subtracts from every element.
  variance = sum((values - mean)**2) / (n - 1)

  print "(a, 1x, i0)", "count", n
  call say("sum", total)
  call say("mean", mean)
  call say("median", median(values))
  call say("min", values(1))     ! Fortran arrays count from 1
  call say("max", values(n))
  call say("range", values(n) - values(1))
  call say("variance", variance)
  call say("std_dev", sqrt(variance))

contains

  function read_column(name) result(column_values)
    character(len=*), intent(in) :: name
    real(dp), allocatable :: column_values(:)
    character(len=256) :: path, line, text
    integer :: unit, status, col
    real(dp) :: value

    call get_command_argument(1, path)
    open(newunit=unit, file=path, status="old", action="read")
    read(unit, "(a)") line
    col = field_index(line, name)
    allocate(column_values(0))
    do
      read(unit, "(a)", iostat=status) line
      if (status /= 0) exit
      text = field(line, col)
      read(text, *) value   ! an internal read parses the number
      column_values = [column_values, value]   ! grows the array by one
    end do
    close(unit)
  end function read_column

  ! 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

  subroutine insertion_sort(a)
    real(dp), intent(inout) :: a(:)
    real(dp) :: key
    integer :: i, j

    do i = 2, size(a)
      key = a(i)
      j = i - 1
      do while (j >= 1)
        if (a(j) <= key) exit
        a(j + 1) = a(j)
        j = j - 1
      end do
      a(j + 1) = key
    end do
  end subroutine insertion_sort

  real(dp) function median(sorted)
    real(dp), intent(in) :: sorted(:)
    integer :: m

    m = size(sorted) / 2   ! integer division
    if (mod(size(sorted), 2) == 1) then
      median = sorted(m + 1)
    else
      median = (sorted(m) + sorted(m + 1)) / 2
    end if
  end function median

  ! 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 descriptive_statistics

// Descriptive statistics on one CSV column, using only the standard library.
//
// Go has no statistics package, so the sum, mean, median, and variance are
// written out by hand.
package main

import (
	"encoding/csv"
	"fmt"
	"math"
	"os"
	"slices"
	"strconv"
)

const column = "minutes"

// readColumn returns one column of a CSV file, found by its header name.
func readColumn(path, name string) ([]float64, error) {
	f, err := os.Open(path)
	if err != nil {
		return nil, err
	}
	defer f.Close()

	records, err := csv.NewReader(f).ReadAll()
	if err != nil {
		return nil, err
	}
	index := slices.Index(records[0], name)
	if index < 0 {
		return nil, fmt.Errorf("no column %q", name)
	}

	values := make([]float64, 0, len(records)-1)
	for _, record := range records[1:] {
		v, err := strconv.ParseFloat(record[index], 64)
		if err != nil {
			return nil, err
		}
		values = append(values, v)
	}
	return values, nil
}

// median expects the values already sorted.
func median(sorted []float64) float64 {
	n := len(sorted)
	middle := n / 2 // integer division: 120 / 2 is 60
	if n%2 == 1 {
		return sorted[middle]
	}
	return (sorted[middle-1] + sorted[middle]) / 2
}

func main() {
	values, err := readColumn(os.Args[1], column)
	if err != nil {
		fmt.Fprintln(os.Stderr, err)
		os.Exit(1)
	}
	slices.Sort(values)

	n := len(values)
	total := 0.0
	for _, v := range values {
		total += v
	}
	// n is an int; Go never converts between number types implicitly.
	mean := total / float64(n)

	// Sample variance divides by n - 1, not n (Bessel's correction).
	squares := 0.0
	for _, v := range values {
		squares += (v - mean) * (v - mean)
	}
	variance := squares / float64(n-1)

	fmt.Printf("count %d\n", n)
	fmt.Printf("sum %.6f\n", total)
	fmt.Printf("mean %.6f\n", mean)
	fmt.Printf("median %.6f\n", median(values))
	fmt.Printf("min %.6f\n", values[0])
	fmt.Printf("max %.6f\n", values[n-1])
	fmt.Printf("range %.6f\n", values[n-1]-values[0])
	fmt.Printf("variance %.6f\n", variance)
	fmt.Printf("std_dev %.6f\n", math.Sqrt(variance))
}

-- Descriptive statistics on one CSV column, using only Haskell's base library.
--
-- Haskell has no loops and no variables that change: every value below is
-- defined once, and totals come from folds such as sum. The base library has
-- no mean, median, or variance, and no CSV reader, so we split lines ourselves.

import Data.List (elemIndex, sort)
import System.Environment (getArgs)
import Text.Printf (printf)

column :: String
column = "minutes"

-- 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 named column, read as Doubles. The header gives its position.
readColumn :: String -> String -> [Double]
readColumn name text = case lines text of
  header : rows
    | Just index <- elemIndex name (splitOnComma header) ->
        [read (splitOnComma row !! index) | row <- rows, not (null row)]
  _ -> error ("no column " ++ name)

-- The middle value of a sorted list, or the mean of the two middle values.
-- Lists count from 0: for n = 120 the middle pair is positions 59 and 60.
median :: [Double] -> Double
median sorted
  | odd n = sorted !! middle
  | otherwise = (sorted !! (middle - 1) + sorted !! middle) / 2
  where
    n = length sorted
    middle = n `div` 2

main :: IO ()
main = do
  [path] <- getArgs
  text <- readFile path
  let values = sort (readColumn column text)
      n = length values
      -- length gives an Int; fromIntegral converts it so it can divide a Double.
      count = fromIntegral n :: Double
      total = sum values
      mean = total / count
      -- Sample variance divides by n - 1, not n (Bessel's correction). The
      -- exponent's type is spelled out, or -Wall warns that it guessed Integer.
      variance = sum [(x - mean) ^ (2 :: Int) | x <- values] / (count - 1)
      results =
        [ ("sum", total)
        , ("mean", mean)
        , ("median", median values)
        , ("min", head values)
        , ("max", last values)
        , ("range", last values - head values)
        , ("variance", variance)
        , ("std_dev", sqrt variance)
        ]
  printf "count %d\n" n
  mapM_ (\(name, value) -> printf "%s %.6f\n" name value) results

// Descriptive statistics on one CSV column, using only the Java standard
// library.
//
// Every Java program lives in a class; the runner calls Main.main. Java's
// streams have average() and summaryStatistics(), which would give the mean,
// min, and max in one call, so we write those out.

import java.io.IOException;
import java.nio.file.Files;
import java.nio.file.Path;
import java.util.Arrays;
import java.util.List;
import java.util.Locale;

class Main {
    static final String COLUMN = "minutes";

    // The named column of a CSV file, as an array of doubles. Reading a file
    // can fail with an IOException, and Java makes us say so with "throws".
    static double[] readColumn(Path path, String name) throws IOException {
        List<String> lines = Files.readAllLines(path);
        int index = Arrays.asList(lines.get(0).split(",")).indexOf(name);
        if (index < 0) {
            throw new IllegalArgumentException("no column " + name);
        }
        double[] values = new double[lines.size() - 1];
        for (int row = 1; row < lines.size(); row++) {
            values[row - 1] = Double.parseDouble(lines.get(row).split(",")[index]);
        }
        return values;
    }

    public static void main(String[] args) throws IOException {
        double[] values = readColumn(Path.of(args[0]), COLUMN);
        Arrays.sort(values);
        int n = values.length;

        double sum = 0;
        for (double value : values) {
            sum += value;
        }
        double mean = sum / n;

        // Arrays count from 0: for n = 120 the middle pair is at 59 and 60.
        // n / 2 divides two ints, so it rounds down to a whole number.
        double median = n % 2 == 1
            ? values[n / 2]
            : (values[n / 2 - 1] + values[n / 2]) / 2;

        // Sample variance divides by n - 1, not n (Bessel's correction).
        double squares = 0;
        for (double value : values) {
            squares += (value - mean) * (value - mean);
        }
        double variance = squares / (n - 1);

        // printf follows the computer's language settings unless told
        // otherwise; Locale.ROOT keeps the decimal point a point.
        Locale root = Locale.ROOT;
        System.out.printf(root, "count %d%n", n);
        System.out.printf(root, "sum %.6f%n", sum);
        System.out.printf(root, "mean %.6f%n", mean);
        System.out.printf(root, "median %.6f%n", median);
        System.out.printf(root, "min %.6f%n", values[0]);
        System.out.printf(root, "max %.6f%n", values[n - 1]);
        System.out.printf(root, "range %.6f%n", values[n - 1] - values[0]);
        System.out.printf(root, "variance %.6f%n", variance);
        System.out.printf(root, "std_dev %.6f%n", Math.sqrt(variance));
    }
}

// Descriptive statistics on one CSV column, using only Node's built-in modules.
// JavaScript has no statistics library, so every formula is written out.

import { readFileSync } from "node:fs";

const COLUMN = "minutes";

// Node has no CSV parser. This file has no quoted fields, so splitting each
// line on commas is enough; the header row tells us which field to read.
function readColumn(path, column) {
  const [header, ...lines] = readFileSync(path, "utf8").trim().split("\n");
  const index = header.split(",").indexOf(column);
  return lines.map((line) => Number(line.split(",")[index]));
}

function median(sortedValues) {
  const n = sortedValues.length;
  const middle = Math.floor(n / 2);
  if (n % 2 === 1) return sortedValues[middle];
  return (sortedValues[middle - 1] + sortedValues[middle]) / 2;
}

const values = readColumn(process.argv[2], COLUMN);
// sort() compares as strings unless given a comparator: [10, 9].sort() is [10, 9].
values.sort((a, b) => a - b);

const n = values.length;
const total = values.reduce((acc, x) => acc + x, 0);
const mean = total / n;
// Sample variance divides by n - 1, not n (Bessel's correction).
const variance = values.reduce((acc, x) => acc + (x - mean) ** 2, 0) / (n - 1);

const results = [
  ["sum", total],
  ["mean", mean],
  ["median", median(values)],
  ["min", values[0]],
  ["max", values.at(-1)],
  ["range", values.at(-1) - values[0]],
  ["variance", variance],
  ["std_dev", Math.sqrt(variance)],
];

console.log(`count ${n}`);
for (const [name, value] of results) {
  console.log(`${name} ${value.toFixed(6)}`);
}

# Descriptive statistics on one CSV column, using only Base Julia.
#
# The Statistics standard library has mean, median, var, and std; they are
# avoided on purpose so the arithmetic stays visible.

using Printf

const COLUMN = "minutes"

function read_column(path, column)
    lines = readlines(path)
    index = findfirst(==(column), split(lines[1], ","))
    return [parse(Float64, split(line, ",")[index]) for line in lines[2:end]]
end

function median_of_sorted(values)
    n = length(values)
    middle = n ÷ 2
    # Julia counts from 1, so for even n the middle pair is middle and middle + 1.
    isodd(n) ? values[middle + 1] : (values[middle] + values[middle + 1]) / 2
end

function main(path)
    values = sort(read_column(path, COLUMN))

    n = length(values)
    total = sum(values)
    mean = total / n
    # Sample variance divides by n - 1, not n (Bessel's correction).
    variance = sum((values .- mean) .^ 2) / (n - 1)

    results = [
        "sum" => total,
        "mean" => mean,
        "median" => median_of_sorted(values),
        "min" => first(values),
        "max" => last(values),
        "range" => last(values) - first(values),
        "variance" => variance,
        "std_dev" => sqrt(variance),
    ]

    println("count $n")
    for (name, value) in results
        @printf("%s %.6f\n", name, value)
    end
end

main(ARGS[1])

# Descriptive statistics on one CSV column, with the Statistics standard library.
#
# The built-in version of ../julia/main.jl: mean, median, var, and std do the
# arithmetic. var and std are the sample versions, dividing by n - 1, unless
# you pass corrected=false.

using Printf
using Statistics

const COLUMN = "minutes"

function read_column(path, column)
    lines = readlines(path)
    index = findfirst(==(column), split(lines[1], ","))
    return [parse(Float64, split(line, ",")[index]) for line in lines[2:end]]
end

function main(path)
    values = read_column(path, COLUMN)
    low, high = extrema(values)

    results = [
        "sum" => sum(values),
        "mean" => mean(values),
        "median" => median(values),
        "min" => low,
        "max" => high,
        "range" => high - low,
        "variance" => var(values),
        "std_dev" => std(values),
    ]

    println("count $(length(values))")
    for (name, value) in results
        @printf("%s %.6f\n", name, value)
    end
end

main(ARGS[1])

-- Descriptive statistics on one CSV column, using only Lua's standard
-- library.
--
-- Lua has no mean, median, or variance, so each is written out. Its one data
-- structure is the table, used here as an array: indexes start at 1, and #t
-- is the length.

local function read_column(path, name)
  local f = assert(io.open(path, "r"))
  local header, col = f:read("l"), nil
  local i = 0
  for field in header:gmatch("[^,]+") do
    i = i + 1
    if field == name then col = i end
  end
  assert(col, "column not found: " .. name)
  local values = {}
  for line in f:lines() do
    local j = 0
    for field in line:gmatch("[^,]+") do
      j = j + 1
      if j == col then values[#values + 1] = tonumber(field) end
    end
  end
  f:close()
  return values
end

local function median(sorted)
  local n = #sorted
  local mid = n // 2               -- // is integer (floor) division
  if n % 2 == 1 then
    return sorted[mid + 1]
  end
  return (sorted[mid] + sorted[mid + 1]) / 2
end

local function say(name, value)
  print(string.format("%s %.6f", name, value))
end

local values = read_column(arg[1], "minutes")
table.sort(values)                 -- sorts in place, smallest first

local n = #values
local total = 0
for _, v in ipairs(values) do
  total = total + v
end
local mean = total / n             -- / always gives a float

-- Sample variance divides by n - 1, not n (Bessel's correction).
local sumsq = 0
for _, v in ipairs(values) do
  sumsq = sumsq + (v - mean) ^ 2
end
local variance = sumsq / (n - 1)

local min, max = values[1], values[n]

print("count " .. n)
say("sum", total)
say("mean", mean)
say("median", median(values))
say("min", min)
say("max", max)
say("range", max - min)
say("variance", variance)
say("std_dev", math.sqrt(variance))