Линейная регрессия теория и практика

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

Определение одномерной модели и её решение

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

Итак, давайте рассмотрим простой пример. Все мы знакомы с законом Ома:

V=I*R,

где V – это напряжение, I – сила тока, а R – сопротивление. Между прочим, это уравнение прямой, которое можно представить в виде y=mx+b. Действительно, если представить, что y=V, x=I, m=R, a b=0, то мы получим наш закон Ома.

Линейная регрессия теория и практика

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

К счастью, для этого у нас есть несколько батарей с разным напряжением, и у нас есть амперметр, измеряющий силу тока. Будем поочередно подключать к цепи эти батареи и записывать значения напряжения и силы тока. У нас получится ряд чисел:

V1, I1; V2, I2; …; Vn, In.

Далее нанесём их на плоскость координат, откладывая по оси x силу тока и напряжение по оси y, и без особого удивления обнаружим, что они выстраиваются в почти идеальную прямую.

Почему мы говорим «почти»? Как вы знаете, в любых экспериментах есть определённая погрешность. Она может быть вызвана погрешностью измерения или, скажем, тем, что используемые медные провода тоже имеют собственное сопротивление. Не забывайте, что физические уравнения – это лишь идеальные модели реального мира. По этому поводу ещё часто говорят: «Карта – это ещё не территория».

Линейная регрессия теория и практика

Что же дальше? Далее мы берём линейку и проводим прямую – линию наилучшую соответствия так, чтобы она проходила через максимальное количество точек. При этом примерно одинаковое количество точек окажется как выше, так и ниже проведённой прямой.

Получив прямую, мы можем рассчитать её наклон, который определяется, как вы знаете, её угловым коэффициентом. Он показывает, насколько изменится напряжение при данном изменении силы тока: ΔV/ΔI, или, в общем случае,

наклон =  .

А поскольку V=I*R, то найденная величина наклона и есть искомое сопротивление:

R = наклон = ΔV/ΔI

Линейная регрессия теория и практика

Данный метод применим, конечно же, не только к закону Ома. По оси y не обязательно откладывать именно напряжение, как и по оси x – не обязательно силу тока. Они, как мы уже говорили, могут обозначать что угодно. К примеру, это может быть зависимость артериального давления от возраста или, скажем, зависимость давления от веса тела. Или биржевые котировки в зависимости от настроений в твиттере.

Учёные знают, что не имеет никакого значения, что именно мы измеряем, так как методы обработки данных остаются прежними. Мне нравится формулировать это обстоятельство одной фразой – «все данные – одинаковы!». Будь вы биологом или финансистом, вы пользуетесь одной и той же линейной регрессией. Нет никакой «линейной регрессии для биологов» или «линейной регрессии для финансистов». Есть просто линейная регрессия.

Поэтому будет лучше, если по оси y мы будем откладывать просто некоторую величину Y, а по оси x – некоторую величину X. Переформулируем теперь нашу задачу в более общем виде с X и Y.

Линейная регрессия теория и практика

Пусть, как и прежде, у нас есть ряд точек: {(x1, y1), (x2, y2), …, (xN, yN)}. Нанесём их на плоскость координат. Теперь нам необходимо найти линию наилучшего соответствия. Конечно же, мы живём в эпоху компьютеров и информатики, и есть куда более точный способ её нахождения, чем проводить её на бумаге с линейкой в руке. Как же её найти?

Начнём с того факта, что искомая линия является прямой и имеет вид y=ax+b, так что и предсказываемое значение величины y (обозначим его ŷ) имеет вид:

ŷi=axi+b.

Линейная регрессия теория и практика

Нам нужно, чтобы все yi, полученные из эксперимента, были как можно ближе к предсказываемому ŷi, рассчитанному из заданных xi.

Как же это сделать?

Первое, что приходит в голову, – посчитать сумму получающихся несоответствий между фактическим значением yi  и предсказываемым значением ŷi при всех значениях i от 1 до N:

\sum_{i=1}^{N} = (y_i - \widehat y_i)

Линейная регрессия теория и практика

Почему так делать нельзя?

Предположим, в одном случае у нас отклонение +5 единиц, а в другом -5 единиц. Если мы их просто сложим, то получим ноль, и будет казаться, что отклонение равно нулю, хотя на самом деле это вовсе не так. На самом же деле нам нужно, чтобы любое отклонение от предсказываемого ŷ выражалось положительным числом. Стандартный способ этого достичь – возвести разницу в квадрат:

\sum_{i=1}^{N} = (y_i - \widehat y_i)^2

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

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

Например, пусть у нас есть формула  E=0,5t2-t,

и мы хочем узнать, при каком значении t значение E становится наименьшим. Тогда мы берём производную по t и приравниваем её к нулю:   dE/dt=t-1=0 ⇒t=1

Линейная регрессия теория и практика

В нашем случае всё не так просто, потому что у нас гораздо больше независимых переменных. Но давайте всё же перепишем формулу и посмотрим, что это нам даст:

E= \sum_{i=1}^{N} (y_i - \widehat y_i)^2

Поскольку ŷi=axi+b, то мы можем записать

E= \sum_{i=1}^{N} (y_i - (ax_i + b))^2

Помните, что yi и xi – это данные, полученные в эксперименте, а найти нам нужно a и b, которые характеризуют наклон прямой и точку её пересечения с осью у. В отличие от первого примера с Е, нам нужно найти минимумы по a и b одновременно, поэтому придётся использовать частные производные.

Итак, первое, что нам нужно, – найти производную по а. Вот она:

\frac{\partial E}{\partial a}=\sum_{i=1}^{N} 2 (y_i-(ax_i+b)) (-x_i).

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

\sum_{i=1}^{N}2(y_i-(ax_i+b)) (-x_i)=0;

- \sum_{i=1}^{N} y_ix_i + a \sum_{i=1}^{N} x_i^2 +b \sum_{i=1}^{N} =0.

Остановитесь на минутку, чтобы понять, как мы получили такой результат, и попробуйте на бумаге проделать расчёты. Перенесём сумму yixi в другую часть уравнения, чтобы избавиться от отрицательных величин:

a \sum_{i=1}^{N} x_i^2 + b \sum_{i=1}^{N} x_i = \sum_{i=1}^{N} y_ix_i.

Итак, у нас есть одно уравнение с двумя неизвестными. Не забывайте, что xi и yi – это наши данные, полученные из эксперимента, а N также известно – это количество опытов.

Но нам ещё нужно взять производную по b. Давайте этим и займёмся. Сделаем всё то же, что и раньше, – найдём производную, приравняем её к нулю и упростим полученное выражение:

\frac{\partial E}{\partial a} =\sum_{i=1}^{N} 2 (y_i - (ax_i + b))(-1)=0;

- \sum_{i=1}^{N} y_i +a\sum_{i=1}^{N} x_i +b \sum_{i=1}^{N}1 =0,

a\sum_{i=1}^{N} x_i+bN=\sum_{i=1}^{N}y_i.

Попытайтесь самостоятельно на бумаге получить данный результат. Итак, у нас получается вот такая интересная ситуация:

a\sum_{i=1}^{N}x_i^2 + b\sum_{i=1}^{N} x_i=\sum_{i=1}^{N}y_ix_i,

a \sum_{i=1}^{N} x_i + bN = \sum_{i=1}^{N} y_i.

Надеюсь, вы помните, что нам необходимо найти a и b, поскольку xi, yi и N нам известны. Таким образом, имеем два уравнения с двумя неизвестными. Это должно дать нам единственное решение для a и b. Упростим наши уравнения, введя буквенные обозначения:

C=\sum_{i=1}^{N} x_i^2, D=\sum_{i=1}^{N}x_i, E = \sum_{i=1}^{N}x_i y_i, F=\sum_{i=1}^{N} y_i.

