Sveučilište u Zagrebu
Fakultet elektrotehnike i računarstva

Strojno učenje 2018/2019

http://www.fer.unizg.hr/predmet/su


Laboratorijska vježba 1: Regresija

Verzija: 1.1
Zadnji put ažurirano: 12. listopada 2018.

(c) 2015-2018 Jan Šnajder, Domagoj Alagić, Mladen Karan

Objavljeno: 12. listopada 2018.
Rok za predaju: 22. listopada 2018. u 07:00h


Upute

Prva laboratorijska vježba sastoji se od osam zadataka. U nastavku slijedite upute navedene u ćelijama s tekstom. Rješavanje vježbe svodi se na dopunjavanje ove bilježnice: umetanja ćelije ili više njih ispod teksta zadatka, pisanja odgovarajućeg kôda te evaluiranja ćelija.

Osigurajte da u potpunosti razumijete kôd koji ste napisali. Kod predaje vježbe, morate biti u stanju na zahtjev asistenta (ili demonstratora) preinačiti i ponovno evaluirati Vaš kôd. Nadalje, morate razumjeti teorijske osnove onoga što radite, u okvirima onoga što smo obradili na predavanju. Ispod nekih zadataka možete naći i pitanja koja služe kao smjernice za bolje razumijevanje gradiva (nemojte pisati odgovore na pitanja u bilježnicu). Stoga se nemojte ograničiti samo na to da riješite zadatak, nego slobodno eksperimentirajte. To upravo i jest svrha ovih vježbi.

Vježbe trebate raditi samostalno. Možete se konzultirati s drugima o načelnom načinu rješavanja, ali u konačnici morate sami odraditi vježbu. U protivnome vježba nema smisla.


In [1]:
# Učitaj osnovne biblioteke...
import numpy as np
import sklearn
import matplotlib.pyplot as plt
import scipy as sp
%pylab inline


Populating the interactive namespace from numpy and matplotlib

Zadatci

1. Jednostavna regresija

Zadan je skup primjera $\mathcal{D}=\{(x^{(i)},y^{(i)})\}_{i=1}^4 = \{(0,4),(1,1),(2,2),(4,5)\}$. Primjere predstavite matrixom $\mathbf{X}$ dimenzija $N\times n$ (u ovom slučaju $4\times 1$) i vektorom oznaka $\textbf{y}$, dimenzija $N\times 1$ (u ovom slučaju $4\times 1$), na sljedeći način:


In [2]:
X = np.array([[0],[1],[2],[4]])
y = np.array([4,1,2,5])

(a)

Proučite funkciju PolynomialFeatures iz biblioteke sklearn i upotrijebite je za generiranje matrice dizajna $\mathbf{\Phi}$ koja ne koristi preslikavanje u prostor više dimenzije (samo će svakom primjeru biti dodane dummy jedinice; $m=n+1$).


In [3]:
from sklearn.preprocessing import PolynomialFeatures
Phi = PolynomialFeatures(1, False, True).fit_transform(X)
print(Phi)


[[1. 0.]
 [1. 1.]
 [1. 2.]
 [1. 4.]]

(b)

Upoznajte se s modulom linalg. Izračunajte težine $\mathbf{w}$ modela linearne regresije kao $\mathbf{w}=(\mathbf{\Phi}^\intercal\mathbf{\Phi})^{-1}\mathbf{\Phi}^\intercal\mathbf{y}$. Zatim se uvjerite da isti rezultat možete dobiti izračunom pseudoinverza $\mathbf{\Phi}^+$ matrice dizajna, tj. $\mathbf{w}=\mathbf{\Phi}^+\mathbf{y}$, korištenjem funkcije pinv.


In [4]:
from numpy import linalg

In [5]:
w = np.dot(np.dot(np.linalg.inv(np.dot(np.transpose(Phi), Phi)), np.transpose(Phi)), y)
print(w)

w2 = np.dot(np.linalg.pinv(Phi), y)
print(w2)


[2.2        0.45714286]
[2.2        0.45714286]

Radi jasnoće, u nastavku je vektor $\mathbf{x}$ s dodanom dummy jedinicom $x_0=1$ označen kao $\tilde{\mathbf{x}}$.

(c)

