2015-09-01

Июль 2015

ТОМ 30 ВЫПУСК 7

Тесты - Линейная регрессия на C#

James McCaffrey Цель решения задачи линейной регрессии — предсказать значение числовой переменной на основе значений одной или более числовых переменных-предикторов. Например, вам может понадобиться спрогнозировать годовой доход некоей персоны на основе его уровня образования, стажа работы и пола (male = 0, female = 1).

Прогнозируемая переменная обычно называется зависимой (dependent variable), а переменные-предикторы — независимыми (independent variables). При наличии только одной переменной-предиктора этот метод иногда называют простой линейной регрессией. Когда имеются две или более переменных-предикторов этот метод называют множественной линейной регрессией (multiple, или multivariate, linear regression).

Хороший способ понять, куда я клоню в этой статье, — взглянуть на демонстрационную программу на рис. 1. Демонстрационная программа на C# прогнозирует годовой доход на основе уровня образования, стажа работы и пола. Демонстрация начинается с генерации 10 синтетических элементов данных. Уровень образование — это значение между 12 и 16, стаж работы — значение между 10 и 30. Пол является индикаторной переменной, где male — эталонное значение, кодируемое как 0, а female кодируется как 1. Доход (в тысячах долларов) записывается в последний столбец. В реальном сценарии вы, вероятно, считывали бы данные из текстового файла, используя какой-то метод с именем наподобие MatrixLoad.

Рис. 1. Линейная регрессия на C#

После генерации синтетических данных демонстрационная программа создает на их основе так называемую матрицу плана (design matrix). Это просто матрица с данными, в которой все значения в ведущем столбце содержат 1.0. Существует несколько алгоритмов, которые можно применять для линейной регрессии; некоторые из них могут использовать матрицу с исходными данными (raw data matrix), тогда как другие работают с матрицей плана. В демонстрационной программе применен метод, требующий матрицу плана.

После создания такой матрицы демонстрационная программа находит значения для четырех коэффициентов (12.0157, 1.0180, 0.5489, –2.9566). Эти коэффициенты иногда называют b-значениями, или бета-значениями. Первое значение, 12.0157, обычно называется в статистике свободным членом (intercept). Это константа, не связанная ни с какими переменными-предикторами. Значения остальных коэффициентов (1.0180, 0.5489, –2.9566) относятся к уровню образования, стажу работы и полу соответственно.

В самой последней части вывода на рис. 1 используются эти значения коэффициентов, чтобы предсказать доход гипотетической персоны с уровнем образования 14, стажем работы 12 лет и полом 0 (мужчина). Прогнозируемый доход составляет 32.86, который вычисляется так:

income = 12.0157 + (1.0180)(14) + (0.5489)(12) + (-2.9566)(0)
                = 12.0157 + 14.2520 + 6.5868 + 0
                = 32.86

Иначе говоря, для прогнозирования с помощью линейной регрессии значения предикторов умножаются на соответствующие коэффициенты и суммируются. Это очень просто. Заметьте, что ведущее значение свободного члена (в этом примере — 12.0157) можно считать коэффициентом, связанным с переменной-предиктором, значение которой всегда равно 1. Отчасти этот факт объясняет наличие столбца со значениями 1.0 в матрице плана.

Суть задачи линейной регрессии заключается в вычислении значений коэффициентов, используя исходные данные или, что эквивалентно, матрицу плана. Это не так легко. В демонстрации применяется метод, называемый обращением матрицы в замкнутой форме (closed form matrix inversion), также известный как обычный метод наименьших квадратов (ordinary least squares method). Альтернативные методы нахождения значений коэффициентов включают итерационный взвешенный метод наименьших квадратов (iteratively reweighted least squares), оценку по методу максимального правдоподобия (maximum likelihood estimation), гребневую регрессию (ridge regression), градиентный спуск (gradient descent) и несколько других.

На рис. 1 демонстрационная программа, прежде чем приступить к прогнозированию, вычисляет показатель, называемый значением R-квадрата (R-squared value), или коэффициентом детерминации (coefficient of determination). R-квадрат — это значение между 0 и 1, которое описывает, насколько хорошо модель прогнозирования подходит для исходных данных. Иногда это выражают в виде «процентной доли отклонения, объясняемой моделью». В свободном толковании можно сказать так: чем ближе R-квадрат к 1, тем лучше модель прогнозирования. В демонстрационной программе это значение составляет 0.7207, или 72%, и для реальных данных оно считалось бы сравнительно высоким (хорошим).

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

Что такое линейная регрессия