Получим:

aC+bD=E,

aD+bN=F.

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

Умножим первое уравнение на D, а второе на C:

 aDC+bD2=ED,

aDC+bNC=FC,

вычтем второе уравнение из первого, благодаря чему сокращается a, и мы можем найти b:

b(D2-NC)=(ED-FC),

Теперь сделаем то же самое для а. Умножим первое уравнение на N, а второе на D и вычтем второе уравнение из первого:

aNC+bND=EN,

aD2+bND=FD.

Получим решение для а:

a(NC-D2)=(EN-FD),

Итак, мы нашли а и b. Перепишем теперь выражения для них в исходном виде, заменив буквы выражениями с xi и yi, поскольку именно они представляют наши данные эксперимента:

 a =\frac{N \sum_{i=1}^{N} y_i x_i - \sum_{i=1}^{N} x_i \sum_{i=1}^{N} y_i} {N \sum_{i=1}^{N} x_i^2 - (\sum_{i=1}^{N} x_i)^2},

b =\frac{\sum_{i=1}^{N} x_i\sum_{i=1}^{N} y_i x_i - \sum_{i=1}^{N} y_i \sum_{i=1}^{N} x_i^2} {(\sum_{i=1}^{N} x_i)^2 - N \sum_{i=1}^{N} x_i^2}.

Обратите внимание, что если в выражении для b поменять знак, то в обоих выражениях получатся одинаковые знаменатели:

a =\frac{N \sum_{i=1}^{N} y_i x_i - \sum_{i=1}^{N} x_i \sum_{i=1}^{N} y_i} {N \sum_{i=1}^{N} x_i^2 - (\sum_{i=1}^{N} x_i)^2},

b =\frac{\sum_{i=1}^{N} y_i\sum_{i=1}^{N} x_i^2 - \sum_{i=1}^{N} x_i \sum_{i=1}^{N} y_i x_i} {N \sum_{i=1}^{N} x_i^2 - (\sum_{i=1}^{N} x_i)^2}.

Следующее, что мы можем сделать, чтобы упростить полученные выражения, – найти среднее значение. Напомним, что средняя величина переменной x определяется как сумма всех x, делённая на количество переменных N:

\overline x =\frac{1}{N} \sum_{i=1} ^ {N}x_i.

То же самое можно сделать и в том случае, если в выражение входят и x, и y:

\overline {xy}=\frac{1}{N} \sum_{i=1}^{N} x_i y_i.

Переписав наши выражения для a и b, мы упростим наши выражения:

a = \frac{\overline {xy} - \overline x \overline y}{\overline {x^2} - \overline x^2};\; b = \frac{\overline y \overline {x^2} - \overline x \overline {xy}}{\overline {x^2} - \overline x^2}.

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

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

total = 0

for i in xrange(N):

  total +=a[i]*b[i].

Вместо этого воспользуйтесь специальной функцией dot:

total = a.dot(b) или total = np.dot(a,b),

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

Программирование одномерной модели на языке Python

В этой разделе статьи мы начнём программировать. Мы напишем код для одномерной линейной регрессии.

Начнем всё с нуля, чтобы вы могли писать код вместе с нами

Первое, что мы сделаем, – создадим новый файл назовём lr_1d.py. У нас уже есть скрипт, который сгенерировал все необходимые данные, так что нам не понадобится их искать, достаточно просто взять csv-файл из репозитория github.

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

NumPy — это библиотека с открытым исходным кодом для языка программирования Python. Возможности:

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

Matplotlib — библиотека на языке программирования Python для визуализации данных двумерной (2D) графикой (3D графика также поддерживается). Получаемые изображения могут быть использованы в качестве иллюстраций в публикациях

import numpy as np

import matplotlib.pyplot as plt

