Regression for yes/no outcomes: model the log-odds as a straight line, fit by maximum likelihood with Newton's method, and read the coefficients as odds ratios.
Many of the questions people most want answered have two outcomes. Will this customer renew? Will this loan default? Is this email spam? Did the treatment work? Logistic regression is the standard model for all of them: it's how credit scores are built, how clinical studies report risk factors, and the simplest classifier in machine learning — a single-layer neural network with a sigmoid output is exactly this model.
It's also where regression stops having a formula. There is no closed-form answer, so the fit is an algorithm, and this module has you write it.
Why not a straight line? Fit least squares to 0/1 outcomes and the line happily predicts probabilities below 0 and above 1, and its errors can't be normal or equally spread. The fix is to model something that can be any number.
Odds and log-odds. If p is the probability of "yes", the odds are p/(1 − p) — 0.8 is odds of 4 to 1 — and the log-odds (or logit) is ln(p/(1 − p)), which runs from −∞ to ∞. Logistic regression says the log-odds is a straight line: ln(p/(1 − p)) = b₀ + b₁x₁ + … Turning that back into a probability is the sigmoid, p = 1/(1 + e^(−η)), an S-curve from 0 to 1.
Reading coefficients. A coefficient b is the change in log-odds per unit of x. Exponentiated, e^b is an odds ratio: each extra unit of x multiplies the odds by e^b. An odds ratio of 1.34 per month means each extra month of membership makes renewal 34% more likely in odds terms — not 34 percentage points.
Maximum likelihood. Instead of minimizing squared error, choose the coefficients that make the observed outcomes most probable. The log-likelihood is Σ[y·ln p + (1 − y)·ln(1 − p)]; higher is better, and its maximum has no closed form.
Newton's method. At each step, the gradient Xᵀ(y − p) says which way is uphill, and the information matrix XᵀWX (with W = p(1 − p) on the diagonal) says how curved the hill is. The step (XᵀWX)⁻¹Xᵀ(y − p) jumps to where a parabola fitted at the current point peaks. It converges in a handful of steps, and the inverse information matrix at the end gives the standard errors for free. Because each step is a weighted least-squares fit, it's also called iteratively reweighted least squares.
Separation. If some line splits the 0s from the 1s perfectly, the likelihood keeps rising as the coefficients grow, forever: there is no maximum. Newton's method doesn't converge, and that's the honest answer. Real software warns about it; yours returns an error.
From probabilities to decisions. A model gives probabilities; acting on them needs a threshold. The confusion matrix counts the four outcomes — true and false positives, true and false negatives — and moving the threshold trades one kind of error for the other.
| Quantity | Formula | In words |
|---|---|---|
| Odds | p / (1 − p) | 0.8 → 4 to 1 |
| Logit | ln(p / (1 − p)) = b₀ + b₁x₁ + … | the linear part |
| Sigmoid | p = 1 / (1 + e^(−η)) | back to a probability |
| Odds ratio | e^b | odds multiplier per unit of x |
| Log-likelihood | Σ [y·η − ln(1 + e^η)] | stable form, η = Xb |
| Newton step | (XᵀWX)⁻¹ Xᵀ(y − p), W = diag(p(1 − p)) | |
| Standard errors | √diag((XᵀWX)⁻¹) at the fit | |
| Wald test | z = b / SE, p = erfc(abs(z) / √2) |
math.Exp(800) is +Inf, and 1/(1+math.Exp(800)) is 0 — fine. But e/(1+e) with e = +Inf is Inf/Inf = NaN. Write the sigmoid in the form that degrades to 0 and 1 rather than NaN.math.Max(η, 0) + math.Log1p(math.Exp(-math.Abs(η))), never rounds to infinity.for iter := 1; iter <= 100; iter++ with a "step smaller than 1e-10" exit is the shape of every numerical solver. The cap turns "runs forever" into a reportable error.var ErrNoConvergence = errors.New(...) lets callers test errors.Is(err, logistic.ErrNoConvergence) and react, instead of parsing a message.(tp, fp, tn, fn int) — document themselves at the call site: tp, fp, tn, fn := Confusion(y, probs, 0.5).package main import ( "fmt" "math" ) func sigmoid(z float64) float64 { return 1 / (1 + math.Exp(-z)) } func main() { // Hours studied and whether each of 12 students passed (made-up data). hours := []float64{0.5, 1, 1.5, 2, 2.5, 3, 3.5, 4, 4.5, 5, 5.5, 6} passed := []float64{0, 0, 0, 0, 1, 0, 1, 0, 1, 1, 1, 1} n := float64(len(hours)) // A straight line through the 0s and 1s (least squares). var mx, my float64 for i := range hours { mx += hours[i] my += passed[i] } mx, my = mx/n, my/n var sxy, sxx float64 for i := range hours { sxy += (hours[i] - mx) * (passed[i] - my) sxx += (hours[i] - mx) * (hours[i] - mx) } slope := sxy / sxx icept := my - slope*mx // Logistic regression by Newton's method: two coefficients, so the // 2×2 information matrix can be inverted by hand. b0, b1 := 0.0, 0.0 for range 25 { var g0, g1, i00, i01, i11 float64 for i, x := range hours { p := sigmoid(b0 + b1*x) w := p * (1 - p) g0 += passed[i] - p g1 += (passed[i] - p) * x i00 += w i01 += w * x i11 += w * x * x } det := i00*i11 - i01*i01 b0 += (i11*g0 - i01*g1) / det b1 += (i00*g1 - i01*g0) / det } fmt.Println("hours straight line logistic") for _, h := range []float64{0, 2, 4, 6, 8} { fmt.Printf("%5.0f %13.2f %8.2f\n", h, icept+slope*h, sigmoid(b0+b1*h)) } fmt.Printf("\nlogistic: log-odds = %.2f + %.2f·hours; each hour multiplies the odds of passing by %.2f\n", b0, b1, math.Exp(b1)) }
It prints:
The straight line says a student who studies zero hours has a −18% chance of passing and one who studies eight hours a 150% chance. The logistic curve bends to stay between 0 and 1, rising fastest in the middle where the outcome is most uncertain. Its slope reads naturally as odds: each extra hour multiplies the odds of passing by about four.
Part A is the logistic function, its inverse, and a numerically stable log-likelihood. Part B fits the model by Newton's method and reports probabilities, odds ratios, and p-values. Part C turns probabilities into decisions with a confusion matrix.
go test ./10-logistic-regression/
The expected coefficients, standard errors, log-likelihood, and p-values come from statsmodels' Logit on 30 made-up subscribers. TestFitSeparable hands Fit perfectly separable data, and passes only if you return ErrNoConvergence rather than looping forever or returning garbage.
Everything below is also in the course download, in 10-logistic-regression/. Work in exercises.go; the tests decide when you’re done.
// Package logistic is Module 10: logistic regression. // // When the outcome is yes/no — renewed or not, defaulted or not, spam or // not — a straight line predicts impossible probabilities. Logistic // regression models the log-odds as a straight line instead, and is fitted // by maximum likelihood rather than least squares. package logistic import ( "errors" "github.com/bozcode/go-statistics-and-regression/internal/linalg" ) // Model is a fitted logistic regression. Index 0 is the intercept. type Model struct { Coef []float64 // b₀, b₁, …: changes in log-odds SE []float64 // standard error of each coefficient LogLik float64 // log-likelihood at the fitted coefficients Iterations int // Newton steps taken } // ErrNoConvergence is returned by Fit when Newton's method doesn't settle, // which usually means the classes are perfectly separable: some line splits // the 0s from the 1s exactly, and the coefficients run off to infinity. var ErrNoConvergence = errors.New("logistic: Newton's method did not converge") // The helpers below are given: a design matrix with a leading 1s column, as // in Module 8, and the matrix tools from internal/linalg. // designMatrix returns the rows [1, cols[0][i], cols[1][i], …]. func designMatrix(cols ...[]float64) [][]float64 { x := make([][]float64, len(cols[0])) for i := range x { x[i] = make([]float64, len(cols)+1) x[i][0] = 1 for j, c := range cols { x[i][j+1] = c[i] } } return x } // inverse returns a⁻¹, or an error if a is singular. func inverse(a [][]float64) ([][]float64, error) { return linalg.Inverse(a) } // mulVec returns the matrix-vector product av. func mulVec(a [][]float64, v []float64) []float64 { return linalg.MulVec(a, v) }
//go:build !solutions package logistic // ============================================================ // Part A — the logistic function and the likelihood // ============================================================ // Sigmoid returns 1 / (1 + e^(−z)), squashing any number into [0, 1]. // Sigmoid(0) = 0.5. Write it in exactly that form: at z = −800, e^800 is // +Inf and 1/Inf is 0, which is right. The look-alike e^z / (1 + e^z) is // not safe — at z = 800 it's Inf/Inf = NaN. func Sigmoid(z float64) float64 { // TODO return 0 } // Logit is the inverse of Sigmoid: the log-odds ln(p / (1 − p)). // Logit(0.5) = 0; Logit(0.9) ≈ 2.197. func Logit(p float64) float64 { // TODO return 0 } // LogLikelihood returns Σ [y·ln p + (1 − y)·ln(1 − p)] with p = Sigmoid(η), // η = row·coef. Computed naively it fails when p rounds to 0 or 1; use the // equivalent, stable form Σ [y·η − softplus(η)], where // softplus(η) = ln(1 + e^η) = max(η, 0) + log1p(e^(−|η|)). func LogLikelihood(x [][]float64, y, coef []float64) float64 { // TODO return 0 } // ============================================================ // Part B — fitting by Newton's method // ============================================================ // Fit finds the coefficients that maximize the log-likelihood, by Newton's // method (also called iteratively reweighted least squares). Start with // every coefficient 0, then repeat up to 100 times: // // p_i = Sigmoid(row_i · coef) // g = Xᵀ(y − p) the gradient // I = Xᵀ W X, W = diag(p(1 − p)) the information matrix // step = I⁻¹ g; coef += step // // Stop when every |step| is below 1e-10. SE[j] = √((I⁻¹)[j][j]) using I at // the final coefficients. Return ErrNoConvergence if it never settles (or // if I can't be inverted). Use the given inverse and mulVec. func Fit(x [][]float64, y []float64) (Model, error) { // TODO return Model{}, nil } // Prob returns the predicted probability for one observation's predictor // values (without the leading 1): Sigmoid(b₀ + b₁x₁ + …). func (m Model) Prob(x []float64) float64 { // TODO return 0 } // OddsRatios returns e^b for every coefficient. For a predictor, it's the // factor the odds are multiplied by for each one-unit increase; for the // intercept, it's the odds when every predictor is 0. func (m Model) OddsRatios() []float64 { // TODO return nil } // PValues returns the two-sided Wald p-value for each coefficient: // z = b / SE, p = erfc(|z| / √2) — the normal distribution's two tails. func (m Model) PValues() []float64 { // TODO return nil } // ============================================================ // Part C — from probabilities to decisions // ============================================================ // Confusion counts outcomes when you predict 1 for every probability >= // threshold: true positives, false positives, true negatives, and false // negatives. y holds the actual 0/1 outcomes. func Confusion(y, probs []float64, threshold float64) (tp, fp, tn, fn int) { // TODO return 0, 0, 0, 0 } // Accuracy returns the share of predictions that were right, // (tp + tn) / total, at the given threshold. func Accuracy(y, probs []float64, threshold float64) float64 { // TODO return 0 }
package logistic import ( "errors" "math" "testing" ) // near reports whether a and b agree to about 9 significant digits. // Expected values come from SciPy and statsmodels' Logit. func near(a, b float64) bool { return math.Abs(a-b) <= 1e-9*math.Max(1, math.Abs(b)) } func slicesNear(t *testing.T, name string, got, want []float64, tol float64) { t.Helper() if len(got) != len(want) { t.Fatalf("%s has %d values, want %d", name, len(got), len(want)) } for i := range want { if math.Abs(got[i]-want[i]) > tol*math.Max(1, math.Abs(want[i])) { t.Errorf("%s[%d] = %v, want %v", name, i, got[i], want[i]) } } } // 30 made-up subscribers: months active, support tickets filed, and whether // they renewed (1) or not (0). var ( months = []float64{23, 6, 20, 4, 13, 4, 17, 20, 11, 23, 20, 9, 14, 18, 20, 22, 4, 18, 4, 22, 6, 21, 8, 23, 13, 8, 7, 15, 9, 17} tickets = []float64{1, 0, 0, 1, 4, 5, 2, 3, 2, 1, 3, 0, 2, 0, 1, 1, 0, 2, 4, 3, 4, 5, 0, 4, 0, 2, 0, 0, 2, 2} renewed = []float64{1, 1, 1, 0, 1, 0, 1, 0, 0, 1, 1, 1, 0, 1, 1, 1, 0, 1, 0, 1, 0, 0, 1, 1, 1, 1, 1, 1, 0, 1} ) // fitted fits the subscriber model, stopping the test (not the whole run) // if Fit isn't written yet. func fitted(t *testing.T) Model { t.Helper() m, err := Fit(designMatrix(months, tickets), renewed) if err != nil { t.Fatalf("Fit: %v", err) } if len(m.Coef) != 3 || len(m.SE) != 3 { t.Fatalf("Fit returned %d coefficients and %d standard errors, want 3 and 3", len(m.Coef), len(m.SE)) } return m } func TestSigmoid(t *testing.T) { for _, tc := range []struct{ z, want float64 }{ {0, 0.5}, {2, 0.8807970779778823}, {-3.5, 0.02931223075135632}, {40, 1}, {-800, 0}, {800, 1}, } { if got := Sigmoid(tc.z); math.IsNaN(got) || !near(got, tc.want) { t.Errorf("Sigmoid(%v) = %v, want %v", tc.z, got, tc.want) } } } func TestLogit(t *testing.T) { if got := Logit(0.9); !near(got, 2.1972245773362196) { t.Errorf("Logit(0.9) = %v, want 2.1972245773362196", got) } if got := Logit(0.2); !near(got, -1.3862943611198906) { t.Errorf("Logit(0.2) = %v, want -1.3862943611198906", got) } for _, p := range []float64{0.01, 0.3, 0.77} { if got := Sigmoid(Logit(p)); !near(got, p) { t.Errorf("Sigmoid(Logit(%v)) = %v", p, got) } } } func TestLogLikelihood(t *testing.T) { x := designMatrix(months, tickets) if got := LogLikelihood(x, renewed, []float64{0, 0, 0}); !near(got, -20.79441541679836) { t.Errorf("LogLikelihood at all-zero coefficients = %v, want 30·ln(0.5) = -20.79441541679836", got) } // p rounds to exactly 1 but y is 0: ln(1 − p) would be −Inf got := LogLikelihood([][]float64{{1, 40}}, []float64{0}, []float64{0, 1}) if math.IsInf(got, 0) || math.IsNaN(got) || math.Abs(got+40) > 1e-9 { t.Errorf("LogLikelihood with η = 40, y = 0 = %v, want about -40 (finite)", got) } } func TestFit(t *testing.T) { m, err := Fit(designMatrix(months, tickets), renewed) if err != nil { t.Fatalf("Fit: %v", err) } slicesNear(t, "Coef", m.Coef, []float64{-0.8887330151170617, 0.2909171504251711, -1.1705918322717863}, 1e-8) slicesNear(t, "SE", m.SE, []float64{1.2070123488098854, 0.12430149023975594, 0.5152566017841241}, 1e-8) if !near(m.LogLik, -11.028737060725605) { t.Errorf("LogLik = %v, want -11.028737060725605", m.LogLik) } if m.Iterations < 3 || m.Iterations > 30 { t.Errorf("Iterations = %d; Newton's method should take a handful", m.Iterations) } } func TestFitSeparable(t *testing.T) { // x = 1, 2 never renew and x = 3, 4 always do: the slope wants to be infinite _, err := Fit(designMatrix([]float64{1, 2, 3, 4}), []float64{0, 0, 1, 1}) if !errors.Is(err, ErrNoConvergence) { t.Errorf("Fit on perfectly separated data: err = %v, want ErrNoConvergence", err) } } func TestProb(t *testing.T) { m := fitted(t) if got := m.Prob([]float64{12, 2}); math.Abs(got-0.564904004087808) > 1e-8 { t.Errorf("Prob(12 months, 2 tickets) = %v, want 0.564904004087808", got) } if got := m.Prob([]float64{3, 4}); math.Abs(got-0.009028058400745104) > 1e-8 { t.Errorf("Prob(3 months, 4 tickets) = %v, want 0.009028058400745104", got) } } func TestOddsAndPValues(t *testing.T) { m := fitted(t) slicesNear(t, "OddsRatios", m.OddsRatios(), []float64{0.411176377125135, 1.3376537553141499, 0.3101833104382333}, 1e-8) slicesNear(t, "PValues", m.PValues(), []float64{0.4615431947538866, 0.019262290105444956, 0.023094853913954324}, 1e-7) } func TestConfusion(t *testing.T) { y := []float64{1, 0, 1, 1, 0, 0, 1} probs := []float64{0.9, 0.6, 0.4, 0.7, 0.2, 0.5, 0.5} if tp, fp, tn, fn := Confusion(y, probs, 0.5); tp != 3 || fp != 2 || tn != 1 || fn != 1 { t.Errorf("Confusion at 0.5 = (%d, %d, %d, %d), want (3, 2, 1, 1)", tp, fp, tn, fn) } if got := Accuracy(y, probs, 0.5); !near(got, 4.0/7) { t.Errorf("Accuracy at 0.5 = %v, want 4/7", got) } // the fitted model on its own data (statsmodels: 18, 3, 7, 2) m := fitted(t) fitted := make([]float64, len(months)) for i := range months { fitted[i] = m.Prob([]float64{months[i], tickets[i]}) } if tp, fp, tn, fn := Confusion(renewed, fitted, 0.5); tp != 18 || fp != 3 || tn != 7 || fn != 2 { t.Errorf("Confusion of the fitted model = (%d, %d, %d, %d), want (18, 3, 7, 2)", tp, fp, tn, fn) } }
Try not to open solution.go until your tests pass, or until you’ve been stuck for a solid half hour. Then compare: different-but-correct is fine; notice where the reference is simpler. go test -tags solutions ./10-logistic-regression/ runs the same tests against it.
solutions · 139 lines›//go:build solutions package logistic import "math" // ============================================================ // Part A — the logistic function and the likelihood // ============================================================ func Sigmoid(z float64) float64 { return 1 / (1 + math.Exp(-z)) } func Logit(p float64) float64 { return math.Log(p / (1 - p)) } func LogLikelihood(x [][]float64, y, coef []float64) float64 { ll := 0.0 for i, row := range x { eta := 0.0 for j, v := range row { eta += v * coef[j] } ll += y[i]*eta - (math.Max(eta, 0) + math.Log1p(math.Exp(-math.Abs(eta)))) } return ll } // ============================================================ // Part B — fitting by Newton's method // ============================================================ func Fit(x [][]float64, y []float64) (Model, error) { p := len(x[0]) coef := make([]float64, p) for iter := 1; iter <= 100; iter++ { info, grad := information(x, y, coef) inv, err := inverse(info) if err != nil { return Model{}, ErrNoConvergence } step := mulVec(inv, grad) biggest := 0.0 for j := range coef { coef[j] += step[j] biggest = math.Max(biggest, math.Abs(step[j])) } if biggest < 1e-10 { info, _ = information(x, y, coef) inv, err = inverse(info) if err != nil { return Model{}, ErrNoConvergence } m := Model{Coef: coef, SE: make([]float64, p), Iterations: iter, LogLik: LogLikelihood(x, y, coef)} for j := range m.SE { m.SE[j] = math.Sqrt(inv[j][j]) } return m, nil } } return Model{}, ErrNoConvergence } // information returns XᵀWX and Xᵀ(y − p) at coef. func information(x [][]float64, y, coef []float64) (info [][]float64, grad []float64) { k := len(coef) info = make([][]float64, k) for a := range info { info[a] = make([]float64, k) } grad = make([]float64, k) for i, row := range x { eta := 0.0 for j, v := range row { eta += v * coef[j] } p := Sigmoid(eta) w := p * (1 - p) for a := range row { grad[a] += row[a] * (y[i] - p) for b := range row { info[a][b] += row[a] * w * row[b] } } } return info, grad } func (m Model) Prob(x []float64) float64 { eta := m.Coef[0] for j, v := range x { eta += m.Coef[j+1] * v } return Sigmoid(eta) } func (m Model) OddsRatios() []float64 { out := make([]float64, len(m.Coef)) for j, b := range m.Coef { out[j] = math.Exp(b) } return out } func (m Model) PValues() []float64 { out := make([]float64, len(m.Coef)) for j := range out { out[j] = math.Erfc(math.Abs(m.Coef[j]/m.SE[j]) / math.Sqrt2) } return out } // ============================================================ // Part C — from probabilities to decisions // ============================================================ func Confusion(y, probs []float64, threshold float64) (tp, fp, tn, fn int) { for i, p := range probs { predicted, actual := p >= threshold, y[i] == 1 switch { case predicted && actual: tp++ case predicted && !actual: fp++ case !predicted && !actual: tn++ default: fn++ } } return tp, fp, tn, fn } func Accuracy(y, probs []float64, threshold float64) float64 { tp, fp, tn, fn := Confusion(y, probs, threshold) return float64(tp+tn) / float64(tp+fp+tn+fn) }
Four questions. Answer each one in your head — or out loud — before you open the answer.
A coefficient for "months active" is 0.29. What is the odds ratio, and what does it mean?
B — Coefficients are changes in log-odds; exponentiate to get the odds multiplier. It's a 34% increase in the odds, not in the probability.
Why is logistic regression fitted with Newton's method rather than a formula like b = (XᵀX)⁻¹Xᵀy?
B — Setting the gradient to zero gives equations with the sigmoid inside them, which can't be solved algebraically. Each Newton step solves a weighted least-squares problem instead.
Every customer with more than 12 months renewed and every one with fewer didn't. What happens when you fit months as a predictor?
B — A steeper and steeper sigmoid keeps raising the likelihood, so there's no finite maximum.
A fraud model has 99% accuracy, but only 1% of transactions are fraud. What should you check?
B — Predicting the majority class every time scores 99% here and catches no fraud. Count the true positives and false negatives.
That’s the course. If every test passes, go run ./cmd/progress says so — and you’ve built every statistic and model in it from scratch.