Линейная регрессия. Простая реализация на Python.

Захотелось глубже разобраться с основными алгоритмами машинного обучения и анализа данных. Для этой цели решил сделать собственные реализации этих алгоритмов на языке Python. Начну с обычной линейной регрессии, а в дальнейшем сделаю аналогичные заметки про нейронные сети, решающие деревья и т.д.

Цитата из Wiki: Линейная регрессия (англ. Linear regression) — используемая в статистике регрессионная модель зависимости одной (объясняемой, зависимой) переменной от другой или нескольких других переменных (факторов, регрессоров, независимых переменных) с линейной функцией зависимости.

Итак. Буду делать упрощенный аналог LinearRegression() из sklearn.

Финальный вариант кода можно скачать или посмотреть на github.

Краткая теория

Если абстрагироваться от умных слов, то основная задача заключается в следующем:

  • нужно реализовать модель, которая будет способна достаточно хорошо описывать зависимость некой результирующей переменной от ряда признаков (регрессоров).

Что значит “хорошо описывать”? Это означает, что предсказание, сделанное моделью, должно минимально отличаться от реальных значений. То есть нам нужно минимизировать ошибку. Делать это я буду с помощью градиентного спуска.

Что это за зверь такой?

Предположим, мы хотим построить модель для предсказания стоимости квартир на основе площади помещения (для простоты пока возьмем всего один признак).

Пусть у нас есть размеченные обучающие данные в формате [площадь, цена]. Будем строить линейную регрессию. Аппроксимируем реальную зависимость цены от площади следующим образом:

Price(w0, w1) = w0 + w1*Square,
где w0, w1 – это веса в нашей модели

Каждой комбинации весов соответствует определенное значение ошибки. Если изобразить все возможные значения на графике, то получится какая-то поверхность. Например, такая: 

surface
(если признаков много, то поверхность многомерная)

Нам нужно найти точку минимума среди множества значений функционала ошибки. На приведенном графике это синяя/фиолетовая область.

Точке минимума соответствуют некоторые веса. Если мы их найдем, то получим модель, которая дает нам наименьшую ошибку. В общем нам как-то надо попасть в точку минимума. Здесь на помощь и приходит метод градиентного спуска.

Как он работает:

  1. На первом шаге нужно “встать” в любую точку среди множества значений функционала ошибки. То есть нужно задать веса в нашей модели случайными значениями.
  2. Далее нужно сделать шаг в направлении, в котором происходит наискорейшее убывание функционала ошибки. Это направление противоположно градиенту (вектору, в направлении которого происходит наискорейший рост).
  3. Если ошибка уменьшается, повторяем пункт 2.
  4. Останавливаемся, когда достигнут локальный минимум.

В последнем пункте ключевое слово – “локальный” минимум. Алгоритм градиентного спуска не гарантирует, что мы найдем глобальный минимум.

Реализация

В своей программе я использовал несколько библиотек:

  • numpy  для работы с массивами
  • matplotlib – для построения графиков в демонстрационных целях
  • sklearn – для сравнения их реализации линейной регрессии с моей моделью
import numpy as np
from matplotlib import pyplot as plt
from sklearn.linear_model import LinearRegression
from sklearn.metrics import mean_squared_error
%matplotlib inline

Далее полная реализация класса LinRegression, в котором реализована линейная регрессия. Дополнительные комментарии после блока с кодом.

