Множественная линейная регрессия используется для оценки взаимосвязи между многими независимыми переменными и одной зависимой переменной.

Фон

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

Множественная линейная регрессия

В то время как простая линейная регрессия имеет уравнение Ŷ = w₁X₁ + w₀X₀, множественная линейная регрессия имеет общую формулу для kнезависимые переменные:

Обратите внимание, что X₀ — это столбец единиц для смещения; это позволяет использовать обобщенную формулу, обсуждавшуюся в первой статье. Как показывает формула, множественная линейная регрессия помогает определить взаимосвязь между многими независимыми переменными (X) и одной зависимой переменной (Ŷ). ). Это достигается путем изучения значений каждого веса (w).

Чтобы продемонстрировать это в действии, можно использовать множественную линейную регрессию с 3 весами:

Эта формула создаст «плоскость наилучшего соответствия» для трехмерных данных:

Реализация

Как и в предыдущем примере, можно использовать уравнение «чертежа» и добавить случайность. Затем модель может попытаться узнать веса. При использовании модели на реальных данных обычно данные разбиваются как минимум на два набора: набор поездов и набор тестов. Набор поездов используется для обучения модели и получения весов. Набор тестов используется для оценки производительности модели на данных, которые она никогда раньше не видела. Если он хорошо работает в обоих случаях, он, скорее всего, будет полезен для реальных приложений. Обычно используется разделение 80% поезд, 20% тест.

Для этого примера можно создать 1000 выборок и разделить их на тестовые и обучающие наборы. Кроме того, поскольку есть две независимые переменные и смещение, будет три столбца. Каждый столбец будет представлять независимую переменную или отклонение, а каждая строка будет представлять выборку переменных (X₀, X₁, X₂). Форма общих данных будет представлять собой матрицу размером (1000, 3); общая форма будет такой (n образцов, количество признаков). Помните, что независимые переменные также известны как функции в машинном обучении. Набор поездов будет иметь размер (800, 3), а тестовый набор будет иметь размер (200, 3).

Данные для этой статьи будут основаны на Y = 6X₂ + 3X₁ + 2. Это означает, что w₀ равно 2, w₁ равно 3 и w₂ равно 6. В приведенном ниже примере для X₁и X₂, генерируется 1000 значений от -250 до 250. и для X₀ генерируется 1000 единиц. Они преобразуются в столбцы и складываются горизонтально, чтобы создать матрицу размером (1000, 3). Выходные данные генерируются с использованием вышеупомянутого уравнения, и добавляются значения из нормального распределения со стандартным отклонением 10.

import torch

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

# create ones for the bias | 1000 ones
X0 = torch.ones(1000).reshape(-1,1)

# create values for the first feature | 1000 numbers from -250 to 250
X1 = (500*(torch.rand(1000) - 0.5)).reshape(-1,1) 

# create values for the second feature | 1000 numbers from -250 to 250
X2 = (500*(torch.rand(1000) - 0.5)).reshape(-1,1)

# stack data together, X0 = X[:,0], X1 = X[:,1], X2 = X[:,2]
X = torch.hstack((X0, X1,X2))

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

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

Эти данные можно просмотреть перед разделением.

import plotly.express as px

fig = px.scatter_3d(x=X[:,1].flatten(),
                    y=X[:,2].flatten(),
                    z=Y.flatten())

fig.update_traces(marker_size=3)
fig.update_layout(scene = dict(xaxis_title='X<sub>1</sub>', 
                               yaxis_title='X<sub>2</sub>', 
                               zaxis_title='Y'))

Теперь данные можно разделить на тестовые и обучающие:

# split the data
Xtrain, Xtest = X[:800], X[800:]
Ytrain, Ytest = Y[:800], Y[800:]

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

# 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)

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

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

torch.manual_seed(5)
w = torch.rand(size=(3, 1))
w
tensor([[0.83],
        [0.13],
        [0.91]])

Текущую плоскость наилучшего соответствия и ее MSE можно проанализировать ниже. Самолет выделен оранжевым цветом, набор поездов красным, а тестовый набор зеленым.

