import scipy as sp 
import numpy as np 
import pandas as pd
from sklearn.metrics import r2_score, precision_score, recall_score, log_loss
from sklearn.linear_model import LinearRegression, Ridge
import sklearn

Some data

df = pd.read_csv("https://download.mlcc.google.com/mledu-datasets/california_housing_train.csv")
df = df[['housing_median_age', 'total_rooms', 'total_bedrooms', 'population', 'households', 'median_income', 'median_house_value']]

Linear Regression

scaled_df = (df - df.min()) / (df.max() - df.min())
X = scaled_df[['housing_median_age', 'total_rooms', 'total_bedrooms', 'population', 'households', 'median_income']].values
y = scaled_df['median_house_value'].values

X_with_intercept = np.hstack((np.ones((len(X), 1)),X))
B = np.linalg.inv(X_with_intercept.T @ X_with_intercept) @ (X_with_intercept.T @ y.reshape(-1, 1))

print("Manual weights: ", B.reshape(-1))
print("Manual score: ", r2_score(y, (X_with_intercept @ B).reshape(-1)))

Manual weights:  [-0.07556544  0.19769139 -1.56087573  1.32234017 -2.57610401  1.59516284
  1.43606576]
Manual score:  0.5713482748283873

from sklearn.metrics import r2_score

RSS = (((X_with_intercept @ B).reshape(-1) - y)**2).sum() # Squared distance from our new regression line
TSS = ((y.mean() - y)**2).sum()                           # Squared distance from the mean
r2 = 1 - RSS / TSS                                        # How much distance did we gained ? Did we reduce the errors ? Are we closer to the actual point values ?
r2, r2_score(y, (X_with_intercept @ B).reshape(-1))

(np.float64(0.5713482748283873), 0.5713482748283873)

Let’s compare those results with sklearn linear regression


lr = LinearRegression().fit(X, y)

print("")
print("Sklearn weights: ", [lr.intercept_] + lr.coef_.tolist() )
print("Sklearn score: ", r2_score(y, lr.predict(X)))


Sklearn weights:  [np.float64(-0.07556543642855307), 0.19769138728528812, -1.5608757342094828, 1.322340171543368, -2.576104006535326, 1.5951628411047347, 1.4360657609756633]
Sklearn score:  0.5713482748283873

Linear regression with regularization (Ridge regression)

Regularization is the action of adding to the loss, a term that contains the weight values. That way these terms are forced to stay small. This helps avoiding overfitting.

Let’s look at the ordinary least sqaure loss and then add the square of each weight to build the regularized loss. Adding the square of each weight means we buil the Ridge regression loss. If we add the absolute value of each weight we build the Lasso regression loss.

e = X_with_intercept @ B - y.reshape(-1, 1)
loss = (e.T @ e).item()
regularized_loss = loss + 0.3 * (B.T @ B).item()
loss, regularized_loss

(416.71131319597765, 421.3531681415839)

The way adding this loss impacts the formula is the following

scaled_df = (df - df.min()) / (df.max() - df.min())
X = scaled_df[['housing_median_age', 'total_rooms', 'total_bedrooms', 'population', 'households', 'median_income']].values
y = scaled_df['median_house_value'].values


X_with_intercept = np.hstack((np.ones((len(X), 1)),X))

I = np.identity(X_with_intercept.shape[1])
I[0,0] = 0
B = np.linalg.inv(X_with_intercept.T @ X_with_intercept + 0.3 * I) @ (X_with_intercept.T @ y.reshape(-1, 1))

print("Manual weights: ", B.reshape(-1))
print("Manual score: ", r2_score(y, (X_with_intercept @ B).reshape(-1)))

Manual weights:  [-0.07457501  0.19926227 -1.4614579   1.30386275 -2.31228351  1.40463349
  1.42708759]
Manual score:  0.5710213053584059


lr = Ridge(alpha=0.3).fit(X, y)

print("")
print("Sklearn weights: ", [lr.intercept_] + lr.coef_.tolist() )
print("Sklearn score: ", r2_score(y, lr.predict(X)))