Prikažite primjere iz $\mathcal{D}$ i funkciju $h(\tilde{\mathbf{x}})=\mathbf{w}^\intercal\tilde{\mathbf{x}}$. Izračunajte pogrešku učenja prema izrazu $E(h|\mathcal{D})=\frac{1}{2}\sum_{i=1}^N(\tilde{\mathbf{y}}^{(i)} - h(\tilde{\mathbf{x}}))^2$. Možete koristiti funkciju srednje kvadratne pogreške mean_squared_error iz modula sklearn.metrics.

Q: Gore definirana funkcija pogreške $E(h|\mathcal{D})$ i funkcija srednje kvadratne pogreške nisu posve identične. U čemu je razlika? Koja je "realnija"?


In [6]:
from sklearn.metrics import mean_squared_error

h = np.dot(Phi, w)
print (h)

error = mean_squared_error(y, h)
print (error)

plt.plot(X, y, '+', X, h, linewidth = 1)
plt.axis([-3, 6, -1, 7])


[2.2        2.65714286 3.11428571 4.02857143]
2.042857142857143
Out[6]:
[-3, 6, -1, 7]

(d)

Uvjerite se da za primjere iz $\mathcal{D}$ težine $\mathbf{w}$ ne možemo naći rješavanjem sustava $\mathbf{w}=\mathbf{\Phi}^{-1}\mathbf{y}$, već da nam doista treba pseudoinverz.

Q: Zašto je to slučaj? Bi li se problem mogao riješiti preslikavanjem primjera u višu dimenziju? Ako da, bi li to uvijek funkcioniralo, neovisno o skupu primjera $\mathcal{D}$? Pokažite na primjeru.


In [7]:
try:
    np.dot(np.linalg.inv(Phi), y)
except LinAlgError as err:
    print(err)


Last 2 dimensions of the array must be square

(e)

Proučite klasu LinearRegression iz modula sklearn.linear_model. Uvjerite se da su težine koje izračunava ta funkcija (dostupne pomoću atributa coef_ i intercept_) jednake onima koje ste izračunali gore. Izračunajte predikcije modela (metoda predict) i uvjerite se da je pogreška učenja identična onoj koju ste ranije izračunali.


In [8]:
from sklearn.linear_model import LinearRegression

lr = LinearRegression().fit(Phi, y)

w2 = [lr.intercept_, lr.coef_[1]]
h2 = lr.predict(Phi)
error2 = mean_squared_error(y, h)

print ('staro: ')
print (w)
print (h)
print (error)

print('novo: ')
print (w2)
print (h2)
print (error2)


staro: 
[2.2        0.45714286]
[2.2        2.65714286 3.11428571 4.02857143]
2.042857142857143
novo: 
[2.2, 0.45714285714285713]
[2.2        2.65714286 3.11428571 4.02857143]
2.042857142857143

2. Polinomijalna regresija i utjecaj šuma

(a)

Razmotrimo sada regresiju na većem broju primjera. Koristite funkciju make_labels(X, f, noise=0) koja uzima matricu neoznačenih primjera $\mathbf{X}_{N\times n}$ te generira vektor njihovih oznaka $\mathbf{y}_{N\times 1}$. Oznake se generiraju kao $y^{(i)} = f(x^{(i)})+\mathcal{N}(0,\sigma^2)$, gdje je $f:\mathbb{R}^n\to\mathbb{R}$ stvarna funkcija koja je generirala podatke (koja nam je u stvarnosti nepoznata), a $\sigma$ je standardna devijacija Gaussovog šuma, definirana parametrom noise. Za generiranje šuma koristi se funkcija numpy.random.normal.

Generirajte skup za učenje od $N=50$ primjera uniformno distribuiranih u intervalu $[-5,5]$ pomoću funkcije $f(x) = 5 + x -2 x^2 -5 x^3$ uz šum $\sigma=200$:


In [9]:
from numpy.random import normal
def make_labels(X, f, noise=0) :
    return map(lambda x : f(x) + (normal(0,noise) if noise>0 else 0), X)

In [10]:
def make_instances(x1, x2, N) :
    return sp.array([np.array([x]) for x in np.linspace(x1,x2,N)])

