# Učitaj osnovne biblioteke...
import sklearn
import mlutils
import numpy as np
import scipy as sp
import matplotlib.pyplot as plt
%pylab inline

"""
Explanation: 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 4: Ansambli i procjena parametara
Verzija: 0.2
Zadnji put ažurirano: 7. prosinca 2018.
(c) 2015-2018 Jan Šnajder, Domagoj Alagić  
Objavljeno: 7. prosinca 2018.
Rok za predaju: 17. prosinca 2018. u 07:00h

Upute
Četvrta laboratorijska vježba sastoji se od četiri zadatka. Kako bi kvalitetnije, ali i na manje zamoran način usvojili gradivo ovog kolegija, potrudili smo se uključiti tri vrste zadataka: 1) implementacija manjih algoritama, modela ili postupaka; 2) eksperimenti s raznim modelima te njihovim hiperparametrima, te 3) primjena modela na (stvarnim) podatcima. Ovim zadatcima pokrivamo dvije paradigme učenja: učenje izgradnjom (engl. learning by building) i učenje eksperimentiranjem (engl. learning by experimenting).
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.
End of explanation
"""


from collections import Counter

class VotingClassifierDIY(object):
    
    SCHEME_COUNTING = "counting"
    SCHEME_AVERAGING = "averaging"
    
    def __init__(self, clfs, voting_scheme=SCHEME_COUNTING):
        self.clfs = clfs
        self.voting_scheme = voting_scheme
        
    def fit(self, X, y):
        [ clf.fit(X,y) for clf in self.clfs ]
    
    def predict_proba(self, X):
        if self.voting_scheme == self.SCHEME_AVERAGING:
            pred_avg = []; pred = []; prob = [];
            
            for i in range(0, len(self.clfs)):
                pred.append(self.clfs[i].predict_proba(X))
                
            for i in range(0, shape(pred)[1]):
                prob = [] 
                for j in range(0, len(pred)):
                    prob.append(pred[j][i])
                pred_avg.append(mean(prob, axis = 0))

            return pred_avg
    
        else:
            raise Exception("Nemoze")
    
        
        
    def predict(self, X):
        if self.voting_scheme == self.SCHEME_AVERAGING:
            avg = self.predict_proba(X)
            return [ argmax(a) for a in avg ]
        
        elif self.voting_scheme == self.SCHEME_COUNTING:
            pred_cnt = []; pred = []; predd = [];
            
            predd = [ clf.predict(X) for clf in self.clfs ]
            
            for i in range(0, shape(predd)[1]): # 0 - 1000
                for j in range(0, len(predd)): # 0 - 3
                    pred.append(predd[j][i])
            
            for i in range(0, len(pred), 3):
                pred_cnt.append((Counter(pred[i:i+3]).most_common()[0])[0])
                
            return pred_cnt
        
        else:
            raise Exception("Nemoze")

"""
Explanation: 1. Ansambli (glasovanje)
(a)
Vaš je zadatak napisati razred VotingClassifierDIY koji implementira glasački ansambl. Konstruktor razreda ima dva parametra: clfs koji predstavlja listu klasifikatora (objekata iz paketa sklearn) i voting_scheme koji označava radi li se o glasovanju prebrojavanjem (SCHEME_COUNTING) ili usrednjavanjem (SCHEME_AVERAGING). Glasovanje prebrojavanjem jednostavno vraća najčešću oznaku klase, dok glasovanje usrednjavanjem uprosječuje pouzdanosti klasifikacije u neku klasu (po svim klasifikatorima) te vraća onu s najvećom pouzdanošću. Primijetite da svi klasifikatori imaju jednake težine. O komplementarnosti klasifikatora vodimo računa tako da koristimo jednake klasifikatore s različitim hiperparametrima.
Razred sadržava metode fit(X, y) za učenje ansambla i dvije metode za predikciju: predict(X) i predict_proba(X). Prva vraća predviđene oznake klasa, a druga vjerojatnosti pripadanja svakoj od klasa za svaki od danih primjera iz X.
NB: Jedan od razreda koji bi Vam mogao biti koristan jest collections.Counter. Također vrijedi i za funkcije numpy.argmax i numpy.dstack.
End of explanation
"""


