Uczenie maszynowe i sztuczne sieci neuronowe/Ćwiczenia 5

Regresja

Jednowymiarowa

W tym ćwiczeniu chcemy zastosować sieć do wykonania regresji nieliniowej.

Załóżmy, że mamy generator który dostarcza żądaną ilość par liczb zgodnie z pewnym modelem, np.:

def gen(ile):
    x = np.sort(5*np.random.rand(ile))
    y = (1+10*x+x**2)/(1+2*x**2)
    y+= 0.1*y*np.random.randn(ile)
    return(x,y)

• Niech generator dostarcza nam po 10 par liczb. Proszę narysować kilka realizacji.
• Proszę skonstruować sieć, którą można nauczyć zależności między punktami wejściowymi (x) i wyjściowymi (y). W tym celu najlepiej wykorzystać sieć z nieliniową warstwą ukrytą typu tangensa hiperbolicznego i liniową warstwę wyjściową buildNetwork(N_wej, N_hid, N_wyj, hiddenclass=TanhLayer, outclass=LinearLayer)

Taka kombinacja warstw nieliniowej i liniowej pozwala na reprezentację dowolniej ciągłej funkcji, pod warunkiem użycia właściwej ilości neuronów w warstwie ukrytej. Liniowa warstwa wyjściowa pozwala na rozszerzenie zbioru wartości.

• Podobnie jak w poprzednim zadaniu proszę zbadać ewolucję:
  • wag
  • błędu na zbiorze uczącym
  • błędu na zbiorze monitorującym
• Powyższe zależności proszę zaobserwować dla kilku rozmiarów warstwy ukrytej.
• Proszę wykreślić funkcję reprezentowaną przez sieć na tle punktów zbioru uczącego i prawdziwej (niezaszumionej) relacji y(x).
• Czy dla ustalonej architektury sieci (rozmiarów warstw) i ustalonego zbioru uczącego sieć zbiega do jednego rozwiązania? Dlaczego?
• Dla kilku ustalonych architektur zbadaj zależność błędu od rozmiaru zbioru uczącego.

Szkielet rozwiązania:

# -*- coding: utf-8 -*-
import matplotlib
#matplotlib.use('TkAgg')
import numpy as np
import pylab as py
  
class siec(object):
    def __init__(self, X, Y, N_hid=3):
        self.X = X
        self.Y = Y
        self.N_wej = X.shape[1] 
        self.N_wyj = Y.shape[1]
        self.N_hid = N_hid
                                
        # inicjujemy połączenia
        # wagi ułożone są tak, że w kolejnych wierszach są kolejne neurony 
        # a w kolumnach wagi od konkretnego neuronu 
        # to +1 jest wagą dla obciążenia
        self.w_1 = (2*np.random.random((self.N_hid, self.N_wej+1)) - 1)/self.N_wej # pomiędzy warstwą pierwszą (wejściem) a warstwą ukrytą
        self.w_2 = (2*np.random.random((self.N_wyj, self.N_hid+1)) - 1)/self.N_hid
        self.dw1 = np.zeros((self.N_hid, self.N_wej+1))
        self.dw2 = np.zeros((self.N_wyj, self.N_hid+1)) 
    
    def g1(self, x):
        y = 1./(1+np.exp(-x))
        return y   
    def g1_prim(self, x):
        y = x*(1-x)
        return y
    def g2(self, x):
        y = x
        return y   
    def g2_prim(self, x):
        y = 1
        return y
    def get_params(self):
        return np.concatenate((self.w_1.reshape(-1), self.w_2.reshape(-1)),1)
        
    def predict(self, x):
        # propagacja "w przód"
        self.a_0 = np.vstack((1,x))  # z warstwy wejściowej (zerowej) wychodzi a_0
        z_1 = np.dot( self.w_1, self.a_0 )# na warstwe 1 wchodzą iloczyny skalarne 
        self.a_1 = np.vstack((1,self.g1(z_1))) # dokładamy 1 i dostaję wyjście z warstwy 1
        z_2 = np.dot( self.w_2, self.a_1 ) # na warstwe 3 wchodzą iloczyny skalarne 
        self.a_2 = self.g2(z_2)
        return self.a_2

    def fit_one_step(self, eta1,eta2):
        self.bl = 0
        D_1 = np.zeros((self.N_hid, self.N_wej+1))
        D_2 = np.zeros((self.N_wyj, self.N_hid+1))
        for i in range(0,self.X.shape[0]):
            # weźmy przykład i-ty        
            x = self.X[i,:].reshape(self.N_wej,1)
            y = self.Y[i,:].reshape(self.N_wyj,1)
            self.a_2 = self.predict(x)
            
            # propagacja "wstecz"
            d_2 = (self.a_2 - y)*self.g2_prim(self.a_2)
            d_1 = np.dot(self.w_2.T, d_2) * self.g1_prim(self.a_1)#z_2
    
            # akumulujemy poprawki 
            D_2 +=  np.dot( d_2, self.a_1.T)
            D_1 +=  np.dot( d_1[1:], self.a_0.T)
    
            self.bl += np.dot(d_2.T,d_2)/self.X.shape[0]
        # uaktualniamy wagi
        self.w_1 -=  eta1*D_1 + eta2*self.dw1
        self.w_2 -=  eta1*D_2+  eta2*self.dw2
        self.dw1  =  eta1*D_1 
        self.dw2  =  eta1*D_2   
        return self.bl
def fun(x):
    return (1+10*x+x**2)/(1+2*x**2)
def gen(ile):
    x = np.sort(5*np.random.rand(ile)).reshape((ile,1))
    y = fun(x).reshape((ile,1))
    y+= 0.05*y*np.random.randn(ile).reshape((ile,1))
    return(x,y)
  
def main(argv=None):  
    #zbiór uczący:
    N_przykladow =37
    X,   Y      = gen(N_przykladow) # przykłady do ciągu uczącego
    X_m, Y_m    = gen(N_przykladow) # przykłady do ciągu monitorującego
    py.figure()
    py.plot(X,Y,'.')
    py.show()
                
    # definiujemy obiekt sieci:
    S = siec( X, Y, N_hid= ...)
    
    # liczba epok uczenia
    N_epochs = 1500
    # inicjuję tablice na ewolucje
    err  = np.zeros(N_epochs) #tablica na błąd zbioru uczącego
    err_m  = np.zeros(N_epochs) #tablica na błąd zbioru monitorującego
    wagi = np.zeros((N_epochs,len(S.get_params()))) #tablica na wagi
 
    eta1 = 0.005   
    eta2 = 0.8                                          
    for cykl in range(N_epochs):
        err[cykl] = ... # wykonaj krok uczenia
        for j, m in enumerate(X): # liczę średni błąd kwadratowy na zbiorze monitorującym:
            err_m[cykl] += (Y_m[j] - ... )**2
        err_m[cykl] /=  X_m.shape[0]# normalizuję aby uzyskać średni błąd kwadratowy
        wagi[cykl,:] = ... #pobieram wagi do zapamiętania
        
    #  rysunki
    py.subplot(2,1,1) # błędów
    py.plot(err,'b',label='zb. uczacy')
    py.plot(err_m,'r',label='zb. monitorujacy')
    py.title(u'błąd')
    py.legend()
    py.ylim([0,3])
 
    py.subplot(2,1,2) #wag
    py.plot(wagi)
    py.title('wagi')
    py.ylim([-3,3])
    py.draw()

    # funkcja reprezentowana przez sieć na tle punktów zbioru uczącego i prawdziwej (niezaszumionej) relacji y(x).
    x_testowe = np.linspace(0.1,7,100)
    y_testowe = np.zeros(100)
    for i,x in enumerate(x_testowe):
        y_testowe[i] = S.predict(x)
    # prawdziwa relacja z(x)
    z = fun(x_testowe)
    # rysunki:
    py.figure() 
    py.plot(x_testowe,y_testowe,'r', label='regresja')
    py.plot(x_testowe,z,'b', label='relacja prawdziwa')
    py.plot(X,Y,'bo',label='zb. uczacy')
    py.plot(X_m,Y_m,'mo',label='zb. monitorujacy')
    py.legend()
    py.show()

if __name__ == "__main__":
    main()

szkilet rozwiązania 

def main(argv=None): 
    N_epok = 1500 
    N_rozmiarow = 7
    N_prob = 6
     
    eta1 = 0.005   
    eta2 = 0.8 
    # inicjuję tablice na ewolucje
    bl    = np.zeros((N_rozmiarow,N_prob)) #tablica na błąd zbioru uczącego
    bl_m  = np.zeros((N_rozmiarow,N_prob)) #tablica na błąd zbioru monitorującego
    rozmiary = np.linspace(10,70,N_rozmiarow)
    for i,N_przykladow in enumerate(rozmiary):
        print ' '
        for proba in range(N_prob): 
            print 'test dla ', N_przykladow, 'przykładów; ''próba: ', proba
            #zbiór uczący:
            X,   Y      = gen(N_przykladow) # przykłady do ciągu uczącego
            X_m, Y_m    = gen(N_przykladow) # przykłady do ciągu monitorującego
            
            # definiujemy obiekt sieci:
            S = siec( X, Y, N_hid=1)
            for cykl in range(N_epok):
                bl[i,proba] = S.fit_one_step(eta1,eta2)
                if cykl%100==0:
                    print 'cykl:', cykl, u'błąd: ',bl[i,proba]
            
            for j, m in enumerate(X): # liczę średni błąd kwadratowy na zbiorze monitorującym:
                bl_m[i,proba] += (Y_m[j] - S.predict(X[j]))**2
            bl_m[i,proba] /=  X_m.shape[0]# normalizuję aby uzyskać średni błąd kwadratowy
     
        
    #  rysunki
    sr_bl    = np.mean(bl, axis = 1)
    std_bl   = np.std(bl, axis = 1) 
    sr_bl_m  = np.mean(bl_m, axis = 1)
    std_bl_m = np.std(bl_m, axis = 1) 
    
    py.figure()
    py.fill_between(rozmiary,sr_bl-std_bl,sr_bl+std_bl , alpha=0.1, color="b")
    py.plot(rozmiary,sr_bl, 'b',label='zb. uczacy')
    py.fill_between(rozmiary,sr_bl_m-std_bl_m,sr_bl_m+std_bl_m , alpha=0.1, color="r")
    py.plot(rozmiary,sr_bl_m, 'r',label='zb. monitorujacy')
    py.title(u'błąd')
    py.legend()
    py.ylim([0,1.5])
    py.show()

Dwuwymiarowa

W zadaniu tym chciałbym abyście zbadali na ile złożona (ile neuronów ukrytych) musi być sieć modelująca/interpolująca funkcję:

[math] f(x,y) = 0.1 + \frac{1 + \sin(2x + 3y)}{3.5 + \sin(x-y)}[/math]

dla [math]x,y \in [-2,2][/math]

• Proszę wykreślić tą funkcję.

• Jako ciąg uczący proszę wykorzystać pary wejścia spróbkowane co 0.2 i odpowiadające im wartości funkcji.
• Test proszę przeprowadzić na danych próbkowanych gęściej, ale tak aby w zbiorze testowym nie było punktów ze zbioru uczącego (np.: np.arange(-2+0.05195,2,0.05195)).

Przydatne mogą być funkcje do:

• generowania danych np.:

def gen(vec):
    x = vec
    y = vec
    f = np.zeros((len(x),len(y)))
    for i, x_i in enumerate(x):
        for j, y_j in enumerate(y):
            f[i,j] = 0.1 + (1+np.sin(2*x_i + 3*y_j))/(3.5 + np.sin(x_i - y_j)) 
    return(x,y,f)

(x,y,f)= gen(np.arange(-2,2,0.2))

• rysowania powierzchni 3D

from matplotlib.pyplot import  plot, show, figure, draw
from mpl_toolkits.mplot3d import axes3d

def rysuj3d(x,y,f,fig):
    X,Y = np.meshgrid(x,y)    
    ax = fig.gca(projection='3d')
    ax.plot_surface(X, Y, f, rstride=3, cstride=3, alpha=0.3)
    ax.set_xlabel('x')
    ax.set_ylabel('y')
    ax.set_zlabel('f')

def scatter3d(x,y,f,fig,kolor):
    X,Y = np.meshgrid(x,y)    
    ax = fig.gca(projection='3d')
    ax.scatter(X, Y, f, c=kolor,marker='o')
    ax.set_xlabel('x')
    ax.set_ylabel('y')
    ax.set_zlabel('f')

#przykładowe użycie:
(x,y,f)= gen(np.arange(-2,2,0.2))
fig = figure()
rysuj3d(x,y,f,fig)
scatter3d(x,y,f,fig,'b')
show