Logistic Regression

Logistic regression is the go-to algorithm for binary classification. It's fast, interpretable, and probabilistic. Unlike methods that just give you a hard yes/no decision, logistic regression tells you “I'm 73% confident this is spam” or “There's a 15% chance this customer will churn.” That probability estimate is incredibly valuable in real-world applications where you need to make risk-based decisions.

Here's the issue with using linear regression for classification: linear regression outputs any real number from negative infinity to positive infinity. But for classification, we need probabilities, which must be between 0 and 1.

If we try to use linear regression for binary classification (treating 0 and 1 as our targets), we'd get predictions like -0.3 or 1.7, which don't make sense as probabilities. We need a way to squash our linear output into the [0, 1] range.

The sigmoid function is our magic transformation. It takes any real number and maps it to a value between 0 and 1:

1def sigmoid(z):
2    """Apply sigmoid function with numerical stability."""
3    return 1 / (1 + np.exp(-np.clip(z, -500, 500)))

Notice the np.clip? That prevents numerical overflow when z gets extremely large or small. Without it, np.exp(-500) would cause issues. These little details matter when you're building from scratch.

Instead of our old linear regression model ŷ = Xw + b, we now work with:

We're still computing a linear combination of features, but we're passing it through the sigmoid. The linear part (Xw + b) is called the logit or log-odds, and the sigmoid transforms it into a probability.

The output ŷ now represents P(y = 1 | X) — the probability that the sample belongs to class 1 given the features. Consequently, P(y = 0 | X) = 1 − ŷ.

Mean squared error doesn't work well for logistic regression. Why? Because when we're predicting probabilities and our true labels are 0 or 1, MSE creates a non-convex cost surface with local minima. We need a cost function specifically designed for probability predictions.

Cross-entropy comes from information theory and measures the difference between two probability distributions. It quantifies how “surprised” we are by the true label given our predictions. If we predict high probability for the correct class, cross-entropy is low. If we predict low probability, cross-entropy is high (we're surprised).

For binary classification, this becomes binary cross-entropy or log loss:

This looks intimidating, but it's actually simple to follow:

1# In the fit method
2z_full = X @ self.coefficients_
3predictions_full = sigmoid(z_full)
4
5# Add small epsilon to prevent log(0)
6cost = -np.mean(
7    y * np.log(predictions_full + 1e-15) + 
8    (1 - y) * np.log(1 - predictions_full + 1e-15)
9)

Something notable about cross-entropy with sigmoid is that the gradient has an incredibly clean form:

Compare this to linear regression's gradient: ∂MSE/∂w = (2/n) X^T(Xw − y). They're almost identical! The difference is we're computing errors on the probability scale (after sigmoid) rather than the raw prediction scale.

1for i in range(0, n_samples, batch_size):
2    batch_indices = indices[i:i+batch_size]
3    X_batch = X[batch_indices]
4    y_batch = y[batch_indices]
5    
6    # Forward pass: compute predictions
7    z = X_batch @ self.coefficients_
8    predictions = sigmoid(z)
9    
10    # Compute gradient
11    gradient = (1 / len(batch_indices)) * X_batch.T @ (predictions - y_batch)
12    
13    # Update weights
14    self.coefficients_ -= self.learning_rate * gradient

Once we've trained the model, making predictions is a two-step process:

def predict_proba(self, X):
    z = X @ self.coefficients_
    return sigmoid(z)

def predict(self, X, threshold=0.5):
    probs = self.predict_proba(X)
    return (probs >= threshold).astype(int)

The default threshold is 0.5 (predict class 1 if probability ≥ 0.5), but this is tunable. If false positives are more costly than false negatives, you might use a higher threshold like 0.7. If you want to catch all possible positive cases, use a lower threshold like 0.3.

For regression, we used MSE and R². For classification, we need different metrics. Accuracy is the obvious one, but it's not always enough.

1def precision_score(y_true, y_pred):
2    """Of predicted positives, how many were actually positive?"""
3    tp = np.sum((y_true == 1) & (y_pred == 1))
4    fp = np.sum((y_true == 0) & (y_pred == 1))
5    return tp / (tp + fp) if (tp + fp) > 0 else 0.0
6
7def recall_score(y_true, y_pred):
8    """Of actual positives, how many did we catch?"""
9    tp = np.sum((y_true == 1) & (y_pred == 1))
10    fn = np.sum((y_true == 1) & (y_pred == 0))
11    return tp / (tp + fn) if (tp + fn) > 0 else 0.0
12
13def f1_score(y_true, y_pred):
14    """Harmonic mean of precision and recall."""
15    p = precision_score(y_true, y_pred)
16    r = recall_score(y_true, y_pred)
17    return 2 * (p * r) / (p + r) if (p + r) > 0 else 0.0

Though we've built logistic regression from scratch, here are various ways to extend and optimize this algorithm: