J'ai essayé de faire ma propre régression linéaire mais j'ai vu qu'il y a des différences avec LinearRegression de sklearn. En effet, c'est ce qui donne ma propre régression linéaire :

Code Python : Sélectionner tout - Visualiser dans une fenêtre à part
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
 
   >>> 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

Ce qui est très différent de ce que sklearn renvoie :

Code Python : Sélectionner tout - Visualiser dans une fenêtre à part
1
2
3
4
5
6
7
8
 
    array([[3.24519438],
           [2.12164625],
           [0.98814432],
           ...,
           [1.13965629],
           [2.09961495],
           [1.84152796]])

Et aussi, bien sûr, des vraies valeurs:

Code Python : Sélectionner tout - Visualiser dans une fenêtre à part
1
2
3
4
5
6
7
8
 
    array([[2.],
           [1.],
           [1.],
           ...,
           [1.],
           [4.],
           [2.]])


Comment se fait-il que j'obtienne des résultats aussi éloignés des valeurs réelles avec mon propre modèle ?

Voici mon modèle :

Code Python : Sélectionner tout - Visualiser dans une fenêtre à part
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
 
   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

J'ai essayé avec un code minimal et reproductible :

Code : Sélectionner tout - Visualiser dans une fenêtre à part
1
2
3
4
5
    >>>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)
Qui retourne:

Code Python : Sélectionner tout - Visualiser dans une fenêtre à part
1
2
3
4
    [0.8444218515250481,
     0.7579544029403025,
     0.420571580830845,
     0.25891675029296335]

Ce qui est assez étrange ...

Avec un taux d'apprentissage provenant du hessian de la perte.

Le problème pourrait être dû à mon taux d'apprentissage, qui aurait pu être trop important. Le meilleur taux d'apprentissage d'une fonction de coût est strictement inférieur à 2/λ, où λ est la plus grande valeur propre du hessien. Je veux obtenir ce taux d'apprentissage de mon algorithme de descente de gradient. J'ai donc essayé de le calculer en utilisant cette réponse :

Formule mathématique

Code Python : Sélectionner tout - Visualiser dans une fenêtre à part
1
2
3
4
5
    # 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)
renvoie:

(74119951381184.84+0j)

J'ai donc changé alpha_0 pour quelque chose de moins que 2/λ Mais les prédictions semblent toujours très différentes des résultats réels :

Code Python : Sélectionner tout - Visualiser dans une fenêtre à part
1
2
3
4
5
6
7
8
9
10
11
12
    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