0
$\begingroup$

I tried to do my own linear regression but I have seen that there are differences with LinearRegression from sklearn. Indeed, this is what gives my own linear regression:

>>> beta = estimate_beta(X_train.values, y_train.values)
>>> (X_test*beta).sum(axis = 1)
1142949     777.876173
696990     1010.574350
1437036     852.938428
1404619     439.875784
1151040    1002.271470
              ...     
848794      396.552199
291118     3137.906208
128267      541.926482
459668      752.088114
1120216     424.453391
Length: 482738, dtype: float64

Which is very different from what sklearn returns :

array([[3.24519438],
       [2.12164625],
       [0.98814432],
       ...,
       [1.13965629],
       [2.09961495],
       [1.84152796]])

And, of course, of the true values:

array([[2.],
       [1.],
       [1.],
       ...,
       [1.],
       [4.],
       [2.]])

How come I have results so far from the real values with my own model?

Here is my model:

import numpy as np

def predict(x_i, beta):
    """assumes that the first element of each x_i is 1"""
    return np.dot(x_i, beta)

def error(x_i, y_i, beta):
    return y_i - predict(x_i, beta)

def squared_error(x_i, y_i, beta):
    return error(x_i, y_i, beta) ** 2

def squared_error_gradient(x_i, y_i, beta):
    """The gradient (with respect to beta)
    correspond to the ith squared error term"""
    return [-2 * x_ij * error(x_i, y_i, beta)
           for x_ij in x_i]

def estimate_beta(x, y):
    beta_initial = [random.random() for x in x[0]]
    return minimize_stochastic(squared_error,
                              squared_error_gradient,
                              x, y,
                              beta_initial,
                              0.01)

def in_random_order(data):
    """generator that returns the elements of data in random order"""
    indexes = [i for i, _ in enumerate(data)] # create a list of indexes
    random.shuffle(indexes)                   # suffle them
    for i in indexes:
        yield data[i]

def minimize_stochastic(target_fn, gradient_fn, x, y, theta_0, alpha_0=0.01):
    data = zip(x, y)
    theta = theta_0                           # the initial guess
    alpha = alpha_0                           # initial step size
    min_theta, min_value = None, float("inf") # the minimum so far
    iterations_with_no_improvement = 0
    
    # if we ever go 100 iterations with no improvement, stop
    while iterations_with_no_improvement <100:
        value = sum(target_fn(x_i, y_i, theta) for x_i, y_i in data)
        
        if value < min_value:
            # if we've found a new minimum, remeber it
            # and go back to the original step size
            min_theta, min_value = theta, value
            iterations_with_no_improvement = 0
            alpha = alpha_0
        else:
            # otherwise we're not improving, so try shinking the step_size
            iterations_with_no_improvement += 1
            alpha *= 0.9

        # and take a gradient step for each of the data points
        for x_i, y_i in in_random_order(data):
            gradient_i = gradient_fn(x_i, y_i, theta)
            theta = np.substract(theta, alpha * gradient_i)
            
    return min_theta

I tried with a minimal and reproducible code:

>>>x = [[1, 49, 4, 2500],[1, 12, 2, 125], [1, 35, 2, 3790], [1, 60, 4, 4500], [1, 10, 4, 5000]]
>>>spendings_customer = [2000, 0, 3600, 0, 3500]
>>>import random
>>>random.seed(0)
>>>beta = estimate_beta(x, spendings_customer)

Which returns to me:

[0.8444218515250481,
 0.7579544029403025,
 0.420571580830845,
 0.25891675029296335]

Which is strange ...

With a learning rate from the hessian of the loss

The problem could have been because of my learning rate, which could have been too big. The best learning rate of a cost function is strictly less than 2/λ, where λ is the largest eigenvalue of the Hessian. I want to get this learning rate of my gradient descent algorithm. So I tried to compute it using this answer:

$$H_L(w) = X^T X$$

# Hessian is X.t*X
h = np.dot(X_test.T,X_test)
from numpy import linalg as LA
w, v = LA.eig(np.array(h))
max(w)

which returns:

(74119951381184.84+0j)

So I changed alpha_0 to something less that 2/λ But the predictions still look very different from the actual results:

1142949    1883.752457
696990     1015.531962
1437036    3342.723244
1404619     397.374124
1151040    3485.172794
              ...     
848794     2366.912144
291118     2037.073368
128267     1853.395186
459668     4395.533697
1120216     786.328257
Length: 482738, dtype: float64
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.