In [11]:
N = 50
sigma = 200
fun = lambda x :5 + x - 2*x**2 - 5*x**3
x = make_instances(-5, 5, N)
y = list(make_labels(x, fun, sigma))

y6a = y
x6a = x

Prikažite taj skup funkcijom scatter.


In [12]:
plt.figure(figsize=(10, 5))
plt.plot(x, fun(x), 'r', linewidth = 1)
plt.scatter(x, y)


Out[12]:
<matplotlib.collections.PathCollection at 0x147c45b0>

(b)

Trenirajte model polinomijalne regresije stupnja $d=3$. Na istom grafikonu prikažite naučeni model $h(\mathbf{x})=\mathbf{w}^\intercal\tilde{\mathbf{x}}$ i primjere za učenje. Izračunajte pogrešku učenja modela.


In [13]:
from sklearn.preprocessing import PolynomialFeatures

Phi = PolynomialFeatures(3).fit_transform(x.reshape(-1, 1))
w = np.dot(np.linalg.pinv(Phi), y)
h = np.dot(Phi, w)
error = mean_squared_error(y, h)
print(error)

plt.figure(figsize=(10,5))
plt.scatter(x, y)
plt.plot(x, h, 'r', linewidth=1)


31571.274119582242
Out[13]:
[<matplotlib.lines.Line2D at 0x14806ad0>]

3. Odabir modela

(a)

Na skupu podataka iz zadatka 2 trenirajte pet modela linearne regresije $\mathcal{H}_d$ različite složenosti, gdje je $d$ stupanj polinoma, $d\in\{1,3,5,10,20\}$. Prikažite na istome grafikonu skup za učenje i funkcije $h_d(\mathbf{x})$ za svih pet modela (preporučujemo koristiti plot unutar for petlje). Izračunajte pogrešku učenja svakog od modela.

Q: Koji model ima najmanju pogrešku učenja i zašto?


In [14]:
Phi_d = []; 
w_d = []; 
h_d = [];
err_d = [];

d = [1, 3, 5, 10, 20]

for i in d:
    Phi_d.append(PolynomialFeatures(i).fit_transform(x.reshape(-1,1)))
    
for i in range(0, len(d)):
    w_d.insert(i, np.dot(np.linalg.pinv(Phi_d[i]), y))
    h_d.insert(i, np.dot(Phi_d[i], w_d[i]))

for i in range(0, len(d)):
    err_d.insert(i, mean_squared_error(y, h_d[i]))
    print (str(d[i]) + ': ' + str(err_d[i]))


fig = plt.figure(figsize=(15, 20))
fig.subplots_adjust(wspace=0.2) 

for i in range(0, len(d)):
    ax = fig.add_subplot(5, 2, i+1)
    ax.scatter(x, y);
    ax.plot(x, h_d[i], 'r', linewidth = 1)


1: 41403.53076562096
3: 31571.274119582242
5: 30963.173055623054
10: 27781.73904463034
20: 21580.616838362606

(b)

Razdvojite skup primjera iz zadatka 2 pomoću funkcije cross_validation.train_test_split na skup za učenja i skup za ispitivanje u omjeru 1:1. Prikažite na jednom grafikonu pogrešku učenja i ispitnu pogrešku za modele polinomijalne regresije $\mathcal{H}_d$, sa stupnjem polinoma $d$ u rasponu $d\in [1,2,\ldots,20]$. Radi preciznosti, funkcije $h(\mathbf{x})$ iscrtajte na cijelom skupu primjera (ali pogrešku generalizacije računajte, naravno, samo na ispitnome skupu). Budući da kvadratna pogreška brzo raste za veće stupnjeve polinoma, umjesto da iscrtate izravno iznose pogrešaka, iscrtajte njihove logaritme.

NB: Podjela na skupa za učenje i skup za ispitivanje mora za svih pet modela biti identična.

Q: Je li rezultat u skladu s očekivanjima? Koji biste model odabrali i zašto?

Q: Pokrenite iscrtavanje više puta. U čemu je problem? Bi li problem bio jednako izražen kad bismo imali više primjera? Zašto?


In [15]:
from sklearn.model_selection import train_test_split

In [16]:
X_train, X_test, y_train, y_test = train_test_split(x, y, test_size = 0.5)
err_train = [];
err_test = [];
d = range(0, 20)

