Logistic Regression From Scratch

Deep dive on logistic regression.

What is Logistic Regression?

• Logistic Regression is a supervised Machine Learning Algorithm used for classification tasks.
• It is suitable for binary (yes/no) and multi-class classification problems, but often you see it on binary.
• In this guide I want to show you the underlying math behind the logistic regression and the implementation of this algorithm using just Numerical Python.

import numpy as np

X = np.array([35, 42, 55, 61, 68, 72, 78, 85, 91, 96])
y_true = np.array([0, 0, 0, 0, 1, 1, 1, 1, 1, 1])

Data

Our data is simple, in one side we have the relative humidity in % as the independent variable (\(X\)), where in this dataset the separation of a rainy day or not, is between 61% and 68%

Important

For the sake of simplicity and understanding, we are not normalizing the data, but it should be normalized before working with logistic regression, later on the code, we normalize the data.

import numpy as np

from plot_helper import plot_class

plot_class(X=X, y_true=y_true, line=59, title='Naive rainy data points').show()

in our data we have a solid separation between a rainy/normal day, pretty much we could draw a vertical line in the middle between 61 and 68, and we have a perfect classification model, if the relative humidity is greater than 64.5% it’s going to rain. let’s have a more real world example where our data is not that clean.

more_rain_points_X = [55, 58, 63, 66, 69, 80, 85, 90, 92, 95, 96, 97, 99]
more_rain_points_y = [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1]
more_normal_points_X = [10, 13, 15, 18, 20, 30, 40, 50, 52, 54, 55, 60, 61]
more_normal_points_y = [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]

X = np.concatenate((more_rain_points_X, X, more_normal_points_X))
y_true = np.concatenate((more_rain_points_y, y_true, more_normal_points_y))

Sigmoid

It’s an \(S\) shaped curved function that can take any value and squash it into a range between \(0\) and \(1\)

\[\sigma(z) = \frac{1}{1 + e^{-z}}\]

can take any real value, from \(\infty\) to \(-\infty\), and the higher the \(z\) we get a value closer to \(1\), and the lower the \(z\) we get a value closer to \(0\), if \(z = 0\) we get exactly \(0.5\)

NoteIntuition

The sigmoid helps us to translate the raw into a calibrated and understandable probability

WarningWhy not a Step Function?

A step function has a derivative of zero everywhere, making gradient descent completely blind no matter how wrongs the prediction, the gradient says “change nothing.” The sigmoid is smooth and differentiable everywhere, providing the learning signal gradient descent needs at every point.

from plot_helper import plot_sigmoid_interactive

plot_sigmoid_interactive(X, y_true).show()

Loss Function

Binary Cross Entropy (Log Loss) the penalty for overconfidence, punishes the model more severely when it is confidently wrong, if the actual label is \(1\) and your model predicts \(0.01\), the loss becomes massive (approaching \(\infty\)).

WarningWhy not MSE?

• if MSE is used will result in a complex surface with many hills and valleys (local minima), where the algorithm might get stuck in a local minimum rather than finding the best possible parameters.
• We get learning stagnation if the model makes a confident but wrong prediction, where the model stops updating its weights even though it is still very inacurate.

\[L(y, \hat{y}) = -[y \log(\hat{y}) + (1 - y) \log(1 - \hat{y})]\] where \(y\) is the actual (true) label and \(\hat{y}\) is the predicted probability.

def sigmoid(z):
    return 1 / (1 + np.exp(-z))


def predict(X, m, b):
    z = (m * X) + b
    return sigmoid(z)


def loss(y_true, y_pred):
    epsilon = 1e-15  # prevents log(0) errors
    y_pred = np.clip(y_pred, epsilon, 1 - epsilon)
    return -np.mean(y_true * np.log(y_pred) + (1 - y_true) * np.log(1 - y_pred))


def update(X, y_true, y_pred, m, b, learning_rate):
    dm = np.mean(X * (y_pred - y_true))
    db = np.mean(y_pred - y_true)
    new_m = m - learning_rate * dm
    new_b = b - learning_rate * db
    return new_m, new_b

The update function uses the gradient descent, the derivative of the Binary Cross Entropy loss with respect to our parameters \(m\) and \(b\) turns out to be the same as the linear regression

Gradient with respect to \(m\):

\((\frac{\partial L}{\partial m}): \text{mean}((\hat{y} - y) \times X)\)

Gradient with respect to \(b\):

\((\frac{\partial L}{\partial b}): \text{mean}(\hat{y} - y)\)

Training

For this example, to make the math more doable and readable, we are going to use 3 data points, and we are going to do the first three epochs.

rel hum %	Rain
50	0
60	1
63	1