Линейная регрессия обычно лучше всего поясняется графиком. Взгляните на график на рис. 2. Данные на этом графике представляют прогнозируемый ежегодный доход на основе всего одной переменной — стажа работы. Каждая из темно-серых точек соответствует точке данных. Например, первый слева элемент данных соответствует стажу = 10 и доходу = 32.06. Линейная регрессия находит два коэффициента: один свободный и один для стажа. Как оказалось, значения коэффициентов равны 27.00 и 0.43.

Рис. 2. Линейная регрессия с одной независимой переменной

Income ($)	Доход ($)
Actual	Реальное значение
Predicted	Прогнозируемое значение
Work (Years)	Стаж (лет)

Значения коэффициентов определяют уравнение прямой, показанное светло-серым цветом на рис. 2. Эта линия (коэффициенты) минимизирует сумму квадратичных отклонений между точками реальных данных (yi) и точками спрогнозированных данных (fi). Два из десяти отклонений отображаются пунктирными линиями на рис. 2. Первым из этих отклонений является yi – fi = 28.6 – 32.6 = –4.0. Заметьте, что отклонения могут быть как положительными, так и отрицательными. Если бы отклонения не были квадратичными, отрицательные и положительные значения могли бы взаимно уничтожаться.

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

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

Решение уравнения наименьших квадратов

Если в задаче линейной регрессии имеются n переменных-предикторов, тогда нужно найти n+1 значений коэффициентов — по одному на каждый предиктор плюс свободный член (intercept). В демонстрационной программе используется самый базовый метод для нахождения значений коэффициентов. Эти значения часто даются с использованием довольно устрашающего уравнения, приведенного на рис. 3. Но это уравнение не настолько сложное, насколько может показаться с первого взгляда.

Рис. 3. Нахождение коэффициентов линейной регрессии с помощью матриц

Греческая буква «бета» напоминает латинскую «B» и представляет значения коэффициентов. Заметьте, что все буквы в уравнении имеют полужирное начертание; в математике это означает, что они представляют объекты с несколькими значениями (матрицы или массивы/векторы), а не простые скалярные величины (обычные числа). Буква X верхнего регистра представляет матрицу плана. Буква X верхнего регистра в степени T означает транспонирование матрицы плана. Символ * обозначает перемножение матриц. Степень –1 указывает на обращение матрицы. Буква Y верхнего регистра — вектор-столбец (матрица с одним столбцом) значений зависимой переменной. Поэтому нахождение значений коэффициентов на самом деле требует понимания матричных операций.

Схемы на рис. 4 иллюстрируют транспонирование, перемножение и обращение матриц. В транспонированной матрице просто меняются местами строки и столбцы. Например, у вас есть матрица 2×3, т. е. одна матрица с двумя строками и тремя столбцами. При транспонировании она превратится в матрицу 3×2, где строки исходной матрицы стали столбцами в транспонированной матрице.

Рис. 4. Три матричные операции, используемые при нахождении коэффициентов линейной регрессии