for i in d:
    Phi_train = PolynomialFeatures(i).fit_transform(X_train.reshape(-1, 1))
    Phi_test = PolynomialFeatures(i).fit_transform(X_test.reshape(-1, 1))
    w_train = np.dot(np.linalg.pinv(Phi_train), y_train)
    h_train = np.dot(Phi_train, w_train)
    h_test = np.dot(Phi_test, w_train)
    
    err_train.insert(i, np.log(mean_squared_error(y_train, h_train)))
    err_test.insert(i, np.log(mean_squared_error(y_test, h_test)))

plt.figure(figsize=(10,5))
plt.plot(d, err_train, d, err_test)
plt.grid()


(c)

Točnost modela ovisi o (1) njegovoj složenosti (stupanj $d$ polinoma), (2) broju primjera $N$, i (3) količini šuma. Kako biste to analizirali, nacrtajte grafikone pogrešaka kao u 3b, ali za sve kombinacija broja primjera $N\in\{100,200,1000\}$ i količine šuma $\sigma\in\{100,200,500\}$ (ukupno 9 grafikona). Upotrijebite funkciju subplots kako biste pregledno posložili grafikone u tablicu $3\times 3$. Podatci se generiraju na isti način kao u zadatku 2.

NB: Pobrinite se da svi grafikoni budu generirani nad usporedivim skupovima podataka, na sljedeći način. Generirajte najprije svih 1000 primjera, podijelite ih na skupove za učenje i skupove za ispitivanje (dva skupa od po 500 primjera). Zatim i od skupa za učenje i od skupa za ispitivanje načinite tri različite verzije, svaka s drugačijom količinom šuma (ukupno 2x3=6 verzija podataka). Kako bi simulirali veličinu skupa podataka, od tih dobivenih 6 skupova podataka uzorkujte trećinu, dvije trećine i sve podatke. Time ste dobili 18 skupova podataka -- skup za učenje i za testiranje za svaki od devet grafova.

Q: Jesu li rezultati očekivani? Obrazložite.


In [17]:
N2 = [100, 200, 1000];
sigma = [100, 200, 500];

X_train4c_temp = [];
X_test4c_temp = [];
y_train4c_temp = [];
y_test4c_temp = [];

x_tmp = np.linspace(-5, 5, 1000);
X_train, X_test = train_test_split(x_tmp, test_size = 0.5)

for i in range(0, 3):
    
    y_tmp_train = list(make_labels(X_train, fun, sigma[i]))
    y_tmp_test = list(make_labels(X_test, fun, sigma[i]))
    for j in range(0,3):   
        X_train4c_temp.append(X_train[0:int(N2[j]/2)])
        X_test4c_temp.append(X_test[0:int(N2[j]/2)])
        y_train4c_temp.append(y_tmp_train[0:int(N2[j]/2)])
        y_test4c_temp.append(y_tmp_test[0:int(N2[j]/2)])

    
err_tr = [];
err_tst = [];

for i in range(0, 9):
    X_train4c = X_train4c_temp[i]
    X_test4c = X_test4c_temp[i]
    y_train4c = y_train4c_temp[i]
    y_test4c = y_test4c_temp[i]
    
    err_train4c = [];
    err_test4c = [];
    d4c = range(0, 20)

    for j in d4c:
        Phi_train4c = PolynomialFeatures(j).fit_transform(X_train4c.reshape(-1, 1))
        Phi_test4c = PolynomialFeatures(j).fit_transform(X_test4c.reshape(-1, 1))
        w_train4c = np.dot(np.linalg.pinv(Phi_train4c), y_train4c)
        h_train4c = np.dot(Phi_train4c, w_train4c)
        h_test4c = np.dot(Phi_test4c, w_train4c)
        err_train4c.insert(j, np.log(mean_squared_error(y_train4c, h_train4c)))
        err_test4c.insert(j, np.log(mean_squared_error(y_test4c, h_test4c)))
        
    err_tr.append(err_train4c);
    err_tst.append(err_test4c);
    
    
fig = plt.figure(figsize=(15, 10))
fig.subplots_adjust(wspace=0.2, hspace = 0.35) 