Sklearn weights:  [np.float64(-0.07457500943073725), 0.19926227134208804, -1.4614578956147584, 1.3038627486537557, -2.312283513756178, 1.4046334910837726, 1.4270875901070879]
Sklearn score:  0.5710213053584055

Logistic Regression

For the logistic regression, we transform the X values in the same way but we add a sigmoid transform at the end in order to map to values between 0 and 1.

We can not use the normal form anymore for computing the weights. We have to resort to other techniques like gradient descent.

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

def log_likelihood(y_hat, y_true):
  # Being far away from the correct class is penalized heavily. 
  return - np.mean( y_true * np.log(y_hat) + (1-y_true) * np.log(1-y_hat) )

def gradient_sigmoid(x):
  sigmoid(X) * (1 - sigmoid(X))


def gradients(X, y, y_hat):
    # Loss = y * log(h) + (1 - y) * log(1-h)
    # where h = sigmoid(z)
    # and z = Xt @ B

    # deriv_loss_to_h = y / h - (1-y) / (1-h) = (y - h) / (h * (1 - h))
    # deriv_h_to_z = sigmoid(h) * (1 - sigmoid(h))
    # deriv_z_to_b = Xt
    # Though chain rule, final derivative 
    # final_derivative = deriv_loss_to_h * deriv_h_to_z * deriv_z_to_b = x * (y - h) = x * (y - y_hat) 
    dw = (1/len(X)) * (X.T @ (y_hat - y))
    return dw

df['median_house_value_cat'] = (df['median_house_value'] > 150_000).astype(int)

scaled_df = (df - df.min()) / (df.max() - df.min())
X = scaled_df[['housing_median_age', 'total_rooms', 'total_bedrooms', 'population', 'households', 'median_income']].values
y = df['median_house_value_cat'].values 

X_with_intercept = np.hstack((np.ones((len(X), 1)),X))

B = np.random.normal(0, 0.1 ,(7, 1))

for i in range(50_000):
  y_hat = sigmoid(X_with_intercept @ B).reshape(-1)
  if i % 5000 == 0 or i ==0: 
    print("loss: ", log_likelihood(y_hat, y))
  deltas = gradients(X_with_intercept, y, y_hat)
  B -= 0.3 * deltas.reshape(-1, 1)


lr = sklearn.linear_model.LogisticRegression().fit(X, y)

loss:  0.7175633823124513
loss:  0.46247742549766496
loss:  0.45400111162106355
loss:  0.45092915252987215
loss:  0.44860694198860407
loss:  0.44654040375802884
loss:  0.44465189016022644
loss:  0.4429190526397757
loss:  0.44132801358504486
loss:  0.43986676128612356

print("Manual weights: ", B.reshape(-1))
print("Manual score: ", 
        precision_score(y, (sigmoid(X_with_intercept @ B).reshape(-1) > 0.5).astype(int) ),
        recall_score(y, (sigmoid(X_with_intercept @ B).reshape(-1) > 0.5).astype(int) ),
      )
print()
print("Sklearn log loss: ", log_loss(y, (sigmoid(X_with_intercept @ B).reshape(-1))))
print("Sklearn weights: ", lr.intercept_.tolist() + lr.coef_.reshape(-1).tolist())
print("Sklearn score", 
      precision_score(y, lr.predict(X)),
      recall_score(y, lr.predict(X))
      )

Manual weights:  [ -4.74777447   2.09463637 -11.67292203   7.04060378  -3.12359728
   8.17776188  18.8596743 ]
Manual score:  0.8215827338129497 0.8513652036156929

Sklearn log loss:  0.43852424076056
Sklearn weights:  [-4.382530400697039, 1.9562829511733213, -10.781475873180312, 6.397256726894048, -2.6864271237826665, 7.662584658270458, 17.40432279277666]
Sklearn score 0.8187667560321715 0.8537880905786972

The weights are not exactly the same but the performances are very similar. This is due to the randomness aspect of training through gradient descent.