Далее мы загружаем данные. Создадим два пустых массива, или списка. Вы можете использовать csv-ридер, но мы предпочитаем сделать всё проще. Поэтому сначала мы преобразуем данные x в вещественные числа и перемещаем в массив X. То же самое проделываем с y. Для удобства  превращаем отступы в коде в пробелы, но это сугубо на ваше усмотрение.

X=[]

Y=[]

for line in open(‘data_1d.csv’):

x, y = line.split(‘,’)

X.append(float(x))

Y.append(float(y))

Теперь преобразуем X и Y в массивы NumPy – так будет удобнее с ними работать.

X=np.array(X)

Y=np.array(Y)

А теперь отобразим наши данные в графическом виде, чтобы понять, как это выглядит, и запускаем.

plt.scatter(X,Y)

plt.show()

При первом прогоне прорисовка может занять некоторое время. Итак, мы видим, что точки располагаются более-менее вдоль прямой.

Теперь используем полученные нами ранее уравнения для вычисления коэффициентов a и b. Не забывайте, что в найденных выражениях коэффициенты имеют общий знаменатель, поэтому нет необходимости вычислять его дважды – мы просто присвоим его значение переменной denominator. Сначала мы вычислим сумму . Конечно, мы могли бы использовать для этого оператор цикла for, но такой способ неэффективен. Куда лучше использовать библиотеку NumPy, которая скрыто работает с массивами, а все вычисления происходят намного быстрее. Воспользуемся специальными функциями библиотеки.

denominator = X.dot(X) – X.mean() * X.sum()

И это всё. Теперь вы понимаете, почему так хорошо использовать библиотеку NumPy.

Теперь вычислим числители коэффициентов a и b.

a = ( X.dot(Y) – Y.mean()*X.sum() ) / denominator

b = ( Y.mean() * X.dot(X) – X.mean() * X.dot(Y) ) / denominator

Теперь вычислим наше предсказываемое Y.

Yhat = a*X + b

И отобразим всё это в графическом виде.

plt.scatter(X, Y)

plt.plot(X, Yhat)

plt.show()

Запустим ещё раз.

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

Определение качественности  модели  – R-квадрат в одномерной линейной регрессии

Предлагаю кратко  поговорить о том, как определить эффективность нашей модели.

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

Теперь мы хотим найти какой-то количественный показатель того, насколько хороша наша модель. Обычно для этого в линейной регрессии – и не только в линейной, но и других – используется коэффициент детерминации R2. Коэффициент детерминации R2 определяется следующим образом:

 R^2 = 1 - \frac{SSres}{SStot},

где , то есть сумма квадратов остатков – это просто сумма квадратов ошибок, а

 SStot = \sum_i (y_i - y)^2.

Итак, что же показывает коэффициент детерминации? Пусть у нас сумма квадратов остатков является величиной, близкой к нулю: SSres≈0. Тогда значение коэффициента детерминации будет близким к единице:

R2 ≈ 1 – 0 ≈ 1.

Если значение коэффициента детерминации равно единице, это значит, что у нас идеальная модель.

Коэффициента детерминации, равное нулю

А что означает значение коэффициента детерминации, равное нулю? Если R2=0, это значит, что  Это значит, что мы рассчитали просто среднее значение y, из чего, в свою очередь, следует, что наша модель получилась не очень удачной.

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

Коэффициент детерминации отрицательный

А если значение коэффициента детерминации оказалось отрицательным?

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

Коэффициент детерминации в коде

Итак, сейчас мы напишем код для расчёта коэффициента детерминации R2.

Числитель включает выражение , и для этой разности мы введём специальную переменную. Помните, что наше Y является вектором, и наша переменная Yhat также является вектором размерности n, поэтому и их разность также является вектором. Но вычитая одно и то же число от каждого элемента вектора Y, мы получим скалярную величину, что очень удобно.

d1 = Y – Yhat

d2 = Y – Y.mean()