Nn = [100, 200, 1000, 100, 200, 1000, 100, 200, 1000]
sgm = [100, 100, 100, 200, 200, 200, 500, 500, 500]

for i in range(0, 9):    
    ax = fig.add_subplot(3, 3, i+1)
    plt.plot(d, err_tr[i], d, err_tst[i]); grid;
    ax.grid();


4. Regularizirana regresija

(a)

U gornjim eksperimentima nismo koristili regularizaciju. Vratimo se najprije na primjer iz zadatka 1. Na primjerima iz tog zadatka izračunajte težine $\mathbf{w}$ za polinomijalni regresijski model stupnja $d=3$ uz L2-regularizaciju (tzv. ridge regression), prema izrazu $\mathbf{w}=(\mathbf{\Phi}^\intercal\mathbf{\Phi}+\lambda\mathbf{I})^{-1}\mathbf{\Phi}^\intercal\mathbf{y}$. Napravite izračun težina za regularizacijske faktore $\lambda=0$, $\lambda=1$ i $\lambda=10$ te usporedite dobivene težine.

Q: Kojih je dimenzija matrica koju treba invertirati?

Q: Po čemu se razlikuju dobivene težine i je li ta razlika očekivana? Obrazložite.


In [18]:
lam = [0, 1, 10]
y = np.array([4,1,2,5])

Phi4a = PolynomialFeatures(3).fit_transform(X)
w_L2 = [];

def w_reg(lam): 
    t1 = np.dot(Phi4a.T, Phi4a) + np.dot(lam, np.eye(4))
    t2 = np.dot(np.linalg.inv(t1), Phi4a.T)
    return np.dot(t2, y)

for i in range(0, 3):
    w_L2.insert(i, w_reg(lam[i]))
    print (w_reg(lam[i]))


[ 4.         -5.91666667  3.375      -0.45833333]
[ 1.79567372 -0.24729075 -0.0175289   0.07014758]
[0.43312265 0.11060671 0.13827839 0.03093411]

(b)

Proučite klasu Ridge iz modula sklearn.linear_model, koja implementira L2-regularizirani regresijski model. Parametar $\alpha$ odgovara parametru $\lambda$. Primijenite model na istim primjerima kao u prethodnom zadatku i ispišite težine $\mathbf{w}$ (atributi coef_ i intercept_).

Q: Jesu li težine identične onima iz zadatka 4a? Ako nisu, objasnite zašto je to tako i kako biste to popravili.


In [19]:
from sklearn.linear_model import Ridge

In [20]:
for i in lam:
    w = [];
    w_L22 = Ridge(alpha = i).fit(Phi4a, y)
    
    w.append(w_L22.intercept_)
    for i in range(0, len(w_L22.coef_[1:])):
        w.append(w_L22.coef_[i])
    
    print (w)


[4.000000000000025, 0.0, -5.916666666666732, 3.375000000000028]
[3.0569614512471652, 0.0, -0.6907936507936514, -0.283174603174602]
[2.494441843122973, 0.0, -0.1589729487341473, -0.13423066536848305]

5. Regularizirana polinomijalna regresija

(a)

Vratimo se na slučaj $N=50$ slučajno generiranih primjera iz zadatka 2. Trenirajte modele polinomijalne regresije $\mathcal{H}_{\lambda,d}$ za $\lambda\in\{0,100\}$ i $d\in\{2,10\}$ (ukupno četiri modela). Skicirajte pripadne funkcije $h(\mathbf{x})$ i primjere (na jednom grafikonu; preporučujemo koristiti plot unutar for petlje).

Q: Jesu li rezultati očekivani? Obrazložite.


In [21]:
x5a = linspace(-5, 5, 50);
f = (5 + x5a - 2*(x5a**2) - 5*(x5a**3));
y5a = f + normal(0, 200, 50);

lamd = [0, 100]
dd = [2, 10]
h5a = []

for i in lamd:
    for j in dd:
        Phi5a = PolynomialFeatures(j).fit_transform(x5a.reshape(-1,1))
        w_5a = np.dot(np.dot(np.linalg.inv(np.dot(Phi5a.T, Phi5a) + np.dot(i, np.eye(j+1))), Phi5a.T), y5a);
        h_5a = np.dot(Phi5a, w_5a)
        h5a.append(h_5a)
       

