Building Logistic Regression in NumPy

In the previous two parts, we separated logistic regression into its main pieces.

The model calculates a weighted score, passes it through the sigmoid function, measures the error with log loss, and updates its parameters using gradient descent.

Now we can place those pieces inside a training loop and build the model ourselves.

We will use scikit-learn to load and prepare the data, but the logistic regression model itself will be written with NumPy.

The dataset

We will use the Iris dataset included in scikit-learn.

The full dataset contains measurements from three iris flower species. Since logistic regression in this example is a binary classifier, we will keep only two:

0 = setosa
1 = versicolor

Each flower has four numerical features:

• sepal length
• sepal width
• petal length
• petal width

import numpy as np
import matplotlib.pyplot as plt
from sklearn.datasets import load_iris
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import StandardScaler
from sklearn.metrics import accuracy_score, log_loss
data = load_iris()
# Keep only setosa and versicolor
binary_mask = data.target < 2
X = data.data[binary_mask]
y = data.target[binary_mask]
print(X.shape)
print(y.shape)

Output:

(100, 4)
(100,)

Each row represents one flower. Each column represents one measurement.

Since the model receives four features, it will learn four weights and one bias.

Preparing the data

We divide the data into training and test sets.

The training set is used to learn the weights and bias. The test set remains separate so that we can later evaluate the model on flowers it has not seen before.

X_train, X_test, y_train, y_test = train_test_split(
    X,
    y,
    test_size=0.20,
    random_state=42,
    stratify=y,
)

The stratify argument keeps the class proportions similar in both sets.

Why scale the features?

The four measurements are not on exactly the same scale.

Gradient descent usually behaves more predictably when the features are numerically comparable. We therefore standardize them so that each training feature has a mean close to 0 and a standard deviation close to 1.

scaler = StandardScaler()
X_train = scaler.fit_transform(X_train)
X_test = scaler.transform(X_test)

The scaler is fitted only on the training data.

We then transform the test set using the means and standard deviations learned from the training set. Fitting the scaler on the complete dataset would allow information from the test set to influence the training process.

This is called data leakage.

The full training process

Before looking at the NumPy code, the process can be written as pseudocode:

initialize the weights and bias
repeat for a fixed number of epochs:
    calculate weighted scores
    apply the sigmoid function
    calculate log loss
    calculate the gradients
    update the weights and bias

Each epoch repeats the same cycle with slightly updated parameters.

Building the model

We first define the sigmoid function:

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

Then we initialize the parameters:

weights = np.zeros(X_train.shape[1])
bias = 0.0
learning_rate = 0.1
epochs = 1000
losses = []

There is one weight for each input feature.

At the beginning, all weights and the bias are zero. This means every weighted score is initially 0.

Since:

sigmoid(0) = 0.5

the model initially assigns every flower a probability of 0.5.

It has not learned anything yet.

The training loop

The full learning process fits inside one loop:

for epoch in range(epochs):
    # Weighted scores
    scores = X_train @ weights + bias
    # Probabilities
    probabilities = sigmoid(scores)
    # Log loss
    epsilon = 1e-12
    safe_probabilities = np.clip(
        probabilities,
        epsilon,
        1 - epsilon,
    )
    loss = -np.mean(
        y_train * np.log(safe_probabilities)
        + (1 - y_train)
        * np.log(1 - safe_probabilities)
    )
    losses.append(loss)
    # Gradients
    errors = probabilities - y_train
    dw = X_train.T @ errors / len(y_train)
    db = np.mean(errors)
    # Parameter updates
    weights -= learning_rate * dw
    bias -= learning_rate * db

This loop contains the entire training process.

What happens inside one epoch?

The first line calculates the weighted scores:

scores = X_train @ weights + bias

This is the matrix version of:$$z=w_1x_1+w_2x_2+w_3x_3+w_4x_4+b$$

The model produces one score for each flower.

The sigmoid function then converts these scores into probabilities:

probabilities = sigmoid(scores)

Each probability represents the model’s estimate that the flower belongs to class 1, which in this example is versicolor.

Next, log loss compares these probabilities with the actual classes.

loss = -np.mean(
    y_train * np.log(safe_probabilities)
    + (1 - y_train)
    * np.log(1 - safe_probabilities)
)

The clipping step prevents the code from calculating log(0), which is undefined.

The model then calculates its prediction errors:

errors = probabilities - y_train

These errors are used to calculate the gradients:

dw = X_train.T @ errors / len(y_train)
db = np.mean(errors)

Finally, the parameters are updated:

weights -= learning_rate * dw
bias -= learning_rate * db

After the update, the next epoch begins with slightly different weights and bias.

Watching the loss decrease

We saved the loss after every epoch.

plt.plot(losses)
plt.xlabel("Epoch")
plt.ylabel("Log loss")
plt.title("Training loss")
plt.show()

At the beginning, the loss should be close to:$$-\log(0.5)\approx0.693$$

This is expected because the initial probabilities are all 0.5.

The loss usually drops quickly during the first epochs. Later, the improvements become smaller as the model approaches a better set of parameters.

A decreasing training loss shows that the model is learning from the training data.

Making predictions

After training, we use the final weights and bias to calculate probabilities for the test set.

test_scores = X_test @ weights + bias
test_probabilities = sigmoid(test_scores)

These outputs are still continuous values between 0 and 1.

To produce class labels, we apply a threshold:

test_predictions = (
    test_probabilities >= 0.5
).astype(int)

A probability of at least 0.5 becomes class 1. A lower probability becomes class 0.

Evaluating the model

We can evaluate the predictions using both accuracy and log loss.

accuracy = accuracy_score(
    y_test,
    test_predictions,
)
test_loss = log_loss(
    y_test,
    test_probabilities,
)
print(f"Accuracy: {accuracy:.3f}")
print(f"Log loss: {test_loss:.3f}")

Accuracy tells us how many final class labels were correct.

Log loss also considers the probabilities behind those decisions.

Two models may classify the same number of flowers correctly, but one may assign much more probability to the correct classes. That model would usually have a lower log loss.

The Iris classes used here are relatively easy to separate, so a high accuracy would not be surprising. The purpose of this example is not to build a competitive classifier. It is to make the training process visible.

The complete model in a few lines

The main calculations can be reduced to this:

scores = X_train @ weights + bias
probabilities = sigmoid(scores)
errors = probabilities - y_train
dw = X_train.T @ errors / len(y_train)
db = np.mean(errors)
weights -= learning_rate * dw
bias -= learning_rate * db

Everything else mainly prepares the data, repeats these steps, tracks the loss, and evaluates the result.

This is the part that is usually hidden behind a method such as:

model.fit(X_train, y_train)

Logistic regression as a single neuron

Our model now has a familiar structure:$$X \rightarrow Xw+b \rightarrow \sigma(z) \rightarrow \hat{y}$$

It receives input features, combines them using weights and a bias, applies an activation function, and produces an output.

That is also the basic structure of an artificial neuron.

Logistic regression can therefore be viewed as a single neuron with a sigmoid activation and no hidden layer.

A neural network extends the same idea by adding more neurons and arranging them in layers.

The model becomes larger, but the learning process remains familiar:

1. Produce an output.
2. Measure the error.
3. Calculate gradients.
4. Update the parameters.
5. Repeat.

Logistic regression gives us the smallest complete version of this process.