Uczenie maszynowe i sztuczne sieci neuronowe/Ćwiczenia 2

Regresja liniowa jako filtr

Poniższy przykład ma nam uświadomić możliwość wykorzystania regresji do filtrowania sygnałów. Przy okazji nauczymy się pracy z ciągiem uczącym

# -*- coding: utf-8 -*-
#from sklearn.neural_network import MLPRegressor #http://scikit-learn.org/dev/modules/generated/sklearn.neural_network.MLPRegressor.html#sklearn.neural_network.MLPRegressor
from sklearn import linear_model
from scipy.signal import firwin, lfilter
import numpy as np
import pylab as py

def spectrum(x,Fs):
    N = len(x)
    P = np.abs(np.fft.fft(x))**2/len(x)
    f = np.fft.fftfreq(len(x),1.0/Fs)
    P = P[0:N/2]
    P[1:-1] *=2
    f = f[0:N/2]
    return f,P

sampling = 128.0
rzad = 30
N_wej = rzad # rozmiar bufora wejściowego
demo = False # czy rysować tworzenie przykładów

T = 10 # ile skund sygnału
t  = np.arange(0,T,1/sampling)
#Losowy sygnał wejściowy
x = np.random.randn(len(t))

# przygotowujemy dane tak aby nauczyć filtrowania dolnoprzepustowego 
b= firwin(rzad,10.0/sampling)
y = lfilter(b,[1],x)


# tworzymy ciąg uczący 
N_przykladow = len(x)-N_wej-1 
X = np.zeros((N_przykladow,N_wej))
Y = np.zeros((N_przykladow,1))

for i in range(N_przykladow):
    X[i,:] = x[i:i+N_wej]
    Y[i] = y[i+N_wej-1]
    if demo:
        py.clf()
        py.plot(t,x,'r')
        py.plot(t,y,'b')
        py.plot(t[i:i+N_wej],X[i,:],'ro')
        py.plot(t[i+N_wej-1],Y[i],'bo')
        py.pause(1)
   
# tworzmy obiekt do regersji liniowej     
model = linear_model.LinearRegression(fit_intercept=False)
# dopasowujemy parametry modelu
model.fit(X, Y) 

# test na zbiorze uczącym: dla kolejnych buforów wejściowych obliczamy wartości wyjściowe
y_est = np.zeros(x.shape)
for i in range(len(x)-N_wej-1):
    bufor_wejsciowy = x[i:i+N_wej]
    y_est[i+N_wej-1] = model.predict(bufor_wejsciowy)
 
# ilustracja dopasowania modelu: 
# - górny panel sygnał wejściowy
# - dolny panel sygnał otrzymany jako predykcje modelu (czerwony) na tle filtrowanego sygnału  (zielony) 
py.figure(1)
py.subplot(2,1,1)
py.plot(t,x)
py.title('wejscie')

py.subplot(2,1,2)
py.plot(t,y,'g')
py.plot(t,y_est,'r')
py.title('wyjscie')


# test na innym zbiorze - dajemy zestaw sinusów o losowej fazie i czestosciach 0-64Hz 
x_test  = np.zeros(len(t))
faza = np.zeros(len(range(1,64,2)))
for i,f in enumerate(range(1,64,2)):
    faza[i] = np.random.rand(1)*2*np.pi
    x_test += np.sin(2*np.pi*f*t + faza[i])



y_test = lfilter(b,[1],x_test)
# test na zbiorze testowym: dla kolejnych buforów wejściowych obliczamy wartości wyjściowe
y_est = np.zeros(x.shape)
for i in range(len(x)-N_wej-1):
    bufor_wejsciowy = x_test[i:i+N_wej]
    y_est[i+N_wej-1] = model.predict(bufor_wejsciowy)
 
#ilustracja działania dopasowanego modelu na sygnał testowy       
py.figure(2)
py.subplot(2,2,1)
py.plot(t,x_test)
py.title(u'wejście')

py.subplot(2,2,3)
py.plot(t,y_est,'r')
py.plot(t,y_test,'g')
py.title(u'wyjście')

py.subplot(2,2,2)
f,P = spectrum(x_test[-sampling:],sampling)
py.stem(f,P)
py.title(u'widmo wejścia')
py.subplot(2,2,4)
f,P = spectrum(y_est[-sampling:],sampling)
py.stem(f,P)
py.title(u'widmo wyjścia')
py.show()


