Concept · intermediate
Logistic regression
Predict a yes-or-no outcome by fitting probabilities with gradient descent.
Logistic regression predicts a yes-or-no outcome. Instead of fitting a straight line through the answers, it fits a probability between 0 and 1, and calls “yes” when that probability is at least 0.5. It’s the standard first model for classification, and a neural network answering a yes-or-no question ends with the same calculation.
The dataset
data.csv holds 200 synthetic students with weekly study hours, nightly
sleep hours, and whether they passed the exam (1) or not (0). We drew it
once with make_data.py (fixed seed 20261101): each student passed with
probability , where . The split column marks 150 rows as train and 50
as test. 92 of the 150 training students passed.
Why not use linear regression? Fit a line to the 0s and 1s here and it predicts a 1.11 “probability” for one student and −0.07 for another: 23 of the 150 land outside 0 to 1.
The model
Logistic regression computes a score from the features, then squeezes it into a probability with the sigmoid function:
A score of 0 gives exactly 0.5. Large positive scores approach 1 and large negative ones approach 0, but never reach them.
To judge a set of weights, we use the log-loss: the average of for students who passed and for those who didn’t. Confident wrong answers cost a lot; confident right ones cost almost nothing. Unlike least squares, there is no formula that solves for the best weights directly, so we search for them.
The method
- Scale each feature to mean 0 and standard deviation 1, using the training rows’ mean and sample standard deviation. Study hours (sd 2.71) and sleep (sd 1.40) then move on the same scale, so one learning rate suits both.
- Start with every weight at 0. Every student gets , and the log-loss is .
- Step downhill 1,000 times. The gradient of the log-loss for weight is the average of over the training rows, and each step subtracts 0.5 times it from every weight at once.
The first step already moves the hours weight furthest, to 0.156, because the students who passed mostly studied more than average. The loss falls to 0.641 after one step, 0.471 after 10, and 0.419409 after 100; by step 200 it has settled at 0.419364. Newton’s method, which homes in on the best weights in a few large steps, lands on the same six decimals, so 1,000 steps is plenty.
Reading the results
The fitted weights are an intercept of 0.753594, 1.924812 for hours, and 0.367271 for sleep. On the scaled features, one standard deviation more study multiplies the odds of passing by ; one more of sleep multiplies them by 1.44.
Take a test student who studied 3.7 hours and slept 6.4. Scaled, that’s −0.7284 for hours and 0.0109 for sleep, so and . The model says “fail”, and the student did fail.
The model gets 37 of the 50 test students right, an accuracy of 0.74. Predicting “pass” for everyone would score 0.68, since 34 of the 50 passed. Nine of the 13 misses are students who passed despite a low predicted probability. Converted back to hours, the fit finds 0.71 per study hour and 0.26 per hour of sleep, against the 0.6 and 0.5 we generated with: 150 students is not much data, and sleep’s smaller effect is the harder one to pin down.
What to watch for across languages
None of these implementations calls a ready-made logistic regression. Of the
twenty-nine languages, only R ships one, glm(..., family = binomial), which finds
the same weights to six decimals in five steps of a cleverer method. The R
tab’s Built-in switch shows it. The
rest use general tools: exp, log, square roots, sums, and loops.
Watch how each takes one step. R, Julia, Octave, Fortran, and APL do it
with matrix products (%*%, * in Julia and Octave, matmul, and +.×),
with no loop over students. In Octave, the dot in the sigmoid’s ./ matters:
with a plain /, the first step quietly builds a 150 by 150 matrix. Python, Ruby, JavaScript, Swift, Haskell, OCaml, Elixir, Common Lisp, Prolog, Tcl, and Perl build each new list of
weights from the old one, and Rust sums
the gradient with iterator chains. awk, BASIC, COBOL, Go, C, C++, Java, Pascal, Lua, Ada, and PHP add it up
in an explicit loop, and PostgreSQL, which has no loops at all, makes each step
one row of a WITH RECURSIVE query.
Watch, too, that every weight moves from the same gradient. Go’s notes show
what happens if you update inside the loop over students: the intercept
comes out 0.758788 instead of 0.753594. And watch COBOL, which does all of
this in fixed-point decimals. With 18 decimal places it matches the
floating-point languages to 14, but its 150,000 calls to EXP make it the
slowest version here.
Data preview
| hours | sleep | passed | split |
|---|---|---|---|
| 1.0 | 7.8 | 0 | train |
| 7.3 | 5.2 | 0 | train |
| 6.7 | 6.4 | 0 | train |
| 3.6 | 4.7 | 0 | train |
| 5.7 | 6.2 | 1 | train |
| 0.6 | 5.1 | 0 | train |
| 8.1 | 7.3 | 1 | train |
| 8.3 | 8.0 | 1 | train |
First 8 of 200 rows · Download data.csv
Expected results
Every implementation must print these, to within a relative tolerance of 10⁻⁴.
train_count | 150 |
test_count | 50 |
intercept | 0.753594 |
weight_hours | 1.924812 |
weight_sleep | 0.367271 |
log_loss | 0.419364 |
correct | 37 |
accuracy | 0.740000 |
Implementations
"""Logistic regression by gradient descent, using only the standard library.
Predicts whether a student passed from study hours and sleep. Each feature is
rescaled with the training rows' mean and standard deviation, then the weights
start at zero and take STEPS fixed steps of batch gradient descent on the
average log-loss.
"""
import csv
import math
import sys
FEATURES = ("hours", "sleep")
LABEL = "passed"
LEARNING_RATE = 0.5
STEPS = 1000
def read(path):
train, test = [], []
with open(path, newline="") as f:
for row in csv.DictReader(f):
example = ([float(row[name]) for name in FEATURES], int(row[LABEL]))
(train if row["split"] == "train" else test).append(example)
return train, test
def sigmoid(z):
return 1 / (1 + math.exp(-z))
def main():
train, test = read(sys.argv[1])
# Scale each feature to mean 0 and (sample) standard deviation 1, using
# the training rows only. The leading 1 multiplies the intercept.
columns = list(zip(*(x for x, _ in train)))
means = [sum(c) / len(c) for c in columns]
sds = [math.sqrt(sum((v - m) ** 2 for v in c) / (len(c) - 1)) for c, m in zip(columns, means)]
def scale(x):
return [1.0] + [(v - m) / s for v, m, s in zip(x, means, sds)]
xs = [scale(x) for x, _ in train]
ys = [y for _, y in train]
n = len(xs)
weights = [0.0] * len(xs[0])
for _ in range(STEPS):
# Gradient of the average log-loss: the mean of (p - y) times each input.
errors = [sigmoid(sum(w * v for w, v in zip(weights, x))) - y for x, y in zip(xs, ys)]
gradient = [sum(e * x[j] for e, x in zip(errors, xs)) / n for j in range(len(weights))]
weights = [w - LEARNING_RATE * g for w, g in zip(weights, gradient)]
probabilities = [sigmoid(sum(w * v for w, v in zip(weights, x))) for x in xs]
log_loss = -sum(y * math.log(p) + (1 - y) * math.log(1 - p) for p, y in zip(probabilities, ys)) / n
def predict(x):
return 1 if sigmoid(sum(w * v for w, v in zip(weights, scale(x)))) >= 0.5 else 0
correct = sum(predict(x) == y for x, y in test)
print(f"train_count {n}")
print(f"test_count {len(test)}")
print(f"intercept {weights[0]:.6f}")
print(f"weight_hours {weights[1]:.6f}")
print(f"weight_sleep {weights[2]:.6f}")
print(f"log_loss {log_loss:.6f}")
print(f"correct {correct}")
print(f"accuracy {correct / len(test):.6f}")
if __name__ == "__main__":
main()Python’s standard library has statistics.linear_regression but nothing for
logistic regression, and scikit-learn’s LogisticRegression isn’t
installed, so there was no built-in to skip. Each step of gradient descent
is two list comprehensions: one for every student’s error, p - y, and one
that averages those errors times each input. Building a new weights list
from the old one means every weight updates from the same gradient.
math.exp raises an OverflowError above about 709, so sigmoid would
fail for a score below −709. Scores here stay between −3.4 and 4.4, but
production code guards against it.
-- Logistic regression by gradient descent, using only Ada's standard
-- library.
--
-- Predicts whether a student passed from study hours and sleep. Each
-- feature is rescaled with the training rows' mean and standard deviation,
-- then the weights start at zero and take Steps fixed steps of batch
-- gradient descent on the average log-loss.
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
Learning_Rate : constant := 0.5;
Steps : constant := 1000;
-- Inputs (0) is always 1, so weight 0 acts as the intercept; 1 and 2
-- are hours and sleep.
type Inputs is array (0 .. 2) of Long_Float;
type Example is record
X : Inputs;
Y : Natural; -- 1 if the student passed, else 0
end record;
package Example_Vectors is new Ada.Containers.Vectors (Positive, Example);
use Example_Vectors;
-- The text of field number Col (counting from 1) in a CSV line.
function Field (Line : String; Col : Positive) return String is
Start : Positive := Line'First;
begin
for I in 1 .. Col - 1 loop
Start := Index (Line (Start .. Line'Last), ",") + 1;
end loop;
declare
Stop : constant Natural := Index (Line (Start .. Line'Last), ",");
begin
return (if Stop = 0 then Line (Start .. Line'Last)
else Line (Start .. Stop - 1));
end;
end Field;
-- The number (counting from 1) of the column called Name.
function Column (Header, Name : String) return Positive is
begin
for C in 1 .. Ada.Strings.Fixed.Count (Header, ",") + 1 loop
if Field (Header, C) = Name then
return C;
end if;
end loop;
raise Constraint_Error with "column not found: " & Name;
end Column;
-- Read the CSV into training and test rows, sorted by the "split" column.
procedure Read_Rows (Path : String; Train, Test : out Vector) is
File : File_Type;
begin
Open (File, In_File, Path);
declare
Header : constant String := Get_Line (File);
Hours : constant Positive := Column (Header, "hours");
Sleep : constant Positive := Column (Header, "sleep");
Passed : constant Positive := Column (Header, "passed");
Split : constant Positive := Column (Header, "split");
begin
while not End_Of_File (File) loop
declare
Line : constant String := Get_Line (File);
Row : constant Example :=
(X => [1.0,
Long_Float'Value (Field (Line, Hours)),
Long_Float'Value (Field (Line, Sleep))],
Y => Natural'Value (Field (Line, Passed)));
begin
if Field (Line, Split) = "train" then
Train.Append (Row);
else
Test.Append (Row);
end if;
end;
end loop;
end;
Close (File);
end Read_Rows;
-- The probability of passing: the sigmoid of the weighted sum of inputs.
function Predict (W, X : Inputs) return Long_Float is
Z : Long_Float := 0.0;
begin
for J in Inputs'Range loop
Z := Z + W (J) * X (J);
end loop;
return 1.0 / (1.0 + Exp (-Z));
end Predict;
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;
Train, Test : Vector;
W : Inputs := [others => 0.0];
Correct : Natural := 0;
Loss : Long_Float := 0.0;
begin
Read_Rows (Ada.Command_Line.Argument (1), Train, Test);
declare
N : constant Long_Float := Long_Float (Train.Length);
begin
-- Scale each feature to mean 0 and (sample) standard deviation 1,
-- using the training rows only, then apply the same shift to the
-- test rows. Input 0 holds the ones, so the features are 1 and 2.
for J in 1 .. 2 loop
declare
Total, Squares : Long_Float := 0.0;
Mean, Sd : Long_Float;
begin
for R of Train loop
Total := Total + R.X (J);
end loop;
Mean := Total / N;
for R of Train loop
Squares := Squares + (R.X (J) - Mean) ** 2;
end loop;
Sd := Sqrt (Squares / (N - 1.0));
-- "for R of" over a variable vector lets us change R in place.
for R of Train loop
R.X (J) := (R.X (J) - Mean) / Sd;
end loop;
for R of Test loop
R.X (J) := (R.X (J) - Mean) / Sd;
end loop;
end;
end loop;
for Step in 1 .. Steps loop
-- Gradient of the average log-loss: the mean of (p - y) times each
-- input. It is summed over every student before any weight moves.
declare
Gradient : Inputs := [others => 0.0];
Error : Long_Float;
begin
for R of Train loop
Error := Predict (W, R.X) - Long_Float (R.Y);
for J in Inputs'Range loop
Gradient (J) := Gradient (J) + Error * R.X (J);
end loop;
end loop;
for J in Inputs'Range loop
W (J) := W (J) - Learning_Rate * Gradient (J) / N;
end loop;
end;
end loop;
for R of Train loop
Loss := Loss - (if R.Y = 1 then Log (Predict (W, R.X))
else Log (1.0 - Predict (W, R.X)));
end loop;
Loss := Loss / N;
end;
for R of Test loop
-- Boolean'Pos turns False and True into 0 and 1.
if Boolean'Pos (Predict (W, R.X) >= 0.5) = R.Y then
Correct := Correct + 1;
end if;
end loop;
Put_Line ("train_count" & Train.Length'Image);
Put_Line ("test_count" & Test.Length'Image);
Say ("intercept", W (0));
Say ("weight_hours", W (1));
Say ("weight_sleep", W (2));
Say ("log_loss", Loss);
Put_Line ("correct" & Correct'Image);
Say ("accuracy", Long_Float (Correct) / Long_Float (Test.Length));
end Main;We hold each student as a record: an Inputs array whose entry 0 is
always 1, so weight 0 acts as the intercept, plus the 0 or 1 label. Ada’s
standard library has no mean, standard deviation, or logistic regression,
so the scaling and every gradient step are written out as loops, and the
gradient is summed over all 150 training students before any weight
moves. Writing for R of Train loop over a vector variable lets you
change each row in place. Boolean'Pos turns a prediction into 0 or 1, so
it can be compared with the label. Alire crates aren’t installed.
⍝ Logistic regression by gradient descent, from scratch in GNU APL.
⍝ Predicts whether a student passed from study hours and sleep. Each feature
⍝ is rescaled with the training rows' mean and sample standard deviation, then
⍝ the weights start at zero and take STEPS steps of batch gradient descent on
⍝ the average log-loss. 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 ¯
∇
RATE←0.5 ⋄ STEPS←1000
⍝ Read the file. ⎕ARG is the command line; ↑¯1↑ takes its last item, the
⍝ data path. ⎕FIO[26] reads the whole file as one string.
TEXT←⎕FIO[26]↑¯1↑⎕ARG
LINES←(TEXT≠⎕UCS 10)⊂TEXT ⍝ ⊂ splits at the 0s of the mask: one item per line
⍝ Split each (¨) line at commas; ⊃ stacks them into a table.
ROWS←⊃{(⍵≠',')⊂⍵}¨LINES
⍝ A column by header name (⍳), minus the header (1↓); ⍎¨ makes numbers.
TEXTS←{1↓ROWS[;ROWS[1;]⍳⊂⍵]}
COLUMN←{⍎¨TEXTS ⍵}
F←⍉⊃COLUMN¨'hours' 'sleep' ⍝ one row per student, one column per feature
Y←COLUMN 'passed'
TRAIN←(TEXTS 'split')∊⊂'train' ⍝ 1 for a training row, 0 for a test row
⍝ Training mean and sample standard deviation of each column. ⌿ works down
⍝ the columns: +⌿ sums each one, and TRAIN⌿ keeps the training rows.
FT←TRAIN⌿F
N←≢FT ⍝ ≢ counts the rows
MEAN←(+⌿FT)÷N
⍝ (⍴FT)⍴MEAN repeats the means on every row; *0.5 is the square root.
SD←((+⌿(FT-(⍴FT)⍴MEAN)*2)÷N-1)*0.5
⍝ Scale every row with the training statistics, and put a 1 in front (1,)
⍝ to multiply the intercept.
X←1,(F-(⍴F)⍴MEAN)÷(⍴F)⍴SD
XT←TRAIN⌿X ⋄ YT←TRAIN/Y
SIGMOID←{÷1+*-⍵} ⍝ monadic ÷ is 1÷; monadic * is e to the power
⍝ One step from weights ⍵. XT+.×⍵ is a matrix product: every row's w·x.
⍝ Minus YT gives each row's error p - y. A second product, errors +.×XT,
⍝ sums error × input for each weight, and ÷N turns the sums into means.
GRADIENT←{(((SIGMOID XT+.×⍵)-YT)+.×XT)÷N}
STEP←{⍵-RATE×GRADIENT ⍵} ⍝ all three weights move together
⍝ The power operator ⍣ applies STEP STEPS times, starting from 3⍴0: three
⍝ weights of 0.
W←(STEP⍣STEPS)3⍴0
P←SIGMOID XT+.×W
LOSS←-(+/(YT×⍟P)+(1-YT)×⍟1-P)÷N ⍝ ⍟ is the natural log
⍝ ~ is not: the test rows. Predict 1 when p ≥ 0.5, and count the matches.
PREDICTED←0.5≤SIGMOID((~TRAIN)⌿X)+.×W
CORRECT←+/PREDICTED=(~TRAIN)/Y
'train_count' SAY N
'test_count' SAY +/~TRAIN
'intercept' SAY6 W[1]
'weight_hours' SAY6 W[2]
'weight_sleep' SAY6 W[3]
'log_loss' SAY6 LOSS
'correct' SAY CORRECT
'accuracy' SAY6 CORRECT÷+/~TRAIN
)OFFBatch gradient descent is two matrix products. XT+.×⍵ gives every
training row’s w·x at once, and after the sigmoid and subtracting the
labels, errors+.×XT sums error × input for all three weights. So a
whole step is one line, ⍵-RATE×GRADIENT ⍵, with no loop over rows.
The loop over steps needs no ∇ function either. The power operator ⍣
applies a function a fixed number of times, so (STEP⍣STEPS)3⍴0 runs all
1000 steps from zero weights. Prepending the intercept’s 1 is 1,, which
adds a column to the left of a table.
# Logistic regression by gradient descent, in POSIX awk.
#
# Predicts whether a student passed from study hours and sleep. Each feature is
# rescaled with the training rows' mean and standard deviation, then the
# weights start at zero and take STEPS fixed steps of batch gradient descent
# on the average log-loss. Rows are stored as they are read; the fitting
# happens in END.
BEGIN { FS = ","; RATE = 0.5; STEPS = 1000 }
# The header line: find each column by name.
NR == 1 {
for (i = 1; i <= NF; i++) {
if ($i == "hours") c_hours = i
if ($i == "sleep") c_sleep = i
if ($i == "passed") c_lab = i
if ($i == "split") c_split = i
}
if (!c_hours || !c_sleep || !c_lab || !c_split) {
print "missing a column (hours, sleep, passed, split)" > "/dev/stderr"
failed = 1
exit 1
}
next
}
# Every other line goes into the training or the test arrays. x[i, j] is
# feature j of training row i; awk fakes 2-D arrays with a joined key.
$c_split == "train" {
n++
x[n, 1] = $c_hours + 0; x[n, 2] = $c_sleep + 0; y[n] = $c_lab + 0
next
}
$c_split == "test" {
m++
q[m, 1] = $c_hours + 0; q[m, 2] = $c_sleep + 0; qy[m] = $c_lab + 0
}
function sigmoid(z) { return 1 / (1 + exp(-z)) }
# w . (1, x1, x2) for already-scaled inputs.
function score(x1, x2) { return w[0] + w[1] * x1 + w[2] * x2 }
END {
if (failed) exit 1
# Training means and sample standard deviations (divide by n - 1), then
# scale the training rows in place. w[0] is the intercept, whose input
# is always 1, so it needs no column.
for (j = 1; j <= 2; j++) {
total = 0
for (i = 1; i <= n; i++) total += x[i, j]
mean[j] = total / n
ss = 0
for (i = 1; i <= n; i++) ss += (x[i, j] - mean[j]) ^ 2
sd[j] = sqrt(ss / (n - 1))
for (i = 1; i <= n; i++) x[i, j] = (x[i, j] - mean[j]) / sd[j]
}
w[0] = w[1] = w[2] = 0
for (step = 1; step <= STEPS; step++) {
# Gradient of the average log-loss: the mean of (p - y) times each
# input, all computed with the old weights before any update.
g0 = g1 = g2 = 0
for (i = 1; i <= n; i++) {
e = sigmoid(score(x[i, 1], x[i, 2])) - y[i]
g0 += e; g1 += e * x[i, 1]; g2 += e * x[i, 2]
}
w[0] -= RATE * g0 / n
w[1] -= RATE * g1 / n
w[2] -= RATE * g2 / n
}
loss = 0
for (i = 1; i <= n; i++) {
p = sigmoid(score(x[i, 1], x[i, 2]))
loss -= y[i] * log(p) + (1 - y[i]) * log(1 - p)
}
# Test rows are scaled with the training means and sds.
correct = 0
for (i = 1; i <= m; i++) {
p = sigmoid(score((q[i, 1] - mean[1]) / sd[1], (q[i, 2] - mean[2]) / sd[2]))
if ((p >= 0.5 ? 1 : 0) == qy[i]) correct++
}
printf "train_count %d\n", n
printf "test_count %d\n", m
printf "intercept %.6f\n", w[0]
printf "weight_hours %.6f\n", w[1]
printf "weight_sleep %.6f\n", w[2]
printf "log_loss %.6f\n", loss / n
printf "correct %d\n", correct
printf "accuracy %.6f\n", correct / m
}awk has no vectors or matrices, so every sum is a written-out loop. That
makes the update rule easy to see: one pass over the training rows adds up
g0, g1, and g2, and only after the pass do the three weights change.
exp, log, and sqrt are built in, which is all the math this needs.
The rows are kept in x[i, j]. awk has no true two-dimensional arrays; it
joins the two indexes with SUBSEP into one key. The fitting runs in END,
after the whole file is read. awk has no regression function, so there was
nothing to skip. It sticks to POSIX awk and runs on mawk and on macOS awk.
10 REM Logistic regression by gradient descent, in classic BASIC.
20 REM
30 REM Predicts whether a student passed from study hours and sleep. Each
40 REM feature is rescaled with the training rows' mean and standard
50 REM deviation, then the weights start at zero and take 1000 steps of
60 REM batch gradient descent on the average log-loss.
70 REM
80 REM run.sh passes the data path in the DATA environment variable,
90 REM because bwbasic can't read command-line arguments.
100 RATE = 0.5 : STEPS = 1000
110 REM A holds training rows and B test rows; column 1 is hours, 2 sleep.
120 DIM A(200, 2), YA(200), B(200, 2), YB(200), FL$(20)
130 REM DEF FN makes a one-line function; it reads the global weights.
140 DEF FNP(X1, X2) = 1 / (1 + EXP(-(W0 + W1 * X1 + W2 * X2)))
150 GOSUB 1000 : REM read the CSV into A, YA, B, and YB
160 J = 1 : GOSUB 2000 : REM scale hours
170 J = 2 : GOSUB 2000 : REM scale sleep
180 W0 = 0 : W1 = 0 : W2 = 0
190 FOR S = 1 TO STEPS
200 REM Gradient of the average log-loss: the mean of (p - y) times each
210 REM input, summed over every student before any weight moves.
220 G0 = 0 : G1 = 0 : G2 = 0
230 FOR I = 1 TO NA
240 E = FNP(A(I, 1), A(I, 2)) - YA(I)
250 G0 = G0 + E : G1 = G1 + E * A(I, 1) : G2 = G2 + E * A(I, 2)
260 NEXT I
270 W0 = W0 - RATE * G0 / NA : W1 = W1 - RATE * G1 / NA : W2 = W2 - RATE * G2 / NA
280 NEXT S
290 LOSS = 0
300 FOR I = 1 TO NA
310 P = FNP(A(I, 1), A(I, 2))
320 IF YA(I) = 1 THEN LOSS = LOSS - LOG(P) ELSE LOSS = LOSS - LOG(1 - P)
330 NEXT I
340 LOSS = LOSS / NA
350 HITS = 0
360 FOR I = 1 TO NB
370 IF FNP(B(I, 1), B(I, 2)) >= 0.5 THEN GUESS = 1 ELSE GUESS = 0
380 IF GUESS = YB(I) THEN HITS = HITS + 1
390 NEXT I
400 PRINT "train_count"; NA
410 PRINT "test_count"; NB
420 L$ = "intercept" : Y = W0 : GOSUB 3000
430 L$ = "weight_hours" : Y = W1 : GOSUB 3000
440 L$ = "weight_sleep" : Y = W2 : GOSUB 3000
450 L$ = "log_loss" : Y = LOSS : GOSUB 3000
460 PRINT "correct"; HITS
470 L$ = "accuracy" : Y = HITS / NB : GOSUB 3000
480 SYSTEM : REM without SYSTEM, bwbasic waits at its prompt
1000 REM Read the CSV into training rows (A, YA) and test rows (B, YB).
1010 F$ = ENVIRON$("DATA")
1020 OPEN "I", #1, F$
1030 LINE INPUT #1, L$
1040 GOSUB 4000 : REM find each column's position in the header
1050 FOR K = 1 TO NF
1060 IF FL$(K) = "hours" THEN CH = K
1070 IF FL$(K) = "sleep" THEN CS = K
1080 IF FL$(K) = "passed" THEN CP = K
1090 IF FL$(K) = "split" THEN CT = K
1100 NEXT K
1110 NA = 0 : NB = 0
1120 IF EOF(1) THEN 1200
1130 LINE INPUT #1, L$
1140 GOSUB 4000
1150 IF FL$(CT) <> "train" THEN 1180
1160 NA = NA + 1 : A(NA, 1) = VAL(FL$(CH)) : A(NA, 2) = VAL(FL$(CS)) : YA(NA) = VAL(FL$(CP))
1170 GOTO 1120
1180 NB = NB + 1 : B(NB, 1) = VAL(FL$(CH)) : B(NB, 2) = VAL(FL$(CS)) : YB(NB) = VAL(FL$(CP))
1190 GOTO 1120
1200 CLOSE #1
1210 RETURN
2000 REM Scale column J to mean 0 and (sample) standard deviation 1,
2010 REM using the training rows only, then shift the test rows the same way.
2020 T = 0
2030 FOR I = 1 TO NA
2040 T = T + A(I, J)
2050 NEXT I
2060 M = T / NA : T = 0
2070 FOR I = 1 TO NA
2080 T = T + (A(I, J) - M) ^ 2
2090 NEXT I
2100 SD = SQR(T / (NA - 1))
2110 FOR I = 1 TO NA
2120 A(I, J) = (A(I, J) - M) / SD
2130 NEXT I
2140 FOR I = 1 TO NB
2150 B(I, J) = (B(I, J) - M) / SD
2160 NEXT I
2170 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 4020The three weights are plain variables, W0, W1, and W2, and the
sigmoid is a one-line function, DEF FNP(X1, X2). It reads the global
weights each time it’s called, so it always predicts with the current
step’s weights. Training rows go in a 2-D array, A(200, 2), and test rows
in B, each sized up front with DIM. There’s no mean, regression, or
matrix function to avoid: every sum is a FOR loop, and the only
built-ins used are EXP, LOG, and SQR. The 1000 steps over 150
students take about 4.5 seconds. Swapping EXP for a square barely
changes that, so the interpreter’s loop, not the exponential, sets the
pace.
/* Logistic regression by gradient descent, using only the C standard library.
*
* Predicts whether a student passed from study hours and sleep. Each feature
* is rescaled with the training rows' mean and standard deviation, then the
* weights start at zero and take STEPS fixed steps of batch gradient descent
* on the average log-loss. */
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#define MAX_LINE 1024
#define MAX_FIELDS 16
#define K 3 /* inputs per row: 1 (for the intercept), hours, sleep */
#define LEARNING_RATE 0.5
#define STEPS 1000
/* One student: x[0] is always 1, so weight 0 acts as the intercept. */
struct example {
double x[K];
int y;
};
/* A growable array of examples; C has no built-in list. */
struct examples {
struct example *rows;
size_t n, capacity;
};
static void append(struct examples *set, struct example row) {
if (set->n == set->capacity) {
set->capacity = set->capacity ? set->capacity * 2 : 64;
set->rows = realloc(set->rows, set->capacity * sizeof *set->rows);
}
set->rows[set->n++] = row;
}
/* Cut a line at its commas into fields[], returning how many there are.
* strtok writes '\0' over each comma, so the fields point into line. */
static int split(char *line, char *fields[]) {
int count = 0;
for (char *f = strtok(line, ",\r\n"); f && count < MAX_FIELDS; f = strtok(NULL, ",\r\n"))
fields[count++] = f;
return count;
}
/* Return the 0-based position of name among the header fields. */
static int column(char *header[], int count, const char *name) {
for (int i = 0; i < count; i++)
if (strcmp(header[i], name) == 0) return i;
fprintf(stderr, "no column %s\n", name);
exit(1);
}
/* Read the CSV into training and test rows, sorted by the "split" column. */
static void read_rows(const char *path, struct examples *train, struct examples *test) {
FILE *file = fopen(path, "r");
if (!file) { perror(path); exit(1); }
char header[MAX_LINE], line[MAX_LINE];
char *fields[MAX_FIELDS];
if (!fgets(header, sizeof header, file)) { fprintf(stderr, "empty file\n"); exit(1); }
int count = split(header, fields);
int c_hours = column(fields, count, "hours"), c_sleep = column(fields, count, "sleep");
int c_passed = column(fields, count, "passed"), c_split = column(fields, count, "split");
while (fgets(line, sizeof line, file)) {
if (split(line, fields) < count) continue;
struct example row = {{1.0, strtod(fields[c_hours], NULL), strtod(fields[c_sleep], NULL)},
atoi(fields[c_passed])};
append(strcmp(fields[c_split], "train") == 0 ? train : test, row);
}
fclose(file);
}
static double sigmoid(double z) { return 1 / (1 + exp(-z)); }
/* The probability of passing: the sigmoid of the weighted sum of inputs. */
static double predict(const double w[K], const double x[K]) {
double z = 0;
for (int j = 0; j < K; j++) z += w[j] * x[j];
return sigmoid(z);
}
int main(int argc, char **argv) {
if (argc < 2) { fprintf(stderr, "usage: %s data.csv\n", argv[0]); return 1; }
struct examples train = {0}, test = {0}; /* {0} sets every field to zero */
read_rows(argv[1], &train, &test);
size_t n = train.n;
/* Scale each feature to mean 0 and (sample) standard deviation 1, using
* the training rows only, then apply the same shift to the test rows.
* Column 0 holds the ones, so the features are columns 1 and 2. */
for (int j = 1; j < K; j++) {
double total = 0, squares = 0;
for (size_t i = 0; i < n; i++) total += train.rows[i].x[j];
double mean = total / n;
for (size_t i = 0; i < n; i++) {
double d = train.rows[i].x[j] - mean;
squares += d * d;
}
double sd = sqrt(squares / (n - 1));
for (size_t i = 0; i < n; i++) train.rows[i].x[j] = (train.rows[i].x[j] - mean) / sd;
for (size_t i = 0; i < test.n; i++) test.rows[i].x[j] = (test.rows[i].x[j] - mean) / sd;
}
double w[K] = {0};
for (int step = 0; step < STEPS; step++) {
/* Gradient of the average log-loss: the mean of (p - y) times each
* input. It is summed over every student before any weight moves. */
double gradient[K] = {0};
for (size_t i = 0; i < n; i++) {
double error = predict(w, train.rows[i].x) - train.rows[i].y;
for (int j = 0; j < K; j++) gradient[j] += error * train.rows[i].x[j];
}
for (int j = 0; j < K; j++) w[j] -= LEARNING_RATE * gradient[j] / n;
}
double loss = 0;
for (size_t i = 0; i < n; i++) {
double p = predict(w, train.rows[i].x);
loss -= train.rows[i].y ? log(p) : log(1 - p);
}
loss /= n;
int correct = 0;
for (size_t i = 0; i < test.n; i++)
correct += (predict(w, test.rows[i].x) >= 0.5) == test.rows[i].y;
printf("train_count %zu\n", n);
printf("test_count %zu\n", test.n);
printf("intercept %.6f\n", w[0]);
printf("weight_hours %.6f\n", w[1]);
printf("weight_sleep %.6f\n", w[2]);
printf("log_loss %.6f\n", loss);
printf("correct %d\n", correct);
/* Cast before dividing, or 37 / 50 is integer division and gives 0. */
printf("accuracy %.6f\n", (double)correct / test.n);
free(train.rows);
free(test.rows);
return 0;
}Each student is a struct holding a fixed array of three inputs, the first
always 1 so that weight 0 acts as the intercept, plus the label. The
training and test rows go into two growable arrays that double with
realloc, and strtok cuts each CSV line at its commas.
Every step adds up the whole gradient in a nested loop over students and
weights, then moves all three weights at once. C has no vectors or matrix
products, so predict writes the weighted sum out as a loop. Watch the
last line: correct / test.n would be integer division, giving 0, so we
cast correct to double first.
// Logistic regression by gradient descent, using only the C++ standard
// library.
//
// Predicts whether a student passed from study hours and sleep. Each feature
// is rescaled with the training rows' mean and standard deviation, then the
// weights start at zero and take 1,000 steps of batch gradient descent on the
// average log-loss. C++ has no regression routine, so every step is a loop.
#include <algorithm>
#include <array>
#include <cmath>
#include <fstream>
#include <numeric>
#include <print>
#include <sstream>
#include <stdexcept>
#include <string>
#include <vector>
constexpr std::size_t K = 3; // inputs per row: 1 (for the intercept), hours, sleep
constexpr double LEARNING_RATE = 0.5;
constexpr int STEPS = 1000;
// std::array has its length in its type, so Inputs is always exactly three
// doubles, and unlike a C array it can be copied, assigned, and returned.
using Inputs = std::array<double, K>;
// One student: x[0] is always 1, so weight 0 acts as the intercept.
struct Example {
Inputs x;
int y;
};
// 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 0-based position of the named column in the header.
std::size_t column(const std::vector<std::string>& header, const std::string& name) {
auto found = std::find(header.begin(), header.end(), name);
if (found == header.end()) {
throw std::runtime_error("no column " + name);
}
return found - header.begin();
}
// Read the CSV into training and test rows, sorted by the "split" column,
// which holds text rather than a number.
void read_rows(const std::string& path, std::vector<Example>& train, std::vector<Example>& test) {
std::ifstream file(path);
std::string line;
std::getline(file, line);
std::vector<std::string> header = split(line);
std::size_t hours = column(header, "hours"), sleep = column(header, "sleep");
std::size_t passed = column(header, "passed"), which = column(header, "split");
while (std::getline(file, line)) {
if (line.empty()) {
continue;
}
std::vector<std::string> fields = split(line);
Example row{{1.0, std::stod(fields.at(hours)), std::stod(fields.at(sleep))},
std::stoi(fields.at(passed))};
(fields.at(which) == "train" ? train : test).push_back(row);
}
}
double sigmoid(double z) { return 1 / (1 + std::exp(-z)); }
// The probability of passing: the sigmoid of the weighted sum of inputs.
// inner_product multiplies the two arrays entry by entry and adds, starting
// from 0.0 so that it adds doubles.
double predict(const Inputs& w, const Inputs& x) {
return sigmoid(std::inner_product(w.begin(), w.end(), x.begin(), 0.0));
}
int main(int argc, char* argv[]) {
if (argc != 2) {
std::println(stderr, "usage: program data.csv");
return 1;
}
std::vector<Example> train, test;
read_rows(argv[1], train, test);
std::size_t n = train.size();
// Scale each feature to mean 0 and (sample) standard deviation 1, using
// the training rows only, then apply the same shift to the test rows.
// Column 0 holds the ones, so the features are columns 1 and 2. The &
// in "Example& row" makes row a reference, so the loop changes the
// stored rows rather than copies of them.
for (std::size_t j = 1; j < K; j++) {
double total = 0.0;
for (const Example& row : train) total += row.x[j];
double mean = total / n;
double squares = 0.0;
for (const Example& row : train) squares += (row.x[j] - mean) * (row.x[j] - mean);
double sd = std::sqrt(squares / (n - 1));
for (Example& row : train) row.x[j] = (row.x[j] - mean) / sd;
for (Example& row : test) row.x[j] = (row.x[j] - mean) / sd;
}
// Braces with nothing inside set every entry to zero.
Inputs w{};
for (int step = 0; step < STEPS; step++) {
// Gradient of the average log-loss: the mean of (p - y) times each
// input. It is summed over every student before any weight moves.
Inputs gradient{};
for (const Example& row : train) {
double error = predict(w, row.x) - row.y;
for (std::size_t j = 0; j < K; j++) gradient[j] += error * row.x[j];
}
for (std::size_t j = 0; j < K; j++) w[j] -= LEARNING_RATE * gradient[j] / n;
}
double loss = 0.0;
for (const Example& row : train) {
double p = predict(w, row.x);
loss -= row.y == 1 ? std::log(p) : std::log(1 - p);
}
loss /= n;
int correct = 0;
for (const Example& row : test) {
int guess = predict(w, row.x) >= 0.5 ? 1 : 0;
if (guess == row.y) correct++;
}
std::println("train_count {}", n);
std::println("test_count {}", test.size());
std::println("intercept {:.6f}", w[0]);
std::println("weight_hours {:.6f}", w[1]);
std::println("weight_sleep {:.6f}", w[2]);
std::println("log_loss {:.6f}", loss);
std::println("correct {}", correct);
// Convert before dividing, or 37 / 50 is integer division and gives 0.
std::println("accuracy {:.6f}", static_cast<double>(correct) / test.size());
}Each student is a struct holding a std::array of three inputs, the first
always 1 so that weight 0 acts as the intercept, plus the label. Two
std::vectors collect the training and test rows, so there’s no realloc
as in C, and std::inner_product writes the weighted sum in one call.
Each step sums the whole gradient into a fresh array, then moves all three
weights in place. A std::array has no arithmetic, so w - gradient won’t
compile and the update is a loop. Watch the & in for (Example& row : train): without it, the scaling loop would change copies and leave the rows
as they were.
*> Logistic regression by gradient descent, in COBOL.
*>
*> Predicts whether a student passed from study hours and sleep.
*> Each feature is rescaled with the training rows' mean and
*> sample standard deviation, then the weights start at zero and
*> take 1000 steps of batch gradient descent on the average
*> log-loss. Values are fixed-point decimals carried to 18
*> decimal places.
IDENTIFICATION DIVISION.
PROGRAM-ID. LOGISTIC-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-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-HOURS-COLUMN PIC 99 VALUE 0.
01 WS-SLEEP-COLUMN PIC 99 VALUE 0.
01 WS-PASSED-COLUMN PIC 99 VALUE 0.
01 WS-SPLIT-COLUMN PIC 99 VALUE 0.
*> Every row: whether it is a training row, its two features (raw,
*> then scaled in place), and its 0/1 label.
01 WS-N PIC 9(4) VALUE 0.
01 WS-TRAIN-COUNT PIC 9(4) VALUE 0.
01 WS-TEST-COUNT PIC 9(4) VALUE 0.
01 WS-ROWS.
05 WS-ROW OCCURS 5000 TIMES.
10 WS-IS-TRAIN PIC X.
10 WS-X1 PIC S9(6)V9(18).
10 WS-X2 PIC S9(6)V9(18).
10 WS-Y PIC 9.
*> Training means and sample standard deviations.
01 WS-MEAN-1 PIC S9(6)V9(18).
01 WS-MEAN-2 PIC S9(6)V9(18).
01 WS-SS-1 PIC S9(9)V9(18).
01 WS-SS-2 PIC S9(9)V9(18).
01 WS-SD-1 PIC S9(6)V9(18).
01 WS-SD-2 PIC S9(6)V9(18).
*> The weights (W0 is the intercept), the gradient sums, and one
*> row's probability and error.
01 WS-W0 PIC S9(6)V9(18) VALUE 0.
01 WS-W1 PIC S9(6)V9(18) VALUE 0.
01 WS-W2 PIC S9(6)V9(18) VALUE 0.
01 WS-G0 PIC S9(9)V9(18).
01 WS-G1 PIC S9(9)V9(18).
01 WS-G2 PIC S9(9)V9(18).
01 WS-P PIC S9(3)V9(18).
01 WS-E PIC S9(3)V9(18).
01 WS-LOSS PIC S9(9)V9(18).
01 WS-LEARNING-RATE PIC 9V9 VALUE 0.5.
01 WS-STEPS PIC 9(4) VALUE 1000.
01 I PIC 9(4).
01 S PIC 9(4).
01 WS-PREDICTED PIC 9.
01 WS-CORRECT PIC 9(4) VALUE 0.
*> Output: a name, then the value edited to six decimal places.
01 WS-NAME PIC X(20).
01 WS-RESULT PIC S9(9)V9(18).
01 WS-RESULT-OUT PIC -(9)9.9(6).
01 WS-COUNT-OUT PIC Z(8)9.
PROCEDURE DIVISION.
MAIN.
ACCEPT WS-FILE-NAME FROM ARGUMENT-VALUE
OPEN INPUT DATA-FILE
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
EVALUATE WS-FIELD(I)
WHEN "hours" MOVE I TO WS-HOURS-COLUMN
WHEN "sleep" MOVE I TO WS-SLEEP-COLUMN
WHEN "passed" MOVE I TO WS-PASSED-COLUMN
WHEN "split" MOVE I TO WS-SPLIT-COLUMN
END-EVALUATE
END-PERFORM
IF WS-HOURS-COLUMN = 0 OR WS-SLEEP-COLUMN = 0
OR WS-PASSED-COLUMN = 0 OR WS-SPLIT-COLUMN = 0
DISPLAY "missing column" 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-ROW
END-READ
END-PERFORM
CLOSE DATA-FILE
PERFORM SCALE-FEATURES
*> Gradient descent: each step computes the whole gradient with
*> the old weights, then moves all three weights together.
PERFORM VARYING S FROM 1 BY 1 UNTIL S > WS-STEPS
MOVE 0 TO WS-G0 WS-G1 WS-G2
PERFORM VARYING I FROM 1 BY 1 UNTIL I > WS-N
IF WS-IS-TRAIN(I) = "Y"
PERFORM SIGMOID
COMPUTE WS-E = WS-P - WS-Y(I)
ADD WS-E TO WS-G0
COMPUTE WS-G1 = WS-G1 + WS-E * WS-X1(I)
COMPUTE WS-G2 = WS-G2 + WS-E * WS-X2(I)
END-IF
END-PERFORM
COMPUTE WS-W0 ROUNDED = WS-W0
- WS-LEARNING-RATE * WS-G0 / WS-TRAIN-COUNT
COMPUTE WS-W1 ROUNDED = WS-W1
- WS-LEARNING-RATE * WS-G1 / WS-TRAIN-COUNT
COMPUTE WS-W2 ROUNDED = WS-W2
- WS-LEARNING-RATE * WS-G2 / WS-TRAIN-COUNT
END-PERFORM
*> Average training log-loss, and test accuracy, with the final
*> weights.
MOVE 0 TO WS-LOSS
PERFORM VARYING I FROM 1 BY 1 UNTIL I > WS-N
PERFORM SIGMOID
IF WS-IS-TRAIN(I) = "Y"
IF WS-Y(I) = 1
COMPUTE WS-LOSS = WS-LOSS - FUNCTION LOG(WS-P)
ELSE
COMPUTE WS-LOSS =
WS-LOSS - FUNCTION LOG(1 - WS-P)
END-IF
ELSE
IF WS-P >= 0.5
MOVE 1 TO WS-PREDICTED
ELSE
MOVE 0 TO WS-PREDICTED
END-IF
IF WS-PREDICTED = WS-Y(I)
ADD 1 TO WS-CORRECT
END-IF
END-IF
END-PERFORM
MOVE WS-TRAIN-COUNT TO WS-COUNT-OUT
DISPLAY "train_count " FUNCTION TRIM(WS-COUNT-OUT)
MOVE WS-TEST-COUNT TO WS-COUNT-OUT
DISPLAY "test_count " FUNCTION TRIM(WS-COUNT-OUT)
MOVE "intercept" TO WS-NAME
MOVE WS-W0 TO WS-RESULT
PERFORM PRINT-RESULT
MOVE "weight_hours" TO WS-NAME
MOVE WS-W1 TO WS-RESULT
PERFORM PRINT-RESULT
MOVE "weight_sleep" TO WS-NAME
MOVE WS-W2 TO WS-RESULT
PERFORM PRINT-RESULT
MOVE "log_loss" TO WS-NAME
COMPUTE WS-RESULT ROUNDED = WS-LOSS / WS-TRAIN-COUNT
PERFORM PRINT-RESULT
MOVE WS-CORRECT TO WS-COUNT-OUT
DISPLAY "correct " FUNCTION TRIM(WS-COUNT-OUT)
MOVE "accuracy" TO WS-NAME
COMPUTE WS-RESULT ROUNDED = WS-CORRECT / WS-TEST-COUNT
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, and count it as training or test.
ADD-ROW.
PERFORM SPLIT-LINE
ADD 1 TO WS-N
IF WS-FIELD(WS-SPLIT-COLUMN) = "train"
MOVE "Y" TO WS-IS-TRAIN(WS-N)
ADD 1 TO WS-TRAIN-COUNT
ELSE
MOVE "N" TO WS-IS-TRAIN(WS-N)
ADD 1 TO WS-TEST-COUNT
END-IF
COMPUTE WS-X1(WS-N) =
FUNCTION NUMVAL(WS-FIELD(WS-HOURS-COLUMN))
COMPUTE WS-X2(WS-N) =
FUNCTION NUMVAL(WS-FIELD(WS-SLEEP-COLUMN))
COMPUTE WS-Y(WS-N) =
FUNCTION NUMVAL(WS-FIELD(WS-PASSED-COLUMN)).
*> Rescale both features of every row with the training rows'
*> mean and sample standard deviation (dividing by n - 1).
SCALE-FEATURES.
MOVE 0 TO WS-MEAN-1 WS-MEAN-2 WS-SS-1 WS-SS-2
PERFORM VARYING I FROM 1 BY 1 UNTIL I > WS-N
IF WS-IS-TRAIN(I) = "Y"
ADD WS-X1(I) TO WS-MEAN-1
ADD WS-X2(I) TO WS-MEAN-2
END-IF
END-PERFORM
COMPUTE WS-MEAN-1 ROUNDED = WS-MEAN-1 / WS-TRAIN-COUNT
COMPUTE WS-MEAN-2 ROUNDED = WS-MEAN-2 / WS-TRAIN-COUNT
PERFORM VARYING I FROM 1 BY 1 UNTIL I > WS-N
IF WS-IS-TRAIN(I) = "Y"
COMPUTE WS-SS-1 =
WS-SS-1 + (WS-X1(I) - WS-MEAN-1) ** 2
COMPUTE WS-SS-2 =
WS-SS-2 + (WS-X2(I) - WS-MEAN-2) ** 2
END-IF
END-PERFORM
COMPUTE WS-SD-1 ROUNDED =
FUNCTION SQRT(WS-SS-1 / (WS-TRAIN-COUNT - 1))
COMPUTE WS-SD-2 ROUNDED =
FUNCTION SQRT(WS-SS-2 / (WS-TRAIN-COUNT - 1))
PERFORM VARYING I FROM 1 BY 1 UNTIL I > WS-N
COMPUTE WS-X1(I) ROUNDED =
(WS-X1(I) - WS-MEAN-1) / WS-SD-1
COMPUTE WS-X2(I) ROUNDED =
(WS-X2(I) - WS-MEAN-2) / WS-SD-2
END-PERFORM.
*> WS-P = sigmoid(w . x) for row I, with the current weights.
SIGMOID.
COMPUTE WS-P ROUNDED = 1 / (1 + FUNCTION EXP(
-(WS-W0 + WS-W1 * WS-X1(I) + WS-W2 * WS-X2(I)))).
*> 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).We expected fixed-point decimals to drift over 1000 steps, but they don’t.
With 18 decimal places, PIC S9(6)V9(18), and ROUNDED on every weight
update, the final weights and log-loss agree with Python’s doubles to all
14 decimal places we printed, so there is no need for COMP-2
floating point. The cost is speed: the 150,000 calls to FUNCTION EXP
take about 10 seconds. Without them, the
loop runs in under 0.2 seconds.
Each step is a PERFORM VARYING over the rows that adds error × input
into three sums. At 260 lines, it’s three and a half times the Python.
GnuCOBOL has MEAN and STANDARD-DEVIATION functions, but no regression.
;;;; Logistic regression by gradient descent, using only standard Common Lisp.
;;;;
;;;; Predicts whether a student passed from study hours and sleep. Each
;;;; feature is rescaled with the training rows' mean and standard deviation,
;;;; then the weights start at zero and take 1,000 steps of batch gradient
;;;; descent on the average log-loss. The students are a list of rows and the
;;;; weights are a list of three numbers: each step builds a new list.
(defparameter *feature-names* '("hours" "sleep"))
(defparameter *label* "passed")
(defparameter *learning-rate* 0.5)
(defparameter *steps* 1000)
;; By default the reader turns 7.8 into a single float, which keeps only
;; about 7 digits. Read doubles, both from the data and in this file.
;; 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-rows (path)
"The CSV file at PATH as a list of rows. Each row is a function from a
column name to that row's field, as a string: (funcall row \"hours\")."
(with-open-file (in path)
(let ((header (split-on-comma (read-line in))))
(flet ((row-lookup (fields)
(lambda (name)
(nth (or (position name header :test #'string=)
(error "no column ~a" name))
fields))))
(loop for line = (read-line in nil)
while line
unless (string= line "")
collect (row-lookup (split-on-comma line)))))))
(defun features (row)
"The feature values of ROW, as numbers."
(mapcar (lambda (name) (parse-number (funcall row name))) *feature-names*))
(defun sigmoid (z)
(/ 1 (+ 1 (exp (- z)))))
(defun dot (u v)
"The sum of the products of U and V, entry by entry."
(reduce #'+ (mapcar #'* u v)))
(defun mean (xs)
(/ (reduce #'+ xs) (length xs)))
(defun sd (xs)
"Sample standard deviation: divide by n - 1."
(let ((m (mean xs)))
(sqrt (/ (reduce #'+ (mapcar (lambda (x) (expt (- x m) 2)) xs))
(1- (length xs))))))
(defun transpose (rows)
"The columns of ROWS. (mapcar #'list row1 row2 ...) collects the first
entry of every row, then the second, and so on; apply passes the rows as
those arguments."
(apply #'mapcar #'list rows))
(defun step-weights (xs ys weights)
"One step of gradient descent. The gradient for weight j is the mean of
(p - y) times input j, summed over every student before any weight moves,
and the result is a new list of weights."
(let* ((errors (mapcar (lambda (x y) (- (sigmoid (dot weights x)) y)) xs ys))
(gradient (mapcar (lambda (column) (/ (dot errors column) (length xs)))
(transpose xs))))
(mapcar (lambda (w g) (- w (* *learning-rate* g))) weights gradient)))
(let* ((rows (read-rows (second sb-ext:*posix-argv*)))
(train (remove-if-not (lambda (r) (string= (funcall r "split") "train")) rows))
(test (remove-if (lambda (r) (string= (funcall r "split") "train")) rows))
;; Scale with the training rows only.
(columns (transpose (mapcar #'features train)))
(means (mapcar #'mean columns))
(sds (mapcar #'sd columns)))
(flet ((scale (x)
;; The leading 1 multiplies the intercept.
(cons 1.0 (mapcar (lambda (v m s) (/ (- v m) s)) x means sds)))
(label (row)
(parse-number (funcall row *label*))))
(let* ((xs (mapcar (lambda (r) (scale (features r))) train))
(ys (mapcar #'label train))
(weights (list 0.0 0.0 0.0)))
;; Each pass replaces the whole list; no weight is changed in place.
(loop repeat *steps*
do (setf weights (step-weights xs ys weights)))
(let ((log-loss (- (/ (reduce #'+ (mapcar (lambda (x y)
(let ((p (sigmoid (dot weights x))))
(+ (* y (log p))
(* (- 1 y) (log (- 1 p))))))
xs ys))
(length xs))))
(correct (count-if (lambda (r)
(= (if (>= (sigmoid (dot weights (scale (features r)))) 0.5)
1 0)
(label r)))
test)))
;; ~d prints an integer, ~,6f a number with six decimals, ~% a newline.
(format t "train_count ~d~%" (length train))
(format t "test_count ~d~%" (length test))
;; destructuring-bind names the three weights in the list.
(destructuring-bind (intercept hours sleep) weights
(format t "intercept ~,6f~%" intercept)
(format t "weight_hours ~,6f~%" hours)
(format t "weight_sleep ~,6f~%" sleep))
(format t "log_loss ~,6f~%" log-loss)
(format t "correct ~d~%" correct)
(format t "accuracy ~,6f~%" (/ correct (length test)))))))The training rows are a list of lists, and the weights are a list of three
numbers. step-weights works out every student’s error from the old
weights, then builds a new list with mapcar, so every weight moves from
the same gradient and nothing is changed in place. (apply #'mapcar #'list rows) turns rows into columns, so each weight’s gradient is one dot.
Watch the decimals in the file itself. The compiler reads the whole file
before the program runs, so a plain setf of *read-default-float-format*
would leave a bare 0.5 there a single float; eval-when makes it apply
while compiling, too. And watch *features*: it’s already a standard
variable, so naming our list that stops the compile.
# Logistic regression by gradient descent, using only Elixir's standard
# library.
#
# Predicts whether a student passed from study hours and sleep. Each feature
# is rescaled with the training rows' mean and standard deviation, then the
# weights start at zero and take 1,000 steps of batch gradient descent on the
# average log-loss. Data never changes in place, so each step builds a new
# list of weights from the old one.
defmodule Logistic do
@learning_rate 0.5 # @name is a constant, fixed at compile time
@steps 1000
# Each row becomes {[hours, sleep], passed, split}.
def read_rows(path) do
[header | lines] = path |> File.read!() |> String.split("\n", trim: true)
names = String.split(header, ",")
col = fn name -> Enum.find_index(names, &(&1 == name)) || raise "column not found: #{name}" end
[hours, sleep, passed, split] = Enum.map(~w(hours sleep passed split), col)
Enum.map(lines, fn line ->
f = line |> String.split(",") |> List.to_tuple() # elem(f, i) reads field i
{[number(elem(f, hours)), number(elem(f, sleep))], String.to_integer(elem(f, passed)), elem(f, split)}
end)
end
# Float.parse returns {number, rest}; "7" parses as 7.0 too.
defp number(text) do
{value, ""} = Float.parse(text)
value
end
def sigmoid(z), do: 1 / (1 + :math.exp(-z))
# zip_reduce walks two lists side by side, carrying a running total.
def dot(u, v), do: Enum.zip_reduce(u, v, 0.0, fn a, b, acc -> acc + a * b end)
def mean(xs), do: Enum.sum(xs) / length(xs)
# Sample standard deviation: divide by n - 1.
def sd(xs) do
m = mean(xs)
:math.sqrt((xs |> Enum.map(&((&1 - m) ** 2)) |> Enum.sum()) / (length(xs) - 1))
end
# The columns of a list of rows: zip_with gathers the first of every row,
# then the second, and so on.
def columns(rows), do: Enum.zip_with(rows, & &1)
# One step. Every error comes from the old weights before any weight moves,
# so every weight moves from the same gradient: the mean of (p - y) times
# its input.
def step(xs, cols, ys, weights) do
n = length(xs)
errors = Enum.zip_with(xs, ys, fn x, y -> sigmoid(dot(weights, x)) - y end)
Enum.zip_with(weights, cols, fn w, col -> w - @learning_rate * dot(errors, col) / n end)
end
def descend(xs, ys) do
cols = columns(xs)
# reduce threads the weights through 1,000 steps; _ ignores the step number.
Enum.reduce(1..@steps, [0.0, 0.0, 0.0], fn _, weights -> step(xs, cols, ys, weights) end)
end
def say(name, value), do: IO.puts("#{name} #{:erlang.float_to_binary(value, decimals: 6)}")
end
[path] = System.argv()
# split_with divides the rows into those that pass the test and the rest.
{train, test} = path |> Logistic.read_rows() |> Enum.split_with(fn {_, _, split} -> split == "train" end)
# Scale with the training rows only. The leading 1 multiplies the intercept.
features = Enum.map(train, fn {x, _, _} -> x end)
ys = Enum.map(train, fn {_, y, _} -> y end)
cols = Logistic.columns(features)
means = Enum.map(cols, &Logistic.mean/1)
sds = Enum.map(cols, &Logistic.sd/1)
scale = fn x -> [1.0 | Enum.zip_with([x, means, sds], fn [v, m, s] -> (v - m) / s end)] end
xs = Enum.map(features, scale)
weights = Logistic.descend(xs, ys)
[intercept, weight_hours, weight_sleep] = weights
n = length(xs)
losses =
Enum.zip_with(xs, ys, fn x, y ->
p = Logistic.sigmoid(Logistic.dot(weights, x))
y * :math.log(p) + (1 - y) * :math.log(1 - p)
end)
log_loss = -Enum.sum(losses) / n
predict = fn x -> if Logistic.sigmoid(Logistic.dot(weights, scale.(x))) >= 0.5, do: 1, else: 0 end
correct = Enum.count(test, fn {x, y, _} -> predict.(x) == y end)
IO.puts("train_count #{n}")
IO.puts("test_count #{length(test)}")
Logistic.say("intercept", intercept)
Logistic.say("weight_hours", weight_hours)
Logistic.say("weight_sleep", weight_sleep)
Logistic.say("log_loss", log_loss)
IO.puts("correct #{correct}")
Logistic.say("accuracy", correct / length(test))The weights are a list, and no step changes it. step works out every
error from the old weights first, then Enum.zip_with pairs each weight
with its column of inputs to build a new list, so every weight moves from
the same gradient. Enum.reduce over the range 1..1000 threads the
weights through the steps and hands back the last list.
Elixir has no matrix module (Matrix isn’t defined) and no Enum.mean,
so sigmoid, dot, mean, and sd are short functions, with exp and
log from Erlang’s :math. Watch Enum.zip_reduce, which dot uses:
given lists of different lengths, it quietly stops at the shorter one
instead of raising. Enum.split_with divides the training and test rows
in one pass.
! Logistic regression by gradient descent, using only standard Fortran.
!
! Predicts whether a student passed from study hours and sleep. Each feature is
! rescaled with the training rows' mean and standard deviation, then the weights
! start at zero and take STEPS fixed steps of batch gradient descent on the
! average log-loss.
program logistic_regression
implicit none
integer, parameter :: dp = kind(1.0d0) ! double precision
real(dp), parameter :: learning_rate = 0.5_dp
integer, parameter :: steps = 1000
! Each row of x is [1, hours, sleep]; the leading 1 multiplies the intercept.
real(dp), allocatable :: x(:, :), x_test(:, :), y(:), y_test(:), p(:)
real(dp) :: means(2), sds(2), weights(3), log_loss
integer :: n, step, j, correct
call read_rows(x, y, x_test, y_test)
n = size(y)
! Scale each feature to mean 0 and (sample) standard deviation 1, using the
! training rows only. Column j + 1 holds feature j (column 1 is the ones).
do j = 1, 2
means(j) = sum(x(:, j + 1)) / n
sds(j) = sqrt(sum((x(:, j + 1) - means(j))**2) / (n - 1))
x(:, j + 1) = (x(:, j + 1) - means(j)) / sds(j)
x_test(:, j + 1) = (x_test(:, j + 1) - means(j)) / sds(j)
end do
weights = 0
do step = 1, steps
! matmul(x, weights) is every row's score at once. The whole gradient,
! the mean of (p - y) times each input, is found before any weight moves.
p = sigmoid(matmul(x, weights))
weights = weights - learning_rate * matmul(transpose(x), p - y) / n
end do
p = sigmoid(matmul(x, weights))
log_loss = -sum(y * log(p) + (1 - y) * log(1 - p)) / n
! merge(a, b, mask) picks a where mask is true: here, the predicted label.
! nint rounds the 0.0 and 1.0 labels to integers for an exact comparison.
correct = count(merge(1, 0, sigmoid(matmul(x_test, weights)) >= 0.5_dp) &
== nint(y_test))
print "(a, 1x, i0)", "train_count", n
print "(a, 1x, i0)", "test_count", size(y_test)
call say("intercept", weights(1)) ! Fortran arrays count from 1
call say("weight_hours", weights(2))
call say("weight_sleep", weights(3))
call say("log_loss", log_loss)
print "(a, 1x, i0)", "correct", correct
call say("accuracy", real(correct, dp) / size(y_test))
contains
! An elemental function works on one number but also applies itself to
! every element of an array, so sigmoid(matmul(x, weights)) is a vector.
elemental real(dp) function sigmoid(z)
real(dp), intent(in) :: z
sigmoid = 1 / (1 + exp(-z))
end function sigmoid
! Read the CSV into training and test rows, sorted by the "split" column.
subroutine read_rows(x, y, x_test, y_test)
real(dp), allocatable, intent(out) :: x(:, :), y(:), x_test(:, :), y_test(:)
character(len=256) :: path, header, line, text
integer :: unit, status, c_hours, c_sleep, c_passed, c_split
real(dp) :: row(3), label
call get_command_argument(1, path)
open(newunit=unit, file=path, status="old", action="read")
read(unit, "(a)") header
c_hours = field_index(header, "hours")
c_sleep = field_index(header, "sleep")
c_passed = field_index(header, "passed")
c_split = field_index(header, "split")
allocate(x(0, 3), y(0), x_test(0, 3), y_test(0))
do
read(unit, "(a)", iostat=status) line
if (status /= 0) exit
row(1) = 1
text = field(line, c_hours)
read(text, *) row(2) ! an internal read parses the number
text = field(line, c_sleep)
read(text, *) row(3)
text = field(line, c_passed)
read(text, *) label
if (field(line, c_split) == "train") then
call append(x, y, row, label)
else
call append(x_test, y_test, row, label)
end if
end do
close(unit)
end subroutine read_rows
! Add one row to a matrix and its label to a vector. Fortran stores
! matrices column by column, so we rebuild the matrix one row longer.
subroutine append(x, y, row, label)
real(dp), allocatable, intent(inout) :: x(:, :), y(:)
real(dp), intent(in) :: row(:), label
real(dp), allocatable :: grown(:, :)
allocate(grown(size(x, 1) + 1, size(x, 2)))
grown(:size(x, 1), :) = x
grown(size(x, 1) + 1, :) = row
call move_alloc(grown, x) ! x takes grown's storage
y = [y, label] ! grows the vector by one
end subroutine append
! 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 logistic_regressionThe training rows are a matrix x whose first column is all ones, for the
intercept. matmul and transpose are intrinsics, so each step is two
lines: p = sigmoid(matmul(x, weights)), then
weights = weights - learning_rate * matmul(transpose(x), p - y) / n. The
whole gradient is computed before any weight changes.
sigmoid is declared elemental: written for one number, it also applies
itself to every element of an array. Where Python’s math.exp raises an
error on overflow, gfortran’s exp(1000) returns Infinity, so sigmoid
gives 0 for very negative scores. Fortran has no logistic regression
intrinsic, and LAPACK isn’t installed, so there was nothing to skip.
// Logistic regression by gradient descent, using only the standard library.
//
// Predicts whether a student passed from study hours and sleep. Each feature
// is rescaled with the training rows' mean and standard deviation, then the
// weights start at zero and take 1000 fixed steps of batch gradient descent
// on the average log-loss.
package main
import (
"encoding/csv"
"fmt"
"math"
"os"
"strconv"
)
const (
label = "passed"
learningRate = 0.5
steps = 1000
)
var features = []string{"hours", "sleep"}
type example struct {
x []float64
y float64
}
func read(path string) (train, test []example, err error) {
f, err := os.Open(path)
if err != nil {
return nil, nil, err
}
defer f.Close()
records, err := csv.NewReader(f).ReadAll()
if err != nil {
return nil, nil, err
}
// Map each header name to its column number.
column := map[string]int{}
for i, name := range records[0] {
column[name] = i
}
for _, record := range records[1:] {
var e example
for _, name := range features {
v, err := strconv.ParseFloat(record[column[name]], 64)
if err != nil {
return nil, nil, err
}
e.x = append(e.x, v)
}
if e.y, err = strconv.ParseFloat(record[column[label]], 64); err != nil {
return nil, nil, err
}
if record[column["split"]] == "train" {
train = append(train, e)
} else {
test = append(test, e)
}
}
return train, test, nil
}
func sigmoid(z float64) float64 {
return 1 / (1 + math.Exp(-z))
}
func dot(w, x []float64) float64 {
total := 0.0
for j := range w {
total += w[j] * x[j]
}
return total
}
func main() {
train, test, err := read(os.Args[1])
if err != nil {
fmt.Fprintln(os.Stderr, err)
os.Exit(1)
}
n := float64(len(train))
// Scale each feature to mean 0 and (sample) standard deviation 1, using
// the training rows only.
means := make([]float64, len(features))
sds := make([]float64, len(features))
for j := range features {
for _, e := range train {
means[j] += e.x[j]
}
means[j] /= n
for _, e := range train {
d := e.x[j] - means[j]
sds[j] += d * d
}
sds[j] = math.Sqrt(sds[j] / (n - 1))
}
// The leading 1 multiplies the intercept.
scale := func(x []float64) []float64 {
scaled := []float64{1}
for j, v := range x {
scaled = append(scaled, (v-means[j])/sds[j])
}
return scaled
}
xs := make([][]float64, len(train))
for i, e := range train {
xs[i] = scale(e.x)
}
weights := make([]float64, len(features)+1)
for range steps {
// Gradient of the average log-loss: the mean of (p - y) times each
// input, computed in full before any weight changes.
gradient := make([]float64, len(weights))
for i, x := range xs {
e := sigmoid(dot(weights, x)) - train[i].y
for j := range gradient {
gradient[j] += e * x[j]
}
}
for j := range weights {
weights[j] -= learningRate * gradient[j] / n
}
}
logLoss := 0.0
for i, x := range xs {
p, y := sigmoid(dot(weights, x)), train[i].y
logLoss -= y*math.Log(p) + (1-y)*math.Log(1-p)
}
logLoss /= n
correct := 0
for _, e := range test {
predicted := 0.0
if sigmoid(dot(weights, scale(e.x))) >= 0.5 {
predicted = 1
}
if predicted == e.y {
correct++
}
}
fmt.Printf("train_count %d\n", len(train))
fmt.Printf("test_count %d\n", len(test))
fmt.Printf("intercept %.6f\n", weights[0])
fmt.Printf("weight_hours %.6f\n", weights[1])
fmt.Printf("weight_sleep %.6f\n", weights[2])
fmt.Printf("log_loss %.6f\n", logLoss)
fmt.Printf("correct %d\n", correct)
fmt.Printf("accuracy %.6f\n", float64(correct)/float64(len(test)))
}Each step makes two passes. The first goes over the 150 training rows and
adds every row’s error times its inputs into a fresh gradient slice; the
second moves the weights. Updating weights in place is safe only because
the whole gradient exists first. Move the update into the row loop and it
becomes stochastic gradient descent, which gives an intercept of 0.758788
instead of 0.753594 and fails the test.
The training loop is for range steps, since Go can range over an integer.
Test rows go through the same scale closure, so they use the training mean
and standard deviation. The standard library has no logistic regression, so
there was no built-in to skip.
-- Logistic regression by gradient descent, using only Haskell's base library.
--
-- Predicts whether a student passed from study hours and sleep. Each feature
-- is rescaled with the training rows' mean and standard deviation, then the
-- weights start at zero and take 1,000 steps of batch gradient descent on the
-- average log-loss. A matrix here is a list of rows, and the weights are a
-- list too: each step builds a new list from the old one.
import Data.List (elemIndex, transpose)
import System.Environment (getArgs)
import Text.Printf (printf)
learningRate :: Double
learningRate = 0.5
steps :: Int
steps = 1000
-- Split a line at every comma, as in the descriptive statistics example.
splitOnComma :: String -> [String]
splitOnComma line = case break (== ',') line of
(field, ',' : rest) -> field : splitOnComma rest
(field, _) -> [field]
-- A row is a function from a column name to that row's field, so
-- row "hours" gives the hours. The header says where each name sits.
type Row = String -> String
readRows :: String -> [Row]
readRows text = case lines text of
header : rows -> [field (splitOnComma row) | row <- rows, not (null row)]
where
names = splitOnComma header
field values name = case elemIndex name names of
Just index -> values !! index
Nothing -> error ("no column " ++ name)
[] -> error "empty file"
sigmoid :: Double -> Double
sigmoid z = 1 / (1 + exp (-z))
dot :: [Double] -> [Double] -> Double
dot u v = sum (zipWith (*) u v)
mean :: [Double] -> Double
mean xs = sum xs / fromIntegral (length xs)
-- Sample standard deviation: divide by n - 1.
sd :: [Double] -> Double
sd xs = sqrt (sum [(x - m) ^ (2 :: Int) | x <- xs] / fromIntegral (length xs - 1))
where
m = mean xs
-- One step of gradient descent. The gradient for weight j is the mean of
-- (p - y) times input j, summed over every student before any weight moves.
-- transpose turns the rows into columns, so each column gives one weight's sum.
step :: [[Double]] -> [Double] -> [Double] -> [Double]
step xs ys weights = zipWith (\w g -> w - learningRate * g) weights gradient
where
errors = [sigmoid (dot weights x) - y | (x, y) <- zip xs ys]
gradient = [dot errors column / fromIntegral (length xs) | column <- transpose xs]
main :: IO ()
main = do
[path] <- getArgs
text <- readFile path
let rows = readRows text
features r = [read (r "hours"), read (r "sleep")] :: [Double]
label r = read (r "passed") :: Double
train = [r | r <- rows, r "split" == "train"]
test = [r | r <- rows, r "split" /= "train"]
-- Scale with the training rows only. The leading 1 multiplies the intercept.
columns = transpose (map features train)
means = map mean columns
sds = map sd columns
scale x = 1 : zipWith3 (\v m s -> (v - m) / s) x means sds
xs = map (scale . features) train
ys = map label train
-- iterate builds the endless list w0, step w0, step (step w0), ...
-- and !! picks the entry after 1,000 steps. Haskell is lazy, so no
-- sum is worked out until a weight is printed, and then each step
-- calls for the one before. A chain 1,000 steps long is small enough
-- that we don't need to force the weights any sooner.
weights = iterate (step xs ys) [0, 0, 0] !! steps
n = fromIntegral (length xs) :: Double
logLoss =
-sum [y * log p + (1 - y) * log (1 - p) | (x, y) <- zip xs ys, let p = sigmoid (dot weights x)] / n
predict x = if sigmoid (dot weights (scale x)) >= 0.5 then 1 else 0
correct = length [r | r <- test, predict (features r) == label r]
printf "train_count %d\n" (length train)
printf "test_count %d\n" (length test)
-- The weights are a list of three, so a pattern takes them apart by name.
case weights of
[intercept, hours, sleep] -> do
printf "intercept %.6f\n" intercept
printf "weight_hours %.6f\n" hours
printf "weight_sleep %.6f\n" sleep
_ -> error "expected three weights"
printf "log_loss %.6f\n" logLoss
printf "correct %d\n" correct
printf "accuracy %.6f\n" (fromIntegral correct / fromIntegral (length test) :: Double)There’s no loop and no weight that changes. step takes one list of
weights and returns the next, and iterate (step xs ys) [0, 0, 0] !! 1000
applies it 1,000 times: iterate builds the endless list of every step’s
weights, and laziness means only the first 1,001 are ever made. Because
step computes all the errors from the old weights before building the new
list, every weight moves from the same gradient.
The training rows are a list of lists. transpose turns them into columns,
so each weight’s gradient is one dot of the errors with a column. Each
CSV row becomes a function from a column name to its field, so
r "hours" reads like a lookup.
// Logistic regression by gradient descent, using only the Java standard
// library.
//
// Predicts whether a student passed from study hours and sleep. Each feature
// is rescaled with the training rows' mean and standard deviation, then the
// weights start at zero and take 1000 fixed steps of batch gradient descent
// on the average log-loss. The JDK has no regression or optimizer, so every
// step is a loop.
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 double LEARNING_RATE = 0.5;
static final int STEPS = 1000;
// The named column of a CSV file's lines, as strings, since "split" holds
// words; the caller parses the numbers.
static String[] 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);
}
String[] values = new String[lines.size() - 1];
for (int row = 1; row < lines.size(); row++) {
values[row - 1] = lines.get(row).split(",")[index];
}
return values;
}
static double sigmoid(double z) {
return 1 / (1 + Math.exp(-z));
}
static double dot(double[] u, double[] v) {
double total = 0;
for (int j = 0; j < u.length; j++) {
total += u[j] * v[j];
}
return total;
}
// 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]));
String[] hours = readColumn(lines, "hours");
String[] sleep = readColumn(lines, "sleep");
String[] passed = readColumn(lines, "passed");
String[] split = readColumn(lines, "split");
// Count the training rows first, since a Java array's length is
// fixed when it is made.
int n = 0;
for (String s : split) {
// == would compare the String objects, not their letters.
if (s.equals("train")) n++;
}
int nTest = split.length - n;
// Each row of x is {1, hours, sleep}; the leading 1 multiplies the
// intercept. A double[][] is an array of row arrays.
double[][] train = new double[n][];
double[] y = new double[n];
double[][] test = new double[nTest][];
double[] yTest = new double[nTest];
for (int i = 0, a = 0, b = 0; i < split.length; i++) {
double[] x = {1, Double.parseDouble(hours[i]), Double.parseDouble(sleep[i])};
if (split[i].equals("train")) {
train[a] = x;
y[a++] = Double.parseDouble(passed[i]);
} else {
test[b] = x;
yTest[b++] = Double.parseDouble(passed[i]);
}
}
// Scale each feature (columns 1 and 2) to mean 0 and (sample) standard
// deviation 1, using the training rows only, then apply the same
// shift and stretch to the test rows.
for (int j = 1; j <= 2; j++) {
double sum = 0;
for (double[] x : train) sum += x[j];
double mean = sum / n;
double squares = 0;
for (double[] x : train) squares += (x[j] - mean) * (x[j] - mean);
double sd = Math.sqrt(squares / (n - 1));
for (double[] x : train) x[j] = (x[j] - mean) / sd;
for (double[] x : test) x[j] = (x[j] - mean) / sd;
}
// new double[3] starts at {0, 0, 0}: Java fills new arrays with zeros.
double[] weights = new double[3];
for (int step = 0; step < STEPS; step++) {
// Gradient of the average log-loss: the mean of (p - y) times
// each input, summed into its own array before any weight moves.
double[] gradient = new double[3];
for (int i = 0; i < n; i++) {
double error = sigmoid(dot(weights, train[i])) - y[i];
for (int j = 0; j < 3; j++) {
gradient[j] += error * train[i][j];
}
}
for (int j = 0; j < 3; j++) {
weights[j] -= LEARNING_RATE * gradient[j] / n;
}
}
double logLoss = 0;
for (int i = 0; i < n; i++) {
double p = sigmoid(dot(weights, train[i]));
logLoss -= y[i] * Math.log(p) + (1 - y[i]) * Math.log(1 - p);
}
logLoss /= n;
int correct = 0;
for (int i = 0; i < nTest; i++) {
double predicted = sigmoid(dot(weights, test[i])) >= 0.5 ? 1 : 0;
if (predicted == yTest[i]) correct++;
}
// 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, "train_count %d%n", n);
System.out.printf(root, "test_count %d%n", nTest);
System.out.printf(root, "intercept %.6f%n", weights[0]);
System.out.printf(root, "weight_hours %.6f%n", weights[1]);
System.out.printf(root, "weight_sleep %.6f%n", weights[2]);
System.out.printf(root, "log_loss %.6f%n", logLoss);
System.out.printf(root, "correct %d%n", correct);
// (double) converts before dividing; 37 / 50 between ints would be 0.
System.out.printf(root, "accuracy %.6f%n", (double) correct / nTest);
}
}Each step makes two passes over the 150 training rows. The first adds every
row’s error times its inputs into a fresh double[] gradient; the second
moves the weights in place. That’s safe only because the whole gradient
exists first: update inside the row loop and the intercept comes out
0.758788 instead of 0.753594.
readColumn returns strings here, since split holds words, and each row
is checked with .equals("train"), because == would compare the String
objects. The JDK has no regression or optimizer, so there was no built-in
to skip; Math.exp and Math.log do the rest.
// Logistic regression by gradient descent, using only Node's built-in modules.
// Predicts whether a student passed from study hours and sleep. Each feature
// is rescaled with the training rows' mean and standard deviation, then the
// weights start at zero and take 1000 fixed steps of batch gradient descent
// on the average log-loss.
import { readFileSync } from "node:fs";
const FEATURES = ["hours", "sleep"];
const LABEL = "passed";
const LEARNING_RATE = 0.5;
const STEPS = 1000;
// 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 sigmoid = (z) => 1 / (1 + Math.exp(-z));
const dot = (w, x) => sum(w.map((wj, j) => wj * x[j]));
const rows = readRows(process.argv[2]).map((row) => ({
x: FEATURES.map((name) => Number(row[name])),
y: Number(row[LABEL]),
split: row.split,
}));
const train = rows.filter((row) => row.split === "train");
const test = rows.filter((row) => row.split === "test");
const n = train.length;
// Scale each feature to mean 0 and (sample) standard deviation 1, using the
// training rows only. The leading 1 multiplies the intercept.
const means = FEATURES.map((_, j) => sum(train.map((row) => row.x[j])) / n);
const sds = FEATURES.map((_, j) =>
Math.sqrt(sum(train.map((row) => (row.x[j] - means[j]) ** 2)) / (n - 1)),
);
const scale = (x) => [1, ...x.map((v, j) => (v - means[j]) / sds[j])];
const xs = train.map((row) => scale(row.x));
const ys = train.map((row) => row.y);
let weights = [0, 0, 0];
for (let step = 0; step < STEPS; step++) {
// Gradient of the average log-loss: the mean of (p - y) times each input.
// map builds a new array, so every gradient uses the old weights.
const errors = xs.map((x, i) => sigmoid(dot(weights, x)) - ys[i]);
const gradient = weights.map((_, j) => sum(errors.map((e, i) => e * xs[i][j])) / n);
weights = weights.map((w, j) => w - LEARNING_RATE * gradient[j]);
}
const logLoss =
-sum(
xs.map((x, i) => {
const p = sigmoid(dot(weights, x));
return ys[i] * Math.log(p) + (1 - ys[i]) * Math.log(1 - p);
}),
) / n;
const predict = (x) => (sigmoid(dot(weights, scale(x))) >= 0.5 ? 1 : 0);
const correct = test.filter((row) => predict(row.x) === row.y).length;
console.log(`train_count ${n}`);
console.log(`test_count ${test.length}`);
console.log(`intercept ${weights[0].toFixed(6)}`);
console.log(`weight_hours ${weights[1].toFixed(6)}`);
console.log(`weight_sleep ${weights[2].toFixed(6)}`);
console.log(`log_loss ${logLoss.toFixed(6)}`);
console.log(`correct ${correct}`);
console.log(`accuracy ${(correct / test.length).toFixed(6)}`);Each step replaces the weights with a new array, weights.map(...), so
every part of the gradient is computed from the old weights, as the method
requires. A leading 1 joins the scaled features with spread syntax:
[1, ...x.map(...)].
The trap is the label. Number(row[LABEL]) looks optional, because
arithmetic converts strings for you: sigmoid(...) - "1" is a number, and
with string labels the weights train to exactly the same values. But strict
equality doesn’t convert, and 1 === "1" is false, so every prediction
would miss and correct would be 0. The standard library has no logistic
regression, so there was no built-in to skip.
# Logistic regression by gradient descent, using only Base Julia.
#
# Predicts whether a student passed from study hours and sleep. Each feature is
# rescaled with the training rows' mean and standard deviation, then the weights
# start at zero and take STEPS fixed steps of batch gradient descent on the
# average log-loss.
using Printf
const LEARNING_RATE = 0.5
const STEPS = 1000
function read_students(path)
train, test = Vector{Float64}[], Vector{Float64}[]
for line in Iterators.drop(eachline(path), 1) # skip the header
hours, sleep, passed, split_name = split(line, ",")
row = parse.(Float64, [hours, sleep, passed])
push!(split_name == "train" ? train : test, row)
end
# Stack the rows into matrices: one row per student, columns hours, sleep, passed.
return reduce(vcat, permutedims.(train)), reduce(vcat, permutedims.(test))
end
sigmoid(z) = 1 / (1 + exp(-z))
function main(path)
train, test = read_students(path)
features, y = train[:, 1:2], train[:, 3]
n = length(y)
# Scale each feature to mean 0 and (sample) standard deviation 1, using
# the training rows only. The column of ones multiplies the intercept.
means = sum(features, dims = 1) / n
sds = sqrt.(sum((features .- means) .^ 2, dims = 1) / (n - 1))
scale(f) = hcat(ones(size(f, 1)), (f .- means) ./ sds)
X = scale(features)
weights = zeros(3)
for _ in 1:STEPS
# Gradient of the average log-loss: X' (p - y) / n. X' is the transpose.
errors = sigmoid.(X * weights) .- y
weights -= LEARNING_RATE * (X' * errors) / n
end
p = sigmoid.(X * weights)
log_loss = -sum(y .* log.(p) .+ (1 .- y) .* log.(1 .- p)) / n
predicted = sigmoid.(scale(test[:, 1:2]) * weights) .>= 0.5
correct = count(predicted .== (test[:, 3] .== 1))
println("train_count $n")
println("test_count $(size(test, 1))")
@printf("intercept %.6f\n", weights[1])
@printf("weight_hours %.6f\n", weights[2])
@printf("weight_sleep %.6f\n", weights[3])
@printf("log_loss %.6f\n", log_loss)
println("correct $correct")
@printf("accuracy %.6f\n", correct / size(test, 1))
end
main(ARGS[1])This version looks like the math. The training rows are a matrix X with
a column of ones for the intercept, so each step is two lines:
errors = sigmoid.(X * weights) .- y, then
weights -= LEARNING_RATE * (X' * errors) / n, where X' is the
transpose. The whole gradient is computed before the weights change.
The trap is the dot in sigmoid.(...). Without it, sigmoid gets a whole
vector and exp stops with a MethodError. Where Python’s math.exp
raises an error above about 709, Julia’s returns Inf, so sigmoid gives
0.0 for very negative scores. There’s no logistic regression in the standard
library to skip.
-- Logistic regression by gradient descent, using only Lua's standard
-- library.
--
-- Predicts whether a student passed from study hours and sleep. Lua has no
-- regression, classifier, or matrix type, so every step is written out. Each
-- feature is rescaled with the training rows' mean and standard deviation,
-- then the weights start at zero and take 1000 fixed steps of batch gradient
-- descent on the average log-loss. Tables serve as arrays, indexed from 1.
local FEATURES = { "hours", "sleep" }
local LEARNING_RATE = 0.5
local STEPS = 1000
-- Returns two arrays, train and test. Each example is a table with the
-- feature values in x and the 0/1 label in y.
local function read_examples(path)
local f = assert(io.open(path, "r"))
local header = {}
for field in f:read("l"):gmatch("[^,]+") do
header[#header + 1] = field
end
local train, test = {}, {}
for line in f:lines() do
local row, i = {}, 0
for field in line:gmatch("[^,]+") do
i = i + 1
row[header[i]] = field -- a table keyed by column name
end
local x = {}
for j, name in ipairs(FEATURES) do
x[j] = tonumber(row[name])
end
local example = { x = x, y = tonumber(row.passed) }
local split = row.split == "train" and train or test
split[#split + 1] = example
end
f:close()
return train, test
end
local function sigmoid(z)
return 1 / (1 + math.exp(-z))
end
local function dot(u, v)
local total = 0
for i = 1, #u do
total = total + u[i] * v[i]
end
return total
end
local function say(name, value)
print(string.format("%s %.6f", name, value))
end
local train, test = read_examples(arg[1])
local n = #train
-- Scale each feature to mean 0 and (sample) standard deviation 1, using the
-- training rows only.
local means, sds = {}, {}
for j = 1, #FEATURES do
local total = 0
for _, e in ipairs(train) do total = total + e.x[j] end
local mean = total / n
local sumsq = 0
for _, e in ipairs(train) do sumsq = sumsq + (e.x[j] - mean) ^ 2 end
means[j], sds[j] = mean, math.sqrt(sumsq / (n - 1))
end
-- The leading 1 multiplies the intercept.
local function scale(x)
local scaled = { 1.0 }
for j, v in ipairs(x) do
scaled[j + 1] = (v - means[j]) / sds[j]
end
return scaled
end
local xs = {}
for i, e in ipairs(train) do xs[i] = scale(e.x) end
local weights = { 0.0, 0.0, 0.0 }
for _ = 1, STEPS do
-- Gradient of the average log-loss: the mean of (p - y) times each input.
local gradient = { 0.0, 0.0, 0.0 }
for i, x in ipairs(xs) do
local err = sigmoid(dot(weights, x)) - train[i].y
for j = 1, #gradient do
gradient[j] = gradient[j] + err * x[j]
end
end
-- Only now, with the whole gradient summed, does any weight move.
for j = 1, #weights do
weights[j] = weights[j] - LEARNING_RATE * gradient[j] / n
end
end
local loss = 0
for i, x in ipairs(xs) do
local p, y = sigmoid(dot(weights, x)), train[i].y
loss = loss - (y * math.log(p) + (1 - y) * math.log(1 - p))
end
local correct = 0
for _, e in ipairs(test) do
local predicted = sigmoid(dot(weights, scale(e.x))) >= 0.5 and 1 or 0
if predicted == e.y then correct = correct + 1 end
end
print("train_count " .. n)
print("test_count " .. #test)
say("intercept", weights[1])
say("weight_hours", weights[2])
say("weight_sleep", weights[3])
say("log_loss", loss / n)
print("correct " .. correct)
say("accuracy", correct / #test) -- / always gives a floatEach example is a small table, { x = x, y = y }, read as e.x and
e.y. Inside each of the 1000 steps, a loop over the students adds every
error times input to a fresh gradient table, and only after that loop
ends does a second loop update weights in place, so every weight moves by
the same gradient.
Lua’s standard library has no regression, classifier, or matrix type, so
there was no built-in to skip; sigmoid, dot, and scale are short
functions. The prediction uses p >= 0.5 and 1 or 0, Lua’s stand-in for a
conditional expression. The accuracy needs no conversion, since / always
gives a float: 37 / 50 is 0.74.
(* Logistic regression by gradient descent, using only OCaml's standard
library.
Predicts whether a student passed from study hours and sleep. Each
feature is rescaled with the training rows' mean and standard deviation,
then the weights start at zero and take 1,000 steps of batch gradient
descent on the average log-loss. The weights are an array, and each step
builds a new array from the old one. *)
let learning_rate = 0.5
let steps = 1000
(* Each row becomes ([| hours; sleep |], passed, split). *)
let read_rows path =
In_channel.with_open_text path (fun ic ->
let lines = In_channel.input_lines ic in
let header = String.split_on_char ',' (List.hd lines) in
let col name =
match List.find_index (String.equal name) header with
| Some i -> i
| None -> failwith ("column not found: " ^ name)
in
let hours = col "hours" and sleep = col "sleep" in
let passed = col "passed" and split = col "split" in
List.map
(fun line ->
let f = Array.of_list (String.split_on_char ',' line) in
( [| float_of_string f.(hours); float_of_string f.(sleep) |],
float_of_string f.(passed),
f.(split) ))
(List.tl lines))
let sigmoid z = 1. /. (1. +. exp (-.z))
(* Array.fold_left and Array.map2 walk two arrays side by side. *)
let dot u v = Array.fold_left ( +. ) 0. (Array.map2 ( *. ) u v)
let mean xs = Array.fold_left ( +. ) 0. xs /. float_of_int (Array.length xs)
(* Sample standard deviation: divide by n - 1. *)
let sd xs =
let m = mean xs in
let sumsq = Array.fold_left (fun acc x -> acc +. ((x -. m) ** 2.)) 0. xs in
sqrt (sumsq /. float_of_int (Array.length xs - 1))
(* One step of gradient descent. The gradient for weight j is the mean of
(p - y) times input j, summed over every student before any weight
moves, so the old array is never changed. *)
let step xs ys weights =
let n = float_of_int (Array.length xs) in
let errors = Array.map2 (fun x y -> sigmoid (dot weights x) -. y) xs ys in
Array.mapi
(fun j w ->
let column = Array.map (fun x -> x.(j)) xs in
w -. (learning_rate *. dot errors column /. n))
weights
let () =
let rows = read_rows Sys.argv.(1) in
(* partition splits one list into the rows that pass a test and the rest. *)
let train, test = List.partition (fun (_, _, split) -> split = "train") rows in
let train = Array.of_list train in
let features = Array.map (fun (x, _, _) -> x) train in
let ys = Array.map (fun (_, y, _) -> y) train in
(* Scale with the training rows only. The leading 1 multiplies the intercept. *)
let columns = Array.init 2 (fun j -> Array.map (fun x -> x.(j)) features) in
let means = Array.map mean columns and sds = Array.map sd columns in
let scale x =
Array.append [| 1. |] (Array.mapi (fun j v -> (v -. means.(j)) /. sds.(j)) x)
in
let xs = Array.map scale features in
(* let rec lets descend call itself: it takes one step, then counts the
remaining steps down until none are left. *)
let rec descend weights k = if k = 0 then weights else descend (step xs ys weights) (k - 1) in
let weights = descend [| 0.; 0.; 0. |] steps in
let n = Array.length xs in
let log_loss =
-.Array.fold_left ( +. ) 0.
(Array.map2
(fun x y ->
let p = sigmoid (dot weights x) in
(y *. log p) +. ((1. -. y) *. log (1. -. p)))
xs ys)
/. float_of_int n
in
let predict x = if sigmoid (dot weights (scale x)) >= 0.5 then 1. else 0. in
let correct = List.length (List.filter (fun (x, y, _) -> predict x = y) test) in
let test_count = List.length test in
Printf.printf "train_count %d\n" n;
Printf.printf "test_count %d\n" test_count;
Printf.printf "intercept %.6f\n" weights.(0);
Printf.printf "weight_hours %.6f\n" weights.(1);
Printf.printf "weight_sleep %.6f\n" weights.(2);
Printf.printf "log_loss %.6f\n" log_loss;
Printf.printf "correct %d\n" correct;
Printf.printf "accuracy %.6f\n" (float_of_int correct /. float_of_int test_count)The weights are an array, and no step changes it. step computes every
error from the old weights, then Array.mapi builds the new array, so
every weight moves from the same gradient. A recursive descend repeats
step 1,000 times and hands back the last array.
OCaml’s standard library has no regression, classifier, or matrix type,
so sigmoid, dot, and sd are short functions. dot pairs two arrays
with Array.map2, which raises Invalid_argument if their lengths
differ. The label is read as a float, since 1 - y with a float y
won’t compile; it has to be 1. -. y. List.partition splits the rows
into training and test sets in one pass.
% Logistic regression by gradient descent, using only core GNU Octave.
%
% Predicts whether a student passed from study hours and sleep. Each feature
% is rescaled with the training rows' mean and standard deviation, then the
% weights start at zero and take STEPS fixed steps of batch gradient descent
% on the average log-loss. The training rows form a matrix X, so each step is
% two matrix products.
LEARNING_RATE = 0.5;
STEPS = 1000;
% The named column of a CSV file, as a cell array of text fields.
function fields = 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
fields = cell(numel(lines) - 1, 1);
for row = 2:numel(lines)
parts = strsplit(strtrim(lines{row}), ",");
fields{row - 1} = parts{index};
end
end
% The ./ divides element by element, so this works on a whole vector of z.
function p = sigmoid(z)
p = 1 ./ (1 + exp(-z));
end
args = argv();
path = args{1};
% str2double turns a whole cell array of text into a column of numbers.
features = [str2double(read_column(path, "hours")), ...
str2double(read_column(path, "sleep"))];
passed = str2double(read_column(path, "passed"));
is_train = strcmp(read_column(path, "split"), "train");
y = passed(is_train);
n = numel(y);
% Scale each feature to mean 0 and (sample) standard deviation 1, using the
% training rows only. sum(A, 1) adds down each column; subtracting the 1-by-2
% row of means from an n-by-2 matrix subtracts it from every row.
train_features = features(is_train, :);
means = sum(train_features, 1) / n;
sds = sqrt(sum((train_features - means) .^ 2, 1) / (n - 1));
% The column of ones multiplies the intercept.
scale = @(f) [ones(rows(f), 1), (f - means) ./ sds];
X = scale(train_features);
weights = zeros(3, 1);
for step = 1:STEPS
% Gradient of the average log-loss: X' * (p - y) / n, where X' is the
% transpose. The whole gradient is computed before any weight changes.
errors = sigmoid(X * weights) - y;
weights = weights - LEARNING_RATE * (X' * errors) / n;
end
p = sigmoid(X * weights);
log_loss = -sum(y .* log(p) + (1 - y) .* log(1 - p)) / n;
predicted = sigmoid(scale(features(~is_train, :)) * weights) >= 0.5;
actual = passed(~is_train) == 1;
test_count = numel(actual);
correct = sum(predicted == actual);
printf("train_count %d\n", n);
printf("test_count %d\n", test_count);
printf("intercept %.6f\n", weights(1));
printf("weight_hours %.6f\n", weights(2));
printf("weight_sleep %.6f\n", weights(3));
printf("log_loss %.6f\n", log_loss);
printf("correct %d\n", correct);
printf("accuracy %.6f\n", correct / test_count);The training rows are a matrix X with a column of ones for the
intercept, so each step is two lines:
errors = sigmoid(X * weights) - y, then
weights = weights - LEARNING_RATE * (X' * errors) / n, where X' is the
transpose. The whole gradient is computed before any weight changes.
The trap is the dot in 1 ./ (1 + exp(-z)). With a plain /, Octave
treats the vector as a system to solve: the first step quietly builds a 150
by 150 matrix, and the second stops with a “nonconformant arguments” error.
We skip optimizers such as fminunc. Core Octave has no glmfit, and the
statistics package isn’t installed.
{ Logistic regression by gradient descent, using only Free Pascal's
run-time library.
Predicts whether a student passed from study hours and sleep. Each
feature is rescaled with the training rows' mean and standard deviation,
then the weights start at zero and take Steps fixed steps of batch
gradient descent on the average log-loss. The Math unit's Mean and
StdDev are left out, so the scaling is written out too. }
program LogisticRegression;
{$mode objfpc}{$H+} // 32-bit Integer, and strings longer than 255 characters
uses
SysUtils;
const
K = 3; // inputs per row: 1 (for the intercept), hours, sleep
LearningRate = 0.5;
Steps = 1000;
type
TVector = array[0..K - 1] of Double; // a fixed-size array type
TExample = record
X: TVector; // X[0] is always 1, so weight 0 acts as the intercept
Y: Integer;
end;
TExamples = array of TExample;
procedure Append(var Rows: TExamples; const Row: TExample);
begin
SetLength(Rows, Length(Rows) + 1); // grow by one
Rows[High(Rows)] := Row;
end;
function Column(const Header: TStringArray; const Name: String): Integer;
var
I: Integer;
begin
for I := 0 to High(Header) do
if Header[I] = Name then
Exit(I); // Exit(value) sets the result and returns
raise Exception.Create('column not found: ' + Name);
end;
{ Read the CSV into training and test rows, sorted by the "split" column. }
procedure ReadRows(const Path: String; out Train, Test: TExamples);
var
F: TextFile;
Line: String;
Fields: TStringArray;
CHours, CSleep, CPassed, CSplit: Integer;
Row: TExample;
begin
Train := nil;
Test := nil;
AssignFile(F, Path);
Reset(F);
ReadLn(F, Line);
Fields := Line.Split(',');
CHours := Column(Fields, 'hours');
CSleep := Column(Fields, 'sleep');
CPassed := Column(Fields, 'passed');
CSplit := Column(Fields, 'split');
while not Eof(F) do
begin
ReadLn(F, Line);
Fields := Line.Split(',');
Row.X[0] := 1;
Row.X[1] := StrToFloat(Fields[CHours]);
Row.X[2] := StrToFloat(Fields[CSleep]);
Row.Y := StrToInt(Fields[CPassed]);
if Fields[CSplit] = 'train' then
Append(Train, Row)
else
Append(Test, Row);
end;
CloseFile(F);
end;
{ The probability of passing: the sigmoid of the weighted sum of inputs. }
function Predict(const W, X: TVector): Double;
var
J: Integer;
Z: Double;
begin
Z := 0;
for J := 0 to K - 1 do
Z := Z + W[J] * X[J];
Result := 1 / (1 + Exp(-Z));
end;
procedure Say(const Name: String; Value: Double);
begin
WriteLn(Name, ' ', Value:0:6); // width 0, six decimals
end;
var
Train, Test: TExamples;
W, Gradient: TVector;
N, I, J, Step, Correct: Integer;
Total, Squares, Mean, Sd, Error, P, Loss: Double;
begin
ReadRows(ParamStr(1), Train, Test);
N := Length(Train);
// Scale each feature to mean 0 and (sample) standard deviation 1, using
// the training rows only, then apply the same shift to the test rows.
// Column 0 holds the ones, so the features are columns 1 and 2.
for J := 1 to K - 1 do
begin
Total := 0;
for I := 0 to N - 1 do
Total := Total + Train[I].X[J];
Mean := Total / N;
Squares := 0;
for I := 0 to N - 1 do
Squares := Squares + Sqr(Train[I].X[J] - Mean);
Sd := Sqrt(Squares / (N - 1));
for I := 0 to N - 1 do
Train[I].X[J] := (Train[I].X[J] - Mean) / Sd;
for I := 0 to High(Test) do
Test[I].X[J] := (Test[I].X[J] - Mean) / Sd;
end;
for J := 0 to K - 1 do
W[J] := 0;
for Step := 1 to Steps do
begin
// Gradient of the average log-loss: the mean of (p - y) times each
// input. It is summed over every student before any weight moves.
for J := 0 to K - 1 do
Gradient[J] := 0;
for I := 0 to N - 1 do
begin
Error := Predict(W, Train[I].X) - Train[I].Y;
for J := 0 to K - 1 do
Gradient[J] := Gradient[J] + Error * Train[I].X[J];
end;
for J := 0 to K - 1 do
W[J] := W[J] - LearningRate * Gradient[J] / N;
end;
Loss := 0;
for I := 0 to N - 1 do
begin
P := Predict(W, Train[I].X);
if Train[I].Y = 1 then
Loss := Loss - Ln(P)
else
Loss := Loss - Ln(1 - P);
end;
Loss := Loss / N;
Correct := 0;
for I := 0 to High(Test) do
// Ord turns a Boolean into 0 or 1, so it compares with the label.
if Ord(Predict(W, Test[I].X) >= 0.5) = Test[I].Y then
Inc(Correct);
WriteLn('train_count ', N);
WriteLn('test_count ', Length(Test));
Say('intercept', W[0]);
Say('weight_hours', W[1]);
Say('weight_sleep', W[2]);
Say('log_loss', Loss);
WriteLn('correct ', Correct);
Say('accuracy', Correct / Length(Test)); // / gives a real, even for integers
end.Each student is a record holding a fixed-size TVector of three inputs,
the first always 1 so that weight 0 acts as the intercept, plus the label.
ReadRows finds the columns by name and appends each row to a training or
test dynamic array, grown with SetLength. We leave out the Math unit’s
Mean and StdDev, so the scaling is written as loops.
Every step sums the whole gradient in a nested loop over students and
weights, then moves all three weights at once. Ord turns the Boolean
p >= 0.5 into 0 or 1 to compare with the label, and the accuracy needs
no cast, since / always gives a real number.
# Logistic regression by gradient descent, using only Perl and its core
# modules.
#
# Predicts whether a student passed from study hours and sleep. Each feature
# is rescaled with the training rows' mean and standard deviation, then the
# weights start at zero and take STEPS fixed steps of batch gradient descent
# on the average log-loss. List::Util, a core module, supplies sum.
use strict;
use warnings;
use List::Util qw(sum);
my @FEATURES = qw(hours sleep);
my $LEARNING_RATE = 0.5;
my $STEPS = 1000;
# Returns two array refs, train and test. Each example is a hash ref,
# { x => [hours, sleep], y => passed }; [...] and {...} build references.
sub read_examples {
my ($path) = @_;
open my $fh, '<', $path or die "can't open $path: $!";
chomp(my $header = <$fh>);
my @names = split /,/, $header;
my (@train, @test);
while (my $line = <$fh>) {
chomp $line;
my %row;
@row{@names} = split /,/, $line; # a hash slice: names to fields
my $example = { x => [ @row{@FEATURES} ], y => $row{passed} };
push @{ $row{split} eq 'train' ? \@train : \@test }, $example;
}
close $fh;
return (\@train, \@test);
}
sub sigmoid { return 1 / (1 + exp(-$_[0])) }
# The dot product of two array refs; $#$u is the last index of @$u.
sub dot {
my ($u, $v) = @_;
return sum(map { $u->[$_] * $v->[$_] } 0 .. $#$u);
}
my ($train, $test) = read_examples($ARGV[0]);
my $n = @$train;
# Scale each feature to mean 0 and (sample) standard deviation 1, using the
# training rows only.
my (@means, @sds);
for my $j (0 .. $#FEATURES) {
my @column = map { $_->{x}[$j] } @$train;
my $mean = sum(@column) / $n;
push @means, $mean;
push @sds, sqrt(sum(map { ($_ - $mean) ** 2 } @column) / ($n - 1));
}
# The leading 1 multiplies the intercept.
sub scale {
my ($x) = @_;
return [ 1.0, map { ($x->[$_] - $means[$_]) / $sds[$_] } 0 .. $#$x ];
}
my @xs = map { scale($_->{x}) } @$train;
my @ys = map { $_->{y} } @$train;
my @weights = (0.0) x @{ $xs[0] }; # x repeats a list: three zeros
for (1 .. $STEPS) {
# Gradient of the average log-loss: the mean of (p - y) times each input.
my @errors = map { sigmoid(dot(\@weights, $xs[$_])) - $ys[$_] } 0 .. $n - 1;
my @gradient = map {
my $j = $_;
sum(map { $errors[$_] * $xs[$_][$j] } 0 .. $n - 1) / $n;
} 0 .. $#weights;
@weights = map { $weights[$_] - $LEARNING_RATE * $gradient[$_] } 0 .. $#weights;
}
my @p = map { sigmoid(dot(\@weights, $_)) } @xs;
my $log_loss = -sum(map {
$ys[$_] * log($p[$_]) + (1 - $ys[$_]) * log(1 - $p[$_])
} 0 .. $n - 1) / $n;
# grep keeps the test examples predicted right; assigned to a single value,
# it gives their count.
my $correct = grep {
(sigmoid(dot(\@weights, scale($_->{x}))) >= 0.5 ? 1 : 0) == $_->{y}
} @$test;
printf "train_count %d\n", $n;
printf "test_count %d\n", scalar @$test;
printf "intercept %.6f\n", $weights[0];
printf "weight_hours %.6f\n", $weights[1];
printf "weight_sleep %.6f\n", $weights[2];
printf "log_loss %.6f\n", $log_loss;
printf "correct %d\n", $correct;
printf "accuracy %.6f\n", $correct / @$test;Each example is a hash ref, { x => [...], y => 0 }, read as
$_->{x} and $_->{y}. The gradient step has no explicit loop over
students: map builds every error, a nested map passed to sum turns
them into one gradient entry per weight, and a third map replaces
@weights all at once, so every weight moves by the same gradient.
Perl has no regression, classifier, or matrix type, so sigmoid, dot,
and scale are short subs, with exp, log, and sqrt as the only math
built-ins. Counting correct predictions uses grep, which gives the
number of matches when assigned to a single value.
<?php
// Logistic regression by gradient descent, using only PHP and the extensions
// built into it.
//
// Predicts whether a student passed from study hours and sleep. Each feature
// is rescaled with the training rows' mean and standard deviation, then the
// weights start at zero and take STEPS fixed steps of batch gradient descent
// on the average log-loss.
declare(strict_types=1);
const FEATURES = ['hours', 'sleep'];
const LABEL = 'passed';
const LEARNING_RATE = 0.5;
const STEPS = 1000;
// No field is quoted, so explode splits each line on commas. array_combine
// pairs the header names with the fields, giving a row keyed by name.
function read_rows(string $path): array
{
$lines = file($path, FILE_IGNORE_NEW_LINES | FILE_SKIP_EMPTY_LINES);
$header = explode(',', array_shift($lines));
return array_map(fn($line) => array_combine($header, explode(',', $line)), $lines);
}
function sigmoid(float $z): float
{
return 1 / (1 + exp(-$z));
}
function dot(array $a, array $b): float
{
$total = 0.0;
foreach ($a as $i => $value) {
$total += $value * $b[$i];
}
return $total;
}
$train = [];
$test = [];
foreach (read_rows($argv[1]) as $row) {
$x = array_map(fn($f) => (float) $row[$f], FEATURES);
$example = [$x, (int) $row[LABEL]];
if ($row['split'] === 'train') {
$train[] = $example;
} else {
$test[] = $example;
}
}
// Scale each feature to mean 0 and (sample) standard deviation 1, using the
// training rows only.
$means = [];
$sds = [];
foreach (array_keys(FEATURES) as $j) {
$column = array_map(fn($ex) => $ex[0][$j], $train);
$m = array_sum($column) / count($column);
$ss = array_sum(array_map(fn($v) => ($v - $m) ** 2, $column));
$means[] = $m;
$sds[] = sqrt($ss / (count($column) - 1));
}
// The leading 1.0 multiplies the intercept. Called with several arrays,
// array_map walks them side by side.
$scale = fn(array $x): array =>
[1.0, ...array_map(fn($v, $m, $s) => ($v - $m) / $s, $x, $means, $sds)];
$xs = array_map(fn($ex) => $scale($ex[0]), $train);
$ys = array_column($train, 1);
$n = count($xs);
$weights = array_fill(0, count($xs[0]), 0.0);
for ($step = 0; $step < STEPS; $step++) {
// Gradient of the average log-loss: the mean of (p - y) times each input.
$gradient = array_fill(0, count($weights), 0.0);
foreach ($xs as $i => $x) {
$error = sigmoid(dot($weights, $x)) - $ys[$i];
foreach ($x as $j => $v) {
$gradient[$j] += $error * $v / $n;
}
}
// Update every weight only after the whole gradient is known.
foreach ($gradient as $j => $g) {
$weights[$j] -= LEARNING_RATE * $g;
}
}
$loss = 0.0;
foreach ($xs as $i => $x) {
$p = sigmoid(dot($weights, $x));
$loss -= $ys[$i] * log($p) + (1 - $ys[$i]) * log(1 - $p);
}
$loss /= $n;
$correct = 0;
foreach ($test as [$x, $y]) {
$predicted = sigmoid(dot($weights, $scale($x))) >= 0.5 ? 1 : 0;
if ($predicted === $y) {
$correct++;
}
}
echo "train_count $n\n";
echo 'test_count ' . count($test) . "\n";
printf("intercept %.6f\n", $weights[0]);
printf("weight_hours %.6f\n", $weights[1]);
printf("weight_sleep %.6f\n", $weights[2]);
printf("log_loss %.6f\n", $loss);
echo "correct $correct\n";
printf("accuracy %.6f\n", $correct / count($test));Each example is a two-item array, the features and the label, and
foreach ($test as [$x, $y]) unpacks it as the loop runs. The gradient
step uses plain loops: one pass over the students adds each error times
each input to a running gradient, and a second loop moves every weight
only once the whole gradient is known.
PHP has no regression, classifier, or matrix type, so sigmoid, dot,
and the $scale closure are written out, with exp, log, and sqrt
as the only math built-ins. Called with several arrays, array_map
walks them side by side, which is how $scale pairs each value with
its mean and standard deviation.
-- Logistic regression by gradient descent, in PostgreSQL.
--
-- Predicts whether a student passed from study hours and sleep. Each feature
-- is rescaled with the training rows' mean and sample standard deviation, then
-- the weights start at zero and take 1000 steps of batch gradient descent on
-- the average log-loss. The loop is a WITH RECURSIVE query whose state is one
-- row: the step number and the three weights.
WITH RECURSIVE train AS (
SELECT hours, sleep, passed FROM data WHERE split = 'train'
),
means AS (
-- avg() would do this; sum and count keep the arithmetic visible.
SELECT count(*) AS n,
sum(hours) / count(*) AS mean_hours,
sum(sleep) / count(*) AS mean_sleep
FROM train
),
spread AS (
-- Sample standard deviation (divide by n - 1), instead of stddev_samp().
SELECT sqrt(sum((hours - mean_hours) ^ 2) / (n - 1)) AS sd_hours,
sqrt(sum((sleep - mean_sleep) ^ 2) / (n - 1)) AS sd_sleep
FROM train, means
GROUP BY n
),
scaled AS (
-- Every row, scaled with the training statistics. The intercept's
-- input is a constant 1, so it needs no column.
SELECT split,
(hours - mean_hours) / sd_hours AS x1,
(sleep - mean_sleep) / sd_sleep AS x2,
passed AS y
FROM data, means, spread
),
steps (step, w0, w1, w2) AS (
SELECT 0, 0::double precision, 0::double precision, 0::double precision
UNION ALL
-- One step: the gradient is the mean of (p - y) times each input, with p
-- from the old weights; then all three weights move together. Aggregates
-- aren't allowed directly in the recursive step, so they sit in LATERAL.
SELECT s.step + 1,
s.w0 - 0.5 * g.g0,
s.w1 - 0.5 * g.g1,
s.w2 - 0.5 * g.g2
FROM steps s,
LATERAL (
SELECT sum(e) / count(*) AS g0,
sum(e * x1) / count(*) AS g1,
sum(e * x2) / count(*) AS g2
FROM (
SELECT x1, x2,
1 / (1 + exp(-(s.w0 + s.w1 * x1 + s.w2 * x2))) - y AS e
FROM scaled
WHERE split = 'train'
) AS errors
) AS g
WHERE s.step < 1000
),
final AS (
SELECT * FROM steps WHERE step = 1000
),
scored AS (
-- Each row's predicted probability under the final weights.
SELECT c.split, c.y,
1 / (1 + exp(-(f.w0 + f.w1 * c.x1 + f.w2 * c.x2))) AS p
FROM scaled c, final f
),
loss AS (
SELECT -sum(y * ln(p) + (1 - y) * ln(1 - p)) / count(*) AS log_loss
FROM scored WHERE split = 'train'
),
test AS (
-- Predict 1 when p >= 0.5, and count the matches.
SELECT count(*) AS n,
sum(CASE WHEN (p >= 0.5) = (y = 1) THEN 1 ELSE 0 END) AS correct
FROM scored WHERE split = 'test'
)
SELECT name, value
FROM means, final, loss, test,
LATERAL (VALUES
(1, 'train_count', means.n::text),
(2, 'test_count', test.n::text),
(3, 'intercept', round(w0::numeric, 6)::text),
(4, 'weight_hours', round(w1::numeric, 6)::text),
(5, 'weight_sleep', round(w2::numeric, 6)::text),
(6, 'log_loss', round(log_loss::numeric, 6)::text),
(7, 'correct', correct::text),
(8, 'accuracy', round((correct::numeric / test.n), 6)::text)
) AS results (position, name, value)
ORDER BY position;Gradient descent is a loop, and SQL has none, so we use WITH RECURSIVE
again. Each row of the recursion is one step: the step number and the
three weights, 1001 rows from step 0 to step 1000. Inside a step, the
gradient is three sum aggregates over the scaled training rows, divided
by count(*). PostgreSQL rejects aggregates in the recursive term itself
(“aggregate functions are not allowed in a recursive query’s recursive
term”), so they sit in a LATERAL subquery that sees the old weights.
The scaling uses sum and count in place of avg and stddev_samp.
Core PostgreSQL has no logistic regression, sigmoid, or logit function to
skip; exp and ln do the rest.
% Logistic regression by gradient descent, using SWI-Prolog and the
% libraries that come with it.
%
% Predicts whether a student passed from study hours and sleep. Each
% feature is rescaled with the training rows' mean and standard deviation,
% then the weights start at zero and take 1,000 steps of batch gradient
% descent on the average log-loss. A Prolog variable takes one value and
% keeps it, so a step can't change the weights: step/4 relates one list of
% weights to the next.
:- initialization(main, main).
% transpose/2, which turns rows into columns, comes from library(clpfd).
:- use_module(library(clpfd), [transpose/2]).
learning_rate(0.5).
steps(1000).
% read_rows(+Path, -Train, -Test): each list holds one Features-Label pair
% per row, such as [1.0, 7.8]-0, split by the "split" column.
read_rows(Path, Train, Test) :-
read_file_to_string(Path, Text, []),
split_string(Text, "\n", "", [HeaderLine|Lines]),
split_string(HeaderLine, ",", "", Header),
nth0(H, Header, "hours"),
nth0(S, Header, "sleep"),
nth0(P, Header, "passed"),
nth0(G, Header, "split"),
findall(Split-([Hours, Sleep]-Passed),
( member(Line, Lines),
Line \== "",
split_string(Line, ",", "", Fields),
nth0(H, Fields, HoursText), number_string(Hours, HoursText),
nth0(S, Fields, SleepText), number_string(Sleep, SleepText),
nth0(P, Fields, PassedText), number_string(Passed, PassedText),
nth0(G, Fields, Split) ),
Rows),
% "train" is a string; findall keeps only the rows that match it.
findall(Row, member("train"-Row, Rows), Train),
findall(Row, member("test"-Row, Rows), Test).
dot(U, V, Dot) :-
foldl([A, B, Acc0, Acc]>>(Acc is Acc0 + A * B), U, V, 0, Dot).
sigmoid(Z, P) :-
P is 1 / (1 + exp(-Z)).
% mean_sd(+Column, -Mean-SD): the mean and sample standard deviation.
mean_sd(Column, Mean-SD) :-
length(Column, N),
sum_list(Column, Sum),
Mean is Sum / N,
foldl([X, Acc0, Acc]>>(Acc is Acc0 + (X - Mean) ** 2), Column, 0, Squares),
SD is sqrt(Squares / (N - 1)).
% scale(+Stats, +Features, -Inputs): rescale each feature, and put a 1 in
% front to multiply the intercept.
scale(Stats, Features, [1.0|Scaled]) :-
maplist([X, Mean-SD, Z]>>(Z is (X - Mean) / SD), Features, Stats, Scaled).
probability(Weights, Inputs, P) :-
dot(Weights, Inputs, Z),
sigmoid(Z, P).
% step(+Xs, +Ys, +Weights, -Next): one step of gradient descent. Every
% error comes from the old weights before any weight moves. Each weight's
% gradient is the mean of (p - y) times its input, a dot product of the
% errors with one column of inputs.
step(Xs, Ys, Weights, Next) :-
maplist([X, Y, E]>>(probability(Weights, X, P), E is P - Y), Xs, Ys, Errors),
length(Xs, N),
transpose(Xs, Columns),
learning_rate(Rate),
maplist([W, Column, W1]>>(dot(Errors, Column, G), W1 is W - Rate * G / N),
Weights, Columns, Next).
% descend(+K, +Xs, +Ys, +Weights, -Final): Final is Weights after K steps.
% The recursive call comes last, so Prolog reuses the same stack frame.
descend(0, _, _, Weights, Weights) :- !.
descend(K, Xs, Ys, Weights, Final) :-
step(Xs, Ys, Weights, Next),
K1 is K - 1,
descend(K1, Xs, Ys, Next, Final).
log_loss_term(P, Y, Loss) :-
Loss is -(Y * log(P) + (1 - Y) * log(1 - P)).
main :-
current_prolog_flag(argv, [Path|_]),
read_rows(Path, Train, Test),
pairs_keys_values(Train, Features, Ys),
transpose(Features, FeatureColumns),
maplist(mean_sd, FeatureColumns, Stats),
maplist(scale(Stats), Features, Xs),
length(Xs, N),
steps(Steps),
descend(Steps, Xs, Ys, [0.0, 0.0, 0.0], Weights),
maplist(probability(Weights), Xs, Ps),
maplist(log_loss_term, Ps, Ys, Losses),
sum_list(Losses, TotalLoss),
LogLoss is TotalLoss / N,
% Count the test rows whose prediction (pass when p >= 0.5) is right.
aggregate_all(count,
( member(F-Y, Test),
scale(Stats, F, X),
probability(Weights, X, P),
( P >= 0.5 -> Y =:= 1 ; Y =:= 0 ) ),
Correct),
length(Test, TestCount),
Accuracy is Correct / TestCount,
Weights = [Intercept, Hours, Sleep],
format("train_count ~d~n", [N]),
format("test_count ~d~n", [TestCount]),
format("intercept ~6f~n", [Intercept]),
format("weight_hours ~6f~n", [Hours]),
format("weight_sleep ~6f~n", [Sleep]),
format("log_loss ~6f~n", [LogLoss]),
format("correct ~d~n", [Correct]),
format("accuracy ~6f~n", [Accuracy]).A Prolog variable gets one value and keeps it, so nothing here updates
the weights. step(Xs, Ys, Weights, Next) relates one list of weights to
the next, and descend counts down from 1,000, passing Next along.
Because the recursive call comes last, Prolog reuses its stack frame: we
ran 20,000 steps inside a 2 MB stack. Every error comes from the old
weights before Next is built, so all three weights move from the same
gradient.
There’s no logistic regression or matrix library in SWI-Prolog to skip.
The [X, Y]>>Goal lambdas come from library(yall), and transpose,
which turns rows into columns, from library(clpfd). Both ship with
SWI-Prolog.
# Logistic regression by gradient descent, using only base R.
#
# Predicts whether a student passed from study hours and sleep.
# glm(passed ~ hours + sleep, family = binomial) would fit it in one call (by
# a different method, iteratively reweighted least squares); here the weights
# start at zero and take a fixed number of batch gradient-descent steps.
learning_rate <- 0.5
steps <- 1000
args <- commandArgs(trailingOnly = TRUE)
data <- read.csv(args[1])
train <- data[data$split == "train", ]
test <- data[data$split == "test", ]
features <- c("hours", "sleep")
# Training means and sample standard deviations (divide by n - 1), written out
# instead of mean() and sd(). sapply runs the function once per column.
means <- sapply(train[features], function(v) sum(v) / length(v))
sds <- sapply(features, function(f) {
v <- train[[f]]
sqrt(sum((v - means[f])^2) / (length(v) - 1))
})
# A matrix with a column of 1s for the intercept, then the scaled features.
# sweep() subtracts each column's mean, then divides by its sd. unname()
# drops the column names, which would otherwise ride along into the weights.
design <- function(rows) {
x <- as.matrix(rows[features])
unname(cbind(1, sweep(sweep(x, 2, means), 2, sds, "/")))
}
sigmoid <- function(z) 1 / (1 + exp(-z))
x <- design(train)
y <- train$passed
n <- nrow(x)
w <- rep(0, ncol(x))
for (step in seq_len(steps)) {
# x %*% w is every row's w . x at once; t(x) %*% errors sums error times
# input for each weight, so dividing by n gives the gradient.
errors <- sigmoid(drop(x %*% w)) - y
gradient <- drop(t(x) %*% errors) / n
w <- w - learning_rate * gradient
}
p <- sigmoid(drop(x %*% w))
log_loss <- -sum(y * log(p) + (1 - y) * log(1 - p)) / n
predicted <- ifelse(sigmoid(drop(design(test) %*% w)) >= 0.5, 1, 0)
correct <- sum(predicted == test$passed)
results <- c(intercept = w[1], weight_hours = w[2], weight_sleep = w[3], log_loss = log_loss)
cat(sprintf("train_count %d\n", n))
cat(sprintf("test_count %d\n", nrow(test)))
cat(sprintf("%s %.6f\n", names(results), results), sep = "")
cat(sprintf("correct %d\n", correct))
cat(sprintf("accuracy %.6f\n", correct / nrow(test)))R can fit this in one call: glm(passed ~ hours + sleep, family = binomial).
We avoided it, along with mean() and sd(). On this data, glm reaches
the same three weights to six decimals in 5 iterations of iteratively
reweighted least squares, where we take 1000 plain gradient steps.
Matrix arithmetic does the looping for us. With a column of 1s for the
intercept, x %*% w gives every row’s score at once, and t(x) %*% errors
sums error times input for each weight. One catch: the matrix keeps its
column names, and %*% passes them on to the weights, so design() calls
unname() to keep the printed labels clean.
# Logistic regression with glm(), R's generalized linear model function.
#
# The built-in version of ../r/main.R. family = binomial makes glm() fit a
# logistic regression. It finds the weights by iteratively reweighted least
# squares, not gradient descent, but reaches the same answer.
args <- commandArgs(trailingOnly = TRUE)
data <- read.csv(args[1])
train <- data[data$split == "train", ]
test <- data[data$split == "test", ]
features <- c("hours", "sleep")
# scale() standardizes each column with its mean and sample standard
# deviation, and remembers both, so the test rows can use the training ones.
scaled_train <- scale(train[features])
center <- attr(scaled_train, "scaled:center")
spread <- attr(scaled_train, "scaled:scale")
scaled_test <- scale(test[features], center = center, scale = spread)
model <- glm(train$passed ~ scaled_train, family = binomial)
w <- unname(coef(model)) # intercept, then hours and sleep
p <- fitted(model) # each training student's probability of passing
y <- train$passed
log_loss <- -mean(y * log(p) + (1 - y) * log(1 - p))
# plogis() is the logistic (sigmoid) function.
predicted <- ifelse(plogis(drop(cbind(1, scaled_test) %*% w)) >= 0.5, 1, 0)
correct <- sum(predicted == test$passed)
results <- c(intercept = w[1], weight_hours = w[2], weight_sleep = w[3], log_loss = log_loss)
cat(sprintf("train_count %d\n", nrow(train)))
cat(sprintf("test_count %d\n", nrow(test)))
cat(sprintf("%s %.6f\n", names(results), results), sep = "")
cat(sprintf("correct %d\n", correct))
cat(sprintf("accuracy %.6f\n", correct / nrow(test)))glm() fits generalized linear models, and family = binomial makes it a
logistic regression. It doesn’t take 1,000 small gradient steps. It uses
iteratively reweighted least squares, which solves a weighted regression
again and again, and here it settles in five rounds on the same weights,
to six decimals, as the from-scratch version.
scale() standardizes the training columns and remembers each mean and
standard deviation, so the test rows are scaled with the training ones.
fitted() gives each training student’s probability for the log-loss, and
plogis(), R’s logistic function, turns the test rows’ weighted sums into
probabilities.
# Logistic regression by gradient descent, using only the standard library.
# Predicts whether a student passed from study hours and sleep. Each feature is
# rescaled with the training rows' mean and standard deviation, then the weights
# start at zero and take STEPS fixed steps of batch gradient descent on the
# average log-loss.
FEATURES = %w[hours sleep].freeze
LABEL = "passed"
LEARNING_RATE = 0.5
STEPS = 1000
# The csv library is a bundled gem since Ruby 3.4, so split lines by hand
# (no quoted fields here) and turn each row into a hash keyed by header name.
def read_rows(path)
header, *lines = File.readlines(path, chomp: true)
names = header.split(",")
lines.map { |line| names.zip(line.split(",")).to_h }
end
def sigmoid(z)
1 / (1 + Math.exp(-z))
end
def dot(a, b)
a.zip(b).sum { |u, v| u * v }
end
rows = read_rows(ARGV[0])
examples = rows.map { |row| [FEATURES.map { |f| Float(row[f]) }, Integer(row[LABEL]), row["split"]] }
train, test = examples.partition { |_, _, split| split == "train" }
# Scale each feature to mean 0 and (sample) standard deviation 1, using the
# training rows only. transpose turns the rows into one array per feature.
columns = train.map(&:first).transpose
means = columns.map { |c| c.sum / c.length }
sds = columns.zip(means).map { |c, m| Math.sqrt(c.sum { |v| (v - m)**2 } / (c.length - 1)) }
# The leading 1.0 multiplies the intercept.
scale = ->(x) { [1.0] + x.zip(means, sds).map { |v, m, s| (v - m) / s } }
xs = train.map { |x, _, _| scale.(x) }
ys = train.map { |_, y, _| y }
n = xs.length
weights = Array.new(xs.first.length, 0.0)
STEPS.times do
# Gradient of the average log-loss: the mean of (p - y) times each input.
errors = xs.zip(ys).map { |x, y| sigmoid(dot(weights, x)) - y }
gradient = weights.each_index.map { |j| errors.zip(xs).sum { |e, x| e * x[j] } / n }
weights = weights.zip(gradient).map { |w, g| w - LEARNING_RATE * g }
end
log_loss = -xs.zip(ys).sum do |x, y|
p = sigmoid(dot(weights, x))
y * Math.log(p) + (1 - y) * Math.log(1 - p)
end / n
correct = test.count { |x, y, _| (sigmoid(dot(weights, scale.(x))) >= 0.5 ? 1 : 0) == y }
puts "train_count #{n}"
puts "test_count #{test.length}"
{ intercept: weights[0], weight_hours: weights[1], weight_sleep: weights[2],
log_loss: log_loss }.each do |name, value|
puts format("%s %.6f", name, value)
end
puts "correct #{correct}"
puts format("accuracy %.6f", correct.to_f / test.length)Ruby keeps each step of gradient descent to one line. transpose turns the
training rows into one array per feature, so the means and standard
deviations are a map over columns. Inside the loop, errors.zip(xs).sum
with a block is the sum of error times input for one weight, and building
the new weights array from the old one means every weight updates from the
same gradient.
The scaling is a lambda, scale, called as scale.(x), so the test rows
reuse the training means and standard deviations. Ruby’s standard library
has no regression or classifier, so there was no built-in to skip.
//! Logistic regression by gradient descent, using only the standard library.
//!
//! Predicts whether a student passed from study hours and sleep. Each feature
//! is rescaled with the training rows' mean and standard deviation, then the
//! weights start at zero and take 1000 fixed steps of batch gradient descent
//! on the average log-loss.
use std::env;
use std::fs;
const FEATURES: [&str; 2] = ["hours", "sleep"];
const LABEL: &str = "passed";
const LEARNING_RATE: f64 = 0.5;
const STEPS: usize = 1000;
/// One row: the two features and the 0/1 label.
type Example = ([f64; 2], f64);
/// Splits the file into training and test rows. There is no CSV parser in
/// the standard library; the file has no quoted fields, so commas suffice.
fn read(path: &str) -> (Vec<Example>, Vec<Example>) {
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 (feature_ix, label_ix, split_ix) = (FEATURES.map(column), column(LABEL), column("split"));
let (mut train, mut test) = (Vec::new(), Vec::new());
for line in lines {
let fields: Vec<&str> = line.split(',').collect();
let example = (
feature_ix.map(|i| fields[i].parse().unwrap()),
fields[label_ix].parse().unwrap(),
);
if fields[split_ix] == "train" {
train.push(example);
} else {
test.push(example);
}
}
(train, test)
}
fn sigmoid(z: f64) -> f64 {
1.0 / (1.0 + (-z).exp())
}
fn dot(w: &[f64; 3], x: &[f64; 3]) -> f64 {
w.iter().zip(x).map(|(a, b)| a * b).sum()
}
fn main() {
let path = env::args().nth(1).expect("usage: main <data.csv>");
let (train, test) = read(&path);
let n = train.len() as f64;
// Scale each feature to mean 0 and (sample) standard deviation 1, using
// the training rows only.
let means: [f64; 2] = [0, 1].map(|j| train.iter().map(|(x, _)| x[j]).sum::<f64>() / n);
let sds: [f64; 2] = [0, 1].map(|j| {
let squares: f64 = train.iter().map(|(x, _)| (x[j] - means[j]).powi(2)).sum();
(squares / (n - 1.0)).sqrt()
});
// The leading 1 multiplies the intercept.
let scale = |x: &[f64; 2]| -> [f64; 3] {
[1.0, (x[0] - means[0]) / sds[0], (x[1] - means[1]) / sds[1]]
};
let xs: Vec<[f64; 3]> = train.iter().map(|(x, _)| scale(x)).collect();
let mut weights = [0.0; 3];
for _ in 0..STEPS {
// Gradient of the average log-loss: the mean of (p - y) times each
// input, computed in full before any weight changes.
let errors: Vec<f64> = xs
.iter()
.zip(&train)
.map(|(x, (_, y))| sigmoid(dot(&weights, x)) - y)
.collect();
let gradient: [f64; 3] = [0, 1, 2].map(|j| {
errors.iter().zip(&xs).map(|(e, x)| e * x[j]).sum::<f64>() / n
});
for (w, g) in weights.iter_mut().zip(gradient) {
*w -= LEARNING_RATE * g;
}
}
let log_loss = -xs
.iter()
.zip(&train)
.map(|(x, (_, y))| {
let p = sigmoid(dot(&weights, x));
y * p.ln() + (1.0 - y) * (1.0 - p).ln()
})
.sum::<f64>()
/ n;
let correct = test
.iter()
.filter(|(x, y)| {
let predicted = if sigmoid(dot(&weights, &scale(x))) >= 0.5 { 1.0 } else { 0.0 };
predicted == *y
})
.count();
println!("train_count {}", train.len());
println!("test_count {}", test.len());
println!("intercept {:.6}", weights[0]);
println!("weight_hours {:.6}", weights[1]);
println!("weight_sleep {:.6}", weights[2]);
println!("log_loss {log_loss:.6}");
println!("correct {correct}");
println!("accuracy {:.6}", correct as f64 / test.len() as f64);
}Fixed-size arrays let the compiler check the shapes. A raw row is
[f64; 2], a scaled row with its leading 1 is [f64; 3], and dot takes
only the second, so forgetting to scale a test row before predicting is a
type error: “expected an array with a size of 3, found one with a size of
2”. [0, 1].map(...) builds the means and standard deviations as arrays of
the same kind.
One precedence trap sits in sigmoid: -z.exp() means -(z.exp()), so
for z = 2 it gives -7.389056 instead of 0.135335. We write (-z).exp().
The standard library has no logistic regression, so there was no built-in
to skip.
// Logistic regression by gradient descent, using only Swift's standard
// library and Foundation, which ships with it.
//
// Predicts whether a student passed from study hours and sleep. Each feature
// is rescaled with the training rows' mean and standard deviation, then the
// weights start at zero and take 1000 fixed steps of batch gradient descent
// on the average log-loss.
import Foundation
let learningRate = 0.5
let steps = 1000
// The named column of a CSV file, as text. A function that can fail is
// marked throws, and every call to it is marked try.
func readColumn(_ path: String, _ name: String) throws -> [String] {
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)")
}
return lines.map { line in
String(line.split(separator: ",", omittingEmptySubsequences: false)[index])
}
}
// Double(_:) gives nil for text that isn't a number, and ! stops the
// program there instead of carrying on with a missing value.
func readNumbers(_ path: String, _ name: String) throws -> [Double] {
try readColumn(path, name).map { Double($0)! }
}
func sigmoid(_ z: Double) -> Double {
1 / (1 + exp(-z))
}
func dot(_ u: [Double], _ v: [Double]) -> Double {
zip(u, v).map { $0 * $1 }.reduce(0, +)
}
func mean(_ xs: [Double]) -> Double {
xs.reduce(0, +) / Double(xs.count)
}
// Sample standard deviation: divide by n - 1.
func sd(_ xs: [Double]) -> Double {
let m = mean(xs)
let squares = xs.map { ($0 - m) * ($0 - m) }.reduce(0, +)
return (squares / Double(xs.count - 1)).squareRoot()
}
let path = CommandLine.arguments[1]
let hours = try readNumbers(path, "hours")
let sleep = try readNumbers(path, "sleep")
let passed = try readNumbers(path, "passed")
let split = try readColumn(path, "split")
// indices is 0..<200; filter keeps the row numbers of each split.
let trainRows = split.indices.filter { split[$0] == "train" }
let testRows = split.indices.filter { split[$0] == "test" }
// Scale with the training rows only. The means and sds are [hours, sleep].
let features = [hours, sleep]
let means = features.map { column in mean(trainRows.map { column[$0] }) }
let sds = features.map { column in sd(trainRows.map { column[$0] }) }
// One scaled row: a leading 1 for the intercept, then each feature.
func scaled(_ row: Int) -> [Double] {
[1] + (0..<2).map { j in (features[j][row] - means[j]) / sds[j] }
}
let xs = trainRows.map(scaled)
let ys = trainRows.map { passed[$0] }
let n = Double(xs.count)
var weights = [0.0, 0.0, 0.0]
for _ in 0..<steps {
// Gradient of the average log-loss: the mean of (p - y) times each
// input, computed in full before any weight changes.
let errors = zip(xs, ys).map { x, y in sigmoid(dot(weights, x)) - y }
let gradient = (0..<3).map { j in
zip(errors, xs).map { e, x in e * x[j] }.reduce(0, +) / n
}
// A new array of weights replaces the old one.
weights = zip(weights, gradient).map { w, g in w - learningRate * g }
}
let logLoss = -zip(xs, ys).map { x, y in
let p = sigmoid(dot(weights, x))
return y * log(p) + (1 - y) * log(1 - p)
}.reduce(0, +) / n
let correct = testRows.filter { row in
let predicted = sigmoid(dot(weights, scaled(row))) >= 0.5 ? 1.0 : 0.0
return predicted == passed[row]
}.count
func show(_ name: String, _ value: Double) {
print(name, String(format: "%.6f", value))
}
print("train_count", trainRows.count)
print("test_count", testRows.count)
show("intercept", weights[0])
show("weight_hours", weights[1])
show("weight_sleep", weights[2])
show("log_loss", logLoss)
print("correct", correct)
show("accuracy", Double(correct) / Double(testRows.count))Each row is a [Double] with a leading 1 for the intercept, so dot is
one zip, map, and reduce(0, +). A step builds the whole gradient
first, then replaces the weights with a new array,
zip(weights, gradient).map { w, g in w - learningRate * g }, so every
weight moves from the same gradient.
Watch the lengths. zip stops at the shorter array, so
dot([1, 2, 3], [10, 20]) quietly gives 50 instead of an error: forget the
leading 1 on a test row and the prediction is just wrong. Foundation
supplies exp, log, and String(format:). Swift has no logistic
regression, so there was no built-in to skip.
# Logistic regression by gradient descent, using only Tcl's built-in
# commands.
#
# Predicts whether a student passed from study hours and sleep. Tcl has no
# regression or classifier, so every step is written out. Each feature is
# rescaled with the training rows' mean and standard deviation, then the
# weights start at zero and take 1000 fixed steps of batch gradient descent
# on the average log-loss. expr does the arithmetic; its {braces} keep Tcl
# from substituting twice.
set features {hours sleep}
set learning_rate 0.5
set steps 1000
# Returns two lists, train and test. Each example in them is a two-element
# list: the feature values, then the 0/1 label.
proc read_examples {path features} {
set f [open $path r]
set header [split [gets $f] ,]
set cols {}
foreach name [concat $features passed split] {
set col [lsearch -exact $header $name]
if {$col < 0} { error "column not found: $name" }
lappend cols $col
}
set train {}
set test {}
while {[gets $f line] >= 0} {
set fields [split $line ,]
set x {}
foreach col [lrange $cols 0 end-2] {
lappend x [lindex $fields $col]
}
set example [list $x [lindex $fields [lindex $cols end-1]]]
if {[lindex $fields [lindex $cols end]] eq "train"} {
lappend train $example
} else {
lappend test $example
}
}
close $f
return [list $train $test]
}
proc sigmoid {z} {
return [expr {1.0 / (1.0 + exp(-$z))}]
}
# foreach can walk two lists side by side, one item from each per pass.
proc dot {u v} {
set total 0.0
foreach a $u b $v {
set total [expr {$total + $a * $b}]
}
return $total
}
proc say {name value} {
puts "$name [format %.6f $value]"
}
# lassign unpacks a list into separate variables.
lassign [read_examples [lindex $argv 0] $features] train test
set n [llength $train]
# Scale each feature to mean 0 and (sample) standard deviation 1, using the
# training rows only. lmap builds a new list from what its body returns.
set means {}
set sds {}
for {set j 0} {$j < [llength $features]} {incr j} {
set column [lmap example $train { lindex $example 0 $j }]
set total 0.0
foreach v $column { set total [expr {$total + $v}] }
set mean [expr {$total / $n}]
set sumsq 0.0
foreach v $column { set sumsq [expr {$sumsq + ($v - $mean) ** 2}] }
lappend means $mean
lappend sds [expr {sqrt($sumsq / ($n - 1))}]
}
# The leading 1.0 multiplies the intercept.
proc scale {x means sds} {
set scaled {1.0}
foreach v $x m $means s $sds {
lappend scaled [expr {($v - $m) / $s}]
}
return $scaled
}
set xs [lmap example $train { scale [lindex $example 0] $means $sds }]
set ys [lmap example $train { lindex $example 1 }]
set weights {0.0 0.0 0.0}
for {set step 0} {$step < $steps} {incr step} {
# Gradient of the average log-loss: the mean of (p - y) times each input.
set gradient {0.0 0.0 0.0}
foreach x $xs y $ys {
set error [expr {[sigmoid [dot $weights $x]] - $y}]
set gradient [lmap g $gradient v $x { expr {$g + $error * $v} }]
}
# Build the new weights from the old ones, so every weight moves by the
# same gradient.
set weights [lmap w $weights g $gradient {
expr {$w - $learning_rate * $g / $n}
}]
}
set loss 0.0
foreach x $xs y $ys {
set p [sigmoid [dot $weights $x]]
set loss [expr {$loss - ($y * log($p) + (1 - $y) * log(1 - $p))}]
}
set correct 0
foreach example $test {
lassign $example x y
set p [sigmoid [dot $weights [scale $x $means $sds]]]
set predicted [expr {$p >= 0.5 ? 1 : 0}]
if {$predicted == $y} { incr correct }
}
puts "train_count $n"
puts "test_count [llength $test]"
say intercept [lindex $weights 0]
say weight_hours [lindex $weights 1]
say weight_sleep [lindex $weights 2]
say log_loss [expr {$loss / $n}]
puts "correct $correct"
say accuracy [expr {double($correct) / [llength $test]}]Each example is a two-element list, the feature values and the label, and
lassign unpacks it into separate variables. A for loop takes the 1000
steps. Inside it, foreach x $xs y $ys walks the inputs and labels side by
side, adding each error times input to a running gradient, and then lmap
builds the new weights from the old ones, so every weight moves by the same
gradient.
Tcl’s built-in commands include no regression, classifier, or matrix type,
so there was no built-in to skip; sigmoid, dot, and scale are short
procs. The accuracy uses double($correct), since 37 / 50 is integer
division in Tcl and gives 0.