We initialize \(w = 0\) (weight) and \(b = 0\) (bias). For logistic regression, zero initialization is perfectly fine (unlike deep neural network); there’s only one layer and no symmetry problem where neurons would all learn the same thing. The model will break symmetry through the data itself since different features pull the weights in different directions from epoch 1.

First epoch

Forward Pass

\(x_1 = 50 : z_1 = (0.1 \times 50) + 0 = 5 \to \hat{y_1} = \sigma(5) \approx 0.993\)

\(x_2 = 60 : z_2 = (0.1 \times 60) + 0 = 6 \to \hat{y_2} = \sigma(6) \approx 0.997\)

\(x_3 = 63 : z_3 = (0.1 \times 63) + 0 = 6.3 \to \hat{y_3} = \sigma(6.3) \approx 0.998\)

Backward Pass

We calculate the error for every data point = \((\hat{y} - y) \times x\)

\(g_{w1} = (0.993 - 0) \times 50 = 49.65\)

\(g_{w2} = (0.997 - 1) \times 60 = -0.18\)

\(g_{w3} = (0.998 - 1) \times 63 = -0.126\)

Mean of \(w\):

\(\frac{49.65 - 0.18 - 0.126}{3} = 16.448\)

Mean of \(b\):

\(\frac{(0.993 - 0) + (0.997 - 1) + (0.998 - 1)}{3} = 0.329\)

Update

We use our learning rate defined at 0.01, with our first (often random) values.

\(w = 0.1 - (0.01 \times 16.448) = -0.064\)

\(b = 0 - (0.01 \times 0.329) = -0.003\)

NoteView Detailed Math for Epochs 2 & 3

Second epoch

Forward Pass

\(x_1 = 50 : z_1 = (-0.064 \times 50) - 0.003 = -3.203 \to \hat{y_1} = \sigma(-3.203) \approx 0.040\)

\(x_2 = 60 : z_2 = (-0.064 \times 60) - 0.003 = -3.843 \to \hat{y_2} = \sigma(-3.843) \approx 0.021\)

\(x_3 = 63 : z_3 = (-0.064 \times 63) - 0.003 = -4.035 \to \hat{y_3} = \sigma(-4.035) \approx 0.017\)

Backward Pass

We calculate the error for every data point = \((\hat{y} - y) \times x\)

\(g_{w1} = (0.040 - 0) \times 50 = 1.95\)

\(g_{w2} = (0.021 - 1) \times 60 = -58.74\)

\(g_{w3} = (0.017 - 1) \times 63 = -61.92\)

Mean of \(w\):

\(\frac{1.95 - 58.74 - 61.92}{3} = -39.57\)

Mean of \(b\):

\(\frac{(0.040 - 0) + (0.021 - 1) + (0.017 - 1)}{3} = -0.64\)

Update

We use our learning rate defined at 0.01, with our first (often random) values.

\(w = -0.064 - (0.01 \times -39.57) = 0.332\)

\(b = -0.003 - (0.01 \times -0.64) = 0.003\)

Third epoch

Forward Pass

\(x_1 = 50 : z_1 = (0.332 \times 50) + 0.003 = 16.55 \to \hat{y_1} = \sigma(16.55) \approx 1\)

\(x_2 = 60 : z_2 = (0.332 \times 60) + 0.003 = 19.86 \to \hat{y_2} = \sigma(19.86) \approx 1\)

\(x_3 = 63 : z_3 = (0.332 \times 63) + 0.003 = 20.85 \to \hat{y_3} = \sigma(20.85) \approx 1\)

Backward Pass

We calculate the error for every data point = \((\hat{y} - y) \times x\)

\(g_{w1} = (1 - 0) \times 50 = 50\)

\(g_{w2} = (1 - 1) \times 60 = 0\)

\(g_{w3} = (1 - 1) \times 63 = 0\)

Mean of \(w\):

\(\frac{50 + 0 + 0}{3} = 16.66\)

Mean of \(b\):

\(\frac{(1 - 0) + (1 - 1) + (1 - 1)}{3} = 0.33\)

Update

We use our learning rate defined at 0.01, with our first (often random) values.

\(w = 0.331 - (0.01 \times 16.66) = 0.164\)

\(b = 0.003 - (0.01 \times 0.33) = 0\)

m = 0
b = 0
epochs = 5000
learning_rate = 0.01
tolerance = 1e-10
patience = 3
X_norm = (X - X.mean()) / X.std()
verbose = False

loss_history = []
slope_history = []

best_loss = float("inf")
bad_epochs = 0