lamdd = [0, 0, 100, 100]
ddd = [2, 10, 2, 10]

fig = plt.figure(figsize=(15, 10))
fig.subplots_adjust(wspace=0.2, hspace = 0.2) 

for i in range(0, len(lamdd)):    
    ax = fig.add_subplot(2, 2, i+1)
    plt.plot(x5a, h5a[i], 'r', linewidth = 2)
    plt.scatter(x5a, y5a);


(b)

Kao u zadataku 3b, razdvojite primjere na skup za učenje i skup za ispitivanje u omjeru 1:1. Prikažite krivulje logaritama pogreške učenja i ispitne pogreške u ovisnosti za model $\mathcal{H}_{d=20,\lambda}$, podešavajući faktor regularizacije $\lambda$ u rasponu $\lambda\in\{0,1,\dots,50\}$.

Q: Kojoj strani na grafikonu odgovara područje prenaučenosti, a kojoj podnaučenosti? Zašto?

Q: Koju biste vrijednosti za $\lambda$ izabrali na temelju ovih grafikona i zašto?


In [22]:
X5a_train, X5a_test, y5a_train, y5a_test = train_test_split(x5a, y5a, test_size = 0.5)
err5a_train = [];
err5a_test = [];
d = 20;
lambda5a = range(0, 50)

for i in lambda5a:
    Phi5a_train = PolynomialFeatures(d).fit_transform(X5a_train.reshape(-1, 1))
    Phi5a_test = PolynomialFeatures(d).fit_transform(X5a_test.reshape(-1, 1))
    w5a_train = np.dot(np.dot(np.linalg.inv(np.dot(Phi5a_train.T, Phi5a_train) + np.dot(i, np.eye(d+1))), Phi5a_train.T), y5a_train);
    h5a_train = np.dot(Phi5a_train, w5a_train)
    h5a_test = np.dot(Phi5a_test, w5a_train)
    
    err5a_train.insert(i, np.log(mean_squared_error(y5a_train, h5a_train)))
    err5a_test.insert(i, np.log(mean_squared_error(y5a_test, h5a_test)))

plt.figure(figsize=(8,4))
plt.plot(lambda5a, err5a_train, lambda5a, err5a_test);
plt.grid(), plt.xlabel('$\lambda$'), plt.ylabel('err');
plt.legend(['ucenje', 'ispitna'], loc='best');


6. L1-regularizacija i L2-regularizacija

Svrha regularizacije jest potiskivanje težina modela $\mathbf{w}$ prema nuli, kako bi model bio što jednostavniji. Složenost modela može se okarakterizirati normom pripadnog vektora težina $\mathbf{w}$, i to tipično L2-normom ili L1-normom. Za jednom trenirani model možemo izračunati i broj ne-nul značajki, ili L0-normu, pomoću sljedeće funkcije:


In [23]:
def nonzeroes(coef, tol=1e-6): 
    return len(coef) - len(coef[sp.isclose(0, coef, atol=tol)])

(a)

Za ovaj zadatak upotrijebite skup za učenje i skup za testiranje iz zadatka 3b. Trenirajte modele L2-regularizirane polinomijalne regresije stupnja $d=20$, mijenjajući hiperparametar $\lambda$ u rasponu $\{1,2,\dots,100\}$. Za svaki od treniranih modela izračunajte L{0,1,2}-norme vektora težina $\mathbf{w}$ te ih prikažite kao funkciju od $\lambda$.

Q: Objasnite oblik obiju krivulja. Hoće li krivulja za $\|\mathbf{w}\|_2$ doseći nulu? Zašto? Je li to problem? Zašto?

Q: Za $\lambda=100$, koliki je postotak težina modela jednak nuli, odnosno koliko je model rijedak?


In [24]:
from sklearn.linear_model import Ridge
from sklearn.linear_model import Lasso

lambda6a = range(1,100)
d6a = 20
X6a_train, X6a_test, y6a_train, y6a_test = train_test_split(x6a, y6a, test_size = 0.5)
Phi6a_train = PolynomialFeatures(d6a).fit_transform(X6a_train.reshape(-1,1))

