1

I am trying to manually predict a logistic regression model using the coefficient and intercept outputs from a scikit-learn model. However, I can't match up my probability predictions with the predict_proba method from the classifier.

I have tried:

from sklearn.datasets import load_iris
from scipy.special import expit
import numpy as np

X, y = load_iris(return_X_y=True)
clf = LogisticRegression(random_state=0).fit(X, y)
# use sklearn's predict_proba function
sk_probas = clf.predict_proba(X[:1, :])
# and attempting manually (using scipy's inverse logit)
manual_probas = expit(np.dot(X[:1], clf.coef_.T)+clf.intercept_)
# with a completely manual inverse logit
full_manual_probas = 1/(1+np.exp(-(np.dot(iris_test, iris_coef.T)+clf.intercept_)))

outputs:

>>> sk_probas
array([[9.81815067e-01, 1.81849190e-02, 1.44120963e-08]])
>>> manual_probas
array([[9.99352591e-01, 9.66205386e-01, 2.26583306e-05]])
>>> full_manual_probas
array([[9.99352591e-01, 9.66205386e-01, 2.26583306e-05]])

I do seem to get the classes to match (using np.argmax), but the probabilities are different. What am I missing?

I've looked at this and this but haven't managed to figure it out yet.

s_pike
  • 1,710
  • 1
  • 10
  • 22
  • What result do you get with `manual_probas = expit(np.dot(X[:1] - np.max(X[:1], axis=1), clf.coef_.T) + clf.intercept_)`? – Sanjar Adilov May 04 '22 at 17:42
  • The outputs to `manual_probas` is shown under outputs in the question. – s_pike May 04 '22 at 18:51
  • In my comment, I **subtract** the maximum of logits from every logit: `X[:1] - np.max(X[:1], axis=1)`. Not the same as you do. I expect this modification to return the values equal to `sk_probas`. Most softmax implementations use tricks alike to prevent overflow. – Sanjar Adilov May 05 '22 at 04:35
  • 1
    Sorry, I missed you'd changed the definition. Your formulation gives something different again. The accepted answer works, and explains what I was missing. – s_pike May 05 '22 at 08:12

1 Answers1

2

The documentation states that

For a multi_class problem, if multi_class is set to be “multinomial” the softmax function is used to find the predicted probability of each class

That is, in order to get the same values as sklearn you have to normalize using softmax, like this:

from sklearn.datasets import load_iris
from sklearn.linear_model import LogisticRegression
import numpy as np

X, y = load_iris(return_X_y=True)
clf = LogisticRegression(random_state=0, max_iter=1000).fit(X, y)
decision = np.dot(X[:1], clf.coef_.T)+clf.intercept_
print(clf.predict_proba(X[:1]))
print(np.exp(decision) / np.exp(decision).sum())

To use sigmoids instead you can do it like this:

from sklearn.datasets import load_iris
from sklearn.linear_model import LogisticRegression
import numpy as np

X, y = load_iris(return_X_y=True)
clf = LogisticRegression(random_state=0, max_iter=1000, multi_class='ovr').fit(X, y) # Notice the extra argument
full_manual_probas = 1/(1+np.exp(-(np.dot(X[:1], clf.coef_.T)+clf.intercept_)))
print(clf.predict_proba(X[:1]))
print(full_manual_probas / full_manual_probas.sum())
Naphat Amundsen
  • 1,519
  • 1
  • 6
  • 17