import plotly.graph_objects as go
def plot_model(x1_range, x2_range):
  """
    Inputs:
      x1_range: x1-axis range [low, high]
      x2_range: x2-axis range [low, high]

    Global Variables:
      Xtrain: array of inputs | (n train samples, num features)
      Ytrain: array of expected outputs | (n train samples, 1)
      Xtest:  array of inputs | (n test samples, num features)
      Xtrain: array of expected outputs | (n test samples, 1)
      
    Output:
      prints plane of best fit
  """ 

  # meshgrid of possible combinations of (X1, X2)
  X1_plot, X2_plot = torch.meshgrid(torch.arange(x1_range[0], x1_range[1], 5),
                                    torch.arange(x2_range[0], x2_range[1], 5))
  X0_plot = torch.ones(X1_plot.shape)
  
  # stack together each point (X1, X2) = (X, Y)
  X_plot = torch.hstack((X0_plot.reshape(-1,1),
                         X1_plot.reshape(-1,1), 
                         X2_plot.reshape(-1,1)))
  
  # all possible model predictions (Yhat = Z)
  Yhat = model(w, X_plot)

  # model's plane of best fit
  fig = go.Figure(data=[go.Mesh3d(x=X_plot[:,1].flatten(), 
                                  y=X_plot[:,2].flatten(), 
                                  z=Yhat.flatten(), 
                                  color='orange', 
                                  opacity=0.50)])
  
  # training data
  fig.add_scatter3d(x=Xtrain[:,1].flatten(),
                    y=Xtrain[:,2].flatten(),
                    z=Ytrain.flatten(), 
                    mode="markers",
                    marker=dict(size=3),
                    name="train")
  
  # test data
  fig.add_scatter3d(x=Xtest[:,1].flatten(),
                    y=Xtest[:,2].flatten(),
                    z=Ytest.flatten(), 
                    mode="markers",
                    marker=dict(size=3),
                    name="test")
  
  # name axes
  fig.update_layout(scene = dict(xaxis_title='X<sub>1</sub>', 
                                 yaxis_title='X<sub>2</sub>', 
                                 zaxis_title='Y'))

  fig.show()

plot_model([-250,250], [-250,250])

MSE(model(w,Xtrain), Ytrain)
tensor(653812.81)

Теперь можно создать обучающий цикл для минимизации MSE. Используя 50 000 эпох и скорость обучения 0,00004, вывод становится чрезвычайно точным. Эти значения были выбраны опытным путем. Также можно увидеть как поезд, так и тестовый MSE.

torch.manual_seed(5)
w = torch.rand(size=(3, 1))

lr = 0.00004
epochs = 50000

# 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) % 10000 == 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_model([-250,250], [-250,250])
epoch: 10000
weights: tensor([[1.51],
        [2.98],
        [6.00]])
Train MSE: tensor(87.52)
Test MSE: tensor(70.03)
==========
epoch: 20000
weights: tensor([[1.82],
        [2.98],
        [6.00]])
Train MSE: tensor(87.20)
Test MSE: tensor(70.04)
==========
epoch: 30000
weights: tensor([[1.96],
        [2.98],
        [6.00]])
Train MSE: tensor(87.12)
Test MSE: tensor(70.11)
==========
epoch: 40000
weights: tensor([[2.02],
        [2.98],
        [6.00]])
Train MSE: tensor(87.09)
Test MSE: tensor(70.16)
==========
epoch: 50000
weights: tensor([[2.05],
        [2.98],
        [6.00]])
Train MSE: tensor(87.08)
Test MSE: tensor(70.18)
==========

Модель предсказала, что плоскость наилучшего соответствия будет Ŷ = 6X₂ + 2,98X₁ + 2,05вместо Y = 6X₂ + 3X₁ + 2. СКО теста и обучения находятся в пределах 17 точек друг от друга, что указывает на то, что модель обобщает невидимые тестовые данные. Единственным ограничением этого подхода является количество эпох, необходимых для минимизации функции потерь. Альтернативным подходом может быть использование решения в закрытой форме, которое было рассмотрено в предыдущей статье Введение в машинное обучение в Python: нормальное уравнение для регрессии в Python. Решение в закрытой форме не требует скорости обучения или эпох для получения весов для минимизации.

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

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

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

w = NormalEquation(Xtrain,Ytrain)
tensor([[2.19],
        [2.98],
        [6.00]])

MSE также можно рассчитать:

MSE(model(w, Xtrain), Ytrain), MSE(model(w, Xtest), Ytest)
(tensor(87.08), tensor(70.20))

Веса и MSE из подходов нормального уравнения и градиентного спуска почти идентичны. В данном случае оба варианта равноправны.

Заключение

Множественная линейная регрессия полезна для определения взаимосвязи между двумя или более независимыми переменными или признаками и одной зависимой переменной. Это также основа полиномиальной регрессии, которая будет рассмотрена в следующей статье: Введение в машинное обучение в Python: полиномиальная регрессия.

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

Рекомендации

  1. Пример регрессии Sklearn