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.

Proszę uzupełnić brakujace fragmenty kodu (w razie wątpliwości pytać prowadzącego):

# -*- coding: utf-8 -*-

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

# po początkowej obserwacji budowy ciągu uczącego warto przełżczyć wartość demo na False
demo = True # czy rysować tworzenie przykładów

T = 1 # 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,20.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: # poniższy fragment kodu ilustruje tworzenie ciągu uczącego
        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 instancję obiektu do regersji liniowej     
model = ... # tu wywołaj odpowiednią funkcję z modułu  linear_model
# dopasowujemy parametry modelu
model.fit(...) # uzupełnij argumenty 

# 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) 
# uzupełnij argumenty zgodnie z powyższym opisem
py.figure(1)
py.subplot(2,1,1)
py.plot(t,...) 
py.title('wejscie')

py.subplot(2,1,2)
py.plot(t,...,'g')
py.plot(t,...,'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(...) #uzupełnij
 
#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..... # uzupełnij

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ć regresję liniową do oszacowania współczynników [math]a_i[/math]?

Proszę uzupełnić poniższy kod, a następnie:

1. porównać wariancje procesu wygenerowanego i wyestymowanego.
2. do jakiego stopnia kolejna próbka w procesie AR jest przewidywalna?

# -*- coding: utf-8 -*-
from sklearn import linear_model
import numpy as np
import pylab as py


N = 20000
a = [0.1, -0.1, 0.5, 0.4,0.1]  #np.array([ -0.5, 0.2]) # inny zestaw parametrów do sprawdzenia
rzad = len(a)
epsilon =0.1

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

for i in range(rzad+1, N):
    for p in range(len(a)):
        x[i] += a[p]*x[i-(p+1)]
    x[i] += epsilon*np.random.randn()
    X[i,:] = x[i-1:i-rzad-1:-1]
    Y[i] = x[i] # chcemy wykonywać predykcję tej właśnie próbki


# wytwarzamy pusty model 
model = ...
# dopasowujemy parametry modelu
model...

# test na zbiorze uczącym
y_est = np.zeros(x.shape)
for i in range(rzad+1,len(x)-1):
    bufor_wejsciowy = x[i:i-rzad:-1]
    y_est[i+1] = model....

py.subplot(2,1,1)
py.plot(x,'b')
py.plot(...,'r')
py.title('procesy wygenerowany i estymowany')
py.subplot(2,1,2)
py.plot(x-y_est,'g')
py.title('residuum')

print 'wariancja reziduum, ', ...
print 'wariancja procesu AR, ', epsilon**2
print 'parametry procesu: ', str(a)
print 'parametry wyestymowane: ', ...
py.show()