class LinRegression:
    '''Linear Regression.'''

    def __init__(self):
        self.about = "Linear Regression by Sergei Bernadsky"
        self.W = [] # model's weights
        self.fscaling = False # is feature scaling used
        
    def cost(self, y_real, y_pred): 
        # cost function for gradient descent algorithm
        return np.sum((y_pred-y_real)**2)/(len(y_real))
    
    def gradient_descent_step(self, learning_rate, dy, m, n, X_tr):
        # one gradient descent step
        s = (np.dot(dy.T, X_tr)).reshape(n, 1)
        dW = 2*(learning_rate*s/m).reshape(n, 1)
        return self.W - dW
    
    def normalize(self, X):
        # normilize X table
        for j in range(X.shape[1]):
            X[:,j] = X[:,j]/np.max(X[:,j])
        return X
    
    def fit(self, X, y, learning_rate = 0.99, nsteps = 3000, e = 0.000000001,
            weight_low = 0, weight_high = 1,
            fscaling = False, kweigths = 1, random_state = 0):
        # train our Linear Regression model
        
        np.random.seed(random_state)
        X = X.astype(float)
        
        # Normilize process
        if fscaling == True:
            X = self.normalize(X)
            self.fscaling = True

        m = X.shape[0]
        # add one's column to X
        X = np.hstack( (np.ones(m).reshape(m, 1), X) )
        n = X.shape[1]
        
        # Weights: random initialization
        self.W = np.random.randint(low = weight_low, high = weight_high, size=(n, 1))
            
        y_pred = np.dot(X, self.W)
        cost0 = self.cost(y, y_pred)

        y = y.reshape(m, 1)
        k = 0
        
        ########## Gradient descent's steps #########
        while True:
            dy = y_pred - y
            W_tmp = self.W
            self.W = self.gradient_descent_step(learning_rate, dy, m, n, X)
            y_pred = np.dot(X, self.W)

            cost1 = self.cost(y, y_pred)
            k += 1

            if (cost1 > cost0):
                self.W = W_tmp
                break    
                
            if ((cost0 - cost1) < e) or (k == nsteps):
                break
                
            cost0 = cost1
        #############################################

        return self.W # return model's weights
    
    def predict(self, X):
        m = X.shape[0]
        if self.fscaling == False:
            return np.dot( np.hstack( (np.ones(m).reshape(m, 1),
                                       X.astype(float)) ) ,
                          self.W)
        else:
            return np.dot( np.hstack( (np.ones(m).reshape(m, 1),
                                       self.normalize(X.astype(float))) ),
                          self.W)

Веса на каждом шаге градиентного спуска обновляются по формуле:

    \[ W_j = W_j - LearningRate\cdot\frac{\partial}{\partial W}J(W) \]

где LearningRate – это коэффициент, определяющий скорость градиентного спуска, J(W) – функционал ошибки (в формуле используется частная производная от J).

Для расчета ошибки J используется такая формула:

    \[ J = \frac{1}{m}\sum_{i=0}^{m}{(y_i-h_i)^2} \]

где y – фактические значения целевой переменной, а h – предсказания нашей модели. То есть я использую в качестве метрики mse – mean squared error (среднеквадратическая ошибка).

Таким образом исходная формула для обновленных весов записывается в таком виде:

    \[W_j = W_j - LearningRate\cdot\frac{2}{m}\sum_{i=0}^{m}{(y_i-h_i)}\cdot X_i_j\]

Демонстрация работы и сравнение с реализацией в sklearn:

Для примера используется датасет diabets из библиотеки sklearn. В нем 10 признаков и 442 строки данных.

from sklearn import datasets

# Load the diabetes dataset
diabetes = datasets.load_diabetes()
X = diabetes.data
y = diabetes.target

Проверка работы моей модели:

lr = LinRegression()
lr.fit(X, y, learning_rate = 0.997, random_state = 0, weight_low = -900, weight_high = 900, nsteps=3000)

xx = [i for i in range(X.shape[0])]
y1 = lr.predict(X)
print 'MSE1 (My LR model):', mean_squared_error(y, y1)

f=0
t=40
plt.plot(xx[f:t], y[f:t], color='r', linewidth=4, label='y')
plt.plot(xx[f:t], y1[f:t], color='b', linewidth=2, label='predicted y')
plt.ylabel('Target label')
plt.xlabel('Line number in dataset')
plt.legend(loc=4)
plt.show()

Проверка линейной регрессии из sklearn:

%%time

lr2 = LinearRegression()
lr2.fit(X,y)
y2 = lr2.predict(X)

plt.figure()
plt.plot(xx[f:t], y[f:t], color='r', linewidth=4, label='y')
plt.plot(xx[f:t], y2[f:t], color='b', linewidth=2, label='predicted y')
plt.ylabel('Target label')
plt.xlabel('Line number in dataset')
plt.legend(loc=4)
plt.show()

Вот фрагмент графиков, демонстрирующих работу моделей (на первом моя реализация, на втором sklearn):

my-model
sklearn-model

По оси X – номер строки в датасете. По Y – значение целевой переменной.

Модели очень похожи. Реализация из sklearn чуть-чуть более точная, но совсем не значительно. Подозреваю, что там более оптимально рассчитывается шаг градиентного спуска. И еще он, скорее всего, не фиксированный, а меняется по мере приближения к точке минимума. В моей реализации этого нет, но и цели пока такой не стояло.

Leave a Reply

Your email address will not be published. Required fields are marked *