for epoch in range(epochs):
    y_pred = predict(X_norm, m, b)
    current_loss = loss(y_true, y_pred)

    loss_history.append(current_loss)
    slope_history.append((m, b))

    if verbose:
        print(f"epoch {epoch}. current loss: {current_loss:.6f}")

    improvement = best_loss - current_loss
    if improvement < tolerance:
        bad_epochs += 1
        if bad_epochs >= patience:
            print(f"Converged at epoch {epoch}. Best loss: {best_loss:.6f}")
            break
    else:
        bad_epochs = 0
        best_loss = current_loss

    m, b = update(X_norm, y_true, y_pred, m, b, learning_rate)

final_loss = loss_history[-1]
final_m, final_b = slope_history[-1]
print(f"Final Loss: {final_loss:.4f}")
print(f"Final Slope: {final_m:.4f}, Final Intercept: {final_b:.4f}")

Final Loss: 0.2203
Final Slope: 3.5544, Final Intercept: 0.0774

from plot_helper import plot_lr_fit

plot_lr_fit(X, y_true, slope_history, loss_history).show()

Evaluation

A logistic regression model outputs a probability, not a class label. To make a final prediction, we apply a decision threshold to any prediction \(\hat{y} \geq 0.5\) is classified as \(1\) (rainy), and anything below as \(0\) (not rainy).

\[\hat{y}_{class} = \begin{cases} 1 & \text{if } \hat{y} \geq 0.5 \\ 0 & \text{if } \hat{y} < 0.5 \end{cases}\]

NoteThreshold is a choice, not a rule

0.5 is the default, but in practice you tune it. In a medical diagnosis model you might lower it to 0.3 so the model catches more positive cases (fewer false negatives), accepting more false alarms in exchange.

# Apply threshold to get class predictions
threshold = 0.5
y_pred_final = predict(X_norm, final_m, final_b)
y_pred_class = (y_pred_final >= threshold).astype(int)

# Accuracy
accuracy = np.mean(y_pred_class == y_true)
print(f"Accuracy: {accuracy * 100:.1f}%")

# Confusion Matrix (manual, no sklearn)
TP = np.sum((y_pred_class == 1) & (y_true == 1))
TN = np.sum((y_pred_class == 0) & (y_true == 0))
FP = np.sum((y_pred_class == 1) & (y_true == 0))
FN = np.sum((y_pred_class == 0) & (y_true == 1))

print("\nConfusion Matrix:")
print(f"                Predicted 0   Predicted 1")
print(f"  Actual 0         {TN:>3}           {FP:>3}")
print(f"  Actual 1         {FN:>3}           {TP:>3}")
print(f"\nTrue Positives:  {TP} (rainy, correctly predicted rainy)")
print(f"True Negatives:  {TN} (normal, correctly predicted normal)")
print(f"False Positives: {FP} (normal, wrongly predicted rainy)")
print(f"False Negatives: {FN} (rainy, wrongly predicted normal)")

def predict_new(x_raw, m, b, mu, sigma):
    z = m * ((x_raw - mu) / sigma) + b
    return sigmoid(z)


mu, sigma = X.mean(), X.std()
for humidity in [30, 55, 65, 73, 90]:
    p = predict_new(humidity, final_m, final_b, mu, sigma)
    label = "rainy" if p >= 0.5 else "normal"
    print(f"humidity = {humidity}% -> P(rain) = {p:.3f} ({label})")

Accuracy: 88.9%

Confusion Matrix:
                Predicted 0   Predicted 1
  Actual 0          15             2
  Actual 1           2            17

True Positives:  17 (rainy, correctly predicted rainy)
True Negatives:  15 (normal, correctly predicted normal)
False Positives: 2 (normal, wrongly predicted rainy)
False Negatives: 2 (rainy, wrongly predicted normal)
humidity = 30% -> P(rain) = 0.015 (normal)
humidity = 55% -> P(rain) = 0.315 (normal)
humidity = 65% -> P(rain) = 0.642 (rainy)
humidity = 73% -> P(rain) = 0.842 (rainy)
humidity = 90% -> P(rain) = 0.982 (rainy)

Takeaways

• Logistic Regression is about modeling the probability of a binary outcome, instead of predicting continuous values, it assigns probabilities between 0 and 1.

• The Sigmoid (or logistic) function plays a central role by transforming raw scores into probabilities. This makes the output interpretable as the likelihood of a class.

• The Loss Function is typically cross entropy (log loss). This function measures the difference between the predicted probabilities and the actual binary outcomes, guiding how well the model is performing.

• Gradient Descent is used iteratively to minimize the loss, and automation is key in practice; techniques like early stopping, converge thresholds and adaptative optimizers help to streamline the training process, while avoiding overfitting and ensuring efficient convergence.

• However, logistic regression isn’t limited only to binary problems; there are extensions that allow it to handle multi class classification.