L0 = [];
L1 = [];
L2 = [];

L1_norm = lambda w: sum(abs(w));
L2_norm = lambda w: math.sqrt(np.dot(w.T, w));

for i in lambda6a:
    w6a = np.dot(np.dot(np.linalg.inv(np.dot(Phi6a_train.T, Phi6a_train) + np.dot(i, np.eye(d6a+1))), Phi6a_train.T), y6a_train);
    
    L0.append(nonzeroes(w6a))
    L1.append(L1_norm(w6a))
    L2.append(L2_norm(w6a))
    
plot(lambda6a, L0, lambda6a, L1, lambda6a, L2, linewidth = 1)
grid()


(b)

Glavna prednost L1-regularizirane regresije (ili LASSO regression) nad L2-regulariziranom regresijom jest u tome što L1-regularizirana regresija rezultira rijetkim modelima (engl. sparse models), odnosno modelima kod kojih su mnoge težine pritegnute na nulu. Pokažite da je to doista tako, ponovivši gornji eksperiment s L1-regulariziranom regresijom, implementiranom u klasi Lasso u modulu sklearn.linear_model.


In [76]:
L0 = [];
L1 = [];
L2 = [];

for i in lambda6a:
    lass = Lasso(alpha = i, tol = 0.115).fit(Phi6a_train, y6a_train)
    w6b = lass.coef_
    
    L0.append(nonzeroes(w6b))
    L1.append(L1_norm(w6b))
    L2.append(L2_norm(w6b))
    
plot(lambda6a, L0, lambda6a, L1, lambda6a, L2, linewidth = 1)
legend(['L0', 'L1', 'L2'], loc = 'best')
grid()


7. Značajke različitih skala

Često se u praksi možemo susreti sa podatcima u kojima sve značajke nisu jednakih magnituda. Primjer jednog takvog skupa je regresijski skup podataka grades u kojem se predviđa prosjek ocjena studenta na studiju (1--5) na temelju dvije značajke: bodova na prijamnom ispitu (1--3000) i prosjeka ocjena u srednjoj školi. Prosjek ocjena na studiju izračunat je kao težinska suma ove dvije značajke uz dodani šum.

Koristite sljedeći kôd kako biste generirali ovaj skup podataka.


In [26]:
n_data_points = 500
np.random.seed(69)

# Generiraj podatke o bodovima na prijamnom ispitu koristeći normalnu razdiobu i ograniči ih na interval [1, 3000].
exam_score = np.random.normal(loc=1500.0, scale = 500.0, size = n_data_points) 
exam_score = np.round(exam_score)
exam_score[exam_score > 3000] = 3000
exam_score[exam_score < 0] = 0

# Generiraj podatke o ocjenama iz srednje škole koristeći normalnu razdiobu i ograniči ih na interval [1, 5].
grade_in_highschool = np.random.normal(loc=3, scale = 2.0, size = n_data_points)
grade_in_highschool[grade_in_highschool > 5] = 5
grade_in_highschool[grade_in_highschool < 1] = 1

# Matrica dizajna.
grades_X = np.array([exam_score,grade_in_highschool]).T

# Završno, generiraj izlazne vrijednosti.
rand_noise = np.random.normal(loc=0.0, scale = 0.5, size = n_data_points)
exam_influence = 0.9
grades_y = ((exam_score / 3000.0) * (exam_influence) + (grade_in_highschool / 5.0) \
            * (1.0 - exam_influence)) * 5.0 + rand_noise
grades_y[grades_y < 1] = 1
grades_y[grades_y > 5] = 5

a)

Iscrtajte ovisnost ciljne vrijednosti (y-os) o prvoj i o drugoj značajki (x-os). Iscrtajte dva odvojena grafa.


In [27]:
plt.figure()
plot(exam_score, grades_y, 'r+')
grid()

plt.figure()
plot(grade_in_highschool, grades_y, 'g+')
grid()


b)

Naučite model L2-regularizirane regresije ($\lambda = 0.01$), na podacima grades_X i grades_y:


In [28]:
w = Ridge(alpha = 0.01).fit(grades_X, grades_y).coef_
    print(w)


[0.00141497 0.09477276]

Sada ponovite gornji eksperiment, ali prvo skalirajte podatke grades_X i grades_y i spremite ih u varijable grades_X_fixed i grades_y_fixed. Za tu svrhu, koristite StandardScaler.


In [29]:
from sklearn.preprocessing import StandardScaler

#grades_y.reshape(-1, 1)

scaler = StandardScaler()
scaler.fit(grades_X)
grades_X_fixed = scaler.transform(grades_X)

scaler2 = StandardScaler()
scaler2.fit(grades_y.reshape(-1, 1))
grades_y_fixed = scaler2.transform(grades_y.reshape(-1, 1))

Q: Gledajući grafikone iz podzadatka (a), koja značajka bi trebala imati veću magnitudu, odnosno važnost pri predikciji prosjeka na studiju? Odgovaraju li težine Vašoj intuiciji? Objasnite.

8. Multikolinearnost i kondicija matrice

a)

Izradite skup podataka grades_X_fixed_colinear tako što ćete u skupu grades_X_fixed iz zadatka 7b duplicirati zadnji stupac (ocjenu iz srednje škole). Time smo efektivno uveli savršenu multikolinearnost.


In [30]:
grades_X_fixed_colinear = [[g[0],g[1],g[1]] for g in grades_X_fixed]

Ponovno, naučite na ovom skupu L2-regularizirani model regresije ($\lambda = 0.01$).


In [31]:
w = Ridge(alpha = 0.01).fit(grades_X_fixed_colinear, grades_y_fixed).coef_
print(w)


[[0.81630364 0.07583957 0.07583957]]

Q: Usporedite iznose težina s onima koje ste dobili u zadatku 7b. Što se dogodilo?

b)

Slučajno uzorkujte 50% elemenata iz skupa grades_X_fixed_colinear i naučite dva modela L2-regularizirane regresije, jedan s $\lambda=0.01$, a jedan s $\lambda=1000$. Ponovite ovaj pokus 10 puta (svaki put s drugim podskupom od 50% elemenata). Za svaki model, ispišite dobiveni vektor težina u svih 10 ponavljanja te ispišite standardnu devijaciju vrijednosti svake od težina (ukupno šest standardnih devijacija, svaka dobivena nad 10 vrijednosti).


In [32]:
w_001s = []
w_1000s = []

for i in range(10):
    X_001, X_1000, y_001, y_1000 = train_test_split(grades_X_fixed_colinear, grades_y_fixed, test_size = 0.5)
    w_001 = Ridge(alpha = 0.01).fit(X_001, y_001).coef_
    w_1000 = Ridge(alpha = 0.01).fit(X_1000, y_1000).coef_
    
    w_001s.append(w_001[0])
    w_1000s.append(w_1000[0])
    #print(w_001)
    #print(w_1000)
    #print()
    
from sklearn.preprocessing import StandardScaler
scaler = StandardScaler()
scaler.fit(w_001s)
print(scaler.mean_)

scaler = StandardScaler()
scaler.fit(w_1000s)
print(scaler.mean_)


[0.81111677 0.06649322 0.06649322]
[0.82191022 0.08521883 0.08521883]

Q: Kako regularizacija utječe na stabilnost težina?
Q: Jesu li koeficijenti jednakih magnituda kao u prethodnom pokusu? Objasnite zašto.

c)

Koristeći numpy.linalg.cond izračunajte kondicijski broj matrice $\mathbf{\Phi}^\intercal\mathbf{\Phi}+\lambda\mathbf{I}$, gdje je $\mathbf{\Phi}$ matrica dizajna (grades_fixed_X_colinear). Ponovite i za $\lambda=0.01$ i za $\lambda=10$.


In [50]:
lam = 0.01

phi = grades_X_fixed_colinear
s = np.add(np.dot(np.transpose(phi), phi), lam * np.identity(len(a)))
print(np.linalg.cond(s))

lam = 10

phi = grades_X_fixed_colinear
s = np.add(np.dot(np.transpose(phi), phi), lam * np.identity(len(a)))
print(np.linalg.cond(s))


100542.86653662546
101.54186653499605

Q: Kako regularizacija utječe na kondicijski broj matrice $\mathbf{\Phi}^\intercal\mathbf{\Phi}+\lambda\mathbf{I}$?