Concept · beginner
Linear regression
Fit a straight line to data by least squares.
Linear regression fits a straight line through a cloud of points, so you can describe how one quantity changes with another and predict one from the other. It is the simplest model that learns from data, and it is built from the pieces of descriptive statistics: means, squared deviations, and sums.
The dataset
data.csv holds 100 students’ weekly study hours and exam score. The
data is synthetic. make_data.py drew each student’s hours uniformly
between 0.5 and 10, then set the score to
plus normal noise with a standard deviation of 6, using the fixed seed
20261030. Scores were kept between 0 and 100, and 3 students hit the cap at
exactly 100. We ran it once and committed the result.
Because we made the data, you can check how close least squares gets to the true line.
The idea: least squares
Any line makes a prediction for each student. The gap between a student’s real score and the line’s prediction is a residual. Least squares picks the intercept and slope that make the sum of the squared residuals as small as possible. As with variance, squaring stops misses above and below the line from canceling out, and it punishes big misses more than small ones.
That minimization has a closed-form answer. With means and , define the sums of squares and cross-products
Then the least-squares slope and intercept, and the correlation , are
Here is the intuition. asks whether students who study more than average also score more than average. Each term is positive when both deviations have the same sign and negative when they differ. In this data, 93 of the 100 terms are positive, and comes to 4720.63. Dividing by , the spread in hours, turns that into points per hour: . The intercept then follows from one fact: the fitted line always passes through the point of means, here 5.027 hours and 68.381 points.
The correlation uses the same , but scales it by the spread in both variables, so it has no units and always lies between −1 and 1.
Reading the results
The fitted slope is 5.202961: each extra hour of study goes with about 5.2 more points. The intercept is 42.225716, the predicted score for a student who studies zero hours. Read it with care, since nobody in the data studied less than half an hour. Both land very close to the 5.2 and 42 we used to generate the data.
The correlation is , a strong positive relationship. Its square, , is the share of the variation in scores that the line explains. The squared deviations of the scores from their mean add up to . The squared residuals around the line add up to only 3980.07, and is exactly . The other 14% is the noise. The root-mean-square residual is 6.31 points, close to the standard deviation of 6 we put in.
What to watch for across languages
None of the from-scratch implementations calls a ready-made regression or
correlation function, even where the language has one: no lm() or cor()
in R, no statistics.linear_regression in Python, no regr_slope or
corr in PostgreSQL, no backslash solve, X \ y, in Julia or Octave,
and no polyfit in Octave. The R, Python,
PostgreSQL, and Octave tabs have a Built-in switch that shows those
calls, with Octave’s polyfit, and they print the same results. The from-scratch versions don’t call a
built-in mean either. They use everyday building blocks: arithmetic, a
square root, and in most languages a built-in sum (PostgreSQL adds
count, C, BASIC, Lua, and Ada, which have no sum, write every total as a loop, OCaml, with no sum either, writes one as a fold, and C++ uses
std::accumulate and std::inner_product).
Watch for how each language makes the two passes the formulas need. The
means come first, then the deviations from them. Watch how it pairs each
hours value with its score: Haskell, for one, reads two lists and
zips them back together, as Swift and Elixir do, and Common Lisp’s mapcar, Prolog’s
maplist, OCaml’s List.map2, PHP’s array_map, and Tcl’s two-list foreach walk two lists side by side. And watch how much code goes into reading
two columns by name compared with the regression itself.
Data preview
| hours | score |
|---|---|
| 5.7 | 60.4 |
| 0.9 | 42.0 |
| 8.0 | 79.9 |
| 3.1 | 55.4 |
| 3.7 | 47.1 |
| 1.9 | 60.7 |
| 1.3 | 46.8 |
| 9.3 | 98.9 |
First 8 of 100 rows · Download data.csv
Expected results
Every implementation must print these, to within a relative tolerance of 10⁻⁴.
count | 100 |
slope | 5.202961 |
intercept | 42.225716 |
r | 0.927659 |
r_squared | 0.860551 |
Implementations
"""Simple linear regression by least squares, using only the standard library.
Fits score = intercept + slope * hours. (statistics.linear_regression would do
it in one call; the sums are written out so the formulas stay visible.)
"""
import csv
import math
import sys
X, Y = "hours", "score"
def main():
with open(sys.argv[1], newline="") as f:
rows = list(csv.DictReader(f))
xs = [float(row[X]) for row in rows]
ys = [float(row[Y]) for row in rows]
n = len(xs)
mean_x = sum(xs) / n
mean_y = sum(ys) / n
# Sums of squared deviations and of cross-products.
sxx = sum((x - mean_x) ** 2 for x in xs)
syy = sum((y - mean_y) ** 2 for y in ys)
sxy = sum((x - mean_x) * (y - mean_y) for x, y in zip(xs, ys))
slope = sxy / sxx
intercept = mean_y - slope * mean_x
r = sxy / math.sqrt(sxx * syy)
print(f"count {n}")
print(f"slope {slope:.6f}")
print(f"intercept {intercept:.6f}")
print(f"r {r:.6f}")
print(f"r_squared {r * r:.6f}")
if __name__ == "__main__":
main()Python stays close to the formulas. csv.DictReader reads both columns by
their header names, and each of the three sums is a single generator
expression. zip(xs, ys) pairs every hours value with its score, so the
cross-product sum, sum((x - mean_x) * (y - mean_y) for x, y in zip(xs, ys)), reads like the formula.
What we left out is the built-in. The standard library’s statistics module
has linear_regression, which returns the slope and intercept, and
correlation, which returns r. Both give the same numbers as this file. We
skipped them so the sums stay visible. In real code you’d reach for them.
"""Simple linear regression with the statistics module.
The built-in version of ../python/main.py. statistics.linear_regression fits
score = intercept + slope * hours by least squares, and
statistics.correlation gives Pearson's r. Both arrived in Python 3.10.
"""
import csv
import statistics
import sys
X, Y = "hours", "score"
def main():
with open(sys.argv[1], newline="") as f:
rows = list(csv.DictReader(f))
xs = [float(row[X]) for row in rows]
ys = [float(row[Y]) for row in rows]
slope, intercept = statistics.linear_regression(xs, ys)
r = statistics.correlation(xs, ys)
print(f"count {len(xs)}")
print(f"slope {slope:.6f}")
print(f"intercept {intercept:.6f}")
print(f"r {r:.6f}")
print(f"r_squared {r * r:.6f}")
if __name__ == "__main__":
main()The statistics module gained linear_regression and correlation in
Python 3.10, so a straight-line fit is now two calls with no imports beyond
the standard library. linear_regression(xs, ys) returns the slope and
intercept as a named tuple, which we unpack into two variables.
Neither function reports r squared, but for a straight line with an
intercept it’s just r times r, the same shortcut the from-scratch version
uses. The order of the arguments matters: xs is the predictor (hours)
and ys is the thing being predicted (score).
-- Simple linear regression by least squares, using only Ada's standard
-- library.
--
-- Fits score = intercept + slope * hours. Ada has no mean, sum, or
-- regression function, so every total is a loop over two arrays.
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;
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;
-- 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.
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;
-- Each column is read on its own; row I of the file is Xs (I) and Ys (I).
Xs : constant Float_Array := Read_Column (Ada.Command_Line.Argument (1), "hours");
Ys : constant Float_Array := Read_Column (Ada.Command_Line.Argument (1), "score");
N : constant Positive := Xs'Length;
Sum_X, Sum_Y, Sxx, Syy, Sxy : Long_Float := 0.0;
begin
if Ys'Length /= N then -- /= is "not equal"
raise Constraint_Error with "columns differ in length";
end if;
-- First pass: the means.
for I in Xs'Range loop
Sum_X := Sum_X + Xs (I);
Sum_Y := Sum_Y + Ys (I);
end loop;
declare
Mean_X : constant Long_Float := Sum_X / Long_Float (N);
Mean_Y : constant Long_Float := Sum_Y / Long_Float (N);
begin
-- Second pass: sums of squared deviations and of cross-products.
for I in Xs'Range loop
declare
Dx : constant Long_Float := Xs (I) - Mean_X;
Dy : constant Long_Float := Ys (I) - Mean_Y;
begin
Sxx := Sxx + Dx * Dx;
Syy := Syy + Dy * Dy;
Sxy := Sxy + Dx * Dy;
end;
end loop;
declare
Slope : constant Long_Float := Sxy / Sxx;
Intercept : constant Long_Float := Mean_Y - Slope * Mean_X;
R : constant Long_Float := Sxy / Sqrt (Sxx * Syy);
begin
Put_Line ("count" & Integer'Image (N));
Say ("slope", Slope);
Say ("intercept", Intercept);
Say ("r", R);
Say ("r_squared", R * R);
end;
end;
end Main;Ada’s standard library has no sum, mean, covariance, or regression
function, so every total is a for loop over two arrays. Ada 2022’s
'Reduce attribute could fold a sum, but loops keep the two passes
clear.
We read each column with the same Read_Column and Field helpers as the
descriptive-statistics page, once for hours and once for score, and
pair them by position: row I of the file is Xs (I) and Ys (I). The
first loop adds up both means, and the second collects all three sums of
deviations at once. The means, slope, and R are constants, declared once
their values are known. Sqrt is the only call from outside the program.
⍝ Simple linear regression by least squares, from scratch in GNU APL.
⍝ Fits score = intercept + slope × hours. APL reads right to left, with no
⍝ operator precedence: 2×3+4 is 2×(3+4). ⍝ starts a comment; ← 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
⍝ Find a column by header name (⍳), drop the header (1↓), and turn each
⍝ text into a number (⍎¨).
COLUMN←{⍎¨1↓ROWS[;ROWS[1;]⍳⊂⍵]}
X←COLUMN 'hours'
Y←COLUMN 'score'
N←≢X ⍝ ≢ counts the items
DX←X-(+/X)÷N ⍝ deviations from the mean; +/ is the sum
DY←Y-(+/Y)÷N
SXX←+/DX*2 ⍝ sum of squared deviations (* is power)
SYY←+/DY*2
SXY←+/DX×DY ⍝ sum of cross-products, pair by pair
SLOPE←SXY÷SXX
INTERCEPT←((+/Y)÷N)-SLOPE×(+/X)÷N
R←SXY÷(SXX×SYY)*0.5 ⍝ *0.5 is the square root
'count' SAY N
'slope' SAY6 SLOPE
'intercept' SAY6 INTERCEPT
'r' SAY6 R
'r_squared' SAY6 R*2
)OFFEach sum is a whole-array expression: DX×DY multiplies the deviations pair
by pair, and +/ adds them up, so the cross-product sum is +/DX×DY with no
loop. The slope, intercept, and correlation follow directly from those
three sums.
APL does have a built-in that would do all of this: ⌹ (matrix divide)
solves least squares, and Y⌹1,[1.5]X returns the intercept and slope in
one step. We checked that it gives the same 42.225716 and 5.202961, and then
left it out, because it performs the method itself. *0.5, a power, is the
square root.
# Simple linear regression by least squares, in POSIX awk.
#
# Fits score = intercept + slope * hours. The means are needed before the
# sums of deviations, so the values are stored in arrays and the work is
# done in END, after the whole file has been read.
BEGIN { FS = ","; xname = "hours"; yname = "score" }
# The header line: find which fields hold the two columns.
NR == 1 {
for (i = 1; i <= NF; i++) {
if ($i == xname) xcol = i
if ($i == yname) ycol = i
}
if (!xcol || !ycol) {
print "missing column " xname " or " yname > "/dev/stderr"
failed = 1
exit 1
}
next
}
# Every other line: keep the pair and add it to the running totals.
{
n++
x[n] = $xcol + 0
y[n] = $ycol + 0
sum_x += x[n]
sum_y += y[n]
}
END {
if (failed) exit 1
mean_x = sum_x / n
mean_y = sum_y / n
# Sums of squared deviations and of cross-products.
for (i = 1; i <= n; i++) {
dx = x[i] - mean_x
dy = y[i] - mean_y
sxx += dx * dx
syy += dy * dy
sxy += dx * dy
}
slope = sxy / sxx
intercept = mean_y - slope * mean_x
r = sxy / sqrt(sxx * syy)
printf "count %d\n", n
printf "slope %.6f\n", slope
printf "intercept %.6f\n", intercept
printf "r %.6f\n", r
printf "r_squared %.6f\n", r * r
}awk splits each line on FS = "," for you, so finding the hours and
score columns in the header takes one short loop, and the running totals
need no setup because every variable starts at zero. sqrt and printf
are built in, so the end of the program reads like the formulas.
The formulas need the means before the deviations, and awk reads a file
once, line by line. So we keep every pair in two arrays, x[n] and y[n],
and do the fit in the END block after the last line. awk has no
regression function, so there was nothing to skip. The program sticks to
POSIX awk and runs unchanged on mawk, which the tests use, and on the awk
that ships with macOS.
10 REM Simple linear regression by least squares, in classic BASIC.
20 REM
30 REM Fits score = intercept + slope * hours. There is no regression or
40 REM correlation function, so the sums are written out as loops.
50 REM
60 REM run.sh passes the data path in the DATA environment variable,
70 REM because bwbasic can't read command-line arguments.
100 DIM X(1000), Y(1000), FL$(20) : REM arrays have a fixed size, set by DIM
110 GOSUB 1000 : REM read hours into X(1..N) and score into Y(1..N)
120 REM First pass: the means.
130 SX = 0 : SY = 0
140 FOR I = 1 TO N
150 SX = SX + X(I) : SY = SY + Y(I)
160 NEXT I
170 MX = SX / N : MY = SY / N
180 REM Second pass: sums of squared deviations and of cross-products.
190 XX = 0 : YY = 0 : XY = 0
200 FOR I = 1 TO N
210 DX = X(I) - MX : DY = Y(I) - MY
220 XX = XX + DX * DX : YY = YY + DY * DY : XY = XY + DX * DY
230 NEXT I
240 B = XY / XX
250 A = MY - B * MX
260 R = XY / SQR(XX * YY)
270 PRINT "count"; N
280 L$ = "slope" : V = B : GOSUB 3000
290 L$ = "intercept" : V = A : GOSUB 3000
300 L$ = "r" : V = R : GOSUB 3000
310 L$ = "r_squared" : V = R * R : GOSUB 3000
320 SYSTEM : REM without SYSTEM, bwbasic waits at its prompt
1000 REM Read the hours and score columns into X(1..N) and Y(1..N).
1010 F$ = ENVIRON$("DATA")
1020 OPEN "I", #1, F$
1030 LINE INPUT #1, L$
1040 GOSUB 4000 : REM find which fields of the header hold each column
1050 FOR K = 1 TO NF
1060 IF FL$(K) = "hours" THEN CX = K
1070 IF FL$(K) = "score" THEN CY = K
1080 NEXT K
1090 N = 0
1100 IF EOF(1) THEN 1150
1110 LINE INPUT #1, L$
1120 GOSUB 4000
1130 N = N + 1 : X(N) = VAL(FL$(CX)) : Y(N) = VAL(FL$(CY))
1140 GOTO 1100
1150 CLOSE #1
1160 RETURN
3000 REM Print L$ and V with six decimals.
3010 PRINT L$; " ";
3020 PRINT USING "#.######"; V
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 4020BASIC has no regression, correlation, or even sum function, so we write the
two passes out. The first FOR loop adds up hours and scores for the means.
The second adds the squared deviations and the cross-products, three running
totals updated on one line separated by colons. The slope, intercept, and
r then come from those sums, with SQR for the square root.
Both columns live in arrays sized by DIM, so X and Y hold at most
1,000 rows each. Names here are short (XX, XY, MX) because every
variable is global and shared with the subroutines. PRINT USING "#.######" prints six decimals, since a plain PRINT shows only about
seven digits.
/* Simple linear regression by least squares, using only the C standard library.
*
* Fits score = intercept + slope * hours. C has no regression or correlation
* function to avoid; the sums are plain loops over two arrays of doubles. */
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#define X "hours"
#define Y "score"
#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;
}
int main(int argc, char **argv) {
if (argc < 2) { fprintf(stderr, "usage: %s data.csv\n", argv[0]); return 1; }
/* Each column is read on its own; row i of the file is xs[i] and ys[i]. */
size_t n, ny;
double *xs = read_column(argv[1], X, &n);
double *ys = read_column(argv[1], Y, &ny);
if (n != ny) { fprintf(stderr, "columns differ in length\n"); return 1; }
/* First pass: the means. */
double sum_x = 0, sum_y = 0;
for (size_t i = 0; i < n; i++) { sum_x += xs[i]; sum_y += ys[i]; }
double mean_x = sum_x / n, mean_y = sum_y / n;
/* Second pass: sums of squared deviations and of cross-products. */
double sxx = 0, syy = 0, sxy = 0;
for (size_t i = 0; i < n; i++) {
double dx = xs[i] - mean_x, dy = ys[i] - mean_y;
sxx += dx * dx;
syy += dy * dy;
sxy += dx * dy;
}
double slope = sxy / sxx;
double intercept = mean_y - slope * mean_x;
double r = sxy / sqrt(sxx * syy);
printf("count %zu\n", n);
printf("slope %.6f\n", slope);
printf("intercept %.6f\n", intercept);
printf("r %.6f\n", r);
printf("r_squared %.6f\n", r * r);
free(xs);
free(ys);
return 0;
}C’s standard library has no regression, correlation, mean, or even sum
function, so there was nothing to skip: every total is a for loop. The
program reads each column with the same read_column helper as the
descriptive-statistics page, once for hours and once for score, and
pairs them by position: row i of the file is xs[i] and ys[i].
The two passes are two loops. The first adds up the means; the second
collects all three sums of deviations at once. The count n is a size_t,
an unsigned integer, but dividing a double by it converts it to double
first, so sum_x / n needs no cast.
// Simple linear regression by least squares, using only the C++ standard
// library.
//
// Fits score = intercept + slope * hours. C++ has no regression or
// correlation function to avoid. std::accumulate and std::inner_product do
// the sums.
#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;
}
// Each column is read on its own; row i of the file is xs[i] and ys[i].
std::vector<double> xs = read_column(argv[1], "hours");
std::vector<double> ys = read_column(argv[1], "score");
if (xs.size() != ys.size()) {
throw std::runtime_error("columns differ in length");
}
std::size_t n = xs.size();
// First pass: the means. The 0.0 makes accumulate add doubles.
double mean_x = std::accumulate(xs.begin(), xs.end(), 0.0) / n;
double mean_y = std::accumulate(ys.begin(), ys.end(), 0.0) / n;
// Second pass: deviations from the means. A lambda subtracts the mean
// from every value, and std::ranges::transform writes the results.
auto deviations = [](const std::vector<double>& values, double mean) {
std::vector<double> result(values.size());
std::ranges::transform(values, result.begin(),
[mean](double value) { return value - mean; });
return result;
};
std::vector<double> dx = deviations(xs, mean_x);
std::vector<double> dy = deviations(ys, mean_y);
// inner_product walks two ranges side by side and adds up a[i] * b[i]:
// the sums of squares and of cross-products.
double sxx = std::inner_product(dx.begin(), dx.end(), dx.begin(), 0.0);
double syy = std::inner_product(dy.begin(), dy.end(), dy.begin(), 0.0);
double sxy = std::inner_product(dx.begin(), dx.end(), dy.begin(), 0.0);
double slope = sxy / sxx;
double intercept = mean_y - slope * mean_x;
double r = sxy / std::sqrt(sxx * syy);
std::println("count {}", n);
std::println("slope {:.6f}", slope);
std::println("intercept {:.6f}", intercept);
std::println("r {:.6f}", r);
std::println("r_squared {:.6f}", r * r);
}Like the C version, this one reads hours and score with the same
read_column helper, once each, and pairs them by position: row i of the
file is xs[i] and ys[i]. C++ has no regression or correlation function,
so there’s nothing to skip, but its <numeric> header replaces C’s loops.
std::accumulate adds up each column for the means, a lambda subtracts the
mean from every value, and std::inner_product, which walks two ranges side
by side, gives all three sums of squares and cross-products.
Like accumulate, inner_product adds in the type of its starting value.
On 28.6, 43.8, and 11.5 dotted with themselves it returns 2867 from 0 and
2868.6499999999996 from 0.0.
*> Simple linear regression by least squares, in COBOL.
*>
*> Fits score = intercept + slope * hours from the sums of squared
*> deviations and cross-products, all in fixed-point decimals.
IDENTIFICATION DIVISION.
PROGRAM-ID. LINEAR-REGRESSION.
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-X-NAME PIC X(40) VALUE "hours".
01 WS-Y-NAME PIC X(40) VALUE "score".
01 WS-EOF PIC X VALUE "N".
88 END-OF-FILE VALUE "Y".
*> One CSV line split into fields, and where the columns are.
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-X-COLUMN PIC 99 VALUE 0.
01 WS-Y-COLUMN PIC 99 VALUE 0.
*> The (x, y) pairs.
01 WS-N PIC 9(4) VALUE 0.
01 WS-POINTS.
05 WS-POINT OCCURS 5000 TIMES.
10 WS-X PIC S9(9)V9(6).
10 WS-Y PIC S9(9)V9(6).
*> Intermediate results carry 12 decimal places.
01 I PIC 9(4).
01 WS-SUM-X PIC S9(12)V9(6) VALUE 0.
01 WS-SUM-Y PIC S9(12)V9(6) VALUE 0.
01 WS-MEAN-X PIC S9(9)V9(12).
01 WS-MEAN-Y PIC S9(9)V9(12).
01 WS-SXX PIC S9(12)V9(12) VALUE 0.
01 WS-SYY PIC S9(12)V9(12) VALUE 0.
01 WS-SXY PIC S9(12)V9(12) VALUE 0.
01 WS-SLOPE PIC S9(9)V9(12).
01 WS-INTERCEPT PIC S9(9)V9(12).
01 WS-R PIC S9(3)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-X-NAME
MOVE I TO WS-X-COLUMN
END-IF
IF WS-FIELD(I) = WS-Y-NAME
MOVE I TO WS-Y-COLUMN
END-IF
END-PERFORM
IF WS-X-COLUMN = 0 OR WS-Y-COLUMN = 0
DISPLAY "missing column hours or score" 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-POINT
END-READ
END-PERFORM
CLOSE DATA-FILE
COMPUTE WS-MEAN-X ROUNDED = WS-SUM-X / WS-N
COMPUTE WS-MEAN-Y ROUNDED = WS-SUM-Y / WS-N
*> Sums of squared deviations and of cross-products.
PERFORM VARYING I FROM 1 BY 1 UNTIL I > WS-N
COMPUTE WS-SXX = WS-SXX + (WS-X(I) - WS-MEAN-X) ** 2
COMPUTE WS-SYY = WS-SYY + (WS-Y(I) - WS-MEAN-Y) ** 2
COMPUTE WS-SXY = WS-SXY
+ (WS-X(I) - WS-MEAN-X) * (WS-Y(I) - WS-MEAN-Y)
END-PERFORM
COMPUTE WS-SLOPE ROUNDED = WS-SXY / WS-SXX
COMPUTE WS-INTERCEPT ROUNDED =
WS-MEAN-Y - WS-SLOPE * WS-MEAN-X
COMPUTE WS-R ROUNDED =
WS-SXY / FUNCTION SQRT(WS-SXX * WS-SYY)
MOVE WS-N TO WS-COUNT-OUT
DISPLAY "count " FUNCTION TRIM(WS-COUNT-OUT)
MOVE "slope" TO WS-NAME
MOVE WS-SLOPE TO WS-RESULT
PERFORM PRINT-RESULT
MOVE "intercept" TO WS-NAME
MOVE WS-INTERCEPT TO WS-RESULT
PERFORM PRINT-RESULT
MOVE "r" TO WS-NAME
MOVE WS-R TO WS-RESULT
PERFORM PRINT-RESULT
MOVE "r_squared" TO WS-NAME
COMPUTE WS-RESULT ROUNDED = WS-R * WS-R
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 (x, y) and add them to the running sums.
ADD-POINT.
PERFORM SPLIT-LINE
ADD 1 TO WS-N
COMPUTE WS-X(WS-N) = FUNCTION NUMVAL(WS-FIELD(WS-X-COLUMN))
COMPUTE WS-Y(WS-N) = FUNCTION NUMVAL(WS-FIELD(WS-Y-COLUMN))
ADD WS-X(WS-N) TO WS-SUM-X
ADD WS-Y(WS-N) TO WS-SUM-Y.
*> 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).COBOL does the arithmetic exactly. Every value is a fixed-point decimal
with a declared size, such as PIC S9(12)V9(12). The inputs have one
decimal place, so the means and all three sums of squares come out exact,
with no floating-point rounding. The first rounding happens at the slope.
Each sum is one COMPUTE inside a PERFORM VARYING loop.
The catch is that precision is declared, not given. A COMPUTE truncates
to the target’s picture unless you write ROUNDED. Everything around the
arithmetic is long too: UNSTRING splits each line, and finding two
columns by name and storing the pairs takes as much code as the regression.
At 154 lines, it’s almost four times the Python. GnuCOBOL has MEAN and
VARIANCE functions but no regression or correlation function.
;;;; Simple linear regression by least squares, using only standard Common
;;;; Lisp.
;;;;
;;;; Fits score = intercept + slope * hours. Common Lisp has no regression,
;;;; correlation, or mean function, and no CSV reader, so we write the sums
;;;; out; reduce, mapcar, and sqrt are built in.
(defparameter *x-column* "hours")
(defparameter *y-column* "score")
;; By default the reader turns 5.7 into a single float, which keeps only
;; about 7 digits. 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)))))))
(defun mean (numbers)
"The arithmetic mean of NUMBERS."
(/ (reduce #'+ numbers) (length numbers)))
;; With --script, the data path is the second item of *posix-argv*.
(let* ((path (second sb-ext:*posix-argv*))
(xs (read-column path *x-column*))
(ys (read-column path *y-column*))
(n (length xs))
;; First pass: the means. Second pass: deviations from them.
(mean-x (mean xs))
(mean-y (mean ys))
(dxs (mapcar (lambda (x) (- x mean-x)) xs))
(dys (mapcar (lambda (y) (- y mean-y)) ys))
;; Given two lists, mapcar walks them side by side, so #'* multiplies
;; each hours deviation by the score deviation from the same row.
(sxx (reduce #'+ (mapcar #'* dxs dxs)))
(syy (reduce #'+ (mapcar #'* dys dys)))
(sxy (reduce #'+ (mapcar #'* dxs dys)))
(slope (/ sxy sxx))
(intercept (- mean-y (* slope mean-x)))
(r (/ sxy (sqrt (* sxx syy)))))
;; ~d prints an integer, ~,6f a number with six decimals, ~% a newline.
(format t "count ~d~%" n)
(format t "slope ~,6f~%" slope)
(format t "intercept ~,6f~%" intercept)
(format t "r ~,6f~%" r)
(format t "r_squared ~,6f~%" (* r r)))We read hours and score as two separate lists of doubles, each found by
its header name, and pair them back up with mapcar. Given two lists,
mapcar walks them side by side, so (mapcar #'* dxs dys) multiplies each
hours deviation by the score deviation from the same row, and
(reduce #'+ ...) adds the products into Sxy.
Common Lisp has no mean, so a small mean function divides the reduce
total by length. The two passes are the order of the let*: the means
come first, then the deviation lists built from them. No regression or
correlation function is called, and there’s none in the standard to skip.
# Simple linear regression by least squares, using only Elixir's standard
# library.
#
# Fits score = intercept + slope * hours. Elixir has Enum.sum but no mean,
# correlation, or regression, so the sums are written out and the formulas
# stay visible.
defmodule Regression do
# One column, found by its header name and read as floats.
def column([header | rows], name) do
col = Enum.find_index(header, &(&1 == name))
if col == nil, do: raise("column not found: #{name}")
Enum.map(rows, fn fields ->
# Float.parse returns {number, rest}; "12" parses as 12.0 too.
{value, ""} = fields |> Enum.at(col) |> Float.parse()
value
end)
end
def mean(values), do: Enum.sum(values) / length(values)
def say(name, value) do
IO.puts("#{name} #{:erlang.float_to_binary(value, decimals: 6)}")
end
end
[path] = System.argv()
# Every line split at its commas; the first is the header.
csv =
path
|> File.read!()
|> String.split("\n", trim: true)
|> Enum.map(&String.split(&1, ","))
xs = Regression.column(csv, "hours")
ys = Regression.column(csv, "score")
# First pass: the means.
mean_x = Regression.mean(xs)
mean_y = Regression.mean(ys)
# Second pass: deviations from the means. Enum.zip pairs each hours value
# with the score on the same row, as a tuple {x, y}.
sxx = xs |> Enum.map(fn x -> (x - mean_x) ** 2 end) |> Enum.sum()
syy = ys |> Enum.map(fn y -> (y - mean_y) ** 2 end) |> Enum.sum()
sxy =
Enum.zip(xs, ys)
|> Enum.map(fn {x, y} -> (x - mean_x) * (y - mean_y) end)
|> Enum.sum()
slope = sxy / sxx
intercept = mean_y - slope * mean_x
r = sxy / :math.sqrt(sxx * syy)
IO.puts("count #{length(xs)}")
Regression.say("slope", slope)
Regression.say("intercept", intercept)
Regression.say("r", r)
Regression.say("r_squared", r * r)We split every line at its commas, find each column by its header name, and
read the values with Float.parse. Then we make the two passes the formulas
need: the means first, then the squared deviations and cross-products.
Elixir has Enum.sum, a general tool, but no mean, correlation, or
regression, so we write those out. To pair each hours value with its
score, Enum.zip turns the two lists into one list of {x, y} tuples,
and the function we map over it pulls both numbers apart by pattern
matching. The square root comes from Erlang’s :math.sqrt. Hex packages
such as Nx aren’t installed.
! Simple linear regression by least squares, using only standard Fortran.
!
! Fits score = intercept + slope * hours. Fortran has no regression or
! correlation intrinsic, so the sums are written out.
program linear_regression
implicit none
integer, parameter :: dp = kind(1.0d0) ! double precision
real(dp), allocatable :: x(:), y(:)
real(dp) :: mean_x, mean_y, sxx, syy, sxy, slope, intercept, r
integer :: n
x = read_column("hours")
y = read_column("score")
n = size(x)
mean_x = sum(x) / n
mean_y = sum(y) / n
! Sums of squared deviations and of cross-products. Whole-array
! arithmetic: x - mean_x is the array of every deviation.
sxx = sum((x - mean_x)**2)
syy = sum((y - mean_y)**2)
sxy = dot_product(x - mean_x, y - mean_y) ! sum of element-wise products
slope = sxy / sxx
intercept = mean_y - slope * mean_x
r = sxy / sqrt(sxx * syy)
print "(a, 1x, i0)", "count", n
call say("slope", slope)
call say("intercept", intercept)
call say("r", r)
call say("r_squared", r * r)
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
! 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 linear_regressionFortran has no regression or correlation intrinsic: asking gfortran for an
intrinsic named mean, cor, or regress is a compile error, so there was
nothing to skip. The sums are written out instead.
Whole-array arithmetic keeps them short. x - mean_x is the array of every
deviation, so sum((x - mean_x)**2) is one line, and
dot_product(x - mean_x, y - mean_y) adds the cross-products. Each column is
found by name in the header, and an internal read turns its text into a
number. The count n is an integer, but dividing a real sum by it promotes it
automatically, so sum(x) / n needs no conversion.
// Simple linear regression by least squares, using only the standard library.
//
// Fits score = intercept + slope * hours. Go has no regression function to
// avoid; the sums are written out so the formulas stay visible.
package main
import (
"encoding/csv"
"fmt"
"math"
"os"
"slices"
"strconv"
)
const xName, yName = "hours", "score"
// readColumns returns the named columns of a CSV file, found by header name.
func readColumns(path string, names ...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
}
columns := make([][]float64, len(names))
for c, name := range names {
index := slices.Index(records[0], name)
if index < 0 {
return nil, fmt.Errorf("no column %q", name)
}
for _, record := range records[1:] {
v, err := strconv.ParseFloat(record[index], 64)
if err != nil {
return nil, err
}
columns[c] = append(columns[c], v)
}
}
return columns, nil
}
func mean(values []float64) float64 {
total := 0.0
for _, v := range values {
total += v
}
return total / float64(len(values))
}
func main() {
columns, err := readColumns(os.Args[1], xName, yName)
if err != nil {
fmt.Fprintln(os.Stderr, err)
os.Exit(1)
}
xs, ys := columns[0], columns[1]
meanX, meanY := mean(xs), mean(ys)
// Sums of squared deviations and of cross-products.
var sxx, syy, sxy float64
for i := range xs {
dx, dy := xs[i]-meanX, ys[i]-meanY
sxx += dx * dx
syy += dy * dy
sxy += dx * dy
}
slope := sxy / sxx
intercept := meanY - slope*meanX
r := sxy / math.Sqrt(sxx*syy)
fmt.Printf("count %d\n", len(xs))
fmt.Printf("slope %.6f\n", slope)
fmt.Printf("intercept %.6f\n", intercept)
fmt.Printf("r %.6f\n", r)
fmt.Printf("r_squared %.6f\n", r*r)
}Go has no regression function, so the from-scratch version is the only
version. One loop builds all three sums at once: each pass computes the
deviations dx and dy and adds dx * dx, dy * dy, and dx * dy.
var sxx, syy, sxy float64 declares them together, and Go starts every
number at zero, so they need no initial value.
The reader takes the column names as a variadic ...string argument and
returns one slice per column, found by header name with slices.Index. As
everywhere in Go, the count is an int and must be converted with
float64() before it can divide a sum.
-- Simple linear regression by least squares, using only Haskell's base library.
--
-- Fits score = intercept + slope * hours. The base library has no regression,
-- correlation, or mean, so the sums are written out and the formulas stay
-- visible. There's no CSV reader either, so we split lines ourselves.
import Data.List (elemIndex)
import System.Environment (getArgs)
import Text.Printf (printf)
xColumn, yColumn :: String
xColumn = "hours"
yColumn = "score"
-- 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)
mean :: [Double] -> Double
mean values = sum values / fromIntegral (length values)
main :: IO ()
main = do
[path] <- getArgs
text <- readFile path
let xs = readColumn xColumn text
ys = readColumn yColumn text
-- First pass: the means.
meanX = mean xs
meanY = mean ys
-- Second pass: deviations from the means. zip pairs each hours value
-- with the score on the same row.
sxx = sum [(x - meanX) * (x - meanX) | x <- xs]
syy = sum [(y - meanY) * (y - meanY) | y <- ys]
sxy = sum [(x - meanX) * (y - meanY) | (x, y) <- zip xs ys]
slope = sxy / sxx
intercept = meanY - slope * meanX
r = sxy / sqrt (sxx * syy)
results =
[ ("slope", slope)
, ("intercept", intercept)
, ("r", r)
, ("r_squared", r * r)
]
printf "count %d\n" (length xs)
mapM_ (\(name, value) -> printf "%s %.6f\n" name value) resultsWe read hours and score as two separate lists of Double, each found by
its header name, and zip pairs them back up row by row for the
cross-products. There are no loops: each sum is a list comprehension handed
to sum, such as [(x - meanX) * (y - meanY) | (x, y) <- zip xs ys], which
reads almost like the formula for Sxy.
Base has no mean, so a small mean function divides sum by length, with
fromIntegral turning the Int count into a Double. The two passes are
just the order of the let: the means are defined first, and the
deviations use them. No regression or correlation function is called.
// Simple linear regression by least squares, using only the Java standard
// library.
//
// Fits score = intercept + slope * hours. Java's streams have average(),
// which would give each mean in one call; we write the sums out instead, so
// the formulas stay visible.
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 X = "hours", Y = "score";
// The named column of a CSV file's lines, as an array of doubles.
static double[] readColumn(List<String> lines, String name) {
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;
}
static double mean(double[] values) {
double sum = 0;
for (double value : values) {
sum += value;
}
return sum / values.length;
}
// 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]));
// Both columns come from the same rows, so xs[i] and ys[i] belong to
// the same student.
double[] xs = readColumn(lines, X);
double[] ys = readColumn(lines, Y);
int n = xs.length;
double meanX = mean(xs), meanY = mean(ys);
// Sums of squared deviations and of cross-products, in one pass.
double sxx = 0, syy = 0, sxy = 0;
for (int i = 0; i < n; i++) {
double dx = xs[i] - meanX, dy = ys[i] - meanY;
sxx += dx * dx;
syy += dy * dy;
sxy += dx * dy;
}
double slope = sxy / sxx;
double intercept = meanY - slope * meanX;
double r = sxy / Math.sqrt(sxx * syy);
// Locale.ROOT keeps the decimal point a point on any computer.
Locale root = Locale.ROOT;
System.out.printf(root, "count %d%n", n);
System.out.printf(root, "slope %.6f%n", slope);
System.out.printf(root, "intercept %.6f%n", intercept);
System.out.printf(root, "r %.6f%n", r);
System.out.printf(root, "r_squared %.6f%n", r * r);
}
}Java’s streams have average() and summaryStatistics(), which would give
each mean in one call, so we skip them and write the sums as plain for
loops over double[] arrays. A small mean method does that for both
columns, and one more loop builds sxx, syy, and sxy together from the
deviations dx and dy.
The file is read once with Files.readAllLines, and readColumn pulls each
column out of those lines by header name. Both arrays come from the same
rows in the same order, so xs[i] and ys[i] always belong to the same
student, and no zipping is needed. Every printf passes Locale.ROOT.
// Simple linear regression by least squares, using only Node's built-in modules.
// Fits score = intercept + slope * hours. JavaScript has no regression
// function to skip; the sums are written out so the formulas stay visible.
import { readFileSync } from "node:fs";
const X = "hours";
const Y = "score";
// Node has no CSV parser. No field here is quoted, so we split on commas and
// turn each line into an object keyed by the header names.
function readRows(path) {
const [header, ...lines] = readFileSync(path, "utf8").trim().split("\n");
const names = header.split(",");
return lines.map((line) => {
const fields = line.split(",");
return Object.fromEntries(names.map((name, i) => [name, fields[i]]));
});
}
const sum = (numbers) => numbers.reduce((acc, x) => acc + x, 0);
const rows = readRows(process.argv[2]);
const xs = rows.map((row) => Number(row[X]));
const ys = rows.map((row) => Number(row[Y]));
const n = xs.length;
const meanX = sum(xs) / n;
const meanY = sum(ys) / n;
// Sums of squared deviations and of cross-products.
const sxx = sum(xs.map((x) => (x - meanX) ** 2));
const syy = sum(ys.map((y) => (y - meanY) ** 2));
const sxy = sum(xs.map((x, i) => (x - meanX) * (ys[i] - meanY)));
const slope = sxy / sxx;
const intercept = meanY - slope * meanX;
const r = sxy / Math.sqrt(sxx * syy);
console.log(`count ${n}`);
console.log(`slope ${slope.toFixed(6)}`);
console.log(`intercept ${intercept.toFixed(6)}`);
console.log(`r ${r.toFixed(6)}`);
console.log(`r_squared ${(r * r).toFixed(6)}`);JavaScript stays close to the formulas, with a little more ceremony than
Python. Node has no CSV parser, so readRows splits each line on commas and
uses Object.fromEntries to build an object keyed by header name. Arrays
have no sum and no zip, so a one-line sum helper wraps reduce, and
the cross-product pairs values by index:
xs.map((x, i) => (x - meanX) * (ys[i] - meanY)).
One thing we never had to think about is integer division. JavaScript has a
single number type, a 64-bit float, so 7 / 2 is 3.5 and no ratio here
can truncate. Its standard library has no regression or correlation
function, so there was no built-in to skip.
# Simple linear regression by least squares, using only Base Julia.
#
# Fits score = intercept + slope * hours. The backslash operator, X \ y,
# would solve it in one step; the sums are written out instead.
using Printf
function read_columns(path)
lines = readlines(path)
header = split(lines[1], ",")
rows = [parse.(Float64, split(line, ",")) for line in lines[2:end]]
column(name) = [row[findfirst(==(name), header)] for row in rows]
return column("hours"), column("score")
end
function main(path)
x, y = read_columns(path)
n = length(x)
mean_x = sum(x) / n
mean_y = sum(y) / n
# Sums of squared deviations and of cross-products. The dot in .- and .*
# applies the operator to every element.
sxx = sum((x .- mean_x) .^ 2)
syy = sum((y .- mean_y) .^ 2)
sxy = sum((x .- mean_x) .* (y .- mean_y))
slope = sxy / sxx
intercept = mean_y - slope * mean_x
r = sxy / sqrt(sxx * syy)
println("count $n")
@printf("slope %.6f\n", slope)
@printf("intercept %.6f\n", intercept)
@printf("r %.6f\n", r)
@printf("r_squared %.6f\n", r^2)
end
main(ARGS[1])Julia can fit this line in one expression: [ones(n) x] \ y solves the
least-squares problem and returns the same intercept and slope, 42.225716
and 5.202961, with no using line needed. Statistics also has cor for
r. We skipped both so the sums stay visible.
Broadcasting keeps the sums short. x .- mean_x is the whole vector of
deviations, so the cross-product sum is sum((x .- mean_x) .* (y .- mean_y)), close to R’s version. Leave out a dot and Julia doesn’t guess:
subtracting a number from a vector with plain - is a MethodError.
-- Simple linear regression by least squares, using only Lua's standard
-- library.
--
-- Fits score = intercept + slope * hours. Lua has no regression,
-- correlation, or even sum function, so every total is a loop. Its one data
-- structure is the table, used here as an array: indexes start at 1, and #t
-- is the length.
-- Read the named columns into one array each, in file order.
local function read_columns(path, x_name, y_name)
local f = assert(io.open(path, "r"))
local header, x_col, y_col = f:read("l"), nil, nil
local i = 0
for field in header:gmatch("[^,]+") do
i = i + 1
if field == x_name then x_col = i end
if field == y_name then y_col = i end
end
assert(x_col and y_col, "column not found")
local xs, ys = {}, {}
for line in f:lines() do
local fields = {}
for field in line:gmatch("[^,]+") do
fields[#fields + 1] = field
end
xs[#xs + 1] = tonumber(fields[x_col])
ys[#ys + 1] = tonumber(fields[y_col])
end
f:close()
return xs, ys
end
local function mean(values)
local total = 0
for _, v in ipairs(values) do
total = total + v
end
return total / #values -- / always gives a float
end
local function say(name, value)
print(string.format("%s %.6f", name, value))
end
local xs, ys = read_columns(arg[1], "hours", "score")
local n = #xs
local mean_x, mean_y = mean(xs), mean(ys)
-- Sums of squared deviations and of cross-products. Lua has no zip, so a
-- numeric for loop uses one index i for both arrays: xs[i] and ys[i] come
-- from the same row.
local sxx, syy, sxy = 0, 0, 0
for i = 1, n do
local dx, dy = xs[i] - mean_x, ys[i] - mean_y
sxx = sxx + dx * dx
syy = syy + dy * dy
sxy = sxy + dx * dy
end
local slope = sxy / sxx
local intercept = mean_y - slope * mean_x
local r = sxy / math.sqrt(sxx * syy)
print("count " .. n)
say("slope", slope)
say("intercept", intercept)
say("r", r)
say("r_squared", r * r)Lua reads each column into its own array, then pairs them by position.
There’s no zip, so a numeric for i = 1, n loop uses one index for both
arrays: xs[i] and ys[i] come from the same row. That loop builds all
three sums at once, and math.sqrt turns them into r.
Lua’s standard library has no regression, correlation, mean, or even
sum, so every total is a loop. There’s no CSV reader or split function
either: gmatch("[^,]+"), a Lua pattern, finds each field between
commas, and tonumber turns it into a number. Remember that / always
gives a float, even when both sides are integers.
(* Simple linear regression by least squares, using only OCaml's standard
library.
Fits score = intercept + slope * hours. The standard library has no
regression, correlation, mean, or even sum, so every total is a fold and
the formulas stay visible. *)
(* The header and the rows of a CSV file, each line split at its commas. *)
let read_csv path =
In_channel.with_open_text path (fun ic ->
match List.map (String.split_on_char ',') (In_channel.input_lines ic) with
| header :: rows -> (header, rows)
| [] -> failwith "empty file")
(* One column, found by its header name and read as floats. find_index
returns an option: Some i, or None if the name is missing. *)
let column (header, rows) name =
match List.find_index (String.equal name) header with
| Some i -> List.map (fun fields -> float_of_string (List.nth fields i)) rows
| None -> failwith ("column not found: " ^ name)
let sum values = List.fold_left ( +. ) 0. values
let () =
let csv = read_csv Sys.argv.(1) in
let xs = column csv "hours" and ys = column csv "score" in
let n = float_of_int (List.length xs) in
(* First pass: the means. *)
let mean_x = sum xs /. n and mean_y = sum ys /. n in
(* Second pass: deviations from the means. map2 walks both lists side by
side, so each hours value meets the score on the same row. *)
let sxx = sum (List.map (fun x -> (x -. mean_x) ** 2.) xs) in
let syy = sum (List.map (fun y -> (y -. mean_y) ** 2.) ys) in
let sxy = sum (List.map2 (fun x y -> (x -. mean_x) *. (y -. mean_y)) xs ys) in
let slope = sxy /. sxx in
let intercept = mean_y -. (slope *. mean_x) in
let r = sxy /. sqrt (sxx *. syy) in
Printf.printf "count %d\n" (List.length xs);
List.iter
(fun (name, value) -> Printf.printf "%s %.6f\n" name value)
[ ("slope", slope); ("intercept", intercept); ("r", r); ("r_squared", r *. r) ]OCaml reads each column into its own list, then List.map2 walks the two
lists side by side, so each hours value meets the score on the same
row. If the lists ever differed in length, map2 would raise
Invalid_argument rather than drop the extra values.
OCaml’s standard library has no regression, correlation, mean, or even
sum; List.sum doesn’t exist, so we write a one-line sum with
List.fold_left ( +. ). Floats get their own operators, +., -.,
*., and /., and float_of_int turns the count into a float before we
divide. sqrt turns the sums into r.
% Simple linear regression by least squares, using only core GNU Octave.
%
% Fits score = intercept + slope * hours. Octave could do it in one call,
% with polyfit(x, y, 1) or the backslash solve [ones(n, 1) x] \ y, and corr
% gives r; the sums are written out instead.
% The named column of a CSV file, as a column vector of 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 = zeros(numel(lines) - 1, 1);
for row = 2:numel(lines)
fields = strsplit(strtrim(lines{row}), ",");
values(row - 1) = str2double(fields{index});
end
end
args = argv();
x = read_column(args{1}, "hours");
y = read_column(args{1}, "score");
n = numel(x);
mean_x = sum(x) / n;
mean_y = sum(y) / n;
% Sums of squared deviations and of cross-products. x - mean_x is the whole
% column of deviations; .* and .^ work element by element, while plain * and
% ^ would be matrix operations.
sxx = sum((x - mean_x) .^ 2);
syy = sum((y - mean_y) .^ 2);
sxy = sum((x - mean_x) .* (y - mean_y));
slope = sxy / sxx;
intercept = mean_y - slope * mean_x;
r = sxy / sqrt(sxx * syy);
printf("count %d\n", n);
printf("slope %.6f\n", slope);
printf("intercept %.6f\n", intercept);
printf("r %.6f\n", r);
printf("r_squared %.6f\n", r ^ 2);Octave can fit this line in one call. polyfit(x, y, 1) returns the slope
and intercept, 5.202961 and 42.225716, and the backslash solve
[ones(n, 1) x] \ y gives the same pair in the other order. corr gives r,
0.927659. We skip all three so the sums stay visible.
Whole-array arithmetic keeps those sums to one line each: x - mean_x is
the column of deviations, and sum((x - mean_x) .* (y - mean_y)) is the
cross-product sum. The dots matter, because plain * would ask for a matrix
product of two 100-by-1 columns, which is an error. The two columns come from
the same read_column helper as descriptive statistics.
% Simple linear regression with polyfit, Octave's polynomial fit.
%
% The built-in version of ../octave/main.m. polyfit(x, y, 1) fits a
% straight line, a polynomial of degree 1, and corr gives r.
% The named column of a CSV file, as a column vector of 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 = zeros(numel(lines) - 1, 1);
for row = 2:numel(lines)
fields = strsplit(strtrim(lines{row}), ",");
values(row - 1) = str2double(fields{index});
end
end
args = argv();
x = read_column(args{1}, "hours");
y = read_column(args{1}, "score");
% polyfit returns the coefficients highest power first: slope, then
% intercept.
p = polyfit(x, y, 1);
r = corr(x, y);
printf("count %d\n", numel(x));
printf("slope %.6f\n", p(1));
printf("intercept %.6f\n", p(2));
printf("r %.6f\n", r);
printf("r_squared %.6f\n", r ^ 2);polyfit(x, y, 1) fits a polynomial of degree 1, a straight line, by
least squares. It returns the coefficients with the highest power first,
so p(1) is the slope and p(2) the intercept. polyfit doesn’t report
how well the line fits, so corr, Octave’s correlation function, gives r,
and squaring it gives R squared.
The from-scratch version gets the same four numbers from sums of squared deviations. Reading the CSV by column name is the same code in both.
{ Simple linear regression by least squares, using only Free Pascal's
run-time library.
Fits score = intercept + slope * hours. The Math unit has Mean and
TotalVariance, but they are left out here, so every sum is a loop over
two arrays of doubles. }
program LinearRegression;
{$mode objfpc}{$H+} // 32-bit Integer, and strings longer than 255 characters
uses
SysUtils;
type
TDoubleArray = array of Double;
function ReadColumn(const Path, Name: String): TDoubleArray;
var
F: TextFile;
Line: String;
Fields: TStringArray;
Col, I: Integer;
begin
Result := nil;
AssignFile(F, Path);
Reset(F);
ReadLn(F, Line);
Fields := Line.Split(',');
Col := -1;
for I := 0 to High(Fields) do
if Fields[I] = Name then
Col := I;
if Col < 0 then
raise Exception.Create('column not found: ' + Name);
while not Eof(F) do
begin
ReadLn(F, Line);
Fields := Line.Split(',');
SetLength(Result, Length(Result) + 1); // grow by one
Result[High(Result)] := StrToFloat(Fields[Col]);
end;
CloseFile(F);
end;
procedure Say(const Name: String; Value: Double);
begin
WriteLn(Name, ' ', Value:0:6); // width 0, six decimals
end;
var
Xs, Ys: TDoubleArray;
N, I: Integer;
SumX, SumY, MeanX, MeanY, Dx, Dy, Sxx, Syy, Sxy: Double;
Slope, Intercept, R: Double;
begin
// Each column is read on its own; row I of the file is Xs[I] and Ys[I].
Xs := ReadColumn(ParamStr(1), 'hours');
Ys := ReadColumn(ParamStr(1), 'score');
N := Length(Xs);
if Length(Ys) <> N then // <> is "not equal"
raise Exception.Create('columns differ in length');
// First pass: the means.
SumX := 0;
SumY := 0;
for I := 0 to N - 1 do
begin
SumX := SumX + Xs[I];
SumY := SumY + Ys[I];
end;
MeanX := SumX / N;
MeanY := SumY / N;
// Second pass: sums of squared deviations and of cross-products.
Sxx := 0;
Syy := 0;
Sxy := 0;
for I := 0 to N - 1 do
begin
Dx := Xs[I] - MeanX;
Dy := Ys[I] - MeanY;
Sxx := Sxx + Dx * Dx;
Syy := Syy + Dy * Dy;
Sxy := Sxy + Dx * Dy;
end;
Slope := Sxy / Sxx;
Intercept := MeanY - Slope * MeanX;
R := Sxy / Sqrt(Sxx * Syy);
WriteLn('count ', N);
Say('slope', Slope);
Say('intercept', Intercept);
Say('r', R);
Say('r_squared', R * R);
end.Free Pascal’s Math unit has Mean and TotalVariance, the sum of
squared deviations, but no covariance, correlation, or regression
function, and we leave out the two it does have. Every sum is a for
loop over two arrays of doubles.
The program reads each column with the same ReadColumn helper as the
descriptive-statistics page, once for hours and once for score, and
pairs them by position: row I of the file is Xs[I] and Ys[I]. The
first loop adds up the means, and the second collects all three sums of
deviations at once. Sqrt is the only call from outside the program.
# Simple linear regression by least squares, using only Perl and its core
# modules.
#
# Fits score = intercept + slope * hours. Perl has no regression or
# correlation function, so the sums are written out. List::Util, a core
# module, supplies sum.
use strict;
use warnings;
use List::Util qw(sum);
# Returns two array references, one per named column. \@xs is a reference:
# a single value that points at the array, so a sub can return several.
sub read_columns {
my ($path, $x_name, $y_name) = @_;
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 (@xs, @ys);
while (my $line = <$fh>) {
chomp $line;
my @row = split /,/, $line;
push @xs, $row[ $col{$x_name} ];
push @ys, $row[ $col{$y_name} ];
}
close $fh;
return (\@xs, \@ys);
}
my ($xs, $ys) = read_columns($ARGV[0], 'hours', 'score');
my @xs = @$xs; # @$xs follows the reference to its array
my @ys = @$ys;
my $n = @xs; # an array used as a number is its length
my $mean_x = sum(@xs) / $n;
my $mean_y = sum(@ys) / $n;
# Sums of squared deviations and of cross-products. map runs its block once
# per item, with the item in $_; the cross-products walk the row numbers so
# each x stays paired with its y.
my $sxx = sum(map { ($_ - $mean_x) ** 2 } @xs);
my $syy = sum(map { ($_ - $mean_y) ** 2 } @ys);
my $sxy = sum(map { ($xs[$_] - $mean_x) * ($ys[$_] - $mean_y) } 0 .. $#xs);
my $slope = $sxy / $sxx;
my $intercept = $mean_y - $slope * $mean_x;
my $r = $sxy / sqrt($sxx * $syy);
print "count $n\n";
printf "slope %.6f\n", $slope;
printf "intercept %.6f\n", $intercept;
printf "r %.6f\n", $r;
printf "r_squared %.6f\n", $r * $r;Perl stays close to the formulas once the columns are read. sum comes
from List::Util, a core module, and map runs its block once per value,
with the value in $_. So the sum of squared deviations is sum over a
map that squares $_ - $mean_x. List::Util also has a zip, but the
cross-product sum simply walks the row numbers, 0 .. $#xs, and reads
$xs[$_] and $ys[$_] together.
Perl has no regression or correlation function, so there was no built-in to
skip. There’s no CSV reader either: split /,/ breaks each line on commas,
and a hash maps each header name to its position. Returning both columns
from one sub takes array references, \@xs and \@ys.
<?php
// Simple linear regression by least squares, using only PHP and the
// extensions built into it.
//
// Fits score = intercept + slope * hours. PHP has no regression or
// correlation function, so the sums are written out.
declare(strict_types=1);
// Returns the two named columns as lists of floats. The caller unpacks
// the pair with [$xs, $ys] = ...
function read_columns(string $path, string $x_name, string $y_name): array
{
$lines = file($path, FILE_IGNORE_NEW_LINES | FILE_SKIP_EMPTY_LINES);
$header = explode(',', array_shift($lines));
$xcol = array_search($x_name, $header, true);
$ycol = array_search($y_name, $header, true);
if ($xcol === false || $ycol === false) {
throw new RuntimeException("column not found");
}
$xs = [];
$ys = [];
foreach ($lines as $line) {
$row = explode(',', $line);
$xs[] = (float) $row[$xcol];
$ys[] = (float) $row[$ycol];
}
return [$xs, $ys];
}
[$xs, $ys] = read_columns($argv[1], 'hours', 'score');
$n = count($xs);
$mean_x = array_sum($xs) / $n;
$mean_y = array_sum($ys) / $n;
// Sums of squared deviations and of cross-products. Given two arrays,
// array_map walks them side by side, so each x stays paired with its y.
$sxx = array_sum(array_map(fn($x) => ($x - $mean_x) ** 2, $xs));
$syy = array_sum(array_map(fn($y) => ($y - $mean_y) ** 2, $ys));
$sxy = array_sum(array_map(
fn($x, $y) => ($x - $mean_x) * ($y - $mean_y),
$xs,
$ys,
));
$slope = $sxy / $sxx;
$intercept = $mean_y - $slope * $mean_x;
$r = $sxy / sqrt($sxx * $syy);
echo "count $n\n";
printf("slope %.6f\n", $slope);
printf("intercept %.6f\n", $intercept);
printf("r %.6f\n", $r);
printf("r_squared %.6f\n", $r * $r);PHP has no regression or correlation function, so we write the sums out.
array_sum adds up each list, and array_map with a short fn squares
each deviation. For the cross-products, we hand array_map both columns
at once: it walks them side by side and passes one x and its y to each
call, much as zip does in Python.
PHP does have a CSV reader, str_getcsv, but no field here is quoted, so
explode splits each line on commas, and array_search finds each column
in the header.
-- Simple linear regression by least squares, in PostgreSQL.
--
-- Fits score = intercept + slope * hours. PostgreSQL has regr_slope(),
-- regr_intercept(), corr(), and regr_r2() built in; they are avoided here so
-- the sums stay visible.
WITH means AS (
SELECT count(*) AS n, sum(hours) / count(*) AS mean_x, sum(score) / count(*) AS mean_y
FROM data
),
sums AS (
-- Sums of squared deviations and of cross-products.
SELECT sum((hours - mean_x) ^ 2) AS sxx,
sum((score - mean_y) ^ 2) AS syy,
sum((hours - mean_x) * (score - mean_y)) AS sxy
FROM data, means
),
fit AS (
SELECT sxy / sxx AS slope,
mean_y - (sxy / sxx) * mean_x AS intercept,
sxy / sqrt(sxx * syy) AS r
FROM means, sums
)
SELECT name, value
FROM means, fit,
LATERAL (VALUES
(1, 'count', n::text),
(2, 'slope', round(slope::numeric, 6)::text),
(3, 'intercept', round(intercept::numeric, 6)::text),
(4, 'r', round(r::numeric, 6)::text),
(5, 'r_squared', round((r * r)::numeric, 6)::text)
) AS results (position, name, value)
ORDER BY position;SQL describes results instead of loops, and the fit is just a few sums
over one table. Three common table expressions follow the formulas in
order: means, then sums of squares and cross-products, then fit. The
second pass over the data is a join to the one-row means result.
The means need their own step because aggregates can’t nest: sum(hours - sum(hours) / count(*)) is an error. Printing also takes care, since
round(value, 6) doesn’t exist for double precision, so each result is
cast to numeric first. PostgreSQL could do the whole job with
regr_slope, regr_intercept, corr, and regr_r2, which match this
query’s results. We avoided them here.
-- Simple linear regression with PostgreSQL's regression aggregates.
--
-- The built-in version of ../postgresql/main.sql. Each regr_* aggregate takes
-- (y, x), with y first: regr_slope(score, hours) fits
-- score = intercept + slope * hours.
WITH fit AS (
SELECT count(*) AS n,
regr_slope(score, hours) AS slope,
regr_intercept(score, hours) AS intercept,
corr(score, hours) AS r,
regr_r2(score, hours) AS r_squared
FROM data
)
SELECT name, value
FROM fit,
LATERAL (VALUES
(1, 'count', n::text),
(2, 'slope', round(slope::numeric, 6)::text),
(3, 'intercept', round(intercept::numeric, 6)::text),
(4, 'r', round(r::numeric, 6)::text),
(5, 'r_squared', round(r_squared::numeric, 6)::text)
) AS results (position, name, value)
ORDER BY position;PostgreSQL has a family of regression aggregates, so the whole fit is one
SELECT over the table. regr_slope, regr_intercept, and regr_r2 give
the line and its R squared, and corr gives r.
The trap is the argument order. Every one of them takes the predicted value
first and the predictor second, regr_slope(score, hours), the reverse of
how you’d write a point (x, y). Swap them and you get the slope of hours on
score instead, a different line. The family also includes regr_sxx,
regr_sxy, and regr_syy, the very sums the from-scratch version writes
out.
% Simple linear regression by least squares, using SWI-Prolog and the
% libraries that come with it.
%
% Fits score = intercept + slope * hours. A rule such as
% deviation(Mean, X, D) is true when D is X - Mean, and running the
% program asks Prolog to find values that make main true. There's no
% mean, regression, or correlation built in, so we state each one.
:- initialization(main, main).
% read_column(+Path, +Name, -Values): Values are the numbers in the column
% called Name. findall collects every Value that satisfies the goals in
% the middle, one per line.
read_column(Path, Name, Values) :-
read_file_to_string(Path, Text, []),
split_string(Text, "\n", "", [HeaderLine|Lines]),
split_string(HeaderLine, ",", "", Header),
nth0(Index, Header, Name),
findall(Value,
( member(Line, Lines),
Line \== "",
split_string(Line, ",", "", Fields),
nth0(Index, Fields, Field),
number_string(Value, Field) ),
Values).
% mean(+Values, -Mean): the sum divided by the count.
mean(Values, Mean) :-
sum_list(Values, Sum),
length(Values, N),
Mean is Sum / N.
% deviation(+Mean, +X, -D): D is X - Mean. "is" does the arithmetic;
% D = X - Mean would only build the term X - Mean.
deviation(Mean, X, D) :-
D is X - Mean.
% product(+A, +B, -P): P is A times B.
product(A, B, P) :-
P is A * B.
% sum_of_products(+As, +Bs, -Sum): multiply the lists element by element
% and add up the results. maplist with two input lists walks them side by
% side, so the first hours value meets the first score, and so on.
sum_of_products(As, Bs, Sum) :-
maplist(product, As, Bs, Products),
sum_list(Products, Sum).
main :-
current_prolog_flag(argv, [Path|_]),
read_column(Path, "hours", Xs),
read_column(Path, "score", Ys),
length(Xs, N),
% First pass: the means.
mean(Xs, MeanX),
mean(Ys, MeanY),
% Second pass: each value's deviation from its mean.
maplist(deviation(MeanX), Xs, Dxs),
maplist(deviation(MeanY), Ys, Dys),
sum_of_products(Dxs, Dxs, Sxx),
sum_of_products(Dys, Dys, Syy),
sum_of_products(Dxs, Dys, Sxy),
Slope is Sxy / Sxx,
Intercept is MeanY - Slope * MeanX,
R is Sxy / sqrt(Sxx * Syy),
RSquared is R * R,
format("count ~d~n", [N]),
forall(member(Label-Number,
[slope-Slope, intercept-Intercept, r-R,
r_squared-RSquared]),
format("~w ~6f~n", [Label, Number])).We read hours and score as two separate lists, each found by its header
name. Pairing them back up takes no zip: maplist(product, Dxs, Dys, Products) walks both lists side by side, so the first hours deviation meets
the first score deviation. The same sum_of_products rule gives Sxx, Syy,
and Sxy.
SWI-Prolog has no mean, regression, or correlation predicate to avoid. Its
library(statistics) reports run time and memory, not data. So mean is a
three-line rule, and the two passes are just the order of the goals in
main: the means come first, and the deviations use them.
# Simple linear regression by least squares, using only base R.
#
# Fits score = intercept + slope * hours. lm(score ~ hours) would do it in one
# call; the sums are written out here so the formulas stay visible.
args <- commandArgs(trailingOnly = TRUE)
data <- read.csv(args[1])
x <- data$hours
y <- data$score
n <- length(x)
mean_x <- sum(x) / n
mean_y <- sum(y) / n
# Sums of squared deviations and of cross-products, with no loop:
# x - mean_x is a whole vector of deviations.
sxx <- sum((x - mean_x)^2)
syy <- sum((y - mean_y)^2)
sxy <- sum((x - mean_x) * (y - mean_y))
slope <- sxy / sxx
intercept <- mean_y - slope * mean_x
r <- sxy / sqrt(sxx * syy)
results <- c(slope = slope, intercept = intercept, r = r, r_squared = r * r)
cat(sprintf("count %d\n", n))
cat(sprintf("%s %.6f\n", names(results), results), sep = "")R was built for this, so the from-scratch version is the long way around.
lm(score ~ hours) fits the line in one call, and cor() gives r. We
avoided both, and mean() too, to show the arithmetic.
Vectorization keeps it short anyway. x - mean_x is a whole vector of
deviations, so sum((x - mean_x) * (y - mean_y)) is the cross-product sum with no loop,
and read.csv hands back the columns by name as data$hours and
data$score. Printing is short too: sprintf also works on vectors, so one
call formats all four named results, and this is the shortest
implementation on the page.
# Simple linear regression with lm(), R's linear-model function.
#
# The built-in version of ../r/main.R. The formula score ~ hours asks for
# score = intercept + slope * hours, and summary() reports R squared.
args <- commandArgs(trailingOnly = TRUE)
data <- read.csv(args[1])
model <- lm(score ~ hours, data = data)
coefs <- coef(model) # named: "(Intercept)" and "hours"
results <- c(
slope = unname(coefs["hours"]),
intercept = unname(coefs["(Intercept)"]),
r = cor(data$hours, data$score),
r_squared = summary(model)$r.squared
)
cat(sprintf("count %d\n", nrow(data)))
cat(sprintf("%s %.6f\n", names(results), results), sep = "")lm() is R’s general linear-model function, and this is how almost
everyone fits a line in R. The formula score ~ hours reads “score
explained by hours”, and R adds the intercept on its own.
lm() does far more than this page needs. summary(model) reports a
standard error, a t value, and a p-value for each coefficient (the slope’s
standard error here is 0.212), along with R squared. We take R squared from
it, the two coefficients from coef(), and r from cor(). Note the
awkward name "(Intercept)": that’s how R labels it.
# Simple linear regression by least squares, using only the standard library.
# Fits score = intercept + slope * hours.
X = "hours"
Y = "score"
# The csv library is a bundled gem since Ruby 3.4, not standard library, so
# split lines by hand. Fine for this data: no quoted fields or embedded commas.
def read_columns(path, *columns)
header, *rows = File.readlines(path, chomp: true).map { |line| line.split(",") }
columns.map do |column|
index = header.index(column) or abort "no column #{column}"
rows.map { |row| Float(row[index]) }
end
end
xs, ys = read_columns(ARGV[0], X, Y)
n = xs.length
mean_x = xs.sum / n
mean_y = ys.sum / n
# Sums of squared deviations and of cross-products.
sxx = xs.sum { |x| (x - mean_x)**2 }
syy = ys.sum { |y| (y - mean_y)**2 }
sxy = xs.zip(ys).sum { |x, y| (x - mean_x) * (y - mean_y) }
slope = sxy / sxx
intercept = mean_y - slope * mean_x
r = sxy / Math.sqrt(sxx * syy)
puts "count #{n}"
{ slope: slope, intercept: intercept, r: r, r_squared: r * r }.each do |name, value|
puts format("%s %.6f", name, value)
endRuby reads almost like the math. sum takes a block, so
xs.sum { |x| (x - mean_x)**2 } is the sum of squared deviations in one
expression, and xs.zip(ys).sum { |x, y| ... } pairs the two columns and
unpacks each pair for the cross-products. A hash keeps the four results in
insertion order, so printing them is one loop.
Ruby’s standard library has no regression or correlation function, so there
was no built-in to skip. The awkward part is reading the file. csv is a
bundled gem rather than part of the standard library, so read_columns
splits each line on commas by hand and finds each column with
header.index. That works for this file but would break on quoted fields.
//! Simple linear regression by least squares, using only the standard library.
//!
//! Fits score = intercept + slope * hours. Rust's standard library has no
//! statistics or CSV support, so the sums and the line splitting are written
//! out; the file has no quoted fields.
use std::env;
use std::fs;
const X: &str = "hours";
const Y: &str = "score";
fn main() {
let path = env::args().nth(1).expect("usage: main <data.csv>");
let text = fs::read_to_string(path).expect("cannot read data file");
let mut lines = text.lines();
let header: Vec<&str> = lines.next().expect("empty file").split(',').collect();
let column = |name| {
header
.iter()
.position(|&h| h == name)
.expect("column not found")
};
let (ix, iy) = (column(X), column(Y));
let mut xs = Vec::new();
let mut ys = Vec::new();
for line in lines {
let fields: Vec<&str> = line.split(',').collect();
xs.push(fields[ix].parse::<f64>().unwrap());
ys.push(fields[iy].parse::<f64>().unwrap());
}
let n = xs.len() as f64;
let mean_x = xs.iter().sum::<f64>() / n;
let mean_y = ys.iter().sum::<f64>() / n;
// Sums of squared deviations and of cross-products.
let sxx: f64 = xs.iter().map(|x| (x - mean_x).powi(2)).sum();
let syy: f64 = ys.iter().map(|y| (y - mean_y).powi(2)).sum();
let sxy: f64 = xs
.iter()
.zip(&ys)
.map(|(x, y)| (x - mean_x) * (y - mean_y))
.sum();
let slope = sxy / sxx;
let intercept = mean_y - slope * mean_x;
let r = sxy / (sxx * syy).sqrt();
println!("count {}", xs.len());
println!("slope {slope:.6}");
println!("intercept {intercept:.6}");
println!("r {r:.6}");
println!("r_squared {:.6}", r * r);
}Rust’s iterators carry the formulas. Each sum is a map followed by sum,
and xs.iter().zip(&ys) pairs every hours value with its score for the
cross-product sum, much as Python’s zip does. The column lookup is a small
closure that finds a header name’s position, called once for hours and
once for score.
Most of the extra length is reading the file. The standard library has no
CSV parser, so we split lines on commas and parse each field into an f64
ourselves. There was no built-in to skip, either: unlike Python and R, Rust
ships no regression or correlation function.
// Simple linear regression by least squares, using only Swift's standard
// library and Foundation, which ships with it.
//
// Fits score = intercept + slope * hours. Swift has no regression or
// correlation function, so the sums are written out; Foundation reads the
// file and formats the numbers.
import Foundation
let xColumn = "hours"
let yColumn = "score"
// The named column of a CSV file, as numbers. A function that can fail is
// marked throws, and every call to it is marked try.
func readColumn(_ path: String, _ name: String) throws -> [Double] {
let text = try String(contentsOfFile: path, encoding: .utf8)
// split skips empty pieces, so a blank last line doesn't become a row.
var lines = text.split(separator: "\n")
let header = lines.removeFirst().split(separator: ",").map(String.init)
guard let index = header.firstIndex(of: name) else {
fatalError("no column \(name)")
}
// Double(_:) gives nil for text that isn't a number, and ! stops the
// program there instead of carrying on with a missing value.
return lines.map { line in
Double(line.split(separator: ",", omittingEmptySubsequences: false)[index])!
}
}
let path = CommandLine.arguments[1]
let xs = try readColumn(path, xColumn)
let ys = try readColumn(path, yColumn)
// First pass: the means. Double(xs.count) converts the Int count, since
// Swift never mixes Int and Double on its own.
let n = Double(xs.count)
let meanX = xs.reduce(0, +) / n
let meanY = ys.reduce(0, +) / n
// Second pass: sums of squared deviations and of cross-products. zip pairs
// each hours value with the score from the same row.
let sxx = xs.map { ($0 - meanX) * ($0 - meanX) }.reduce(0, +)
let syy = ys.map { ($0 - meanY) * ($0 - meanY) }.reduce(0, +)
let sxy = zip(xs, ys).map { x, y in (x - meanX) * (y - meanY) }.reduce(0, +)
let slope = sxy / sxx
let intercept = meanY - slope * meanX
let r = sxy / (sxx * syy).squareRoot()
func show(_ name: String, _ value: Double) {
print(name, String(format: "%.6f", value))
}
print("count", xs.count)
show("slope", slope)
show("intercept", intercept)
show("r", r)
show("r_squared", r * r)Swift has no regression or correlation function, so we write out the two
passes the formulas need. The first pass takes the means with
reduce(0, +), and the second sums the squared deviations and
cross-products. zip(xs, ys) pairs each hours value with its score,
and squareRoot(), a method on Double, finishes the correlation.
The same readColumn helper as descriptive statistics reads each column
by name, so the file is read twice. Swift won’t divide a Double by an
Int, so the count becomes Double(xs.count) before it divides the sums.
Foundation, which ships with Swift, reads the file and formats six
decimals with String(format:).
# Simple linear regression by least squares, using only Tcl's built-in
# commands.
#
# Fits score = intercept + slope * hours. Tcl has no regression or
# correlation command (Tcllib's math::statistics has one, but it isn't
# installed), so the sums are written out. expr does the arithmetic; its
# {braces} keep Tcl from substituting twice.
proc read_column {path name} {
set f [open $path r]
set header [split [gets $f] ,]
set col [lsearch -exact $header $name]
if {$col < 0} { error "column not found: $name" }
set values {}
while {[gets $f line] >= 0} {
lappend values [lindex [split $line ,] $col]
}
close $f
return $values
}
proc mean {values} {
set total 0.0
foreach v $values {
set total [expr {$total + $v}]
}
return [expr {$total / [llength $values]}]
}
proc say {name value} {
puts "$name [format %.6f $value]"
}
set path [lindex $argv 0]
set xs [read_column $path hours]
set ys [read_column $path score]
set n [llength $xs]
set mean_x [mean $xs]
set mean_y [mean $ys]
# Sums of squared deviations and of cross-products. foreach can walk two
# lists side by side: each pass takes the next x from xs and the next y
# from ys, so every hours value meets its own score.
set sxx 0.0
set syy 0.0
set sxy 0.0
foreach x $xs y $ys {
set dx [expr {$x - $mean_x}]
set dy [expr {$y - $mean_y}]
set sxx [expr {$sxx + $dx * $dx}]
set syy [expr {$syy + $dy * $dy}]
set sxy [expr {$sxy + $dx * $dy}]
}
set slope [expr {$sxy / $sxx}]
set intercept [expr {$mean_y - $slope * $mean_x}]
set r [expr {$sxy / sqrt($sxx * $syy)}]
puts "count $n"
say slope $slope
say intercept $intercept
say r $r
say r_squared [expr {$r * $r}]Tcl reads each column on its own with read_column, then pairs them in
one loop. foreach x $xs y $ys walks both lists side by side, so each pass
gets one hours value and the score from the same row. That loop builds
all three sums at once, and expr turns them into the slope, intercept,
and r.
Tcl’s own commands have no regression, correlation, or even mean. Tcllib
has them in math::statistics, but it isn’t installed, so each sum is
written out. One thing to know about that two-list foreach: if one list
runs short, Tcl pads it with empty strings instead of stopping, so both
columns need the same length.