Введение в машинное обучение в Python: множественная линейная регрессия

Полиномиальная регрессия может идентифицировать нелинейную связь между независимой переменной и зависимой переменной.

Фон

Эта статья является третьей в серии статей о регрессии, градиентном спуске и MSE. Предыдущие статьи посвящены Простой линейной регрессии, Нормальному уравнению регрессии и Множественной линейной регрессии.

Полиномиальная регрессия

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

Обратите внимание, что X₀ – это столбец с единицами смещения; это позволяет использовать обобщенную формулу, обсуждавшуюся в первой статье. Используя приведенное выше уравнение, каждую «независимую» переменную можно рассматривать как возведенную в степень версию X₁.

Это позволяет использовать одну и ту же модель из множественной линейной регрессии, поскольку необходимо идентифицировать только коэффициенты каждой переменной. В качестве примера можно создать простую полиномиальную модель третьей степени. Его уравнение следующее:

Обобщенные функции для модели, градиентного спуска и MSE можно использовать из предыдущих статей:

# line of best fit
def model(w, X):
  """
    Inputs:
      w: array of weights | (num features, 1)
      X: array of inputs  | (n samples, num features)

    Output:
      returns the output of X@w | (n samples, 1)
  """

  return torch.matmul(X, w)
# mean squared error (MSE)
def MSE(Yhat, Y):
  """
    Inputs:
      Yhat: array of predictions | (n samples, 1)
      Y: array of expected outputs | (n samples, 1)
    Output:
      returns the loss of the model, which is a scalar
  """
  return torch.mean((Yhat-Y)**2) # mean((error)^2)
# optimizer
def gradient_descent(w):
  """
    Inputs:
      w: array of weights | (num features, 1)

    Global Variables / Constants:
      X: array of inputs  | (n samples, num features)
      Y: array of expected outputs | (n samples, 1)
      lr: learning rate to scale the gradient

    Output:
      returns the updated weights
  """ 

  n = X.shape[0]

  return w - (lr * 2/n) * (torch.matmul(-Y.T, X) + torch.matmul(torch.matmul(w.T, X.T), X)).reshape(w.shape)

Создание данных

Теперь все, что требуется, — это некоторые данные для обучения модели. Можно использовать функцию «чертежа» и добавить случайность. Здесь используется тот же подход, что и в предыдущих статьях. Схему можно увидеть ниже:

Можно создать набор поездов размером (800, 4) и тестовый набор размером (200, 4). Обратите внимание, что каждая функция, кроме смещения, является экспоненциальной версией первой.

import torch

torch.manual_seed(5)
torch.set_printoptions(precision=2)

# features
X0 = torch.ones((1000,1))
X1 = (100*(torch.rand(1000) - 0.5)).reshape(-1,1) # generates 1000 random numbers from -50 to 50
X2, X3 = X1**2, X1**3
X = torch.hstack((X0,X1,X2,X3))

# normal distribution with a mean of 0 and std of 8
normal = torch.distributions.Normal(loc=0, scale=8)

# targets
Y = (3*X[:,3] + 2*X[:,2] + 1*X[:,1] + 5 + normal.sample(torch.ones(1000).shape)).reshape(-1,1)

# train, test
Xtrain, Xtest = X[:800], X[800:]
Ytrain, Ytest = Y[:800], Y[800:]

После определения начальных весов данные могут быть нанесены на график с линией наилучшего соответствия.

torch.manual_seed(5)
w = torch.rand(size=(4, 1))
w
tensor([[0.83],
        [0.13],
        [0.91],
        [0.82]])
import matplotlib.pyplot as plt

def plot_lbf():
  """
    Output:
      prints the line of best fit in comparison to the train and test data
  """

  # plot the train and test sets
  plt.scatter(Xtrain[:,1],Ytrain,label="train")
  plt.scatter(Xtest[:,1],Ytest,label="test")

  # plot the line of best fit
  X1_plot = torch.arange(-50, 50.1,.1).reshape(-1,1) 
  X2_plot, X3_plot = X1_plot**2, X1_plot**3
  X0_plot = torch.ones(X1_plot.shape)
  X_plot = torch.hstack((X0_plot,X1_plot,X2_plot,X3_plot))

  plt.plot(X1_plot.flatten(), model(w, X_plot).flatten(), color="red", zorder=4)

  plt.xlim(-50, 50)
  plt.xlabel("$X$")
  plt.ylabel("$Y$")
  plt.legend()
  plt.show()

plot_lbf()

Обучение модели

Чтобы частично минимизировать функцию стоимости, можно использовать скорость обучения 5e-11 и 500 000 эпох с градиентным спуском.

lr = 5e-11
epochs = 500000

# update the weights 1000 times
for i in range(0, epochs):
  # update the weights
  w = gradient_descent(w)

  # print the new values every 10 iterations
  if (i+1) % 100000 == 0:
    print("epoch:", i+1)
    print("weights:", w)
    print("Train MSE:", MSE(model(w,Xtrain), Ytrain))
    print("Test MSE:", MSE(model(w,Xtest), Ytest))
    print("="*10)

plot_lbf()
epoch: 100000
weights: tensor([[0.83],
        [0.13],
        [2.00],
        [3.00]])
Train MSE: tensor(163.87)
Test MSE: tensor(162.55)
==========
epoch: 200000
weights: tensor([[0.83],
        [0.13],
        [2.00],
        [3.00]])
Train MSE: tensor(163.52)
Test MSE: tensor(162.22)
==========
epoch: 300000
weights: tensor([[0.83],
        [0.13],
        [2.00],
        [3.00]])
Train MSE: tensor(163.19)
Test MSE: tensor(161.89)
==========
epoch: 400000
weights: tensor([[0.83],
        [0.13],
        [2.00],
        [3.00]])
Train MSE: tensor(162.85)
Test MSE: tensor(161.57)
==========
epoch: 500000
weights: tensor([[0.83],
        [0.13],
        [2.00],
        [3.00]])
Train MSE: tensor(162.51)
Test MSE: tensor(161.24)
==========

Даже с 500 000 эпох и чрезвычайно малой скоростью обучения модель не может определить первые два веса. Хотя текущее решение является очень точным с MSE 161,24, для его полной минимизации, вероятно, потребуются миллионы эпох. Это одно из ограничений градиентного спуска для полиномиальной регрессии.

Нормальное уравнение

В качестве альтернативы можно использовать нормальное уравнение из второй статьи для прямого вычисления оптимизированных весов:

def NormalEquation(X, Y):
  """
    Inputs:
      X: array of input values | (n samples, num features)
      Y: array of expected outputs | (n samples, 1)
      
    Output:
      returns the optimized weights | (num features, 1)
  """
  
  return torch.inverse(X.T @ X) @ X.T @ Y

w = NormalEquation(Xtrain, Ytrain)
w
tensor([[4.57],
        [0.98],
        [2.00],
        [3.00]])

Нормальное уравнение способно сразу определить правильные значения для каждого веса, а MSE для каждого набора примерно на 100 пунктов ниже, чем при градиентном спуске:

MSE(model(w,Xtrain), Ytrain), MSE(model(w,Xtest), Ytest)
(tensor(60.64), tensor(63.84))

Заключение

При реализации простой линейной, множественной линейной и полиномиальной регрессии следующие две статьи будут посвящены регрессии Лассо и Риджа. Эти типы регрессии вводят два важных понятия в машинном обучении: переоснащение и регуляризация.

Пожалуйста, не забудьте поставить лайк и подписаться! :)