from sklearn.datasets import make_classification
from sklearn.ensemble import VotingClassifier
from sklearn.linear_model import LogisticRegression

X_voting, y_voting = make_classification(n_samples=1000, n_features=4, n_redundant=0, n_informative=3, n_classes=3, n_clusters_per_class=2)

Logisticka1 = LogisticRegression(multi_class = 'auto', solver = 'lbfgs', C = 1, max_iter = 1000).fit(X_voting, y_voting)
Logisticka2 = LogisticRegression(multi_class = 'auto', solver = 'lbfgs', C = 100, max_iter = 1000).fit(X_voting, y_voting)
Logisticka3 = LogisticRegression(multi_class = 'auto', solver = 'lbfgs', C = 1000, max_iter = 1000).fit(X_voting, y_voting)

clfs = [Logisticka1, Logisticka2, Logisticka3]

moj = VotingClassifierDIY(clfs = clfs)
njihov = VotingClassifier(estimators= [('lr1', Logisticka1), ('lr2', Logisticka2), ('lr3', Logisticka3)])

a = moj.fit(X_voting, y_voting)
a = njihov.fit(X_voting, y_voting)

moj.voting_scheme = moj.SCHEME_COUNTING
njihov.voting = 'hard'

print(all(moj.predict(X_voting) == njihov.predict(X_voting)))

moj.voting_scheme = moj.SCHEME_AVERAGING
njihov.voting = 'soft'

print(all(moj.predict_proba(X_voting) == njihov.predict_proba(X_voting)))


"""
Explanation: (b)
Uvjerite se da Vaša implementacija radi jednako onoj u razredu ensemble.VotingClassifier, i to pri oba načina glasovanja (parametar voting). Parametar weights ostavite na pretpostavljenoj vrijednosti. Za ovu provjeru koristite tri klasifikatora logističke regresije s različitom stopom regularizacije i brojem iteracija. Koristite skup podataka dan u nastavku. Ekvivalentnost implementacije najlakše je provjeriti usporedbom izlaza funkcije predict (kod prebrojavanja) i funkcije predict_proba (kod usrednjavanja).
NB: Ne koristimo SVM jer njegova ugrađena (probabilistička) implementacija nije posve deterministička, što bi onemogućilo robusnu provjeru Vaše implementacije.
End of explanation
"""


from sklearn.model_selection import train_test_split

X_bag, y_bag = make_classification(n_samples=1000, n_features=20, n_redundant=1, n_informative=17, n_classes=3, n_clusters_per_class=2)
X_bag_train, X_bag_test, y_bag_train, y_bag_test = train_test_split(X_bag, y_bag, train_size=0.7, test_size = 0.3, random_state=69)

"""
Explanation: Q: Kada je prebrojavanje bolje od usrednjavanja? Zašto? A obratno?
Q: Bi li se ovakav algoritam mogao primijeniti na regresiju? Kako?
2. Ansambli (bagging)
U ovom zadatku ćete isprobati tipičnog predstavnika bagging-algoritma, algoritam slučajnih šuma. Pitanje na koje želimo odgovoriti jest kako se ovakvi algoritmi nose s prenaučenošću, odnosno, smanjuje li bagging varijancu modela.
Eksperiment ćete provesti na danom skupu podataka:
End of explanation
"""


from sklearn.tree import DecisionTreeClassifier
from sklearn.metrics import zero_one_loss

tree1 = DecisionTreeClassifier(max_depth = None).fit(X_bag_train, y_bag_train)
h1 = tree1.predict(X_bag_test)

print(zero_one_loss(y_bag_test, h1))


tree2 = DecisionTreeClassifier(max_depth= 10000000000).fit(X_bag_train, y_bag_train)
h2 = tree2.predict(X_bag_test)

