<jupyter_start><jupyter_text>Logistic Regression on ProbabilitiesTypical uses of logistic regression have the dependent variable's values encoded as 1 or 0. However, sometimes, the dependent variable may come as probabilities. Scikit-Learn will not work when the dependent variable is not encoded as binary values. What can we do? We may choose an arbitrary threshold (e.g. 0.5) and set all values greater than or equal to such threshold as 1 and the others as 0. But the threshold is arbitrary. Of course, we can do other types of analyses to see which threshold is the best one, but nevertheless, our learned model (the learned coefficients) will change.You can actually perform logistic regression on a continuous dependent variable, however, you will have to code up the algorithm to learn the coefficients yourself. You may be wondering, when would I want to perform logistic regression on a continuous output? Take the [Bradley-Terry Model](https://en.wikipedia.org/wiki/Bradley%E2%80%93Terry_model), which estimates paired ranking preferences. In these models, the design matrix `X` is regressed against continuous values to estimate rankings.   DataThis data comes from the [<NAME>](https://www.amazon.com/Introduction-Categorical-Analysis-Probability-Statistics/dp/1119405262/) book `Introduction to Categorical Analysis`. The data is composed of head-to-head matches between 5 professional men's tennis players. For each row, $n_{ij}$ is the number of wins and $n_{ji}$ is the number of losses for the player marked by 1 against the player marked by -1. The $y$ value is derived as $y = \dfrac{n_{ij}}{n_{ij} + n_{ji}}$. The design matrix `X` is just the data without $n_{ij}$, $n_{ji}$ and $y$. The dependent variable $y$ is continuous and each value is a probability. If we perform logistic regression $y \sim X$, the coefficients give us the rankings of each player.<jupyter_code>import autograd.numpy as np
import pandas as pd
from autograd import grad

df = pd.DataFrame([
    [1, -1, 0, 0, 0, 9, 6],
    [1, 0, -1, 0, 0, 14, 3],
    [1, 0, 0, -1, 0, 9, 2],
    [1, 0, 0, 0, -1, 4, 3],
    [0, 1, -1, 0, 0, 5, 0],
    [0, 1, 0, -1, 0, 5, 1],
    [0, 1, 0, 0, -1, 7, 2],
    [0, 0, 1, -1, 0, 2, 4],
    [0, 0, 1, 0, -1, 2, 2],
    [0, 0, 0, 1, -1, 4, 3]
], columns=['Djokovic', 'Federer', 'Murray', 'Nadal', 'Wawrinka', 'n_ij', 'n_ji'])
df['y'] = df.n_ij / (df.n_ij + df.n_ji)

df
X = df[[c for c in df.columns if c not in ['n_ij', 'n_ji', 'y']]].values
y = df.y.values<jupyter_output><empty_output><jupyter_text>Scikit failureWe can see that Scikit-Learn cannot perform logistic regression when $y$ is continuous.<jupyter_code>from sklearn.linear_model import LogisticRegression

m = LogisticRegression()
m.fit(X, y)<jupyter_output><empty_output><jupyter_text>Log lossWe can use gradient descent to learn the model (estimate the coefficients) on our own. The loss function we are trying to optimize is given as follows.- $\hat{y} = Xw$- $\dfrac{\sum -y \hat{y} + \log(1 + \exp(\hat{y}))}{n}$<jupyter_code>from autograd.numpy import exp, log, sqrt

def loss(w, X, y):
    n = float(len(X))
    y_pred = np.dot(X, w)
    return np.sum(-(y_pred * y) + log(1.0 + exp(y_pred))) / n

loss_grad = grad(loss)
w = np.array([0.0 for _ in range(X.shape[1])])
alpha=0.05

for i in range(10_000):
    loss = loss_grad(w, X, y)
    w = w - (loss * alpha)<jupyter_output><empty_output><jupyter_text>Here are the weights/coefficients learned. Notice there is no intercept estimation in the Bradley-Terry model. If we needed an intercept, we would add a column vector of 1 to `X`. These weights corresponding to rankings of each player (sort them descendingly to see who is the top; it should be Federer, Djokovic, Nadal, Wawrinka and then Murray).<jupyter_code>w<jupyter_output><empty_output><jupyter_text>Do the following to calculate the probability of Federer beating Djokovic.<jupyter_code>np.exp(w[1] - w[0]) / (1 + np.exp(w[1] - w[0]))<jupyter_output><empty_output><jupyter_text>Do the following the calculate the probability of Djokovic beating Federer.<jupyter_code>np.exp(w[0] - w[1]) / (1 + np.exp(w[0] - w[1]))<jupyter_output><empty_output><jupyter_text>Comparing the true $y$ values to the predicted ones $\hat{y}$. Notice that Federer has beaten Murray 5-0 or 100%,  but the model estimates 87% chance of win for Federer.<jupyter_code>import pandas as pd

pd.DataFrame({
    'y_true': y,
    'y_pred': 1 / (1 + np.exp(-X.dot(w)))
})<jupyter_output><empty_output><jupyter_text>MSE lossHere's another interesting thing to do. Let's define the loss function using Mean Squared Error (MSE) and learn the model using this loss function.- $\hat{y} = \dfrac{1}{1 + \exp(Xw)}$- $\dfrac{\sum (\hat{y} - y)^2}{n}$<jupyter_code>def loss(w, X, y):
    n = float(len(X))
    y_pred = np.dot(X, w)
    y_pred = 1 / (1 + np.exp(-y_pred))
    
    loss = ((y_pred - y) ** 2.0)
    return loss.mean(axis=None)

loss_grad = grad(loss)
w = np.array([0.0 for _ in range(X.shape[1])])
alpha=0.05

for i in range(20_000):
    loss = loss_grad(w, X, y)
    w = w - (loss * alpha)<jupyter_output><empty_output><jupyter_text>Here are the weights.<jupyter_code>w<jupyter_output><empty_output><jupyter_text>What's the chance of Federer beating Djokovic?<jupyter_code>np.exp(w[1] - w[0]) / (1 + np.exp(w[1] - w[0]))<jupyter_output><empty_output><jupyter_text>What's the chance of Djokovic beating Federer?<jupyter_code>np.exp(w[0] - w[1]) / (1 + np.exp(w[0] - w[1]))<jupyter_output><empty_output><jupyter_text>Compare the $y$ against $\hat{y}$.<jupyter_code>pd.DataFrame({
    'y_true': y,
    'y_pred': 1 / (1 + np.exp(-X.dot(w)))
})<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1359773371104816, 'fraction_numerical': 0.028151558073654392, 'mean_word_length': 4.416107382550336, 'pattern_counts': {'":': 0, '<': 50, '<?xml version=': 0, '>': 50, 'https://': 2, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '12', 'dir': 'jupyter-structured-clean-dedup', 'id': '467708', 'license': 'CC-BY-4.0', 'n_tokens_mistral': 2124, 'n_tokens_neox': 1984, 'n_words': 700, 'path': 'sphinx/datascience/source/logreg-y-probability.ipynb', 'repo_name': 'oneoffcoder/books', 'size': '5652'}