Examples
import numpy as np
import pandas as pd
rng = np.random.default_rng(24)
def logistic(x):
return 1 / (1 + np.exp(-x))
def data_generator(x, rng=rng):
x = np.asarray(x, dtype=float)
y = np.where(logistic(2 * x - 1) > rng.random(len(x)), "A", "B")
return pd.DataFrame({"x": x, "y": y})
data = data_generator([0., 2., 3.])
data| x | y | |
|---|---|---|
| 0 | 0.0 | B |
| 1 | 2.0 | A |
| 2 | 3.0 | A |
from sklearn.linear_model import LogisticRegression
X = data[["x"]].values
m = LogisticRegression(penalty=None).fit(X, data["y"])
print("classes: ", m.classes_)
print("intercept:", round(-m.intercept_[0], 4))
print("slope: ", round(-m.coef_[0][0], 4))classes_[1], which is B. This page writes everything in terms of \(P(\text{A}|x)\), and flipping the sign of both parameters is what turns one into the other.
classes: ['A' 'B']
intercept: -8.158
slope: 8.5716The model gives a probability for each class.
grid = np.arange(-1, 2.5, 0.5).reshape(-1, 1)
pd.DataFrame({"x": grid[:, 0],
"P(A|x)": m.predict_proba(grid)[:, 0].round(3),
"predicted": m.predict(grid)})| x | P(A|x) | predicted | |
|---|---|---|---|
| 0 | -1.0 | 0.000 | B |
| 1 | -0.5 | 0.000 | B |
| 2 | 0.0 | 0.000 | B |
| 3 | 0.5 | 0.020 | B |
| 4 | 1.0 | 0.602 | A |
| 5 | 1.5 | 0.991 | A |
| 6 | 2.0 | 1.000 | A |
The prediction is the class with the highest probability, which with two classes means a threshold of 0.5. That threshold is a choice, not something the model dictates, and further down this page you can move it and watch the decision change.
The log-likelihood of a parameter pair is the sum, over the training points, of the log probability the model gives to the class that was actually observed. Our three labels are B, A, A, so it is
\[ \ell(\beta_0, \beta_1) = \log P(\text{B}|0) + \log P(\text{A}|2) + \log P(\text{A}|3). \]
from sklearn.metrics import log_loss
def loglik(theta):
b0, b1 = theta
pA = logistic(b0 + b1 * data["x"].values)
observed = np.where(data["y"].values == "A", pA, 1 - pA)
return np.sum(np.log(observed))
theta_hat = [-m.intercept_[0], -m.coef_[0][0]]
by_hand = loglik(theta_hat)
by_sklearn = -log_loss(data["y"], m.predict_proba(X), normalize=False)
print(f"by hand: {by_hand:.6f}")
print(f"sklearn: {by_sklearn:.6f}")LogisticRegression sorts the classes alphabetically and reports the coefficients of the second one, here B. Ours are written for A, and \(s(-z) = 1 - s(z)\), so the two differ by a sign.
by hand: -0.000412
sklearn: -0.000412log_loss with normalize=False is the negative log-likelihood. The two agree.
Look at the size of that fitted slope, and at how close the log-likelihood is to zero. With these three points every label is predicted almost perfectly, because one straight cut at any \(x\) between 0 and 2 separates the B from the two As. Nothing then stops the fit from making the step steeper and the likelihood larger, and what the numbers above report is not a maximum but wherever the optimizer chose to stop. Three points are not enough to pin a curve down.
The hand computation above evaluates \(\ell\) at a single pair of parameters. Fitting is the search for the pair that makes it as large as possible, and that search is easier to see with enough data to have an answer. The panel below draws sixty points from the same generator, \(P(\text{A}|x) = s(2x - 1)\).
Two things are worth doing with it. Set the sliders to the generator, intercept \(-1\) and slope \(2\), and the two curves coincide — but the log-likelihood there, \(-13.01\), is not the best the grid offers. Now raise the slope one notch to \(2.5\) and it improves to \(-12.54\).
That gap is the whole idea. Maximum likelihood aims at the data, not at the generator, and with sixty points the two are close but not equal. The fit is the parameters that best explain the labels we happened to draw, and a different sixty points would move it somewhere else nearby.
Nothing so far used the fact that \(x\) was a single number. With two predictors,
\[ P(Y = \text{A}|x_1, x_2) = s(\beta_0 + \beta_1 x_1 + \beta_2 x_2), \]
the probability is a surface over the input plane, and the threshold is a horizontal plane cutting through it. Where the two meet is where the decision changes.
That meeting is a straight line, and this is the one thing worth taking away from the picture below. The logistic function is monotone, so the probability is the same wherever \(\beta_0 + \beta_1 x_1 + \beta_2 x_2\) is the same, and that is a straight line.
Two things in that last panel are worth separating. The small dots are what the model decides, and they change only when the boundary moves. The large dots are what the generator drew, and near the boundary the same input can come out either way. Even with the parameters that produced the data, some points sit on the wrong side. That is not an error of the fit. It is the noise in the problem, and no classifier can remove it.
The cross-entropy is what we fit. It is not always what we report. Once a prediction is a class rather than a probability, the natural summary is the misclassification rate, the average of the zero-one loss
\[ L(y, \hat y) = \begin{cases} 0 & \hat y = y\\ 1 & \hat y \neq y \end{cases} \]
and its complement, the accuracy.
big = data_generator(rng.normal(0, 2, 500))
m2 = LogisticRegression(penalty=None).fit(big[["x"]].values, big["y"])
pred = m2.predict(big[["x"]].values)
print("accuracy: ", round((pred == big["y"]).mean(), 3))
print("misclassification: ", round((pred != big["y"]).mean(), 3))accuracy: 0.852
misclassification: 0.148