print(zero_one_loss(y_bag_test, h2))

"""
Explanation: Razred koji implementira stablo odluke jest tree.DecisionTreeClassifier. Prvo naučite stablo odluke (engl. decision tree) na skupu za učenje, ali tako da je taj model presložen. To možete postići tako da povećate najveću moguću dubinu stabla (parametar max_depth). Ispišite pogrešku na skupu za ispitivanje (pogrešku 0-1; pogledajte paket metrics).
End of explanation
"""


from sklearn.ensemble import RandomForestClassifier

L = list(range(0, 31))
error_train = []
error_test = []

for l in L:
    rfc = RandomForestClassifier(n_estimators=l+1).fit(X_bag_train, y_bag_train)
    h_train = rfc.predict(X_bag_train)
    h_test = rfc.predict(X_bag_test)
    error_train.append(zero_one_loss(y_bag_train, h_train))
    error_test.append(zero_one_loss(y_bag_test, h_test))
    
minimal_error = min(error_test)
index_of_min = error_test.index(minimal_error)

print(minimal_error)
print(index_of_min)

figure(figsize(10, 5))
plot(L, error_train)
plot(L, error_test)
scatter(L[index_of_min], minimal_error)
grid()

"""
Explanation: Sada isprobajte algoritam slučajnih šuma (dostupan u razredu ensemble.RandomForestClassifier) za različit broj stabala $L \in [1, 30]$. Iscrtajte pogrešku na skupu za učenje i na skupu za ispitivanje u ovisnosti o tom hiperparametru. Ispišite najmanju pogrešku na skupu za ispitivanje.
End of explanation
"""


from sklearn.datasets import make_circles

circ_X, circ_y = make_circles(n_samples=400, noise=0.1, factor=0.4)
mlutils.plot_2d_clf_problem(circ_X, circ_y)

"""
Explanation: Q: Što možete zaključiti iz ovih grafikona?
Q: Kako bagging postiže diverzifikaciju pojedinačnih osnovnih modela?
Q: Koristi li ovaj algoritam složeni ili jednostavni osnovni model? Zašto?
3. Ansambli (boosting)
U ovom zadatku pogledat ćemo klasifikacijski algoritam AdaBoost, koji je implementiran u razredu ensemble.AdaBoostClassifier. Ovaj algoritam tipičan je predstavnik boosting-algoritama.
Najprije ćemo generirati eksperimentalni skup podataka koristeći datasets.make_circles. Ova funkcija stvara dvodimenzijski klasifikacijski problem u kojem su dva razreda podataka raspoređena u obliku kružnica, tako da je jedan razred unutar drugog.
End of explanation
"""


from sklearn.ensemble import AdaBoostClassifier
from sklearn.tree import DecisionTreeClassifier

L = [1, 2, 3, 50]

figure(figsize(12, 8))

for i, l in enumerate(L):
    ada = AdaBoostClassifier(base_estimator = DecisionTreeClassifier(max_depth = 1), n_estimators = l).fit(circ_X, circ_y)
    
    subplot(3,3, i +1)
    mlutils.plot_2d_clf_problem(circ_X, circ_y, ada.predict)


"""
Explanation: (a)
Boosting, kao vrsta ansambla, također se temelji na kombinaciji više klasifikatora s ciljem boljih prediktivnih sposobnosti. Međutim, ono što ovakav tip ansambla čini zanimljivim jest to da za osnovni klasifikator traži slabi klasifikator (engl. weak classifier), odnosno klasifikator koji radi tek malo bolje od nasumičnog pogađanja. Često korišteni klasifikator za tu svrhu jest panj odluke (engl. decision stump), koji radi predikciju na temelju samo jedne značajke ulaznih primjera. Panj odluke specijalizacija je stabla odluke (engl. decision tree) koje smo već spomenuli. Panj odluke stablo je dubine 1. Stabla odluke implementirana su u razredu tree.DecisionTreeClassifier.
Radi ilustracije, naučite ansambl (AdaBoost) koristeći panj odluke kao osnovni klasifikator, ali pritom isprobavajući različit broj klasifikatora u ansamblu iz skupa $L \in {1, 2, 3, 50}$. Prikažite decizijske granice na danom skupu podataka za svaku od vrijednosti korištenjem pomoćne funkcije mlutils.plot_2d_clf_problem.
NB: Još jedan dokaz da hrvatska terminologija zaista može biti smiješna. :)
End of explanation
"""


X_boost, y_boost = make_classification(n_samples=1000, n_features=20, n_redundant=0, n_informative=18, n_classes=3, n_clusters_per_class=1)
X_boost_train, X_boost_test, y_boost_train, y_boost_test = train_test_split(X_boost, y_boost, train_size=0.7, random_state=69)

"""
Explanation: Q: Kako AdaBoost radi? Ovise li izlazi pojedinih osnovnih modela o onima drugih?
Q: Je li AdaBoost linearan klasifikator? Pojasnite.
(b)
Kao što je i za očekivati, broj klasifikatora $L$ u ansamblu predstavlja hiperparametar algoritma AdaBoost. U ovom zadatku proučit ćete kako on utječe na generalizacijsku sposobnost Vašeg ansambla. Ponovno, koristite panj odluke kao osnovni klasifikator.
Poslužite se skupom podataka koji je dan niže.
End of explanation
"""



L = list(range(0, 81))
error_train = []
error_test = []

for l in L:
    rfc = AdaBoostClassifier(base_estimator = DecisionTreeClassifier(max_depth = 1), n_estimators=l+1).fit(X_bag_train, y_bag_train)
    h_train = rfc.predict(X_bag_train)
    h_test = rfc.predict(X_bag_test)
    error_train.append(zero_one_loss(y_bag_train, h_train))
    error_test.append(zero_one_loss(y_bag_test, h_test))
    
minimal_error = min(error_test)
index_of_min = error_test.index(minimal_error)

print(minimal_error)
print(index_of_min)

figure(figsize(10, 5))
plot(L, error_train)
plot(L, error_test)
scatter(L[index_of_min], minimal_error)
grid()

"""
Explanation: Iscrtajte krivulje pogrešaka na skupu za učenje i ispitivanje u ovisnosti o hiperparametru $L \in [1,80]$. Koristite pogrešku 0-1 iz paketa metrics. Ispišite najmanju ostvarenu pogrešku na skupu za ispitivanje, te pripadajuću vrijednost hiperparametra $L$.
End of explanation
"""



D = list(range(1, 21))
error_train = []
error_test = []

for d in D:
    rfc = DecisionTreeClassifier(max_depth = d).fit(X_bag_train, y_bag_train)
    h_train = rfc.predict(X_bag_train)
    h_test = rfc.predict(X_bag_test)
    error_train.append(zero_one_loss(y_bag_train, h_train))
    error_test.append(zero_one_loss(y_bag_test, h_test))
    
minimal_error = min(error_test)
index_of_min = error_test.index(minimal_error)

print(minimal_error)
print(index_of_min)

figure(figsize(10, 5))
plot(D, error_train)
plot(D, error_test)
scatter(D[index_of_min], minimal_error)
grid()

"""
Explanation: Q: Može li uopće doći do prenaučenosti pri korištenju boosting-algoritama?
(c)
Kao što je rečeno na početku, boosting-algoritmi traže slabe klasifikatore kako bi bili najefikasniji što mogu biti. Međutim, kako se takav ansambl mjeri s jednim jakim klasifikatorom (engl. strong classifier)? To ćemo isprobati na istom primjeru, ali korištenjem jednog optimalno naučenog stabla odluke.
Ispišite pogrešku ispitivanja optimalnog stabla odluke. Glavni hiperparametar stabala odluka jest njihova maksimalna dubina $d$ (parametar max_depth). Iscrtajte krivulje pogrešaka na skupu za učenje i ispitivanje u ovisnosti o dubini stabla $d \in [1,20]$.
End of explanation
"""


def L(u, N, m):
    return (u**m) * (1-u)**(N-m)

"""
Explanation: Q: Isplati li se koristiti ansambl u obliku boostinga? Idu li grafikoni tome u prilog?<br>
Q: Koja je prednost boostinga nad korištenjem jednog jakog klasifikatora?
4. Procjena maksimalne izglednosti i procjena maksimalne aposteriorne vjerojatnosti
(a)
Definirajte funkciju izglednosti $\mathcal{L}(\mu|\mathcal{D})$ za skup $\mathcal{D}={x^{(i)}}_{i=1}^N$ Bernoullijevih varijabli. Neka od $N$ varijabli njih $m$ ima vrijednost 1 (npr. od $N$ bacanja novčića, $m$ puta smo dobili glavu). Definirajte funkciju izglednosti tako da je parametrizirana s $N$ i $m$, dakle definirajte funkciju $\mathcal{L}(\mu|N,m)$.
End of explanation
"""


m1 = [1,2,5,9]
N1 = 10
m2 = [1, 10, 50, 90]
N2 = 100
u = linspace(0, 1, 300)

figure(figsize(10,8))

subplot(2,1,1)
grid()
for m in m1:
    plot(u, L(u, N1, m))
    
subplot(2,1,2)
grid()
for m in m2:
    plot(u, L(u, N2, m))

"""
Explanation: (b)
Prikažite funkciju $\mathcal{L}(\mu|N,m)$ za (1) $N=10$ i $m=1,2,5,9$ te za (2) $N=100$ i $m=1,10,50,90$ (dva zasebna grafikona).
End of explanation
"""


N = 10

figure(figsize(10, 8))
plot(u, L(u, N, 0))
plot(u, L(u, N, 9))
grid()

"""
Explanation: Q: Koja vrijednost odgovara ML-procjenama i zašto?
(c)
Prikažite funkciju $\mathcal{L}(\mu|N,m)$ za $N=10$ i $m={0,9}$.
End of explanation
"""


from scipy.stats import beta

alphas = [ 1, 2, 3, 0.5, 1, 2, 3, 0.5 ]
betas = [ 1, 2, 3, 0.5, 0.5, 3, 2, 1 ]

for a, b in zip(alphas, betas):
    plot(u, beta.pdf(u, a, b))
grid()
axis([0, 1, 0, 3])

"""
Explanation: Q: Koja je ML-procjena za $\mu$ i što je problem s takvom procjenom u ovome slučaju?
(d)
Prikažite beta-distribuciju $B(\mu|\alpha,\beta)$ za različite kombinacije parametara $\alpha$ i $\beta$, uključivo $\alpha=\beta=1$ te $\alpha=\beta=2$.
End of explanation
"""


def P(u, N, m, a, b):
    return L(u, N, m) * beta.pdf(u, a, b)

N = 10
m = 9

for a, b in zip(alphas, betas):
    plot(u, P(u, N, m, a, b))
grid()

"""
Explanation: Q: Koje parametere biste odabrali za modeliranje apriornog znanja o parametru $\mu$ za novčić za koji mislite da je "donekle pravedan, ali malo češće pada na glavu"? Koje biste parametre odabrali za novčić za koji držite da je posve pravedan? Zašto uopće koristimo beta-distribuciju, a ne neku drugu?
(e)
Definirajte funkciju za izračun zajedničke vjerojatnosti $P(\mu,\mathcal{D}) = P(\mathcal{D}|\mu) \cdot P(\mu|\alpha,\beta)$ te prikažite tu funkciju za $N=10$ i $m=9$ i nekolicinu kombinacija parametara $\alpha$ i $\beta$.
End of explanation
"""


N = 10
m = 1
a = 4
b = 2

plot(u, beta.pdf(u, a, b))
plot(u, L(u, N, m))
plot(u, P(u, N, m, a, b))
grid()

"""
Explanation: Q: Koje vrijednosti odgovaraju MAP-procjeni za $\mu$? Usporedite ih sa ML-procjenama.
(f)
Za $N=10$ i $m=1$, na jednome grafikonu prikažite sve tri distribucije: $P(\mu,\mathcal{D})$, $P(\mu|\alpha,\beta)$ i $\mathcal{L}(\mu|\mathcal{D})$.
End of explanation
"""


from sklearn.datasets import load_iris
import itertools as it

iris = load_iris()
y = iris.target == 1
X = iris.data[y]

x = it.combinations(range(1, 5), 2)
figure(figsize(12, 28))

for i, (a, b) in enumerate(x):
    subplot(6, 1, i+1)
    grid()
    scatter(X[:, a-1], X[:, b-1])

"""
Explanation: (g)
Pročitajte ove upute o učitavanju oglednih skupova podataka u SciPy. Učitajte skup podataka Iris. Taj skup sadrži $n=4$ značajke i $K=3$ klase. Odaberite jednu klasu i odaberite sve primjere iz te klase, dok ostale primjere zanemarite (u nastavku radite isključivo s primjerima iz te jedne klase). Vizualizirajte podatke tako da načinite 2D-prikaze za svaki par značajki (šest grafikona; za prikaz je najjednostavnije koristiti funkciju scatter).
NB: Mogla bi Vam dobro dući funkcija itertools.combinations.
End of explanation
"""


from scipy.stats import norm

def log_izglednost(D, mi_ML, sigm2_ML):
    N = len(D)
    return -1*((N/2)*log(2*pi) + N*log(sqrt(sigm2_ML)) + (1/(2*sigm2_ML))*sum((D - mi_ML)**2))


"""
Explanation: (h)
Implementirajte funkciju log-izglednosti za parametre $\mu$ i $\sigma^2$ normalne distribucije.
End of explanation
"""


def mi_ML(D): 
    return sum(D)/len(D)

def sigm2_ML(D): 
    return sum((D - mi_ML(D))**2)/len(D)

for i in range(4):
    D = X[:,i]
    print('n=', i+1, ', mi_ML = ', str(mi_ML(D))[:4], ', sigm2_ML = ', str(sigm2_ML(D))[:7], ', izglednost = ', str(log_izglednost(D, mi_ML(D), sigm2_ML(D)))[:7]  )
    

"""
Explanation: (i)
Izračunajte ML-procjene za $(\mu, \sigma^2)$ za svaku od $n=4$ značajki iz skupa Iris. Ispišite log-izglednosti tih ML-procjena.
End of explanation
"""


from scipy.stats import pearsonr


x = it.combinations(range(0,4), 2)

for i,j in x:
    pear = pearsonr(X[:,i], X[:,j])[0]
    print('from ', i, ' ', j, ' pear: ', pear)

"""
Explanation: Q: Možete li, na temelju dobivenih log-izglednosti, zaključiti koja se značajka najbolje pokorava normalnoj distribuciji? 
(j)
Proučite funkciju pearsonr za izračun Pearsonovog koeficijenta korelacije. Izračunajte koeficijente korelacije između svih četiri značajki u skupu Iris.
End of explanation
"""


N = [int(len(X)/4), int(len(X)/2), len(X)]

for br in N:
    nepristran = cov(X[0:br,], rowvar = False, bias= False)
    pristran = cov(X[0:br,], rowvar = False, bias = True)
    
    diff_abs = mean(abs(nepristran - pristran))
    diff_sqr = mean(pow((nepristran - pristran),2))
    
    print('abs diff ', str(diff_abs)[:9], ' square diff ', diff_sqr, ' N ', br)

"""
Explanation: (k)
Proučite funkciju cov te izračunajte ML-procjenu za kovarijacijsku matricu za skup Iris. Usporedite pristranu i nepristranu procjenu. Pokažite da se razlika (srednja apsolutna i kvadratna) smanjuje s brojem primjera (npr. isprobajte za $N/4$ i $N/2$ i $N$ primjera).
End of explanation
"""


