B-preview writes when the script runsDictionary of Applied Machine Learning · logistic regression
Numerical companion to the entry logistic regression: it recomputes what the entry states and prints one line per check
Two binary classification trainsets. First, one year of weather records from the GeoSphere Austria station Krems an der Donau (station id 3805): each day is a data point with two numeric features — the maximum temperature of the day and that of the preceding day — and binary label +1 if the following day stays free of freezing temperatures (minimum temperature at least 0 deg C) and -1 otherwise; the entry opens with the scatter plot of these data points. Second, a synthetic one-feature trainset on which logistic regression is fit by a hand-written gradient descent loop: the average logistic loss decreases monotonically, the gradient vanishes at the learned parameters, the GD update leaves them (approximately) unchanged — its fixed point — and thresholding the learned hypothesis at zero classifies the trainset better than always answering with the majority label. Self-contained (numpy/matplotlib only, stdlib urllib for the download), fixed seed.
Run it without installing anything:uv run https://dictionaryofml.org/terms/logreg.py
uv downloads this script and the pinned NumPy and Matplotlib it needs, then runs it; the script fetches any input file it uses. To keep the output files, download logreg.py into a folder and run uv run logreg.py there. With NumPy and Matplotlib already installed, python3 logreg.py, from any directory — it writes its output files into the current directory. Fixed seeds, so the printed numbers reproduce exactly. Download logreg.py · Notebook · Open in Colab
One cell per block of the script: the code, and what that code printed when it last ran here
"""
logreg.py — numerical companion to the glossary entry 'logistic regression'.
Purpose
-------
Two binary classification trainsets. First, one year of weather records
from the GeoSphere Austria station Krems an der Donau (station id 3805):
each day is a data point with two numeric features — the maximum
temperature of the day and that of the preceding day — and binary label
+1 if the following day stays free of freezing temperatures (minimum
temperature at least 0 deg C) and -1 otherwise; the entry opens with the
scatter plot of these data points. Second, a synthetic one-feature
trainset on which logistic regression is fit by a hand-written gradient
descent loop: the average logistic loss decreases monotonically, the
gradient vanishes at the learned parameters, the GD update leaves them
(approximately) unchanged — its fixed point — and thresholding the
learned hypothesis at zero classifies the trainset better than always
answering with the majority label. Self-contained (numpy/matplotlib
only, stdlib urllib for the download), fixed seed.
Blocks
------
[B-fetch] Download the daily maximum and minimum temperature at Krems
for 2024 from the GeoSphere data hub; write the 366 records
to logreg_weather.csv and pin them to the archive.
[B-days] Build 364 data points (features: the maximum temperature of
the day and of the preceding day; label: +1 if the following
day's minimum temperature is at least 0 deg C, -1 if it
drops below); write every third day of each class to
logreg_nofreeze.csv and logreg_freeze.csv for the entry's
scatter plot (thinned so the markers do not overlap).
[B-fit2d] Fit the two-feature model h(x) = w^T x + b to the weather
data by gradient descent (in standardized coordinates, mapped
back to the original temperature units) and write its decision
boundary w^T x + b = 0 (logreg_boundary.csv) and the normal
vector w (logreg_normal.csv) for the entry's scatter plot. The
fitted classifier beats the constant majority rule.
[B-data] m = 40 data points with a single feature x drawn uniformly
from [-3, 3]; the binary label y in {-1, +1} is +1 with
probability sigmoid(2 x - 1), so the labels are noisy around
the point where 2 x - 1 = 0.
[B-gd] Gradient descent on the average logistic loss
f(w) = (1/m) sum_r log(1 + exp(-y^(r) w^T x^(r))) with the
constant feature 1 absorbing the intercept: the loss never
increases along the run, the gradient norm at the learned
parameters is small, and one further GD update moves them by
(approximately) nothing — the learned parameters are a fixed
point of the update.
[B-clf] The classifier sign(w^T x) obtained by thresholding the
learned hypothesis classifies a larger fraction of the
trainset correctly than the constant rule that always
answers with the majority label.
[B-loss] The logistic, squared-error and hinge loss of a data point
with label y = +1 as a function of the hypothesis value
h(x), on a grid; written to logreg_losses.csv for the
entry's loss-comparison figure. The squared-error loss turns
upward for a large correct h(x) while the logistic and hinge
losses do not.
[B-preview] The matplotlib preview of the scatter and the loss curves.
Outputs
-------
logreg_weather.csv : date, tlmax, tlmin at Krems, 366 days of 2024.
logreg_nofreeze.csv : xtoday,xyest for the thinned days with label +1.
logreg_freeze.csv : xtoday,xyest for the thinned days with label -1.
logreg_boundary.csv : x1,x2 — two endpoints of the decision boundary line.
logreg_normal.csv : x1,x2 — base and tip of the normal-vector w arrow.
logreg_losses.csv : h, logistic, squared, hinge — the three losses of a
data point with label y = +1 on a grid of h(x).
logreg.png : matplotlib preview of the figures (checking only).
"""
import json
import urllib.request
from pathlib import Path
import numpy as np
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
OUT_DIR = Path(__file__).parent
report = []
def check(name, ok):
report.append((name, bool(ok)))
print(f" [{'ok' if ok else 'FAIL'}] {name}")
def sigmoid(z):
return 1.0 / (1.0 + np.exp(-z))
Download the daily maximum and minimum temperature at Krems for 2024 from the GeoSphere data hub; write the 366 records to logreg_weather.csv and pin them to the archive.
PARAMS = ["tlmax", "tlmin"]
URL = ("https://dataset.api.hub.geosphere.at/v1/station/historical/"
f"klima-v2-1d?parameters={','.join(PARAMS)}&station_ids=3805"
"&start=2024-01-01&end=2024-12-31")
with urllib.request.urlopen(URL, timeout=120) as resp:
payload = json.load(resp)
cols = payload["features"][0]["properties"]["parameters"]
days = [stamp[:10] for stamp in payload["timestamps"]]
raw = np.array([cols[p]["data"] for p in PARAMS], dtype=float).T
with open(OUT_DIR / "logreg_weather.csv", "w") as f:
f.write("date," + ",".join(PARAMS) + "\n")
for day, row in zip(days, raw):
f.write(day + "," + ",".join(f"{v:g}" for v in row) + "\n")
check("[B-fetch] 366 days downloaded for 2024", len(days) == 366)
check("[B-fetch] no observation is missing", not np.isnan(raw).any())
check("[B-fetch] the record matches the archive (Jan 1: max 9.6, "
"min 2.6)", raw[0].tolist() == [9.6, 2.6])
[ok] [B-fetch] 366 days downloaded for 2024
[ok] [B-fetch] no observation is missing
[ok] [B-fetch] the record matches the archive (Jan 1: max 9.6, min 2.6)
364 days, 67 followed by a freezing day
Build 364 data points (features: the maximum temperature of the day and of the preceding day; label: +1 if the following day's minimum temperature is at least 0 deg C, -1 if it drops below); write every third day of each class to logreg_nofreeze.csv and logreg_freeze.csv for the entry's scatter plot (thinned so the markers do not overlap).
tmax, tmin = raw[:, 0], raw[:, 1]
xtoday = tmax[1:-1] # max temperature of the day
xyest = tmax[:-2] # max temperature of the preceding day
label = np.where(tmin[2:] >= 0.0, 1.0, -1.0)
n_freeze = int((label < 0).sum())
print(f" {len(label)} days, {n_freeze} followed by a freezing day")
check("[B-days] 364 data points with two features each",
len(xtoday) == 364 and len(xyest) == 364 and len(label) == 364)
check("[B-days] both labels occur", (label > 0).any() and (label < 0).any())
check("[B-days] 67 days are followed by a freezing day", n_freeze == 67)
THIN = 3 # plot every third day so markers show
for fname, keep in [("logreg_nofreeze.csv", label > 0),
("logreg_freeze.csv", label < 0)]:
with open(OUT_DIR / fname, "w") as f:
f.write("xtoday,xyest\n")
for a, b in zip(xtoday[keep][::THIN], xyest[keep][::THIN]):
f.write(f"{a:g},{b:g}\n")
[ok] [B-days] 364 data points with two features each [ok] [B-days] both labels occur [ok] [B-days] 67 days are followed by a freezing day
Fit the two-feature model h(x) = w^T x + b to the weather data by gradient descent (in standardized coordinates, mapped back to the original temperature units) and write its decision boundary w^T x + b = 0 (logreg_boundary.csv) and the normal vector w (logreg_normal.csv) for the entry's scatter plot. The fitted classifier beats the constant majority rule.
Xw = np.c_[xtoday, xyest] # two raw temperature features
muw, sdw = Xw.mean(axis=0), Xw.std(axis=0)
Zb = np.c_[(Xw - muw) / sdw, np.ones(len(Xw))] # standardized + intercept
def grad2(wz):
return -(Zb * (label * sigmoid(-label * (Zb @ wz)))[:, None]).mean(axis=0)
wz = np.zeros(3)
for _ in range(4000):
wz = wz - 0.5 * grad2(wz) # GD in standardized coordinates
w2d = wz[:2] / sdw # weights in the original units
b2d = float(wz[2] - np.sum(wz[:2] * muw / sdw)) # offset in the original units
acc2d = float(np.mean(np.sign(Xw @ w2d + b2d) == label))
maj_acc = float(max((label > 0).mean(), (label < 0).mean()))
check("[B-fit2d] the fitted two-feature classifier beats the majority rule",
acc2d > maj_acc)
def x2_on_boundary(x1): # w2d . x + b2d = 0 as x2(x1)
return -(w2d[0] * x1 + b2d) / w2d[1]
# clip the boundary line to the data bounding box (the two features are highly
# collinear, so the line is steep and leaves the plotted range quickly)
x1lo, x1hi = float(xtoday.min()), float(xtoday.max())
x2lo, x2hi = float(xyest.min()), float(xyest.max())
grid = np.linspace(x1lo, x1hi, 800)
on = (x2_on_boundary(grid) >= x2lo) & (x2_on_boundary(grid) <= x2hi)
xe = (float(grid[on][0]), float(grid[on][-1])) # x1 of the two box crossings
with open(OUT_DIR / "logreg_boundary.csv", "w") as f:
f.write("x1,x2\n")
for x1 in xe:
f.write(f"{x1:.4f},{x2_on_boundary(x1):.4f}\n")
# the weight vector w (unit length, pointing to y=+1), based at the point of
# the boundary nearest the data centroid so the arrow sits in the data cloud
wn = w2d / float(np.linalg.norm(w2d))
cen = Xw.mean(axis=0)
base = cen - (float(w2d @ cen) + b2d) / float(w2d @ w2d) * w2d
NLEN = 7.0 # arrow length in deg C
with open(OUT_DIR / "logreg_normal.csv", "w") as f:
f.write("x1,x2\n")
f.write(f"{base[0]:.4f},{base[1]:.4f}\n")
f.write(f"{base[0] + NLEN * wn[0]:.4f},{base[1] + NLEN * wn[1]:.4f}\n")
[ok] [B-fit2d] the fitted two-feature classifier beats the majority rule
m = 40 data points with a single feature x drawn uniformly from [-3, 3]; the binary label y in {-1, +1} is +1 with probability sigmoid(2 x - 1), so the labels are noisy around the point where 2 x - 1 = 0.
rng = np.random.default_rng(3)
m = 40
x = rng.uniform(-3.0, 3.0, m)
p_true = sigmoid(2.0 * x - 1.0)
y = np.where(rng.uniform(size=m) < p_true, 1.0, -1.0)
X = np.c_[x, np.ones(m)] # constant feature absorbs the intercept
check("[B-data] both labels occur", (y > 0).any() and (y < 0).any())
[ok] [B-data] both labels occur
Gradient descent on the average logistic loss f(w) = (1/m) sum_r log(1 + exp(-y^(r) w^T x^(r))) with the constant feature 1 absorbing the intercept: the loss never increases along the run, the gradient norm at the learned parameters is small, and one further GD update moves them by (approximately) nothing — the learned parameters are a fixed point of the update.
def avg_logloss(w):
return float(np.mean(np.log1p(np.exp(-y * (X @ w)))))
def grad(w):
return -(X * (y * sigmoid(-y * (X @ w)))[:, None]).mean(axis=0)
eta = 0.5
w = np.zeros(2)
losses = [avg_logloss(w)]
for _ in range(2000):
w = w - eta * grad(w) # the GD update w <- T(w)
losses.append(avg_logloss(w))
check("[B-gd] the average logistic loss never increases and ends "
"below its start",
all(b <= a + 1e-12 for a, b in zip(losses, losses[1:]))
and losses[-1] < losses[0])
check("[B-gd] the gradient vanishes at the learned parameters",
float(np.linalg.norm(grad(w))) < 1e-4)
w_next = w - eta * grad(w) # one further update
check("[B-gd] the learned parameters are a fixed point of the update",
float(np.linalg.norm(w_next - w)) < 1e-4)
[ok] [B-gd] the average logistic loss never increases and ends below its start
[ok] [B-gd] the gradient vanishes at the learned parameters
[ok] [B-gd] the learned parameters are a fixed point of the update
correctly classified: logistic regression 0.90, majority label 0.60
The classifier sign(w^T x) obtained by thresholding the learned hypothesis classifies a larger fraction of the trainset correctly than the constant rule that always answers with the majority label.
frac_lr = float(np.mean(np.sign(X @ w) == y))
majority = 1.0 if (y > 0).sum() >= (y < 0).sum() else -1.0
frac_const = float(np.mean(y == majority))
print(f" correctly classified: logistic regression {frac_lr:.2f}, "
f"majority label {frac_const:.2f}")
check("[B-clf] sign(w^T x) beats the constant majority rule",
frac_lr > frac_const)
[ok] [B-clf] sign(w^T x) beats the constant majority rule
The logistic, squared-error and hinge loss of a data point with label y = +1 as a function of the hypothesis value h(x), on a grid; written to logreg_losses.csv for the entry's loss-comparison figure. The squared-error loss turns upward for a large correct h(x) while the logistic and hinge losses do not.
h = np.linspace(-3.0, 3.0, 241) # symmetric about the boundary h(x) = 0
loss_logistic = np.log1p(np.exp(-h)) # log(1 + exp(-y h)) at y = +1
loss_squared = (1.0 - h) ** 2 # (y - h)^2 at y = +1
loss_hinge = np.maximum(0.0, 1.0 - h) # max(0, 1 - y h) at y = +1
# the squared-error loss is the only one that penalizes a large correct h(x)
i_correct = h > 1.5
check("[B-loss] the squared-error loss rises for a large correct h(x), "
"the logistic and hinge losses fall",
loss_squared[i_correct][-1] > loss_squared[i_correct][0]
and loss_logistic[i_correct][-1] < loss_logistic[i_correct][0]
and loss_hinge[i_correct][-1] <= loss_hinge[i_correct][0])
with open(OUT_DIR / "logreg_losses.csv", "w") as f:
f.write("h,logistic,squared,hinge\n")
for hi, a, b, c in zip(h, loss_logistic, loss_squared, loss_hinge):
f.write(f"{hi:.4f},{a:.4f},{b:.4f},{c:.4f}\n")
[ok] [B-loss] the squared-error loss rises for a large correct h(x), the logistic and hinge losses fall 13/13 checks pass wrote logreg_boundary.csv, logreg_normal.csv, logreg_losses.csv, logreg.png in /Users/junga1/dictionaryappliedml/pythondemos
The matplotlib preview of the scatter and the loss curves.
fig, (ax0, ax) = plt.subplots(1, 2, figsize=(10.0, 3.8))
ax0.plot(xtoday[label > 0][::THIN], xyest[label > 0][::THIN], "o",
color="0.6", ms=5, label="no freeze the next day ($y = +1$)")
ax0.plot(xtoday[label < 0][::THIN], xyest[label < 0][::THIN], "^",
color="black", mfc="none", ms=7, label="freeze the next day ($y = -1$)")
ax0.plot([xe[0], xe[1]], [x2_on_boundary(xe[0]), x2_on_boundary(xe[1])], "k--",
lw=1.4, label="decision boundary")
tipx, tipy = base[0] + NLEN * wn[0], base[1] + NLEN * wn[1]
ax0.annotate("", xy=(tipx, tipy), xytext=(base[0], base[1]),
arrowprops=dict(arrowstyle="->", lw=1.6))
ax0.text(tipx, tipy, r" $\mathbf{w}$", fontsize=11)
ax0.set_xlim(x1lo - 2, x1hi + 2)
ax0.set_ylim(x2lo - 2, x2hi + 2)
ax0.set_xlabel("max temperature of the day (deg C)")
ax0.set_ylabel("max temperature of the preceding day (deg C)")
ax0.set_aspect("equal")
ax0.set_title("Krems: will the next day bring frost?")
ax0.legend(frameon=False, fontsize=8, loc="upper left")
ax.plot(h, loss_logistic, "k-", lw=1.6, label="logistic loss")
ax.plot(h, loss_squared, "k--", lw=1.6, label="squared error loss")
ax.plot(h, loss_hinge, "k:", lw=1.8, label="hinge loss")
ax.axvline(0.0, color="0.5", ls="--", lw=1.2)
ax.annotate("", xy=(2.9, -0.55), xytext=(0.15, -0.55),
arrowprops=dict(arrowstyle="->", lw=1.4))
ax.annotate("", xy=(-2.9, -0.55), xytext=(-0.15, -0.55),
arrowprops=dict(arrowstyle="->", lw=1.4))
ax.text(1.5, -0.9, "correct", fontsize=9, ha="center")
ax.text(-1.5, -0.9, "wrong", fontsize=9, ha="center")
ax.set_xlabel("hypothesis value $h(x)$")
ax.set_ylabel("loss at label $y = +1$")
ax.set_xlim(-3, 3)
ax.set_ylim(-1.2, 4.2)
ax.set_title("loss vs $h(x)$ for a data point with $y = +1$")
ax.legend(frameon=False, loc="upper right", fontsize=8)
fig.tight_layout()
fig.savefig(OUT_DIR / "logreg.png", dpi=110)
n_ok = sum(ok for _, ok in report)
print(f"\n{n_ok}/{len(report)} checks pass")
print(f"wrote logreg_boundary.csv, logreg_normal.csv, logreg_losses.csv, "
f"logreg.png in {OUT_DIR}")
if n_ok != len(report):
raise SystemExit(1)
B-preview writes when the script runs