np.set_printoptions(precision =3, suppress = True)
print "Porównanie współczynników"
print "filtr:"
print b
print "model"
print model.coef_

Modelowanie i predykcja szeregu czasowego

Załóżmy, że mamy pewien szereg czasowy. Dalej zakładamy, że jest on wynikiem pewnego procesu autoregresyjnego, tzn. że jego kolejna próbka jest pewną kombinacją liniową poprzednich próbek z dodanym szumem:

[math] x_t = \sum_{i=1}^p a_i x_{t-i} +\epsilon_t[/math]

Czy dałoby się zastosować techniki uczenia sieci neuronowej do oszacowania współczynników [math]a_i[/math]?

1. Przy pomocy poniższego kodu proszę wygenerować realizację procesu AR o długości N = 2000 i współczynnikach a = [-0.5, 0.2], epsilon niech będzie = 0.1:

def generujAR(a, epsilon, N):
    x=np.zeros(N)
    rzad = len(a)
    for i in range(rzad,N):
        for p in range(len(a)):
            x[i] += a[p]*x[i-(p+1)]
        x[i] += epsilon*np.random.randn()
    return x

1. Proszę wykreślić tą realizację
2. Skonstruować ciąg uczący gdzie wejściem jest p poprzednich próbek zaś wyjściem kolejna próbka.
3. Skonstruuj sieć o ilości wejść równej rzędowi modelu p
4. Ucz sieć i obserwuj jak zmieniają się wagi.
5. Zbadaj co dzieje się z estymowanymi parametrami oraz błędem popełnianym przez sieć gdy rozmiar wejścia nie jest zgodny z rzędem modelu.

# -*- coding: utf-8 -*-
from pybrain.structure import FeedForwardNetwork, LinearLayer, FullConnection
from pybrain.supervised.trainers import BackpropTrainer
from pybrain.datasets import SupervisedDataSet
from scipy.signal import firwin, lfilter
import numpy as np
import pylab as py


def generujAR(a, epsilon, N):
    x=np.zeros(N)
    rzad = len(a)
    for i in range(rzad,N):
        for p in range(len(a)):
            x[i] += a[p]*x[i-(p+1)]
        x[i] += epsilon*np.random.randn()
    return x

N = 2000
x = np.zeros(N)
a = np.array([ -0.5, 0.2]) #[0.1 -0.2 0.3 0.1 -0.1]
epsilon =0.1
x = generujAR(a,epsilon,N)
py.plot(x)
py.show()

# tworzymy ciąg uczący CU
N_wej = 2
CU = SupervisedDataSet(N_wej, 1)
for i in range(len(x)-N_wej-1):
    bufor_wejsciowy = x[i:i+N_wej]
    wartosc_wyjsciowa = x[i+N_wej]
    CU.addSample( bufor_wejsciowy[::-1], wartosc_wyjsciowa)#




# wytwarzamy pustą sieć
siec = FeedForwardNetwork()

# tworzymy węzły wejściowe i wyjściowe

warstwaWejsciowa = LinearLayer(N_wej)
warstwaWyjsciowa = LinearLayer(1)

# dodajemy węzły do sieci
siec.addInputModule(warstwaWejsciowa)
siec.addOutputModule(warstwaWyjsciowa)

# łączymy węzły
wej_do_wyj = FullConnection(warstwaWejsciowa, warstwaWyjsciowa)
siec.addConnection(wej_do_wyj)

# inicjujemy strukturę sieci
siec.sortModules()
  
    
trainer = BackpropTrainer(siec, CU,learningrate=0.01, momentum=0.8,verbose=True)
#trainer.trainUntilConvergence(maxEpochs = 500)
for i in range(10):
    trainer.trainEpochs(1)
    print i, siec.params

# test na zbiorze uczącym
y_est = np.zeros(x.shape)
for i in range(len(x)-N_wej-1):
    bufor_wejsciowy = x[i:i+N_wej]
    y_est[i+N_wej] = siec.activate(bufor_wejsciowy[::-1])
#py.plot(t,x)
py.subplot(2,1,1)
py.plot(x,'b')
py.plot(y_est,'r')
py.subplot(2,1,2)
py.plot(x-y_est)
print 'wariancja reziduum, ', np.var(x-y_est)
print 'wariancja procesu AR, ', epsilon**2
print 'parametry procesu: ', str(a)
print 'parametry wyestymowane: ', str(siec.params)
py.show()