Введение в машинное обучение.
Вводятся базовые элементы пакета линейной алгебры (Numpy) относфщейся к питону (Python версии 3.xx). Конкретно речь идет о алгебре и, соотвествено, линейной регресии. Последнее, в частности, используется для ввода ключевых понятий из машинного обучения: выборка, обучающее множество, тестовое/валидационное множество, поиск "оптимального" решения и регуляризация. Уровни значимости (F-статистика) при построении регрессий (statmodels).
Это предварительная версия! Любые замечания приветствуются.
Отрисовка кривой
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
# Для начала случайно сгенерируем точки в которых будет вычислена функция.
np.random.uniform(10, 20, 5) # Создаем матрицу равномерных распределений. Первые два числа казывают диапазон.
x = np.random.uniform( 3, 10, 10 ) # 10 точек на отрезке [3, 10]
x
x.shape # Это действительно массив из 10 чисел.
y = x * 2 - 15 # Вычисляем значение в каждой из точек. Функция линейная.
y
plt.plot( x, y )
yy = y + np.random.randn(10) # Добавим шум к каждому значению (измерению)
plt.plot( x, yy ) # Кривая как-то странно выглядит.
xx = np.sort( x ) # Дело в том, что необходимо отсортировать числа по оси x.
yy = xx * 2 - 15 + np.random.randn(10) # Теперь заново вычисляем.
yy.shape, x.shape # Массивы соответсвуют друг другу по размеру.
plt.plot( xx, yy )
Линейные уравнения.
# Решаем систему уравнений.
a = np.array( [ [2, -1], [1, 3] ] ) # 2x -y = 9
b = np.array( [ 9, 8] ) # 1x +3y = 8
x = np.linalg.solve(a, b)
x
np.allclose(np.dot(a, x), b) # Проверяем на равенство с учетом машинного эпсилон, т.е. числа могуть чуть отличаться.
5 - np.sqrt( 5 ) ** 2 # Например, должны были получить 0, но жизнь иначе распорядилась.
# Можно вычислить и обратную матрицу в явном виде
# a = np.array( [ [2, -1], [1, 3] ] )
ainv = np.linalg.inv( a )
ainv
# Решаем систему. Для этого используем произведение вектора на матрицу:
ainv.dot( b ) # Получили тот же ответ, что и раньше.
# Можно умножить и на матрицу. Тогда произведение матрицы на её обратное должно дать единичную.
np.allclose( np.dot(a, ainv), np.eye(2) )
Линейная регрессия
yy = xx * 2 - 15 + np.random.randn(10)/2 # Тоже самое но с меньшим шумом.
plt.plot( xx, yy, 'g.--') # Цвет зеленый, соединения пунктиром, точки маленькие.
q=xx.reshape(10,1) # Превращяем массив в вектор столбец, т.е. массив массивов (из одного элемента).
q
np.ones(10).reshape(10,1) # Создаем вектор столбец единиц.
A = np.concatenate( [np.ones((10,1)), q ], axis=1 ) # Создаем общую матрицу. Присваеваем её переменной A.
A
xx
AA = np.linalg.pinv(A) # Вычисляем псевдо обратную матрицу от A. Метод решения обычных систем не годится.
AA
d = AA.dot( yy.reshape(10,1) ) # Теперь её используем для решения системы.
d # В общем параметры приблизительно правельные.
p = np.linspace(3, 10, 100) # Отрисуем получившуюся прямую. Используем для этого 100 точек.
yyy = p*d[1] + d[0] # Параметры вычисленной линии хранятся в переменной d.
plt.plot( xx, yy, 'g.-', p, yyy, 'r-')
# Какая точность этого приближения?
xx[0]*d[1] + d[0], yy[0] # В даннойт точк значения вычисленные и табулированые:
(xx[0]*d[1] + d[0])[0], yy[0]
for i in range(xx.shape[0]): # Цикл пробегает по всем нашим точкам
print ( xx[i]*d[1] + d[0], yy[i] ) # Выводит два значения: вычисленное и табулированное.
xx.shape[0], len(xx)
def evalLine(d, x): # Создадим функцию, которая по параметрам прямой (d) считаем значение в точке (x).
return x * d[1] + d[0]
evalLine( d, xx[0]) # Вычислим значение в точке.
# Значение согласуется с предыдущими результатами.
# Посчитаем среднею ошибку по всем точкам.
sum1 = 0. # Линейная ошибка.
sum2 = 0. # Квадратичная ошибка.
for i in range(xx.shape[0]):
err = evalLine( d, xx[i]) - yy[i]
sum1 += err
sum2 += err * err
print ( xx[i]*d[1] + d[0], err )
print( "integral sum = ", sum1/xx.shape[0], sum2/xx.shape[0] )
ДЗ1 Почему сумма ошибка почти равна нулю?
# Обособим данный подсчет в функцию.
def errLine(d, x, y): # По параметрам прямой (d) и точкам (x,y) вычислим общую ошибку.
sum1 = 0.
sum2 = 0.
for i in range(x.shape[0]):
err = evalLine( d, x[i]) - y[i]
sum1 += err
sum2 += err * err
return sum1/x.shape[0], sum2/x.shape[0]
errLine(d, xx, yy)
# Обособим код вычисляющий линейную регрессию в функцию.
def fitLin(x, y): # x, y задают данные. На выходе параметры линейной регрессии.
sz = x.shape[0]
A = np.concatenate( [np.ones((sz,1)), x.reshape( sz, 1) ], axis=1 )
Ai=np.linalg.pinv(A)
return Ai.dot( y.reshape(sz,1) )
d0 = fitLin( xx, yy) # Проверим, что соответсвует редыдущим результатам.
d0
errLine(d0, xx, yy)
Learn and Test
perm = np.random.permutation( xx.shape[0] ) # Создаем перестановку числе от 0 до xx.shape[0] не включительно.
perm
#np.random.randint(0, xx.shape[0], 5)
ii = perm[:5] # Берем первые 5 чисел.
ii
d1 = fitLin( xx[ii], yy[ii]) # Вычисляем параметры линейной регрессии для подмножества числе (ii)
d1
ee = errLine(d1, xx[ii], yy[ii]) # Вычисляем ошибку на подмножестве.
ee
errLine(d1, xx, yy) # Вычисляем ошибку на всем множестве точек.
def experimentLin(n, d, p, k ):
x = np.random.uniform( 3, 10, (n) )
y = d[1] * x + d[0] + np.random.randn( x.shape[0] )/2
d0 = fitLin( x, y)
e0 = errLine(d0, x, y)
ee10 = np.array([0.])
ee11 = np.array([0.])
for j in range(k):
perm = np.random.permutation( xx.shape[0] )
ii = perm[:int(x.shape[0] * p)]#np.random.randint(0, x.shape[0], int(x.shape[0] * p) )
d1 = fitLin( x[ii], y[ii])
e1 = errLine(d1, x, y)
ee10 += e1[0]
ee11 += e1[1]
print( e0, [ee10 / k, ee11 / k] )
experimentLin(10, np.array([-15, 2.5]), 0.85, 1000)
ДЗ2: Построить статистику для ошибки. Вычислить её аналетически.
Polynomial
# Рассмотрим случай полиномов. Тонее на них мы изучим другое важное поянтие машинного обучения. Регулиризация.
import numpy.polynomial.polynomial as poly
poly.polyval(xx[0], d ) # Оказывается есть функция poly.polyval которая умеет вычилсять значение полинома.
xx.shape[0]
x.shape[0]
ee[0]
d = np.array( [-10, -28, 2.2] ) # Параметры задающие параболу.
n = 10
x = np.random.uniform( 3, 10, (n) )
x = np.sort( x )
x
y = poly.polyval(x, d ) + np.random.randn( x.shape[0] )/2
y
plt.plot( x, y, 'g*')
dd = poly.polyfit( x, y, x.shape[0] - 2 ) # Есть библиотечная функция для вычисления полиномиальной регрессии.
dd
p = np.linspace(3, 9, 100) # Построим вычисленую кривую.
plt.plot( x, y, 'g*', p, poly.polyval(p, dd ), 'r-' )
q = x.reshape( n, 1) # Решим эту же задачу вручную.
q
q**2
A = np.concatenate( [np.ones((n,1)), q, q**2, q**3, q**4, q**5, q**6, q**7, q**8], axis=1 )
Ai = np.linalg.pinv( A )
dd = Ai.dot( y.reshape(n, 1) )
dd = dd.reshape( dd.shape[0] )
dd #, dd.shape
#print( A.shape, Ai.shape, n )
AA = A.transpose().dot( A )
AAinv = np.linalg.inv( AA )
AApinv = AAinv.dot( A.transpose() )
#print( A.shape, AAinv.shape, AApinv.shape )
dd = AApinv.dot( y.reshape(n, 1) )
dd = dd.reshape( dd.shape[0] )
dd, dd.shape
AApinv2, residuals, rank, s= np.linalg.lstsq( A, y )
AApinv2
residuals, rank # Сумма ошибок отклонения. И ранг.
s # Собственные значения.
p = np.linspace(3, 9, 100)
plt.plot( x, y, 'g*', p, poly.polyval(p, dd ), 'r-' )
# Получилась очень страная кривая. Ветви параболы направлены вниз. Через точки почти не проходит
AA = A.transpose().dot( A )
AA = AA + 1.5 * np.eye( AA.shape[0] ) # Добавляем регулиризацию.
AAinv = np.linalg.inv( AA )
AApinv = AAinv.dot( A.transpose() )
#print( A.shape, AAinv.shape, AApinv.shape )
dd = AApinv.dot( y.reshape(n, 1) )
dd = dd.reshape( dd.shape[0] )
dd, dd.shape
p = np.linspace(3, 9, 100)
plt.plot( x, y, 'g*', p, poly.polyval(p, dd ), 'r-' )
AA = A.transpose().dot( A )
U, s, V = np.linalg.svd( AA )
s # Собственые значения. Очень большой порядок их изменения.
s[0]/s[-1] # Собственные значения отсортированы. Делим наибольшее на наименьшее.
AA = AA + 1.5 * np.eye( AA.shape[0] )
U, s, V = np.linalg.svd( AA )
s
s[0]/s[-1] # Порядок существенно уменьшен, хотя всеравно большой. np.trace?linalg.norm
a = np.array( [ [1, -2], [1, -2] ] )
U, s, V = np.linalg.svd(a, full_matrices=True)
U.shape, V.shape, s.shape
print( U )
print( V )
print( s )
S = np.diag( s )
np.allclose(a, np.dot(U, np.dot(S, V)))
Произведение векторов
# Скаляроне произведение.
a = np.array([ 1, 2, 3])
b = np.array([ -1, 0, 2])
np.inner( a, b) # Скяларное произведение
a * b # Поэлементное произведение векторов.
np.sum( a * b ) # Но можно и так.
b
# Комплексные числа.
1 + 2j
(1 + 2j) * (3 + 4j) # 1*3 + 1*4j + 2j*3 + 2j*4j = 3 + 4j + 6j + 8jj = -5 + 10j
# Сопряженное скалярное произведение.
a = np.array([ 1 + 2j, 3 + 4j])
b = np.array([ 1 - 2j, 3 + 4j])
np.vdot(a, b)
statmodels
import statsmodels.api as sm
import statsmodels.sandbox.regression.predstd as predstd
nsample = 100
x = np.linspace(-5, 5, 100)
X = np.column_stack((x, x, x**2)) # Отмечу, что переменная x повторена. т.е. есть избыток в данных.
beta = np.array([1, 0.05, 0.05, 0.5])
e = np.random.normal(size=nsample)
X = sm.add_constant(X)
y = np.dot(X, beta) + e
model = sm.OLS(y, X)
results = model.fit()
print(results.summary())
print('Parameters: ', results.params)
print('R2: ', results.rsquared)
yy = np.dot( X, results.params)
plt.plot( x, yy)
plt.plot( x, y, ".")
results = sm.OLS(y, sm.add_constant(x)).fit() # Все внимание на F-statistic, Porb(...), R-squared.
print(results.summary()) # Они показывают, что мы плохо приблизили. "Слуйчайная прямая" дала бы схожий результат.
yy = np.dot( sm.add_constant(x), results.params)
plt.plot( x, yy)
plt.plot( x, y, ".")
Доверительный интервал
model = sm.OLS(y, X) # Возвращаеся к прошлой модели.
results = model.fit()
yy = np.dot( X, results.params)
prstd, iv_l, iv_u = predstd.wls_prediction_std( results )
plt.plot(x, yy, 'g-', label="data")
plt.plot(x, y, 'b.', label="True")
plt.plot(x, results.fittedvalues, 'r--.', label="OLS")
plt.plot(x, iv_u, 'r--')
plt.plot(x, iv_l, 'r--')
plt.legend(loc='best');
R = [[0, 1, 0, 0], [0, 0, 1, 0]]
print(np.array(R))
print(results.f_test(R)) # т.е. нестоит отбрасывать гипотезу о том, что 2 и 3 параметр идентичны.
R = [[1, 0, 0, 0], [0, 0, 1, 0]]
print(np.array(R))
print(results.f_test(R)) # первая переменная и предпоследняя явно разные.