M is (2x3)	M — (2×3)
MT (transpose) is (3x2	MT (транспонированная) — (3×2)
M is (3x3)	M — (3×3)
M-1 (inverse) is (3x3)	M–1 (обращенная) — (3×3)
A is (3x4)	A — (3×4)
B is (4x2)	B — (4×2)
A * B (product) is (3x2)	A * B (произведение) — (3×2)

Перемножение матриц может показаться слегка странным, если вы не встречались раньше с такими операциями. Если умножить матрицу размером (n x m) на матрицу размером (m x p), результатом будет матрица размером (n x p). Например, перемножение матриц 3×4 * 4×2 дает матрицу 3×2. Подробное обсуждение перемножения матриц выходит за рамки этой статьи, но, как только вы увидите несколько примеров, процесс станет понятен, и вы сможете реализовать его в коде.

Третья матричная операция, необходимая для нахождения значений коэффициентов линейной регрессии, — обращение матрицы, что, увы, трудно понять и сложно реализовать. Для целей этой статьи достаточно знать, что обращение матрицы определяется, только когда матрица имеет одинаковое количество строк и столбцов (квадратная матрица). На рис. 4 показана матрица 3×3 и ее обращение.

Существует несколько алгоритмов, применимых для обращения матрицы. В демонстрационной программе используется метод разложения Дулитла (Doolittle’s decomposition).

Подведем итоги. Задача линейной регрессии с n переменными-предикторами включает нахождение значений для n+1 коэффициентов. Это можно сделать, используя матричные операции: транспонирование, перемножение и обращение. Транспонирование и перемножение просты в понимании и реализации, но нахождение обращенной матрицы — задача трудная..

Структура демонстрационной программы

Чтобы создать демонстрационную программу, я запустил Visual Studio и выбрал шаблон C# Console Application. Я назвал проект LinearRegression. В этой программе нет значимых зависимостей от .NET Framework, поэтому подойдет любая версия Visual Studio.

После загрузки кода шаблона в редактор я переименовал в окне Solution Explorer файл Program.cs в более описательный LinearRegressionProgram.cs, и Visual Studio автоматически переименовала класс Program за меня. В начале кода я удалил все лишние выражения using, оставив только ссылку на пространство имен верхнего уровня System.

Общая структура программы с небольшими правками для экономии места представлена на рис. 5. Вся управляющая логика программы находится в Main. Демонстрационная программа использует подход на основе статических методов, а не объектно-ориентированное программирование.

Рис. 5. Структура программы, демонстрирующей линейную регрессию

using System;
namespace LinearRegression
{
  class LinearRegressionProgram
  {
    static void Main(string[] args)
    {
      Console.WriteLine("Begin linear regression demo");
      // Генерируем синтетические данные.
      // Создаем матрицу плана.
      // Находим коэффициенты линейной регрессии.
      // Вычисляем R-квадратичное значение.
      // Делаем прогноз.
      Console.WriteLine("End linear regression demo");
      Console.ReadLine();
    }
    static double Income(double x1, double x2,
      double x3, double[] coef) { . . }
    static double RSquared(double[][] data,
      double[] coef) { . . }
    static double[][] DummyData(int rows,
      int seed) { . . }
    static double[][] Design(double[][] data) { . . }
    static double[] Solve(double[][] design) { . . }
    static void ShowMatrix(double[][] m, int dec) { . . }
    static void ShowVector(double[] v, int dec) { . . }
    // ----------
    static double[][] MatrixTranspose(double[][] matrix)
      { . . }
    static double[][] MatrixProduct(double[][] matrixA,
      double[][] matrixB) { . . }
    static double[][] MatrixInverse(double[][] matrix)
      { . . }
  // Здесь находятся другие матричные процедуры
  }
} // ns

Метод Income возвращает прогнозируемый доход на основе входных параметров со значениями для уровня образования, стажа работы и пола, используя массив значений коэффициентов. Метод RSquared возвращает R-квадратичное значение модели на основе данных и коэффициентов. Метод DummyData генерирует синтетические данные, применяемые в демонстрации.

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

Большую часть трудной работы берет на себя набор статических методов, выполняющих матричные операции. В демонстрации определена предельно простая матрица — массив массивов. Альтернатива заключается в создании класса Matrix, но, на мой взгляд, этот подход здесь неоправданно сложен. Иногда обычные массивы лучше программно определяемых объектов.

Метод MatrixTranspose возвращает транспонированную матрицу, метод MatrixProduct — результат перемножения двух матриц, а метод MatrixInverse — обращенную матрицу. В демонстрации много вспомогательных методов. В частности, метод MatrixInverse вызывает вспомогательные методы MatrixDuplicate, MatrixDecompose и HelperSolve.

Метод Solve

Центральное место в программе, демонстрирующей линейную регрессию, занимает метод Solve. Определение этого метода начинается с:

static double[] Solve(double[][] design)
{
  int rows = design.Length;
  int cols = data[0].Length;
  double[][] X = MatrixCreate(rows, cols - 1);
  double[][] Y = MatrixCreate(rows, 1);
...

Единственным входным параметром является матрица плана (design). Можно подумать об альтернативном подходе: передавать матрицу с исходными данными, а затем заставлять Solve вызывать вспомогательный метод Design, чтобы получить матрицу плана. Вспомогательный метод MatrixCreate выделяет пространство под матрицу с указанным количеством строк и столбцов и возвращает ее. Локальная матрица X содержит значения независимых переменных (с ведущим значением 1.0). Локальная матрица Y имеет только один столбец и хранит значения зависимой переменной (в нашем примере — годовой доход).

Затем ячейки матриц X и Y заполняются, используя значения в матрице плана:

int j;
for (int i = 0; i < rows; ++i)
{
  for (j = 0; j < cols - 1; ++j)
  {
    X[i][j] = design[i][j];
  }
  Y[i][0] = design[i][j]; // последний столбец
}

Заметьте, что переменная-индекс j объявляется вне вложенных циклов for, чтобы ее можно было использовать для заполнения матрицы Y. Располагая матрицами X и Y, вы можете найти коэффициенты линейной регрессии согласно уравнению, показанному на рис. 3:

...
  double[][] Xt = MatrixTranspose(X);
  double[][] XtX = MatrixProduct(Xt, X);
  double[][] inv = MatrixInverse(XtX);
  double[][] invXt = MatrixProduct(inv, Xt);
  double[][] mResult = MatrixProduct(invXt, Y);
  double[] result = MatrixToVector(mResult);
  return result;
} // Solve

В демонстрации матрица X имеет размер 10×4, поэтому ее транспонированная версия, Xt, получает размер 4×10. Результат произведения Xt и X имеет размер 4×4, равно как и обращение матрицы, inv. В целом, для задачи линейной регрессии с n независимыми переменными-предикторами при использовании метода обращения матрицы вам потребуется найти обращение матрицы размером (n+1)×(n+1). То есть метод обращения не годится для задач линейной регрессии с огромным числом переменных-предикторов.

Результат перемножения обращенной матрицы 4×4 и транспонированной матрицы 4×10 (в коде это переменная invXt) имеет размер 4×10. Результат перемножения invXt и матрицы Y размером 10×1 (в коде это mResult) получается размером 4×1. Эти значения являются нужными вам коэффициентами. Для удобства значения в матрице Y с одним столбцом перемещаются в обычный массив вызовом вспомогательного метода MatrixToVector.

Вычисление R-квадрата

Как уже отмечалось, показатель R-squared (R-квадрат, или коэффициент детерминации) является мерой того, насколько хорошо точки реальных данных соответствуют вычисленной линии регрессии. В математических терминах R-квадрат определяется как R2 = 1 – (SSres / SStot). Член SSres обычно называется остаточной суммой квадратов (residual sum of squares). Это сумма квадратов разности (squared differences) между реальными и предсказанными Y-значениями, как иллюстрирует график на рис. 2. Член SStot — это общая сумма квадратов (total sum of squares). Этот член является суммой квадратов разности между каждым реальным Y-значением и средним всех реальных Y-значений.

Показатель R-squared в линейной регрессии также называется коэффициентом детерминации и имеет отношение (но отличается от него) к другому статистическому показателю с названием «r-squared» («малый r-квадрат»). Интерпретация R-квадрата не совсем однозначна и зависит от предметной области конкретной задачи. В естественных и общественных науках, где данные, как правило, запутанные и неполные, значение R-квадрата от 0.6 и выше часто трактуется как довольно хорошее.

Существует альтернативная мера отклонений, поясняемая моделью регрессии, под названием «скорректированный R-квадрат» (adjusted R-squared). Этот показатель учитывает количество переменных-предикторов и число элементов данных. Для большинства целей применение обычного значения R-квадрата вполне достаточно, чтобы получить представление о прогностическом качестве модели линейной регрессии.

Заключение

Если вы поищете в Интернете примеры того, как выполнять линейную регрессию на каком-нибудь языке программирования, то не найдете много ссылок. Думаю, что есть две основные причины такой относительной нехватки информации. Во-первых, нахождение коэффициентов линейной регрессии с помощью матричных операций довольно сложно — в основном из-за операции обращения матрицы. В некоторых отношениях я считаю метод MatrixInverse в демонстрационной программе одной из самых сложных процедур, которые я когда-либо писал. Во-вторых, существует масса автономных средств, способных выполнять линейную регрессию, в частности Excel с надстройкой Data Analysis. Необходимость встраивания решения линейной регрессии прямо в код программной системы возникает относительно редко.

Линейная регрессия изучалась десятилетиями, и есть много способов расширения этого метода. Например, вы можете ввести так называемые эффекты взаимодействия, которые комбинируют две или более переменных-предикторов. Эти расширения иногда называют универсальными линейными моделями, чтобы отличать их от базовой формы линейной регрессии.

По моему мнению, линейная регрессия является своего рода «Hello World» среди методов классической статистики. Четкого и общепризнанного различия между классической статистикой и машинным обучением нет, но я склонен считать методами классической статистики те из них, которые впервые были исследованы математиками в начале 1900-х. В моем понимании методы машинного обучения вроде классификации на основе нейронных сетей появились гораздо позднее — первые упоминания встречаются в 1950-х. Линейная регрессия из классической статистики тесно связана с методом машинного обучения, называемым логистической регрессией, о которой я рассказывал в нескольких выпусках этой рубрики.

Джеймс Маккафри (Dr. James McCaffrey) работает на Microsoft Research в Редмонде (штат Вашингтон). Принимал участие в создании нескольких продуктов Microsoft, в том числе Internet Explorer и Bing. С ним можно связаться по адресу jammc@microsoft.com.

Выражаю благодарность за рецензирование статьи эксперту Microsoft Research Чарльзу Паркеру (Charles Parker).