Коэффициент детерминации равен единице минус сумма квадратов остатков, делённая на сумму общих квадратов. Разность по каждому элементу массива у нас уже есть, и нам нужно возвести их в квадрат:

r2 = 1 – d1.dot(d1) / d2.dot(d2)

Даём команду на печать

print “the r-squared is:”, r2

И запускаем программу.

Коэффициент детерминации оказался равным 0,99, что очень близко к единице. Собственно, это и ожидалось, поскольку я специально подобрал данные для такого результата.

Демонстрация закона Мура в коде

Давайте исследуем закон Мура в коде языка Python.

Первое, с чего мы начнём – импорт библиотек. Начнём с Re, поскольку там есть одна занятная штучка. Также мы импортируем NumPy и Matplotlib.

import re

import numpy as np

import matplotlib.pyplot as plt

Далее мы создаём пустые массивы для X и Y. Кроме того, нам надо убрать дробную часть из наших входных данных. Если вы не знаете, как это сделать, не волнуйтесь, это не играет важной роли в нашем курсе.

X = [ ]

Y = [ ]

non_decimal = re.compile(r ‘ [^\d]+’)

for line in open (‘moore.csv’):

r = line.split(‘\t’)

x = int(non_decimal.sub(‘ ‘, r[2].split(‘[‘)[0]))

y = int(non_decimal.sub(‘ ‘, r[1].split(‘[‘)[0]))

X.append(x)

Y.append(y)

Не забывайте, что года – это наше x – находятся в третьем столбце таблицы, а число транзисторов – наше y – во втором столбце. Кроме того, преобразуем наши данные в np-массивы, поскольку так с ними будет легче работать.

X = np.array(X)

Y = np.array(Y)

И отобразим наши данные в графическом виде.

plt.scatter(X, Y)

plt.show()

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

Y = np.log(Y)

plt.scatter(X, Y)

plt.show()

Теперь мы можем заняться собственно нахождением решения для нашей линейной регрессии. Вычислим знаменатель, a, b и искомую прямую Ŷ.

denominator = X.dot(X) – X.mean() * X.sum()

a – ( X.dot(Y) – Y.mean() * X.sum() ) / denominator

b = ( Y.mean() * X.dot(X) – X.mean() * X.dot(Y) ) / denominator

Yhat = a*X + b

Вычислив искомую линейную регрессию, мы, конечно, захотим воочию увидеть результат.

plt.scatter(X, Y)

plt.plot (X, Yhat)

plt.(show)

Наша следующая задача – определить, насколько хороша наша линейная модель, вычислив коэффициент детерминации. Итак, давайте вычислим, насколько далеки мы от истинных значений. Результаты вычислений отобразим на экране.

d1 = Y – Yhat

d2 = Y – Y.mean()

r2 = 1 – d1.dot(d1) / d2.dot(d2)

print(‘’a:’’, a, ‘’b:’’, b)

print(‘’the r-squared is:’’, r2)

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

Итак, логарифм количества транзисторов равен нашему a, умноженному на количество лет, плюс b: ln(tc) = a*год + b – это, собственно, и есть найденная нами линейная регрессия. Следовательно, само количество транзисторов равно tc = eb*ea*год. Проведя соответствующие вычисления, можем вычислить время удвоения.

print(‘’time to double:’’, np.log(2)/a, ‘’years’’)

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

Мы видим также, что наш коэффициент детерминации равен 0,95. Это очень хороший показатель, означающий, что наша зависимость действительно линейная. Время удвоения равна приблизительно 1,97 года, что очень близко к двум годам. Следовательно, время удвоения количества транзисторов на микросхеме действительно равно примерно двум годам.

Понравилась статья? Поделить с друзьями:
Добавить комментарий

;-) :| :x :twisted: :smile: :shock: :sad: :roll: :razz: :oops: :o :mrgreen: :lol: :idea: :grin: :evil: :cry: :cool: :arrow: :???: :?: :!: