Uczenie maszynowe i sztuczne sieci neuronowe/Ćwiczenia 10
Z Brain-wiki
Uczenie_maszynowe_i_sztuczne_sieci_neuronowe_cw/SVM1
Klasyfikacja za pomocą algorytmu wektorów wspierających (SVM)
Materiały
Na tych ćwiczeniach zapoznamy się z zastosowaniem SVM do klasyfikacji. Poniżej znajduje się moduł dostarczający kilku funkcji, z których dziś będziemy korzystać. Proszę zapisać go w bieżącym katalogu.
Kod modułu svm_modul.py
# -*- coding: utf-8 -*-
import numpy as np
import pylab as py
#####################################################################
# definicje funkcji pomocniczych, których nie musisz modyfikować
#####################################################################
def rysujDaneGrup(X, y, marker, xlabel, ylabel,legend_list):
'''X - macierz wartości wejściowych zorganizowana tak, że kolejne przykłady są
w wierszach, kolumny to kolejne wynmiary wejścia,
y - wektor określający przynależność do grupy, indeksy tego wektora odpowiadają wireszom macierzy X,
marker - zestaw markerów do oznaczania elementów grup, markerów powinno być tyle ile jest grup'''
p=[]
for i,g in enumerate(np.unique(y)):
g = int(g)
tmp =py.plot(X[np.where(y==g),0],X[np.where(y==g),1],marker[i])
p.append(tmp[0])
py.legend(p,legend_list)
# Dodajemy napisy
py.xlabel(xlabel)
py.ylabel(ylabel)
def rysujPodzial(model, X):
# wytworzymy siatkę punktów pokrywających obszar danych:
N = 100 # ilość punktów siatki w jednym wymiarze
os_x = np.linspace(X.min(),X.max(),N)
klasa = np.zeros((N,N))
for ix1, x1 in enumerate(os_x):
for ix2, x2 in enumerate(os_x):
XX = np.array([x1,x2]).reshape(1,2)
klasa[ix1,ix2] = svmPredict(model, XX) # dla każdego punktu siatki obliczamy jego klasę
x1_grid,x2_grid = np.meshgrid(os_x,os_x)
py.contourf(x1_grid, x2_grid, klasa.T,2)
def gaussianKernel(Xi,Xj,sigma=0.1):
z = np.dot((Xi-Xj).T,(Xi-Xj))
S = 2*sigma*sigma
Z= z/S
return np.exp(-Z)
def svmTrain(X, Y, C, kernelFunction, tol = 1e-3, max_passes = 5, sigma=0.1):
'''VMTRAIN Trenuje klasyfikator SVM za pomocą uproszczonego algorytmu SMO.
X - macierz wejściowa przykładów z ciągu uczącego wiersze - przyklady, kolumny - cechy
Y - wektor etykiet klas {-1,1}
C - regularyzacja SVM
tol - tolerancja na odstępstwa od wrunków KTT
max_passes - ile iteracji bez zmian mnożników Lagrangaea wykonać zanim uznamy, że już nie ma co poprawiać
kernelFunction - funkcja jądra, zaimplementowane są:
- gaussianKernel
- linearKernel
sigma - standardowe odchylenie dla jądra gaussowskiego
funkcja zwraca parametry doapsowanego modelu w słowniku model
Uwaga: To jest wersja uproszczona algorytmu wzorowana na
https://www.microsoft.com/en-us/research/wp-content/uploads/2016/02/tr-98-14.pdf
Nie jest ani specjalnie elegancka ani szybka.
Do praktycznych zastosowań zaleca się stosowanie zoptymalizowanych bibliotek:
LIBSVM (http://www.csie.ntu.edu.tw/~cjlin/libsvm/)
SVMLight (http://svmlight.joachims.org/)
lub sklearn
'''
# Pobieramy rozmiary
m,n = X.shape #m - ilość przykładów, n - wymiar wejścia
# Zmienne
alphas = np.zeros(m)
b = 0
E = np.zeros(m)
passes = 0
eta = 0
L = 0
H = 0
# Pre-compute the Kernel Matrix since our dataset is small
# (in practice, optimized SVM packages that handle large datasets
# gracefully will _not_ do this)
#
# We have implemented optimized vectorized version of the Kernels here so
# that the svm training will run faster.
print 'Obliczam macierz jądra'
if kernelFunction =='linearKernel':
# to jądro można policzyć od razu dla wszystkich przykładów
K = np.dot(X,X.T)
else:
# Jak nie możemy wymyśleć wektoryzacji obliczeń to
# obliczamy każdy element macierzy jądra osobno
K = np.zeros((m,m))
for i in range(m):
for j in range(i,m):
K[i,j] = gaussianKernel(X[i,:].T, X[j,:].T,sigma)
K[j,i] = K[i,j] #the matrix is symmetric
print 'Trenuję ...'
dots = 12
while passes < max_passes:
num_changed_alphas = 0
for i in range(m): #dla każdego przykładu z ciągu uczącego
# obliczamy błąd predykcji dla wektora i
E[i] = b + np.sum (alphas*Y*K[:,i]) - Y[i]
# jeśli jest co poprawiać:
if (( (Y[i]*E[i] < -tol) & (alphas[i] < C)) | ((Y[i]*E[i] > tol) & (alphas[i] > 0))):
# In practice, there are many heuristics one can use to select
# the i and j. In this simplified code, we select them randomly.
j = np.floor(m * np.random.rand())
while j == i: # Make sure i \neq j
j = np.floor(m * np.random.rand())
# Obliczamy błąd predykcji dla wektora j.
E[j] = b + np.sum (alphas*Y*K[:,j]) - Y[j]
# Save old alphas
alpha_i_old = alphas[i]
alpha_j_old = alphas[j]
# Oblicz przycięcia do pudełka [0,C]
if (Y[i] == Y[j]):
L = np.max((0, alphas[j] + alphas[i] - C))
H = np.min((C, alphas[j] + alphas[i]))
else:
L = np.max((0, alphas[j] - alphas[i]))
H = np.min((C, C + alphas[j] - alphas[i]))
if (L == H):
# continue to next i.
continue
# Compute eta by (15).
eta = 2 * K[i,j] - K[i,i] - K[j,j]
if (eta >= 0):
# continue to next i.
continue
# Compute and clip new value for alpha j using (16) and (17).
alphas[j] = alphas[j] - (Y[j] * (E[i] - E[j])) / eta
# Clip
alphas[j] = np.min ((H, alphas[j]))
alphas[j] = np.max ((L, alphas[j]))
#Check if change in alpha is significant
if (np.abs(alphas[j] - alpha_j_old) < tol):
# continue to next i.
# replace anyway
alphas[j] = alpha_j_old
continue
# Determine value for alpha i using (16).
alphas[i] = alphas[i] + Y[i]*Y[j]*(alpha_j_old - alphas[j])
# Compute b1 and b2 using (20) and (21) respectively.
b1 = b - E[i] - Y[i] * (alphas[i] - alpha_i_old) * K[i,j] - Y[j] * (alphas[j] - alpha_j_old) * K[i,j].T
b2 = b - E[j] - Y[i] * (alphas[i] - alpha_i_old) * K[i,j] - Y[j] * (alphas[j] - alpha_j_old) * K[j,j].T
# Compute b by (19).
if ( (0 < alphas[i]) & (alphas[i] < C)):
b = b1
elif (0 < alphas[j]) & (alphas[j] < C):
b = b2
else:
b = (b1+b2)/2
num_changed_alphas = num_changed_alphas + 1
if (num_changed_alphas == 0):
passes = passes + 1
else:
passes = 0
print num_changed_alphas
print ' Gotowe! \n\n'
# Save the model
idx = alphas > 0
model = {}
model['X'] = X[idx,:]
model['Y'] = Y[idx]
model['kernelFunction'] = kernelFunction
model['b'] = b
model['alphas']= alphas[idx]
model['w'] = (np.dot((alphas*Y).T, X)).T
model['sigma'] = sigma
print 'ilość wektorów wspierających: ', len(model['alphas'])
return model
def svmPredict(model,X):
'''model - model otrzymany z funkcji svmTrain
X - macierz m x n ,
w której wierszach są przykłady do sklasyfikowania (m)
każdy przykład ma wymiar n
funkcja zwraca wektor pred - i-ty element to predykcja dla i-tego przykładu
'''
# pobieramy rozmiary:
m,n = X.shape
#print 'm,n',m,n
# przygotowujemy tablice:
pred = np.zeros(m) # predyktory
margines = np.zeros(m) # wartości marginesów
if model['kernelFunction'] == 'linearKernel':
margines = np.dot(X , model['w']) + model['b']
elif model['kernelFunction'] =='gaussianKernel':
for i in range(m): #ta pętla iteruje po przykładach z macierzy X
for j in range(len(model['alphas'])): # ta pętlla iteruje po wektorach wspierających
margines[i] += model['alphas'][j]*model['Y'][j]* gaussianKernel(X[i,:],model['X'][j,:],model['sigma'])
margines[i] += model['b']
else:
print 'niezaimplementowane jądro '+ model['kernelFunction']
pred[margines >= 0] = 1
pred[margines < 0] = -1
return pred
Potrzebne będą nam też następujące zestawy danych:
Plik:Dane1.txt, Plik:Dane2.txt, Plik:Dane3.txt.
Proszę pobrać te pliki i zapisać je w bieżącym katalogu.
Ćwiczenie 1: Dane separowalne liniowo
Poniższy kod prezentuje zastosowanie SVM do problemu, który jest separowalny liniowo. Wykonując poniższy kod proszę zwrócić uwagę na punkt należący do klasy1 o współrzędnych (0.09, 4).
Jak pamiętamy z wykładu parametr C to współczynnik regularyzacji SVM, który karze za naruszanie marginesów. Proszę wykonać kod dla C o wartościach {1,2,5,10,20,30,60,120} i zaobserwować wyniki.
# -*- coding: utf-8 -*-
#importujemy potrzebne moduły i klasy
import numpy as np
import pylab as py
from svm_modul import *
#==================================================================
# Program
#==================================================================
# wczytywanie danych
dane = np.loadtxt('Dane1.txt') # dane zorganizowane są w trzech kolumnach
N_przyk, N_wej = dane.shape
X = dane[:,0:2] # pierwsze dwie kolumny to wejście
y = dane[:,2] # trzecia kolumna to etykiety klas
# narysujmy te dane
rysujDaneGrup(X, y, marker=('or','xb'), xlabel='x0', ylabel='x1',legend_list=('klasa0','klasa1'))
py.show()
# trenujemy model
model = svmTrain(X, y, C=100, kernelFunction = 'linearKernel', tol = 1e-3, max_passes = 20,sigma = 10)
# prezentujemy podział przestrzeni wejść reprezentowany przez model
rysujDaneGrup(X, y, marker=('or','xb'), xlabel='x0', ylabel='x1',legend_list=('klasa0','klasa1'))
rysujPodzial(model,X)
py.show()
Ćwiczenie 2: jądro Gaussowskie
W poprzednim programie proszę zmodyfikować wywołanie svmTrain:
- podmienić funkcję jądra na 'gaussianKernel'.
- ustawić C = 10
- zmieniać sigma na wartości: {0.1, 0.2, 0.4, 0.8, 1, 2, 4, 8}
Ćwiczenie 3: skomplikowany podział nieliniowy
Przy pomocy kodu z ćwiczenia 2 proszę dobrać parametry C i sigma aby otrzymać sensownie wyglądający podział przestrzeni dla danych zawartych w pliku dane2.txt.
Ćwiczenie 4: automatyzacja dobierania parametrów C i sigma
W wielu prawdziwych zastosowaniach chcielibyśmy aby nasz wybór parametrów był optymalny a jednocześnie możliwie obiektywny.
Powszechnie stosowaną metodą jest przeszukanie przestrzeni parametrów (C,sigma). Generuje się siatkę wartości (C,sigma) i dla każdego punktu siatki:
- estymuje się model
- ocenia się jakość generalizacji
Do oceny jakości w każdym punkcie siatki można zastosować albo zbiór monitorujący albo metody typu leave-one-out.
Uwaga: podział przestrzeni często wykonuje się w skali logarytmicznej.
Ćwiczenie wykonamy dla zbioru uczącego z pliku dane3.txt.
Musimy ten zbiór podzielić na dane do trenowania i dane do testowania np. w proporcji 3:1, Można to zrobić tak:
# wczytywanie danych
dane = np.loadtxt('dane3.txt') # dane zorganizowane są w trzech kolumnach
N_przyk, N_wej = dane.shape
X = dane[:,0:2] # pierwsze dwie kolumny to wejście
y = dane[:,2] # trzecia kolumna to etykiety klas
#podział na zbiór uczący i testujący
grupa0, = np.where(y==-1)
grupa1, = np.where(y==1)
# mieszamy kolejność indexów
np.random.shuffle(grupa0)
np.random.shuffle(grupa1)
# kopiujemy dane do zbioru uczącego (pierwsze 75% grupy0 i grupy1)
Xu = X[np.concatenate((grupa0[0: int(0.75*len(grupa0))],grupa1[0:int(0.75*len(grupa0))]))]
yu = y[np.concatenate((grupa0[0: int(0.75*len(grupa0))],grupa1[0:int(0.75*len(grupa0))]))]
# kopiujemy dane do zbioru testowego (końcowe 25% grupy0 i grupy1)
Xt = X[np.concatenate((grupa0[int(0.75*len(grupa0)):], grupa1[int(0.75*len(grupa0)):]))]
yt = y[np.concatenate((grupa0[int(0.75*len(grupa0)):], grupa1[int(0.75*len(grupa0)):]))]
# narysujmy te dane
rysujDaneGrup(Xu, yu, marker=('xr','xb'), xlabel='x0', ylabel='x1',legend_list=('klasa0','klasa1'))
rysujDaneGrup(Xt, yt, marker=('or','ob'), xlabel='x0', ylabel='x1',legend_list=('klasa0_test','klasa1_test'))
py.show()
Mając taki podział danych możemy dopasować model SVM do części uczącej:
model = svmTrain(Xu, yu, C=10, kernelFunction = 'gaussianKernel', tol = 1e-3, max_passes = 20,sigma = 0.5)
A następnie ocenić jego jakość na części testowej (funkcja svmPredict dostarczana jest przez moduł svm_modul.py):
TPR = np.sum(yt == svmPredict(model,Xt))/float(len(yt))
Proszę napisać program, który
- skanuje przestrzeń (C,sigma): C w zakresie od 0.1 do 100, sigma w zakresie od 0.1 do 10. Do wygenerowania zakresu ze skalą logarytmiczną można wykorzystać np. takie polecenie: zakresC = np.logspace(np.log2(0.1),np.log2(100),8, base=2)
- znajduje najlepsze parametry
- rysuje podział przestrzeni dla najlepszych parametrów.
Dodatek: implementacja w oparciu o bibliotekę LIBSVM
Do praktycznych zastosowań zaleca się stosowanie zoptymalizowanych bibliotek, np. LIBSVM (http://www.csie.ntu.edu.tw/~cjlin/libsvm/). Poniżej prezentujemy przykłady implementacji ćwiczenia 3 i ćwiczenia 4 w oparciu o tą bibliotekę. Aby móc wykonać te programy potrzebna jest skompilowana biblioteka libsvm i moduły pythonowe svm.py oraz svmutil. Cały zestaw wraz z instrukcją komilacji można pobrać z http://www.csie.ntu.edu.tw/~cjlin/libsvm/#download. Uwaga: w bibliotece tej jest nieco inna konwencja notacji jądra Gaussowskiego: [math]K(x,z) =\exp(-g ||x-z||^2)[/math], tzn. parametr [math]g = \frac{1}{2 \sigma^2}[/math]
cwiczenie3
# -*- coding: utf-8 -*-
#importujemy potrzebne moduły i klasy
import numpy as np
import pylab as py
from svm import *
from svmutil import *
def rysujDaneGrup(X, y, marker, xlabel, ylabel,legend_list):
'''X - macierz wartości wejściowych zorganizowana tak, że kolejne przykłady są
w wierszach, kolumny to kolejne wynmiary wejścia,
y - wektor określający przynależność do grupy, indeksy tego wektora odpowiadają wireszom macierzy X,
marker - zestaw markerów do oznaczania elementów grup, markerów powinno być tyle ile jest grup'''
p=[]
for i,g in enumerate(np.unique(y)):
g = int(g)
tmp =py.plot(X[np.where(y==g),0],X[np.where(y==g),1],marker[i])
p.append(tmp[0])
py.legend(p,legend_list)
# Dodajemy napisy
py.xlabel(xlabel)
py.ylabel(ylabel)
def rysujPodzial(model, X):
# wytworzymy siatkę punktów pokrywających obszar danych:
N = 100 # ilość punktów siatki w jednym wymiarze
os_x = np.linspace(X.min(),X.max(),N)
klasa = np.zeros((N,N))
for ix1, x1 in enumerate(os_x):
for ix2, x2 in enumerate(os_x):
XX = [[x1,x2]]#np.array([x1,x2]).reshape(1,2)
#print XX
p_label, p_acc, p_val =svm_predict([0], XX, model, '-b 1')
klasa[ix1,ix2] = p_label[0]
#svmPredict(model, XX) # dla każdego punktu siatki obliczamy jego klasę
x1_grid,x2_grid = np.meshgrid(os_x,os_x)
py.contourf(x1_grid, x2_grid, klasa.T,2)
#==================================================================
# Program
#==================================================================
# wczytywanie danych
dane = np.loadtxt('dane2.txt') # dane zorganizowane są w trzech kolumnach
N_przyk, N_wej = dane.shape
X = dane[:,0:2].tolist() # pierwsze dwie kolumny to wejście
y = dane[:,2].tolist() # trzecia kolumna to etykiety klas
prob = svm_problem(y, X, isKernel=False)
param = svm_parameter('-t 2 -c 10 -g 50 -b 1 -q')
''''options':
-s svm_type : set type of SVM (default 0)
0 -- C-SVC
1 -- nu-SVC
2 -- one-class SVM
3 -- epsilon-SVR
4 -- nu-SVR
-t kernel_type : set type of kernel function (default 2)
0 -- linear: u'*v
1 -- polynomial: (gamma*u'*v + coef0)^degree
2 -- radial basis function: exp(-gamma*|u-v|^2)
3 -- sigmoid: tanh(gamma*u'*v + coef0)
4 -- precomputed kernel (kernel values in training_set_file)
-d degree : set degree in kernel function (default 3)
-g gamma : set gamma in kernel function (default 1/num_features)
-r coef0 : set coef0 in kernel function (default 0)
-c cost : set the parameter C of C-SVC, epsilon-SVR, and nu-SVR (default 1)
-n nu : set the parameter nu of nu-SVC, one-class SVM, and nu-SVR (default 0.5)
-p epsilon : set the epsilon in loss function of epsilon-SVR (default 0.1)
-m cachesize : set cache memory size in MB (default 100)
-e epsilon : set tolerance of termination criterion (default 0.001)
-h shrinking : whether to use the shrinking heuristics, 0 or 1 (default 1)
-b probability_estimates : whether to train a SVC or SVR model for probability estimates, 0 or 1 (default 0)
-wi weight : set the parameter C of class i to weight*C, for C-SVC (default 1)
-v n: n-fold cross validation mode
-q : quiet mode (no outputs)
'''
m = svm_train(prob, param)
p_label, p_acc, p_val = svm_predict(y, X, m, '-b 1')
ACC, MSE, SCC = evaluations(y, p_label)
print ACC, MSE, SCC
# prezentujemy podział przestrzeni wejść reprezentowany przez model
rysujDaneGrup(np.array(X), np.array(y), marker=('or','xb'), xlabel='x0', ylabel='x1',legend_list=('klasa0','klasa1'))
rysujPodzial(m,np.array(X))
py.show()
cwiczenie4
# -*- coding: utf-8 -*-
#importujemy potrzebne moduły i klasy
import numpy as np
import pylab as py
from svm import *
from svmutil import *
def rysujDaneGrup(X, y, marker, xlabel, ylabel,legend_list):
'''X - macierz wartości wejściowych zorganizowana tak, że kolejne przykłady są
w wierszach, kolumny to kolejne wynmiary wejścia,
y - wektor określający przynależność do grupy, indeksy tego wektora odpowiadają wireszom macierzy X,
marker - zestaw markerów do oznaczania elementów grup, markerów powinno być tyle ile jest grup'''
p=[]
for i,g in enumerate(np.unique(y)):
g = int(g)
tmp =py.plot(X[np.where(y==g),0],X[np.where(y==g),1],marker[i])
p.append(tmp[0])
py.legend(p,legend_list)
# Dodajemy napisy
py.xlabel(xlabel)
py.ylabel(ylabel)
def rysujPodzial(model, X):
# wytworzymy siatkę punktów pokrywających obszar danych:
N = 100 # ilość punktów siatki w jednym wymiarze
os_x = np.linspace(X.min(),X.max(),N)
klasa = np.zeros((N,N))
for ix1, x1 in enumerate(os_x):
for ix2, x2 in enumerate(os_x):
XX = [[x1,x2]]#np.array([x1,x2]).reshape(1,2)
#print XX
p_label, p_acc, p_val =svm_predict([0], XX, model, '-b 1')
klasa[ix1,ix2] = p_label[0]
#svmPredict(model, XX) # dla każdego punktu siatki obliczamy jego klasę
x1_grid,x2_grid = np.meshgrid(os_x,os_x)
py.contourf(x1_grid, x2_grid, klasa.T,2)
#==================================================================
# Program
#==================================================================
# wczytywanie danych
dane = np.loadtxt('dane3.txt') # dane zorganizowane są w trzech kolumnach
N_przyk, N_wej = dane.shape
X = dane[:,0:2] # pierwsze dwie kolumny to wejście
y = dane[:,2] # trzecia kolumna to etykiety klas
#podział na zbiór uczący i testujący
grupa0, = np.where(y==-1)
grupa1, = np.where(y==1)
# mieszamy kolejność indexów
np.random.shuffle(grupa0)
np.random.shuffle(grupa1)
# kopiujemy dane do zbioru uczącego (pierwsze 75% grupy0 i grupy1)
Xu = X[np.concatenate((grupa0[0: int(0.75*len(grupa0))],grupa1[0:int(0.75*len(grupa0))]))]
yu = y[np.concatenate((grupa0[0: int(0.75*len(grupa0))],grupa1[0:int(0.75*len(grupa0))]))]
# kopiujemy dane do zbioru testowego (końcowe 25% grupy0 i grupy1)
Xt = X[np.concatenate((grupa0[int(0.75*len(grupa0)):], grupa1[int(0.75*len(grupa0)):]))]
yt = y[np.concatenate((grupa0[int(0.75*len(grupa0)):], grupa1[int(0.75*len(grupa0)):]))]
# narysujmy te dane
rysujDaneGrup(Xu, yu, marker=('xr','xb'), xlabel='x0', ylabel='x1',legend_list=('klasa0','klasa1'))
rysujDaneGrup(Xt, yt, marker=('or','ob'), xlabel='x0', ylabel='x1',legend_list=('klasa0_test','klasa1_test'))
py.show()
# trenujemy model
yul = yu.tolist()
Xul = Xu.tolist()
ytl = yt.tolist()
Xtl = Xt.tolist()
prob = svm_problem(yul, Xul, isKernel=False)
Ng = 20
zakresC = np.logspace(np.log2(0.1),np.log2(10),20, base=2)
zakresS = np.logspace(np.log2(0.1),np.log2(10),10, base=2)
bestC = 0
bestS = 0
bestTPR =0
TPR=np.zeros((20,10))
for i,C in enumerate(zakresC):
for j,S in enumerate(zakresS):
param_string = '-t 2 -q -b 1 -e 0.00001 -c '+str(C)+' -g '+str(S)
param = svm_parameter(param_string)
m = svm_train(prob, param)
p_label, p_acc, p_val = svm_predict(ytl, Xtl, m, '-b 1')
TPR[i,j] = np.sum(np.array(ytl)==np.array(p_label))/float(len(ytl))
if TPR[i,j]>bestTPR:
bestTPR = TPR[i,j]
bestModel = m
bestC = C
bestS = S
print C,S,TPR[i,j]
# prezentujemy podział przestrzeni wejść reprezentowany przez model
rysujDaneGrup(X, y, marker=('or','xb'), xlabel='x0', ylabel='x1',legend_list=('klasa0','klasa1'))
rysujPodzial(bestModel,X)
py.title('C: '+str(bestC)+' S: '+str(bestS))
py.figure()
py.pcolor(zakresC,zakresS,(TPR.T))
py.xlabel('C')
py.ylabel('S')
py.colorbar()
py.show()
