この記事では、単純 (対応) 回帰直線の数学的方程式を決定するいくつかの方法について説明します。
ここで説明する方程式を解くすべての方法は、最小二乗法に基づいています。 メソッドを次のように表します。
- 分析ソリューション
- 勾配降下法
- 確率的勾配降下法
この記事では、直線の方程式を解く方法ごとにさまざまな関数を提供します。これらの関数は、主にライブラリを使用せずに記述された関数に分けられます。 NumPy および計算に使用するもの NumPy。 上手に利用すると思われます NumPy コンピューティングコストが削減されます。
記事に記載されているすべてのコードは言語で書かれています python 2.7 を使用して ジュピターノート。 ソースコードとサンプルデータを含むファイルは以下に掲載されています。
この記事は、初心者と、人工知能の非常に幅広い分野である機械学習の研究をすでに徐々に習得し始めている人の両方を対象としています。
この内容を説明するために、非常に簡単な例を使用します。
条件例
私たちは依存を特徴づけるXNUMXつの価値観を持っています Y から X (表1):
表No.1「条件例」
値を次のように仮定します。 はその年の月であり、 — 今月の収益。 言い換えれば、収益はその年の月によって決まります。 - 収益が左右される唯一の兆候。
この例は、収益が年の月に条件付きで依存するという観点からも、値の数という観点からも、まあまあです。それらはほとんどありません。 しかし、このような単純化により、彼らが言うように、必ずしも簡単ではないが、初心者が理解できる内容を説明することが可能になります。 また、数値が単純であるため、多大な人件費をかけずに紙の上で例を解くことができます。
この例で示した依存関係は、次の形式の単純 (対) 回帰直線の数学方程式によって非常によく近似できると仮定します。
どこ 収益を受け取った月です。 — その月に対応する収益、 и は推定直線の回帰係数です。
係数に注意してください 多くの場合、推定された直線の傾きまたは勾配と呼ばれます。 の量を表します。 それが変わるとき .
明らかに、この例でのタスクは方程式内のそのような係数を選択することです。 и 、月ごとに計算された収益値の真の答えからの偏差、つまりサンプルに示されている値は最小限になります。
最小二乗法
最小二乗法によれば、偏差は二乗して計算されます。 この手法を使用すると、偏差の符号が反対の場合に、偏差の相互キャンセルを回避できます。 たとえば、あるケースで偏差が次のようになったとします。 +5 (プラス XNUMX)、もう XNUMX つは -5 (マイナス 0)、偏差の合計は互いに打ち消し合い、XNUMX (ゼロ) になります。 偏差を二乗するのではなく、係数の特性を使用することも可能です。そうすると、すべての偏差が正になり、累積されます。 この点については詳しく説明しませんが、計算の便宜上、偏差を XNUMX 乗するのが通例であることだけを示します。
これは、平方偏差 (誤差) の最小和を決定する式がどのようになるかを示しています。
どこ は真の答え (つまり、計算した収益) の近似関数です。
は本当の答え (サンプルで提供される収益)、
サンプルインデックス (偏差が決定される月の番号)
関数を微分し、偏微分方程式を定義して、解析的解決に進む準備をしましょう。 しかしその前に、微分とは何かについて簡単に説明し、導関数の幾何学的意味を思い出してみましょう。
差別化
微分は、関数の導関数を見つける操作です。
導関数は何に使用されますか? 関数の導関数は関数の変化率を特徴づけ、その方向を示します。 特定の点での導関数が正の場合、関数は増加し、そうでない場合、関数は減少します。 また、絶対導関数の値が大きくなるほど、関数値の変化率が高くなり、関数グラフの傾きも急になります。
たとえば、デカルト座標系の条件下では、点 M(0,0) での導関数の値は次と等しくなります。 +25 値がシフトされた特定の時点でのことを意味します 従来の単位、値で右へ 従来の25ユニット増加します。 グラフで見ると、値がかなり急激に上昇しているように見えます ある地点から。
もう一つの例。 微分値は等しい -0,1 位置がずれたときのことを意味します 従来単位XNUMX個あたりの値 従来のユニットのわずか0,1の減少です。 同時に、関数のグラフでは、かろうじて目立つ下向きの傾斜が観察できます。 山に例えると、非常に急な頂上に登らなければならなかった前の例とは異なり、山から緩やかな斜面を非常にゆっくりと下っているようなものです:)
したがって、関数を微分した後、 オッズで и 、1階の偏微分方程式を定義します。 方程式を決定した後、XNUMX つの方程式からなる系を受け取り、これを解くことで係数の値を選択できるようになります。 и 、指定された点での対応する導関数の値は、非常にわずかな量だけ変化しますが、解析解の場合はまったく変化しません。 言い換えれば、これらの点での偏導関数の値がゼロに等しくなるため、見つかった係数での誤差関数は最小値に達します。
したがって、微分の法則に従って、係数に関する 1 次の偏微分方程式は次のようになります。 次の形式になります:
に関する一次偏微分方程式 次の形式になります:
その結果、非常に単純な解析解を持つ連立方程式を受け取りました。
開始{式*}
開始{件}
na + bsumlimits_{i=1}^nx_i — sumlimits_{i=1}^ny_i = 0
sumlimits_{i=1}^nx_i(a +bsumlimits_{i=1}^nx_i — sumlimits_{i=1}^ny_i) = 0
終了{件}
終了{式*}
方程式を解く前に、プリロードを行い、ロードが正しいことを確認し、データをフォーマットしましょう。
データのロードとフォーマット
分析ソリューション、その後の勾配および確率的勾配降下法では、コードを XNUMX つのバリエーションで使用することに注意してください。ライブラリを使用する NumPy これを使用しない場合は、適切なデータ形式が必要になります (コードを参照)。
データのロードと処理のコード
# импортируем все нужные нам библиотеки
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import math
import pylab as pl
import random
# графики отобразим в Jupyter
%matplotlib inline
# укажем размер графиков
from pylab import rcParams
rcParams['figure.figsize'] = 12, 6
# отключим предупреждения Anaconda
import warnings
warnings.simplefilter('ignore')
# загрузим значения
table_zero = pd.read_csv('data_example.txt', header=0, sep='t')
# посмотрим информацию о таблице и на саму таблицу
print table_zero.info()
print '********************************************'
print table_zero
print '********************************************'
# подготовим данные без использования NumPy
x_us = []
[x_us.append(float(i)) for i in table_zero['x']]
print x_us
print type(x_us)
print '********************************************'
y_us = []
[y_us.append(float(i)) for i in table_zero['y']]
print y_us
print type(y_us)
print '********************************************'
# подготовим данные с использованием NumPy
x_np = table_zero[['x']].values
print x_np
print type(x_np)
print x_np.shape
print '********************************************'
y_np = table_zero[['y']].values
print y_np
print type(y_np)
print y_np.shape
print '********************************************'可視化
ここで、最初にデータをロードし、次にロードの正確性をチェックし、最後にデータをフォーマットした後、最初の視覚化を実行します。 これによく使用される方法は、 ペアプロット ライブラリー シーボーン。 この例では、数が限られているため、ライブラリを使用する意味がありません。 シーボーン。 通常の図書館を利用します matplotlib そして散布図を見てください。
散布図コード
print 'График №1 "Зависимость выручки от месяца года"'
plt.plot(x_us,y_us,'o',color='green',markersize=16)
plt.xlabel('$Months$', size=16)
plt.ylabel('$Sales$', size=16)
plt.show()グラフNo.1「収益の月依存性」
分析ソリューション
最も一般的なツールを使用しましょう パイソン そして連立方程式を解きます。
開始{式*}
開始{件}
na + bsumlimits_{i=1}^nx_i — sumlimits_{i=1}^ny_i = 0
sumlimits_{i=1}^nx_i(a +bsumlimits_{i=1}^nx_i — sumlimits_{i=1}^ny_i) = 0
終了{件}
終了{式*}
クレイマーの法則によると 一般行列式と次の行列式を求めます。 とによって 、その後、行列式を次の値で割ります。 一般行列式への係数を求める 、同様に係数を求めます。 .
分析ソリューションコード
# определим функцию для расчета коэффициентов a и b по правилу Крамера
def Kramer_method (x,y):
# сумма значений (все месяца)
sx = sum(x)
# сумма истинных ответов (выручка за весь период)
sy = sum(y)
# сумма произведения значений на истинные ответы
list_xy = []
[list_xy.append(x[i]*y[i]) for i in range(len(x))]
sxy = sum(list_xy)
# сумма квадратов значений
list_x_sq = []
[list_x_sq.append(x[i]**2) for i in range(len(x))]
sx_sq = sum(list_x_sq)
# количество значений
n = len(x)
# общий определитель
det = sx_sq*n - sx*sx
# определитель по a
det_a = sx_sq*sy - sx*sxy
# искомый параметр a
a = (det_a / det)
# определитель по b
det_b = sxy*n - sy*sx
# искомый параметр b
b = (det_b / det)
# контрольные значения (прооверка)
check1 = (n*b + a*sx - sy)
check2 = (b*sx + a*sx_sq - sxy)
return [round(a,4), round(b,4)]
# запустим функцию и запишем правильные ответы
ab_us = Kramer_method(x_us,y_us)
a_us = ab_us[0]
b_us = ab_us[1]
print ' 33[1m' + ' 33[4m' + "Оптимальные значения коэффициентов a и b:" + ' 33[0m'
print 'a =', a_us
print 'b =', b_us
print
# определим функцию для подсчета суммы квадратов ошибок
def errors_sq_Kramer_method(answers,x,y):
list_errors_sq = []
for i in range(len(x)):
err = (answers[0] + answers[1]*x[i] - y[i])**2
list_errors_sq.append(err)
return sum(list_errors_sq)
# запустим функцию и запишем значение ошибки
error_sq = errors_sq_Kramer_method(ab_us,x_us,y_us)
print ' 33[1m' + ' 33[4m' + "Сумма квадратов отклонений" + ' 33[0m'
print error_sq
print
# замерим время расчета
# print ' 33[1m' + ' 33[4m' + "Время выполнения расчета суммы квадратов отклонений:" + ' 33[0m'
# % timeit error_sq = errors_sq_Kramer_method(ab,x_us,y_us)得られたものは次のとおりです。
したがって、係数の値が見つかり、偏差の二乗の合計が確立されました。 求められた係数に従って、散乱ヒストグラム上に直線を描いてみましょう。
回帰直線コード
# определим функцию для формирования массива рассчетных значений выручки
def sales_count(ab,x,y):
line_answers = []
[line_answers.append(ab[0]+ab[1]*x[i]) for i in range(len(x))]
return line_answers
# построим графики
print 'Грфик№2 "Правильные и расчетные ответы"'
plt.plot(x_us,y_us,'o',color='green',markersize=16, label = '$True$ $answers$')
plt.plot(x_us, sales_count(ab_us,x_us,y_us), color='red',lw=4,
label='$Function: a + bx,$ $where$ $a='+str(round(ab_us[0],2))+',$ $b='+str(round(ab_us[1],2))+'$')
plt.xlabel('$Months$', size=16)
plt.ylabel('$Sales$', size=16)
plt.legend(loc=1, prop={'size': 16})
plt.show()図表No.2「計算された正解」
月ごとの偏差グラフがご覧いただけます。 私たちの場合、そこから重要な実用的な価値は得られませんが、単純な線形回帰方程式が収益の月への依存性をどの程度適切に特徴付けるかについての好奇心を満たすことができます。
偏差チャートコード
# определим функцию для формирования массива отклонений в процентах
def error_per_month(ab,x,y):
sales_c = sales_count(ab,x,y)
errors_percent = []
for i in range(len(x)):
errors_percent.append(100*(sales_c[i]-y[i])/y[i])
return errors_percent
# построим график
print 'График№3 "Отклонения по-месячно, %"'
plt.gca().bar(x_us, error_per_month(ab_us,x_us,y_us), color='brown')
plt.xlabel('Months', size=16)
plt.ylabel('Calculation error, %', size=16)
plt.show()チャートNo.3「偏差、%」
完璧ではありませんが、タスクは完了しました。
係数を決定する関数を書いてみましょう и 図書館を利用する NumPyより正確には、XNUMX つの関数を作成します。XNUMX つは擬似逆行列を使用し (プロセスは計算的に複雑で不安定であるため、実際には推奨されません)、もう XNUMX つは行列方程式を使用します。
分析ソリューション コード (NumPy)
# для начала добавим столбец с не изменяющимся значением в 1.
# Данный столбец нужен для того, чтобы не обрабатывать отдельно коэффицент a
vector_1 = np.ones((x_np.shape[0],1))
x_np = table_zero[['x']].values # на всякий случай приведем в первичный формат вектор x_np
x_np = np.hstack((vector_1,x_np))
# проверим то, что все сделали правильно
print vector_1[0:3]
print x_np[0:3]
print '***************************************'
print
# напишем функцию, которая определяет значения коэффициентов a и b с использованием псевдообратной матрицы
def pseudoinverse_matrix(X, y):
# задаем явный формат матрицы признаков
X = np.matrix(X)
# определяем транспонированную матрицу
XT = X.T
# определяем квадратную матрицу
XTX = XT*X
# определяем псевдообратную матрицу
inv = np.linalg.pinv(XTX)
# задаем явный формат матрицы ответов
y = np.matrix(y)
# находим вектор весов
return (inv*XT)*y
# запустим функцию
ab_np = pseudoinverse_matrix(x_np, y_np)
print ab_np
print '***************************************'
print
# напишем функцию, которая использует для решения матричное уравнение
def matrix_equation(X,y):
a = np.dot(X.T, X)
b = np.dot(X.T, y)
return np.linalg.solve(a, b)
# запустим функцию
ab_np = matrix_equation(x_np,y_np)
print ab_np係数の決定に費やした時間を比較してみましょう и 、提示された 3 つの方法に従って。
計算時間を計算するコード
print ' 33[1m' + ' 33[4m' + "Время выполнения расчета коэффициентов без использования библиотеки NumPy:" + ' 33[0m'
% timeit ab_us = Kramer_method(x_us,y_us)
print '***************************************'
print
print ' 33[1m' + ' 33[4m' + "Время выполнения расчета коэффициентов с использованием псевдообратной матрицы:" + ' 33[0m'
%timeit ab_np = pseudoinverse_matrix(x_np, y_np)
print '***************************************'
print
print ' 33[1m' + ' 33[4m' + "Время выполнения расчета коэффициентов с использованием матричного уравнения:" + ' 33[0m'
%timeit ab_np = matrix_equation(x_np, y_np)
少量のデータでは、Cramer の方法を使用して係数を見つける「自己作成」関数が先に登場します。
ここで、係数を見つける他の方法に進むことができます。 и .
勾配降下法
まず、グラデーションとは何かを定義しましょう。 簡単に言えば、勾配は関数の最大成長の方向を示すセグメントです。 山に登ることに例えると、勾配面が山の頂上まで最も急な登りとなる場所です。 山での例を展開すると、実際には、できるだけ早く低地、つまり機能が増加または減少しない最小の場所に到達するには、最も急な下り坂が必要であることを思い出します。 この時点で、導関数はゼロに等しくなります。 したがって、勾配ではなく反勾配が必要です。 反勾配を見つけるには、勾配に次の値を掛けるだけです。 -1 (マイナス XNUMX)。
関数には複数の最小値が存在する可能性があり、以下に提案するアルゴリズムを使用してそのうちの XNUMX つに到達すると、見つかった最小値よりも低い別の最小値を見つけることはできなくなるという事実に注意してください。 リラックスしましょう、これは私たちに対する脅威ではありません! 私たちの場合、私たちの関数は単一の最小値を扱っています。 グラフ上は規則的な放物線です。 そして、学校の数学の授業で私たち全員がよく知っているように、放物線には最小値が XNUMX つだけあります。
勾配が必要な理由と、勾配がセグメント、つまり、まったく同じ係数を持つ指定された座標を持つベクトルであることがわかった後、 и 勾配降下法を実装できます。
始める前に、降下アルゴリズムに関するいくつかの文を読んでおくことをお勧めします。
- 擬似ランダムな方法で係数の座標を決定します。 и 。 この例では、ゼロに近い係数を決定します。 これは一般的な方法ですが、場合によっては独自の方法がある場合があります。
- 座標から 点における 1 次偏導関数の値を減算します。 。 したがって、導関数が正の場合、関数は増加します。 したがって、微分の値を差し引くと、成長とは逆の方向、つまり下降方向に進むことになります。 導関数が負の場合、この時点での関数は減少し、導関数の値を減算することで下降方向に移動します。
- 座標に対しても同様の操作を実行します。 : 点における偏導関数の値を減算します。 .
- 最小値を飛び越えて深宇宙に飛び出さないようにするには、降下方向のステップサイズを設定する必要があります。 一般に、計算コストを削減するために、ステップを正しく設定する方法と、降下プロセス中にステップを変更する方法について、記事全体を書くことができます。 しかし今、私たちは少し異なるタスクを前にしています。「つつく」という科学的方法を使用して、または一般的な用語で言うところの、経験的にステップサイズを確立します。
- 指定された座標から出発したら и 導関数の値を減算すると、新しい座標が得られます и 。 すでに計算された座標から次のステップ (減算) を実行します。 そして、必要な収束が達成されるまで、このサイクルが何度も繰り返されます。
全て！ さあ、マリアナ海溝の最も深い峡谷を探しに行く準備が整いました。 始めましょう。
勾配降下法のコード
# напишем функцию градиентного спуска без использования библиотеки NumPy.
# Функция на вход принимает диапазоны значений x,y, длину шага (по умолчанию=0,1), допустимую погрешность(tolerance)
def gradient_descent_usual(x_us,y_us,l=0.1,tolerance=0.000000000001):
# сумма значений (все месяца)
sx = sum(x_us)
# сумма истинных ответов (выручка за весь период)
sy = sum(y_us)
# сумма произведения значений на истинные ответы
list_xy = []
[list_xy.append(x_us[i]*y_us[i]) for i in range(len(x_us))]
sxy = sum(list_xy)
# сумма квадратов значений
list_x_sq = []
[list_x_sq.append(x_us[i]**2) for i in range(len(x_us))]
sx_sq = sum(list_x_sq)
# количество значений
num = len(x_us)
# начальные значения коэффициентов, определенные псевдослучайным образом
a = float(random.uniform(-0.5, 0.5))
b = float(random.uniform(-0.5, 0.5))
# создаем массив с ошибками, для старта используем значения 1 и 0
# после завершения спуска стартовые значения удалим
errors = [1,0]
# запускаем цикл спуска
# цикл работает до тех пор, пока отклонение последней ошибки суммы квадратов от предыдущей, не будет меньше tolerance
while abs(errors[-1]-errors[-2]) > tolerance:
a_step = a - l*(num*a + b*sx - sy)/num
b_step = b - l*(a*sx + b*sx_sq - sxy)/num
a = a_step
b = b_step
ab = [a,b]
errors.append(errors_sq_Kramer_method(ab,x_us,y_us))
return (ab),(errors[2:])
# запишем массив значений
list_parametres_gradient_descence = gradient_descent_usual(x_us,y_us,l=0.1,tolerance=0.000000000001)
print ' 33[1m' + ' 33[4m' + "Значения коэффициентов a и b:" + ' 33[0m'
print 'a =', round(list_parametres_gradient_descence[0][0],3)
print 'b =', round(list_parametres_gradient_descence[0][1],3)
print
print ' 33[1m' + ' 33[4m' + "Сумма квадратов отклонений:" + ' 33[0m'
print round(list_parametres_gradient_descence[1][-1],3)
print
print ' 33[1m' + ' 33[4m' + "Количество итераций в градиентном спуске:" + ' 33[0m'
print len(list_parametres_gradient_descence[1])
print
私たちはマリアナ海溝の一番底まで潜ったところ、すべて同じ係数値が見つかりました。 и 、まさに予想通りでした。
もう一度潜ってみましょう。今度だけ、私たちの深海船には他のテクノロジー、つまりライブラリが搭載されます。 NumPy.
勾配降下法 (NumPy) のコード
# перед тем определить функцию для градиентного спуска с использованием библиотеки NumPy,
# напишем функцию определения суммы квадратов отклонений также с использованием NumPy
def error_square_numpy(ab,x_np,y_np):
y_pred = np.dot(x_np,ab)
error = y_pred - y_np
return sum((error)**2)
# напишем функцию градиентного спуска с использованием библиотеки NumPy.
# Функция на вход принимает диапазоны значений x,y, длину шага (по умолчанию=0,1), допустимую погрешность(tolerance)
def gradient_descent_numpy(x_np,y_np,l=0.1,tolerance=0.000000000001):
# сумма значений (все месяца)
sx = float(sum(x_np[:,1]))
# сумма истинных ответов (выручка за весь период)
sy = float(sum(y_np))
# сумма произведения значений на истинные ответы
sxy = x_np*y_np
sxy = float(sum(sxy[:,1]))
# сумма квадратов значений
sx_sq = float(sum(x_np[:,1]**2))
# количество значений
num = float(x_np.shape[0])
# начальные значения коэффициентов, определенные псевдослучайным образом
a = float(random.uniform(-0.5, 0.5))
b = float(random.uniform(-0.5, 0.5))
# создаем массив с ошибками, для старта используем значения 1 и 0
# после завершения спуска стартовые значения удалим
errors = [1,0]
# запускаем цикл спуска
# цикл работает до тех пор, пока отклонение последней ошибки суммы квадратов от предыдущей, не будет меньше tolerance
while abs(errors[-1]-errors[-2]) > tolerance:
a_step = a - l*(num*a + b*sx - sy)/num
b_step = b - l*(a*sx + b*sx_sq - sxy)/num
a = a_step
b = b_step
ab = np.array([[a],[b]])
errors.append(error_square_numpy(ab,x_np,y_np))
return (ab),(errors[2:])
# запишем массив значений
list_parametres_gradient_descence = gradient_descent_numpy(x_np,y_np,l=0.1,tolerance=0.000000000001)
print ' 33[1m' + ' 33[4m' + "Значения коэффициентов a и b:" + ' 33[0m'
print 'a =', round(list_parametres_gradient_descence[0][0],3)
print 'b =', round(list_parametres_gradient_descence[0][1],3)
print
print ' 33[1m' + ' 33[4m' + "Сумма квадратов отклонений:" + ' 33[0m'
print round(list_parametres_gradient_descence[1][-1],3)
print
print ' 33[1m' + ' 33[4m' + "Количество итераций в градиентном спуске:" + ' 33[0m'
print len(list_parametres_gradient_descence[1])
print
係数値 и 変わらない。
勾配降下中に誤差がどのように変化したか、つまり偏差二乗和が各ステップでどのように変化したかを見てみましょう。
偏差二乗和をプロットするコード
print 'График№4 "Сумма квадратов отклонений по-шагово"'
plt.plot(range(len(list_parametres_gradient_descence[1])), list_parametres_gradient_descence[1], color='red', lw=3)
plt.xlabel('Steps (Iteration)', size=16)
plt.ylabel('Sum of squared deviations', size=16)
plt.show()グラフNo.4「勾配降下時の偏差の二乗和」
グラフでは、ステップごとに誤差が減少し、一定回数の反復の後、ほぼ水平な線が観察されることがわかります。
最後に、コードの実行時間の違いを見積もってみましょう。
勾配降下の計算時間を決定するコード
print ' 33[1m' + ' 33[4m' + "Время выполнения градиентного спуска без использования библиотеки NumPy:" + ' 33[0m'
%timeit list_parametres_gradient_descence = gradient_descent_usual(x_us,y_us,l=0.1,tolerance=0.000000000001)
print '***************************************'
print
print ' 33[1m' + ' 33[4m' + "Время выполнения градиентного спуска с использованием библиотеки NumPy:" + ' 33[0m'
%timeit list_parametres_gradient_descence = gradient_descent_numpy(x_np,y_np,l=0.1,tolerance=0.000000000001)
おそらく何か間違ったことをしているのかもしれませんが、これもライブラリを使用しない単純な「自作」関数です。 NumPy ライブラリを使用した関数の計算時間を上回る NumPy.
しかし、私たちは立ち止まっているわけではなく、単純な線形回帰方程式を解くための別の刺激的な方法の研究に向かって進んでいます。 会う！
確率的勾配降下法
確率的勾配降下法の動作原理をすぐに理解するには、通常の勾配降下法との違いを理解するのが良いでしょう。 勾配降下の場合、次の導関数の方程式では、 и サンプル内で利用可能なすべての特徴と真の回答の値の合計 (つまり、すべての特徴の合計) を使用しました。 и ）。 確率的勾配降下法では、サンプルに存在するすべての値を使用するのではなく、いわゆるサンプル インデックスを擬似ランダムに選択し、その値を使用します。
たとえば、インデックスが番号 3 (XNUMX) であると判断された場合、次の値を取得します。 и 、次に、その値を微分方程式に代入して、新しい座標を決定します。 次に、座標を決定したら、再度擬似ランダムにサンプルのインデックスを決定し、そのインデックスに対応する値を偏微分方程式に代入し、新しい方法で座標を決定します。 и 等コンバージェンスが緑色に変わるまで。 一見すると、これはまったく機能しないように思えるかもしれませんが、実際には機能します。 確かに、ステップごとに誤差が減少するわけではありませんが、傾向があることは確かです。
確率的勾配降下法の従来のものと比較した利点は何ですか? サンプル サイズが非常に大きく、数万の値で測定される場合は、サンプル全体を処理するよりも、たとえばランダムに数千個を処理する方がはるかに簡単です。 ここで確率的勾配降下法が登場します。 もちろん、私たちの場合、大きな違いはわかりません。
コードを見てみましょう。
確率的勾配降下法のコード
# определим функцию стох.град.шага
def stoch_grad_step_usual(vector_init, x_us, ind, y_us, l):
# выбираем значение икс, которое соответствует случайному значению параметра ind
# (см.ф-цию stoch_grad_descent_usual)
x = x_us[ind]
# рассчитывыаем значение y (выручку), которая соответствует выбранному значению x
y_pred = vector_init[0] + vector_init[1]*x_us[ind]
# вычисляем ошибку расчетной выручки относительно представленной в выборке
error = y_pred - y_us[ind]
# определяем первую координату градиента ab
grad_a = error
# определяем вторую координату ab
grad_b = x_us[ind]*error
# вычисляем новый вектор коэффициентов
vector_new = [vector_init[0]-l*grad_a, vector_init[1]-l*grad_b]
return vector_new
# определим функцию стох.град.спуска
def stoch_grad_descent_usual(x_us, y_us, l=0.1, steps = 800):
# для самого начала работы функции зададим начальные значения коэффициентов
vector_init = [float(random.uniform(-0.5, 0.5)), float(random.uniform(-0.5, 0.5))]
errors = []
# запустим цикл спуска
# цикл расчитан на определенное количество шагов (steps)
for i in range(steps):
ind = random.choice(range(len(x_us)))
new_vector = stoch_grad_step_usual(vector_init, x_us, ind, y_us, l)
vector_init = new_vector
errors.append(errors_sq_Kramer_method(vector_init,x_us,y_us))
return (vector_init),(errors)
# запишем массив значений
list_parametres_stoch_gradient_descence = stoch_grad_descent_usual(x_us, y_us, l=0.1, steps = 800)
print ' 33[1m' + ' 33[4m' + "Значения коэффициентов a и b:" + ' 33[0m'
print 'a =', round(list_parametres_stoch_gradient_descence[0][0],3)
print 'b =', round(list_parametres_stoch_gradient_descence[0][1],3)
print
print ' 33[1m' + ' 33[4m' + "Сумма квадратов отклонений:" + ' 33[0m'
print round(list_parametres_stoch_gradient_descence[1][-1],3)
print
print ' 33[1m' + ' 33[4m' + "Количество итераций в стохастическом градиентном спуске:" + ' 33[0m'
print len(list_parametres_stoch_gradient_descence[1])
私たちは係数を注意深く見て、「どうしてこんなことが起こり得るのか?」という疑問を抱いていることに気づきます。 他の係数値も取得しました и 。 おそらく確率的勾配降下法により、方程式にとってより最適なパラメーターが見つかったのではないでしょうか? 残念だけど違う。 偏差の二乗和を見て、係数の新しい値では誤差が大きくなることを確認するだけで十分です。 私たちは急いで絶望するつもりはありません。 誤差の変化のグラフを作成してみましょう。
確率的勾配降下法で偏差の二乗和をプロットするコード
print 'График №5 "Сумма квадратов отклонений по-шагово"'
plt.plot(range(len(list_parametres_stoch_gradient_descence[1])), list_parametres_stoch_gradient_descence[1], color='red', lw=2)
plt.xlabel('Steps (Iteration)', size=16)
plt.ylabel('Sum of squared deviations', size=16)
plt.show()グラフNo.5「確率的勾配降下法における偏差の二乗和」
スケジュールを見ると、すべてが適切な位置に収まったので、これからすべてを修正します。
どうしたの？ 次のことが起こりました。 月をランダムに選択すると、その選択された月について、アルゴリズムが収益計算の誤差を削減しようとします。 次に、別の月を選択して計算を繰り返しますが、XNUMX 番目に選択した月の誤差は減少します。 ここで、最初の XNUMX か月は単純な線形回帰方程式の直線から大幅に逸脱していることを思い出してください。 これは、これら XNUMX つの月のいずれかが選択されると、それぞれの誤差を減らすことによって、アルゴリズムがサンプル全体の誤差を大幅に増加させることを意味します。 じゃあ何をすればいいの？ 答えは簡単です。下降ステップを減らす必要があるからです。 結局のところ、下降ステップを減らすことで、エラーが上下に「ジャンプ」することもなくなります。 というか、「ジャンプ」エラーは止まらないのですが、そんなにすぐには止まりません:) 確認してみましょう。
SGD をより小さい増分で実行するコード
# запустим функцию, уменьшив шаг в 100 раз и увеличив количество шагов соответсвующе
list_parametres_stoch_gradient_descence = stoch_grad_descent_usual(x_us, y_us, l=0.001, steps = 80000)
print ' 33[1m' + ' 33[4m' + "Значения коэффициентов a и b:" + ' 33[0m'
print 'a =', round(list_parametres_stoch_gradient_descence[0][0],3)
print 'b =', round(list_parametres_stoch_gradient_descence[0][1],3)
print
print ' 33[1m' + ' 33[4m' + "Сумма квадратов отклонений:" + ' 33[0m'
print round(list_parametres_stoch_gradient_descence[1][-1],3)
print
print ' 33[1m' + ' 33[4m' + "Количество итераций в стохастическом градиентном спуске:" + ' 33[0m'
print len(list_parametres_stoch_gradient_descence[1])
print 'График №6 "Сумма квадратов отклонений по-шагово"'
plt.plot(range(len(list_parametres_stoch_gradient_descence[1])), list_parametres_stoch_gradient_descence[1], color='red', lw=2)
plt.xlabel('Steps (Iteration)', size=16)
plt.ylabel('Sum of squared deviations', size=16)
plt.show()
グラフNo.6「確率的勾配降下法（80万ステップ）における偏差の二乗和」
係数は改善されましたが、まだ理想的ではありません。 仮に、これはこの方法で修正できます。 たとえば、最後の 1000 回の反復で最小誤差が生じた係数の値を選択します。 確かに、このためには、係数自体の値も書き留める必要があります。 私たちはこれを行わず、むしろスケジュールに注意を払います。 見た目はスムーズで誤差も均等に減っているように見えます。 実はこれは真実ではありません。 最初の 1000 回の反復を見て、最後の反復と比較してみましょう。
SGD チャートのコード (最初の 1000 ステップ)
print 'График №7 "Сумма квадратов отклонений по-шагово. Первые 1000 итераций"'
plt.plot(range(len(list_parametres_stoch_gradient_descence[1][:1000])),
list_parametres_stoch_gradient_descence[1][:1000], color='red', lw=2)
plt.xlabel('Steps (Iteration)', size=16)
plt.ylabel('Sum of squared deviations', size=16)
plt.show()
print 'График №7 "Сумма квадратов отклонений по-шагово. Последние 1000 итераций"'
plt.plot(range(len(list_parametres_stoch_gradient_descence[1][-1000:])),
list_parametres_stoch_gradient_descence[1][-1000:], color='red', lw=2)
plt.xlabel('Steps (Iteration)', size=16)
plt.ylabel('Sum of squared deviations', size=16)
plt.show()グラフNo.7「偏差二乗和SGD（最初の1000ステップ）」
グラフNo.8「偏差二乗和SGD（最後の1000ステップ）」
降下の最初の段階では、誤差がかなり均一かつ急激に減少していることが観察されます。 最後の反復では、誤差が 1,475 の値を回って、ある瞬間にはこの最適値と等しくなることがわかりますが、それでも誤差は増加します... 繰り返しますが、次の値を書き留めることができます。係数 и を選択し、誤差が最小のものを選択します。 しかし、より深刻な問題がありました。最適値に近い値を得るには 80 万ステップ (コードを参照) を実行する必要がありました。 そして、これはすでに、勾配降下法と比較して確率的勾配降下法で計算時間を節約するという考えと矛盾しています。 何を修正および改善できるでしょうか? 最初の反復では自信を持って下降していることに気づくのは難しくありません。したがって、最初の反復では大きなステップを残し、先に進むにつれてステップを減らす必要があります。 この記事ではこれについては説明しません。すでに長すぎます。 望む人は、これを行う方法を自分で考えることができます。それは難しいことではありません:)
次に、ライブラリを使用して確率的勾配降下法を実行してみましょう NumPy (そして、前に見つけた石につまずかないようにしましょう)
確率的勾配降下法 (NumPy) のコード
# для начала напишем функцию градиентного шага
def stoch_grad_step_numpy(vector_init, X, ind, y, l):
x = X[ind]
y_pred = np.dot(x,vector_init)
err = y_pred - y[ind]
grad_a = err
grad_b = x[1]*err
return vector_init - l*np.array([grad_a, grad_b])
# определим функцию стохастического градиентного спуска
def stoch_grad_descent_numpy(X, y, l=0.1, steps = 800):
vector_init = np.array([[np.random.randint(X.shape[0])], [np.random.randint(X.shape[0])]])
errors = []
for i in range(steps):
ind = np.random.randint(X.shape[0])
new_vector = stoch_grad_step_numpy(vector_init, X, ind, y, l)
vector_init = new_vector
errors.append(error_square_numpy(vector_init,X,y))
return (vector_init), (errors)
# запишем массив значений
list_parametres_stoch_gradient_descence = stoch_grad_descent_numpy(x_np, y_np, l=0.001, steps = 80000)
print ' 33[1m' + ' 33[4m' + "Значения коэффициентов a и b:" + ' 33[0m'
print 'a =', round(list_parametres_stoch_gradient_descence[0][0],3)
print 'b =', round(list_parametres_stoch_gradient_descence[0][1],3)
print
print ' 33[1m' + ' 33[4m' + "Сумма квадратов отклонений:" + ' 33[0m'
print round(list_parametres_stoch_gradient_descence[1][-1],3)
print
print ' 33[1m' + ' 33[4m' + "Количество итераций в стохастическом градиентном спуске:" + ' 33[0m'
print len(list_parametres_stoch_gradient_descence[1])
print
使用せずに下降した場合とほぼ同じ値であることが判明 NumPy。 ただし、これは論理的です。
確率的勾配降下法にどれくらいの時間がかかったかを調べてみましょう。
SGD 計算時間を決定するコード (80 ステップ)
print ' 33[1m' + ' 33[4m' +
"Время выполнения стохастического градиентного спуска без использования библиотеки NumPy:"
+ ' 33[0m'
%timeit list_parametres_stoch_gradient_descence = stoch_grad_descent_usual(x_us, y_us, l=0.001, steps = 80000)
print '***************************************'
print
print ' 33[1m' + ' 33[4m' +
"Время выполнения стохастического градиентного спуска с использованием библиотеки NumPy:"
+ ' 33[0m'
%timeit list_parametres_stoch_gradient_descence = stoch_grad_descent_numpy(x_np, y_np, l=0.001, steps = 80000)
森の奥に行くほど、雲は暗くなります。やはり、「自分で作成した」式が最良の結果を示します。 これらすべては、ライブラリを使用するさらに巧妙な方法があるに違いないことを示唆しています。 NumPyこれにより、計算処理が大幅に高速化されます。 この記事ではそれらについては学びません。 暇なときに考えることがあるでしょう:)
要約しよう
要約する前に、親愛なる読者からおそらく寄せられたであろう質問に答えたいと思います。 実際、なぜそのような下り坂の「拷問」が必要なのか、貴重な低地を見つけるためになぜ山を上り下り（ほとんどは下山）しなければならないのか、もし私たちがそのような強力で単純な装置を手に持っているなら、分析ソリューションの形式で、私たちを正しい場所に瞬時にテレポートします。
この質問に対する答えは表面にあります。 ここでは非常に単純な例を見ていきました。この例の本当の答えは次のとおりです。 一つのサインに依存する 。 人生ではこれを頻繁に見ることはありません。ですから、2、30、50、あるいはそれ以上のサインがあると想像してみましょう。 これに、各属性の値を数千、さらには数万個加えてみましょう。 この場合、分析溶液は試験に耐えられず不合格になる可能性があります。 次に、勾配降下法とそのバリエーションは、ゆっくりと、しかし確実に目標、つまり関数の最小値に近づけます。 速度については心配しないでください。おそらく、歩幅 (つまり速度) を設定および調整できる方法を検討することになります。
そして、実際の簡単な要約です。
まず、この記事で紹介されている内容が、初心者の「データ サイエンティスト」にとって、単純な (それだけではない) 線形回帰方程式の解き方を理解するのに役立つことを願っています。
次に、方程式を解くいくつかの方法を検討しました。 状況に応じて、問題の解決に最適なものを選択できるようになりました。
XNUMX 番目に、追加の設定、つまり勾配降下ステップ長の威力を確認しました。 このパラメータは無視できません。 上で述べたように、計算コストを削減するには、下降中に歩幅を変更する必要があります。
第 XNUMX に、私たちの場合、「自分で作成した」関数が計算に最適な結果を示しました。 これはおそらく、ライブラリの機能があまり専門的に使用されていないことが原因です。 NumPy。 しかし、それはともかく、次の結論はそれ自体を示唆しています。 一方では、確立された意見を疑う価値がある場合もありますが、他方では、すべてを複雑にする価値があるとは限りません。逆に、問題を解決するより単純な方法の方が効果的である場合もあります。 そして、私たちの目標は、単純な線形回帰方程式を解くための XNUMX つのアプローチを分析することであったため、「自作」関数を使用するだけで十分でした。
文学（またはそのようなもの）
1. 線形回帰
2. 最小二乗法
3. デリバティブ
4。 グラジエント
5. 勾配降下法
6. NumPyライブラリ
出所： habr.com