Введение в машинное обучение.
Вводятся базовые элементы пакета линейной алгебры (Numpy) относфщейся к питону (Python версии 3.xx). Конкретно речь идет о алгебре и, соотвествено, линейной регресии. Последнее, в частности, используется для ввода ключевых понятий из машинного обучения: выборка, обучающее множество, тестовое/валидационное множество, поиск "оптимального" решения и регуляризация. Уровни значимости (F-статистика) при построении регрессий (statmodels).
Это предварительная версия! Любые замечания приветствуются.
Введение
Регрессия
Отрисовка кривой
In [1]:
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
In [2]:
# Для начала случайно сгенерируем точки в которых будет вычислена функция.
np.random.uniform(10, 20, 5) # Создаем матрицу равномерных распределений. Первые два числа казывают диапазон.
Out[2]:
In [3]:
x = np.random.uniform( 3, 10, 10 ) # 10 точек на отрезке [3, 10]
In [4]:
x
Out[4]:
In [5]:
x.shape # Это действительно массив из 10 чисел.
Out[5]:
In [6]:
y = x * 2 - 15 # Вычисляем значение в каждой из точек. Функция линейная.
In [7]:
y
Out[7]:
In [8]:
plt.plot( x, y )
Out[8]:
In [9]:
yy = y + np.random.randn(10) # Добавим шум к каждому значению (измерению)
In [10]:
plt.plot( x, yy ) # Кривая как-то странно выглядит.
Out[10]:
In [11]:
xx = np.sort( x ) # Дело в том, что необходимо отсортировать числа по оси x.
yy = xx * 2 - 15 + np.random.randn(10) # Теперь заново вычисляем.
In [12]:
yy.shape, x.shape # Массивы соответсвуют друг другу по размеру.
Out[12]:
In [13]:
plt.plot( xx, yy )
Out[13]:
In [14]:
plt.plot( xx, yy, '*' )
Out[14]:
Линейная регрессия
см заметку по питону: массивы и графики
In [15]:
yy = xx * 2 - 15 + np.random.randn(10)/2 # Тоже самое но с меньшим шумом.
plt.plot( xx, yy, 'g.--') # Цвет зеленый, соединения пунктиром, точки маленькие.
Out[15]:
In [16]:
qq = np.array([1,2,3,4,5,6,7,8,9])
qq
Out[16]:
In [17]:
# Меняем форму.
qq.reshape(3,3) # 9 = 3 * 3
# Такое же количество элементов.
Out[17]:
In [18]:
qq # Исходный массив не изменился.
Out[18]:
In [19]:
# Если оличество элементов не будет совпадать.
qq.reshape(2,5) # 9 != 2 * 5
In [20]:
# Вектор столбец.
qq.reshape(9,1) # 9 - строк,
# 1 - колонка.
Out[20]:
In [21]:
# Превращяем массив в вектор столбец,
q=xx.reshape(10,1)
q # т.е. массив массивов (из одного элемента).
Out[21]:
In [22]:
np.ones(10).reshape(10,1) # Создаем вектор столбец единиц.
Out[22]:
In [23]:
# Можно сразу. ones -- одни единицы.
np.ones( (10, 1) )
Out[23]:
In [24]:
# Создаем общую матрицу. Присваеваем её переменной A.
A = np.concatenate( [np.ones((10,1)), q ], axis=1 )
A
Out[24]:
In [25]:
xx
Out[25]:
In [26]:
# Вычисляем псевдо обратную матрицу от A.
AA = np.linalg.pinv(A)
AA # Метод решения обычных систем не годится.
Out[26]:
In [27]:
# Теперь её используем для решения системы.
d = AA.dot( yy.reshape(10,1) )
d # В общем параметры приблизительно правельные.
Out[27]:
In [28]:
# Отрисуем получившуюся прямую. Используем для этого 100 точек.
p = np.linspace(3, 10, 100)
# Параметры вычисленной линии хранятся в переменной d.
yyy = p*d[1] + d[0]
plt.plot( xx, yy, 'g.-', p, yyy, 'r-')
Out[28]:
Функция потерь -- вычисление ошибки
In [29]:
# Какая точность этого приближения?
xx[0]*d[1] + d[0], yy[0] # В даннойт точк значения вычисленные и табулированые:
Out[29]:
In [30]:
(xx[0]*d[1] + d[0])[0], yy[0]
Out[30]:
In [31]:
for i in range(xx.shape[0]): # Цикл пробегает по всем нашим точкам
print ( xx[i]*d[1] + d[0], yy[i] ) # Выводит два значения: вычисленное и табулированное.
In [32]:
xx.shape[0], len(xx)
Out[32]:
In [33]:
def evalLine(d, x): # Создадим функцию, которая по параметрам прямой (d) считаем значение в точке (x).
return x * d[1] + d[0]
In [34]:
evalLine( d, xx[0]) # Вычислим значение в точке.
# Значение согласуется с предыдущими результатами.
Out[34]:
In [35]:
# Посчитаем среднею ошибку по всем точкам.
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 Почему сумма ошибка почти равна нулю?
In [36]:
# Обособим данный подсчет в функцию.
# По параметрам прямой (d) и точкам (x,y) вычислим общую ошибку.
def errLineCycle(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]
# А можно так. Без циклов.
def errLine(d, x, y):
xy = zip(x, y)
err = list(map( lambda z : evalLine( d, z[0]) - z[1], xy))
sum1 = sum( err )
err = np.array( err )
sum2 = sum( err * err )
return sum1/x.shape[0], sum2/x.shape[0]
Упр. Как ещё сильнее сократить код за счет функции mean?
In [37]:
errLine(d, xx, yy), errLineCycle(d, xx, yy)
Out[37]:
Обособим код вычисляющий линейную регрессию в функцию.
In [38]:
# x, y задают данные. На выходе параметры линейной регрессии.
def fitLin(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) )
In [39]:
d0 = fitLin( xx, yy) # Проверим, что соответсвует редыдущим результатам.
d0
Out[39]:
In [40]:
errLine(d0, xx, yy)
Out[40]:
Обучение и проверка
Learn and TestСлучайная выборка
In [41]:
perm = np.random.permutation( xx.shape[0] ) # Создаем перестановку числе от 0 до xx.shape[0] не включительно.
perm
Out[41]:
In [42]:
#np.random.randint(0, xx.shape[0], 5)
ii = perm[:5] # Берем первые 5 чисел.
ii
Out[42]:
In [43]:
# Вычисляем параметры линейной регрессии для подмножества числе (ii)
d1 = fitLin( xx[ii], yy[ii])
d1
Out[43]:
In [44]:
# Вычисляем ошибку на подмножестве.
ee = errLine(d1, xx[ii], yy[ii])
ee
Out[44]:
In [45]:
jj=perm[5:]
In [46]:
errLine(d1, xx[jj], yy[jj]) # Вычисляем ошибку на всем множестве точек.
Out[46]:
In [47]:
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] )
In [48]:
experimentLin(10, np.array([-15, 2.5]), 0.85, 1000)
ДЗ2: Построить статистику для ошибки. Вычислить её аналетически.
Регуляризация
Свободу модули нужно как-то регулировать.
Количество параметров -- Случай полинома
In [49]:
# Рассмотрим случай полиномов. Тонее на них мы изучим другое важное поянтие машинного обучения. Регулиризация.
import numpy.polynomial.polynomial as poly
In [50]:
d
Out[50]:
In [51]:
# Оказывается есть функция poly.polyval,
poly.polyval(xx[0], d ) # которая умеет вычилсять значение полинома.
Out[51]:
In [52]:
xx.shape[0]
Out[52]:
In [53]:
x.shape[0]
Out[53]:
In [54]:
ee[0]
Out[54]:
In [55]:
d = np.array( [-10, -28, 2.2] ) # Параметры задающие параболу.
In [56]:
n = 10
x = np.random.uniform( 3, 10, (n) )
x = np.sort( x )
x
Out[56]:
In [57]:
y = poly.polyval(x, d ) + np.random.randn( x.shape[0] )/2
In [58]:
y
Out[58]:
In [59]:
plt.plot( x, y, 'g*')
Out[59]:
In [60]:
# Есть библиотечная функция для вычисления полиномиальной регрессии.
dd = poly.polyfit( x, y, x.shape[0] - 2 )
dd
Out[60]:
In [61]:
p = np.linspace(3, 9, 100) # Построим вычисленую кривую.
plt.plot( x, y, 'g*', p, poly.polyval(p, dd ), 'r-' )
Out[61]:
In [62]:
p = np.linspace(3.5, 8.5, 100) # Построим вычисленую кривую.
plt.plot( x, y, 'g*', p, poly.polyval(p, dd ), 'r-' )
Out[62]:
При полной свободе криая прошла через все точки. В данном случае точки из обучающего множества.
Ручной режим
Решим эту же задачу вручную.
In [63]:
q = x.reshape( n, 1)
q
Out[63]:
In [64]:
q**2
Out[64]:
In [65]:
A = np.concatenate( [np.ones((n,1)), q, q**2, q**3, q**4, q**5, q**6, q**7, q**8], axis=1 )
In [66]:
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 )
Out[66]:
In [67]:
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
Out[67]:
Стабильность
In [68]:
AApinv2, residuals, rank, s= np.linalg.lstsq( A, y )
AApinv2
Out[68]:
In [69]:
residuals, rank # Сумма ошибок отклонения. И ранг.
Out[69]:
In [70]:
s # Собственные значения.
Out[70]:
In [71]:
p = np.linspace(3, 9, 100)
plt.plot( x, y, 'g*', p, poly.polyval(p, dd ), 'r-' )
# Получилась очень страная кривая. Ветви параболы направлены вниз. Через точки почти не проходит
Out[71]:
In [72]:
p = np.linspace(3.5, 8.5, 100)
plt.plot( x, y, 'g*', p, poly.polyval(p, dd ), 'r-' )
# Получилась очень страная кривая. Ветви параболы направлены вниз. Через точки почти не проходит
Out[72]:
Чего-то полного совпадения нет...
Магнитуда коэффициентов
Можно ограничить диапазон значений коэффициентов, т.е. мы считаем, что коэффициенты для рассматриваемых данных не должен привывашь тысячи....а может и 100.
Для этого нужно перейти к совсем ручному вычислению псевдо обратной матрицы.
In [73]:
AA = A.transpose().dot( A )
AA = AA + 15 * np.eye( AA.shape[0] ) # Добавляем регулиризацию. 150
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
Out[73]:
In [74]:
p = np.linspace(3, 9, 100)
plt.plot( x, y, 'g*', p, poly.polyval(p, dd ), 'r-' )
Out[74]:
In [75]:
AA = A.transpose().dot( A )
U, s, V = np.linalg.svd( AA )
s # Собственые значения. Очень большой порядок их изменения.
Out[75]:
In [76]:
s[0]/s[-1] # Собственные значения отсортированы. Делим наибольшее на наименьшее.
Out[76]:
In [77]:
AA = AA + 15 * np.eye( AA.shape[0] )
U, s, V = np.linalg.svd( AA )
s
Out[77]:
In [78]:
s[0]/s[-1] # Порядок существенно уменьшен, хотя всеравно большой. np.trace?linalg.norm
Out[78]:
In [79]:
a = np.array( [ [1, -2], [1, -2] ] )
U, s, V = np.linalg.svd(a, full_matrices=True)
U.shape, V.shape, s.shape
Out[79]:
In [80]:
print( U )
print( V )
print( s )
In [81]:
S = np.diag( s )
In [82]:
np.allclose(a, np.dot(U, np.dot(S, V)))
Out[82]:
In [83]:
# optimze gradient
In [ ]: