%%javascript
$.getScript('https://kmahelona.github.io/ipython_notebook_goodies/ipython_notebook_toc.js')

"""
Explanation: Les déterminants du choix de contraception en Indonésie
Présentation littéraire
Selon de récentes estimations, l'Indonésie a une population d'environ 255 millions d'habitants. Petit à petit elle voit grâce à l'intervention de son gouvernement son taux de fertilité diminuait jusqu'à atteindre en 2013 2.3 naissances par femme. Il est donc intéressant de voir les effets en pratique de sa politique et ses traductions dans le choix de contraception des femmes.
Présentation technique
Il s'agit d'un problème de classification multiclasse puisque la variable à prédire a trois modalités. Une spécificté plus particulière de notre jeu de données est de ne contenir quasiment que des variables catégorielles (7/9) ce qui modifie sensiblement notre approche, notamment dans la première partie de l'étude. Le déroulement de l'analyse est organisé en deux grandes parties, l'une reposant sur une description statistiques des données et l'autre sur un 
<h1 id="tocheading">Table des matières</h1>
<div id="toc"></div>
End of explanation
"""


import requests
presentation = "https://archive.ics.uci.edu/ml/machine-learning-databases/cmc/cmc.names"
f = requests.get(presentation)

print (f.text)


"""
Explanation: Récupération des données
Je cherche dans cette partie à avoir les instruments pour pouvoir commencer l'analyse.
Ainsi je récupère les données ainsi que leur présentation afin de ne pas avoir à recopier les indications.
Dans un deuxième temps je les partitionne en trois tables : celle contenant et les données et la variable à prédire (df), celle ne contenant que les données (data), et celle contenant que la variable à prédire (y).
Par ailleurs pour pouvoir documenter au mieux mes résultats, je crée un vecteur avec les intitulés des modalités de la variable à prédire afin de les ajouter sur les graphiques (y_names).
Importation de la présentation des données sur UCI
End of explanation
"""


import pandas, urllib.request
furl = urllib.request.urlopen("https://archive.ics.uci.edu/ml/machine-learning-databases/cmc/cmc.data")
df = pandas.read_csv(furl, names=["Age", "Educ", "Educ mari","Nbr Enfant","Religion","Statut trav","Statut mari","Niveau de vie","Exp Media","Contraception"])
print(df.head())
print(df.shape) 
# nous avons donc 1473 observations et 10 variables comme indiqué en introduction
 

#Création de trois tables et un vecteur : data y df y_names
copy= df.copy() #je ne veux pas écraser ma table d'entrée puisque je pourrais encore en avoir besoin
y= copy.pop('Contraception')# je ne récupère que la variable d'intérêt
print(y.shape)
data = copy# je récupère les autres colonnes qui restent soit les variables de prédiction
print(data.head())
print(data.shape)
y_names=["","No-Use","Long-Term","Short-Term"]# je garde en mémoire la signification des modalités de y pour pouvoir lire aisément la siginification des tableaux en sortie
print(y_names)


"""
Explanation: Importation des données
End of explanation
"""


data.hist(figsize=(12,12), alpha=0.40 )
#comme la plupart des variables sont des variables catégorielles 
# il est d'autant plus pertinent de ne choisir que la représentation en histogramme

"""
Explanation: Première approche: les statistiques descriptives
On cherche dans cette partie à avoir une première appréhension des données grâce à des outils standards de statistiques.
Représentation de la répartition des variables
Les variables explicatives: traitement du dataframe "data"
End of explanation
"""


plt.hist(y, 3, facecolor='red',alpha=0.20, ) 
plt.xlabel('Contraception')
plt.ylabel('Effectif')
plt.annotate('No-use', xy=(1.5, 200), xytext=(1.15, 200))
plt.annotate('Long-Term', xy=(1.5,200), xytext=(1.9, 200))
plt.annotate('Short-Term', xy=(3,200), xytext=(2.6, 200))
plt.figure()


"""
Explanation: La variable à expliquer : traitement de y
End of explanation
"""


print(df.Contraception.value_counts()) 
# Nous sommes rassurés d'autant plus que la modalité la plus importante est le non-use.
# Or d'un point de vue sociétal c'est sur celle-là qu'une politique publique peut vouloir jouer

"""
Explanation: Nous cherchons par ailleurs pour la suite à savoir si une des modalités est bien moins fréquente que les autres ce qui pourrait fausser l'analyse de classification.
End of explanation
"""


# Une vision en une multitude de graphe représentant les variables deux à deux
import seaborn
seaborn.pairplot(df)
# Cependant le fait que nos variables explicatives sont surtout des données catégorielles 
#  réduit le potentiel explicatif d'un tel graphique

"""
Explanation: Étude des interactions entre les variables
Nous sommmes en présence d'un jeu de données contenant majoritairement des variables catégorielles (8/10) ce qui nous empêche de mener une étude des corrélations, pourtant très commode pour avoir une vue synthétique des relations entre nos variables. Nous allons donc procéder ainsi: d'abord visualiser toutes les variables deux à deux pour comprendre leur covariation puis nous allons chercher un subsitut d'une analyse en composante prinicpales et enfin nous allons nous concentrer plus spécifiquement sur notre variable d'intérêt soit la contraception.
Entre toutes les variables : Pairplot
Nous utilisons dans cette partie délibérement des graphes qui permettent de présenter les résultats de façon synthétique vu le grand nombre de variables à comparer.
End of explanation
"""


#de l'art du copier/coller
from scipy.linalg import diagsvd
import numpy as np
import pandas as pd
import functools


def process_df(DF, cols, ncols):
    if cols:  # if you want us to do the dummy coding
        K = len(cols)  # the number of categories
        X = dummy(DF, cols)
    else:  # if you want to dummy code it yourself or do all the cols
        K = ncols
        if ncols is None:  # be sure to pass K if you didn't multi-index
            K = len(DF.columns)  # ... it with mca.dummy()
            if not K:
                raise ValueError("Your DataFrame has no columns.")
        elif not isinstance(ncols, int) or ncols <= 0 or \
                        ncols > len(DF.columns):  # if you dummy coded it yourself
            raise ValueError("You must pass a valid number of columns.")
        X = DF
    J = X.shape[1]
    return X, K, J


def dummy(DF, cols=None):
    """Dummy code select columns of a DataFrame."""
    return pd.concat((pd.get_dummies(DF[col])
        for col in (DF.columns if cols is None else cols)),
        axis=1, keys=DF.columns)


def _mul(*args):
    """An internal method to multiply matrices."""
    return functools.reduce(np.dot, args)


class MCA(object):
    """Run MCA on selected columns of a pd DataFrame.
    
    If the column are specified, assume that they hold
    categorical variables that need to be replaced with
    dummy indicators, otherwise process the DataFrame as is.
    'cols': The columns of the DataFrame to process.
    'ncols': The number of columns before dummy coding. To be passed if cols isn't.
    'benzecri': Perform Benzécri correction (default: True)
    'TOL': value below which to round eigenvalues to zero (default: 1e-4)
    """

    def __init__(self, DF, cols=None, ncols=None, benzecri=True, TOL=1e-4):

        X, self.K, self.J = process_df(DF, cols, ncols)
        S = X.sum().sum()
        Z = X / S  # correspondence matrix
        self.r = Z.sum(axis=1)
        self.c = Z.sum()
        self._numitems = len(DF)
        self.cor = benzecri
        self.D_r = np.diag(1/np.sqrt(self.r))
        Z_c = Z - np.outer(self.r, self.c)  # standardized residuals matrix
        self.D_c = np.diag(1/np.sqrt(self.c))

        # another option, not pursued here, is sklearn.decomposition.TruncatedSVD
        self.P, self.s, self.Q = np.linalg.svd(_mul(self.D_r, Z_c, self.D_c))

        self.E = None
        E = self._benzecri() if self.cor else self.s**2
        self.inertia = sum(E)
        self.rank = np.argmax(E < TOL)
        self.L = E[:self.rank]

    def _benzecri(self):
        if self.E is None:
            self.E = np.array([(self.K/(self.K-1.)*(_ - 1./self.K))**2
                              if _ > 1./self.K else 0 for _ in self.s**2])
        return self.E

    def fs_r(self, percent=0.9, N=None):
        """Get the row factor scores (dimensionality-reduced representation),
        choosing how many factors to retain, directly or based on the explained
        variance.
        'percent': The minimum variance that the retained factors are required
                                to explain (default: 90% = 0.9)
        'N': The number of factors to retain. Overrides 'percent'.
                If the rank is less than N, N is ignored.
        """
        if not 0 <= percent <= 1:
                raise ValueError("Percent should be a real number between 0 and 1.")
        if N:
                if not isinstance(N, (int, np.int64)) or N <= 0:
                        raise ValueError("N should be a positive integer.")
                N = min(N, self.rank)
                # S = np.zeros((self._numitems, N))
        # else:
        self.k = 1 + np.flatnonzero(np.cumsum(self.L) >= sum(self.L)*percent)[0]
        #  S = np.zeros((self._numitems, self.k))
        # the sign of the square root can be either way; singular value vs. eigenvalue
        # np.fill_diagonal(S, -np.sqrt(self.E) if self.cor else self.s)
        num2ret = N if N else self.k
        s = -np.sqrt(self.L) if self.cor else self.s
        S = diagsvd(s[:num2ret], self._numitems, num2ret)
        self.F = _mul(self.D_r, self.P, S)
        return self.F

    def fs_c(self, percent=0.9, N=None):
        """Get the column factor scores (dimensionality-reduced representation),
        choosing how many factors to retain, directly or based on the explained
        variance.
        'percent': The minimum variance that the retained factors are required
                                to explain (default: 90% = 0.9)
        'N': The number of factors to retain. Overrides 'percent'.
                If the rank is less than N, N is ignored.
        """
        if not 0 <= percent <= 1:
                raise ValueError("Percent should be a real number between 0 and 1.")
        if N:
                if not isinstance(N, (int, np.int64)) or N <= 0:
                        raise ValueError("N should be a positive integer.")
                N = min(N, self.rank)  # maybe we should notify the user?
                # S = np.zeros((self._numitems, N))
        # else:
        self.k = 1 + np.flatnonzero(np.cumsum(self.L) >= sum(self.L)*percent)[0]
        #  S = np.zeros((self._numitems, self.k))
        # the sign of the square root can be either way; singular value vs. eigenvalue
        # np.fill_diagonal(S, -np.sqrt(self.E) if self.cor else self.s)
        num2ret = N if N else self.k
        s = -np.sqrt(self.L) if self.cor else self.s
        S = diagsvd(s[:num2ret], len(self.Q), num2ret)
        self.G = _mul(self.D_c, self.Q.T, S)  # important! note the transpose on Q
        return self.G

    def cos_r(self, N=None):  # percent=0.9
        """Return the squared cosines for each row."""

        if not hasattr(self, 'F') or self.F.shape[1] < self.rank:
                self.fs_r(N=self.rank)  # generate F
        self.dr = np.linalg.norm(self.F, axis=1)**2
        # cheaper than np.diag(self.F.dot(self.F.T))?

        return np.apply_along_axis(lambda _: _/self.dr, 0, self.F[:, :N]**2)

    def cos_c(self, N=None):  # percent=0.9,
        """Return the squared cosines for each column."""

        if not hasattr(self, 'G') or self.G.shape[1] < self.rank:
                self.fs_c(N=self.rank)  # generate
        self.dc = np.linalg.norm(self.G, axis=1)**2
        # cheaper than np.diag(self.G.dot(self.G.T))?

        return np.apply_along_axis(lambda _: _/self.dc, 0, self.G[:, :N]**2)

    def cont_r(self, percent=0.9, N=None):
        """Return the contribution of each row."""

        if not hasattr(self, 'F'):
            self.fs_r(N=self.rank)  # generate F
        return np.apply_along_axis(lambda _: _/self.L[:N], 1,
                np.apply_along_axis(lambda _: _*self.r, 0, self.F[:, :N]**2))

    def cont_c(self, percent=0.9, N=None):  # bug? check axis number 0 vs 1 here
        """Return the contribution of each column."""

        if not hasattr(self, 'G'):
            self.fs_c(N=self.rank)  # generate G
        return np.apply_along_axis(lambda _: _/self.L[:N], 1,
                np.apply_along_axis(lambda _: _*self.c, 0, self.G[:, :N]**2))

    def expl_var(self, greenacre=True, N=None):
        """
        Return proportion of explained inertia (variance) for each factor.
        :param greenacre: Perform Greenacre correction (default: True)
        """
        if greenacre:
            greenacre_inertia = (self.K / (self.K - 1.) * (sum(self.s**4)
                                 - (self.J - self.K) / self.K**2.))
            return (self._benzecri() / greenacre_inertia)[:N]
        else:
            E = self._benzecri() if self.cor else self.s**2
            return (E / sum(E))[:N]

    def fs_r_sup(self, DF, N=None):
        """Find the supplementary row factor scores.
        ncols: The number of singular vectors to retain.
        If both are passed, cols is given preference.
        """
        if not hasattr(self, 'G'):
            self.fs_c(N=self.rank)  # generate G

        if N and (not isinstance(N, int) or N <= 0):
                raise ValueError("ncols should be a positive integer.")
        s = -np.sqrt(self.E) if self.cor else self.s
        N = min(N, self.rank) if N else self.rank
        S_inv = diagsvd(-1/s[:N], len(self.G.T), N)
        # S = scipy.linalg.diagsvd(s[:N], len(self.tau), N)
        return _mul(DF.div(DF.sum(axis=1), axis=0), self.G, S_inv)[:, :N]

    def fs_c_sup(self, DF, N=None):
        """Find the supplementary column factor scores.
        ncols: The number of singular vectors to retain.
        If both are passed, cols is given preference.
        """
        if not hasattr(self, 'F'):
            self.fs_r(N=self.rank)  # generate F

        if N and (not isinstance(N, int) or N <= 0):
                raise ValueError("ncols should be a positive integer.")
        s = -np.sqrt(self.E) if self.cor else self.s
        N = min(N, self.rank) if N else self.rank
        S_inv = diagsvd(-1/s[:N], len(self.F.T), N)
        # S = scipy.linalg.diagsvd(s[:N], len(self.tau), N)
        return _mul((DF/DF.sum()).T, self.F, S_inv)[:, :N]
    


# Discrétisation de la table df : création de dfdis
dfdis = df.copy()
dfdis["Nbr Enfant"]=pandas.qcut(dfdis["Nbr Enfant"],3,labels=["Enf1","Enf2","Enf3"])
dfdis["Age"]= pandas.qcut(dfdis["Age"],3,labels=["Age1","Age2","Age3"])
dfdis.head()
 


# Transformation des variables en dummies : création de la table dc
totdis=dfdis.copy()
totdis["Educ"]=pandas.Categorical(totdis["Educ"],ordered=False)
totdis["Educ"]=totdis["Educ"].cat.rename_categories(["Ed1","Ed2","Ed3","Ed4"])
totdis["Educ mari"]=pandas.Categorical(totdis["Educ mari"],ordered=False)
totdis["Educ mari"]=totdis["Educ mari"].cat.rename_categories(["Edm1","Edm2","Edm3","Edm4"])
totdis["Statut mari"]=pandas.Categorical(totdis["Statut mari"],ordered=False)
totdis["Statut mari"]=totdis["Statut mari"].cat.rename_categories(["St1","St2","St3","St4"])
totdis["Niveau de vie"]=pandas.Categorical(totdis["Niveau de vie"],ordered=False)
totdis["Niveau de vie"]=totdis["Niveau de vie"].cat.rename_categories(["Nv1","Nv2","Nv3","Nv4"])
totdis["Contraception"]=pandas.Categorical(totdis["Contraception"],ordered=False)
totdis["Contraception"]=totdis["Contraception"].cat.rename_categories(["Ctr1","Ctr2","Ctr3"])

dc=pandas.DataFrame(pandas.get_dummies(totdis[["Age","Nbr Enfant", "Educ","Educ mari","Statut mari","Niveau de vie","Exp Media","Religion","Statut trav","Contraception"]]))

dc.head()

# Création de la table des variations explicatives discrétisées : datadis

copy= dc.copy() #je ne veux pas écraser ma table d'entrée puisque je pourrais encore en avoir besoin
datadis=copy.drop(["Contraception_Ctr1","Contraception_Ctr2", "Contraception_Ctr3"], axis=1)
print(datadis.shape)
print(dc.shape)

# Exécution de l'ACM
mca_dc=MCA(dc,benzecri=False)
col=[1,2,3,4,4,4,5,5,5,6,6,6,6,7,7,7,7,8,8,8,8,9,9,9,9,10,10,10]
plt.scatter(mca_dc.fs_c()[:, 0],mca_dc.fs_c()[:, 1],c=col, cmap=plt.cm.Reds)
for i, j, nom in zip(mca_dc.fs_c()[:, 0],
    mca_dc.fs_c()[:, 1], dc.columns):
    plt.text(i, j, nom)
   
    
    
show()

"""
Explanation: Entre toutes les variables : ACM
Nous cherchons un graphique qui serait plus parlant. Vu le nombre de variable catégorielles, nous pensons à effectuer une analsye des correspondances multiples. Malheureusement il n'y a pas de modules évidents pour pouvoir en effectuer une. Nous sommes donc aller prendre le code source de Emre Safak sur Github. Il faut dans un second temps aussi transformer toutes nos variables en variables binaires pour pouvoir appliquer une acm
End of explanation
"""


import scipy # on ne l'étudie que sur les variables catégorielles
for i in (list(data.T.index[1:3,])): 
    u=pandas.crosstab(df.Contraception, df[i])
    print(u)
    print(scipy.stats.chi2_contingency(u))

for i in (list(data.T.index[4::,])):
    u=pandas.crosstab(df.Contraception, df[i])
    print(u)
    print(scipy.stats.chi2_contingency(u))


"""
Explanation: Je reconnais que l'ACM n'est pas facilement visible, le seul changement que j'ai réussi à faire : la couleur des points. Je n'ai pas trouvé comment l'agrandir. On remarque surtout les instances qui n'ont pas beaucoup d'impact sur la contraception quelle qu'elle soit : le grand nbr d'enfant, un âge important, éducation faible du mari et de la femme, une exposition nulle au média.
Cela nous encourage à aller regarder de plus près l'intéraction de notre variable d'intérêt avec les autres pour voir quels sont les facteurs d'influence.
Entre notre variable d'intérêt et les autres : test du chi2
On cherche dans un premier temps grâce au test du chi2, à savoir si toutes les variables socio-démographiques à notre disposition ont un impact sur le choix de contraception. On teste donc à partir des tableaux de contingence l'hypothèse nulle que chacune des catégories prises séparément n'influence pas la contraception, i.e on teste l'indépendance entre les variables
End of explanation
"""


#L'exposition au média
seaborn.countplot(x="Exp Media", hue="Contraception", data=df, palette='Greens_d')

# Le nombre d'enfant 
u = df.copy() # Il faut d'abord la discrétiser
u["Nbr Enfant"]= pandas.qcut(u["Nbr Enfant"],6)# le choix du nombre de catégories s'est effectuée selon l'allure de l'histogramme préenté au début
# Nous choisission délibérément qcut parce que proche de la distribution initiale du fait la prise en compte des quantiles

seaborn.countplot(x="Nbr Enfant", hue="Contraception", data=u, palette='Greens_d')

u.head()


"""
Explanation: Entre notre variable d'intérêt et les autres : histogramme conditionnel
Nous utilisons une fonction de seaborn qui permet de représenter sous forme d'histogramme les tableaux de contingence. On peut donc comparer la structure de l'histogramme de chaque modalité des variables explicatives organisé selon la contraception.
Détail technique: J'aurais souhaiter obtenir cette représentation selon chaque variable plotée dans un même gaphe pour avoir une vue panoramique. Cependant que ce soit avec la commande subplot de matplotlib ou Facetgrid de seaborn j'ai pas réussi. Je présente donc ici juste deux variables: l'exposition au média et le nombre d'enfants.
End of explanation
"""


from sklearn.cross_validation import train_test_split
data_train, data_test, y_train, y_test = train_test_split(data, y, test_size=0.33, random_state=42)
df.shape, data_test.shape, y_test.shape


import timeit
tic=timeit.default_timer()
from sklearn import tree
DecTree = tree.DecisionTreeClassifier(min_samples_leaf=10, min_samples_split=10)
DecTree.fit(data_train, y_train)
y_predDec = DecTree.fit(data_train, y_train).predict(data_test)
toc=timeit.default_timer()
toc - tic 

from sklearn.metrics import confusion_matrix, f1_score
def plot_confusion_matrix(cm, title='Matrice de confusion', cmap=plt.cm.Blues):
    plt.imshow(cm, interpolation='nearest', cmap=cmap)
    plt.title(title)
    plt.colorbar()
    plt.tight_layout()
    plt.ylabel('Vraie contraception')
    plt.xlabel('Contraception prédite')
    

# Compute confusion matrix
cm = confusion_matrix(y_test, y_predDec)
np.set_printoptions(precision=2)
print('Matrice de confusion')
print(cm)
plt.figure()
plot_confusion_matrix(cm)
    
plt.show()
f1Dec = f1_score(y_test, y_predDec)
print('f1Dec',f1Dec)

"""
Explanation: Ainsi on peut en déduire que le potentiel explicatif de Exp Media est moins important que celui de Nbr Enf car la distribution selon le mode de contraception a la même structure selon chacune des modalités de Exp médi alors qu'elle différe beaucoup suivant le nbr d'enfants.
Seconde approche : machine learning et supervised classification
Grand balayage : Decision Tree, AdaBoost et GaussianNB
Decisition Tree
End of explanation
"""


#Application avec 800 estimateurs et temps de calcul
#Il est intéressant pour AdaBoost de voir le temps mis car c'est un modèle qui met du temps à tourner sur mon jeu de données
import timeit
tic=timeit.default_timer()

from sklearn.ensemble import AdaBoostClassifier
AdaBoost = AdaBoostClassifier(DecTree,
                         algorithm='SAMME',
                         n_estimators=800,
                         learning_rate=0.5)
y_predAda= AdaBoost.fit(data_train, y_train).predict(data_test)

toc=timeit.default_timer()
toc - tic 

# AdaBoost : performance de prédiction
cm = confusion_matrix(y_test, y_predAda)
np.set_printoptions(precision=2)
print('Matrice de confusion')
print(cm)
plt.figure()
plot_confusion_matrix(cm)

f1Ada = f1_score(y_test, y_predAda)
print('f1Ada',f1Ada)

"""
Explanation: AdaBoost
End of explanation
"""


tic=timeit.default_timer()
from sklearn.naive_bayes import GaussianNB
Gaussian = GaussianNB()
y_predGaus= Gaussian.fit(data_train, y_train).predict(data_test)

toc=timeit.default_timer()
toc - tic 

cm = confusion_matrix(y_test, y_predGaus)
np.set_printoptions(precision=2)
print('Matrice de confusion')
print(cm)
plt.figure()
plot_confusion_matrix(cm)
f1Gaus = f1_score(y_test, y_predGaus)
print('f1Gaus',f1Gaus)


l = [ {"Méthode":"DécisionTree", "f1":f1Dec},
         {"Méthode":"AdaBoost", "f1":f1Ada},
         {"Méthode":"Gaussian", "f1":f1Gaus}]
f1=pd.DataFrame(l)
f1

"""
Explanation: Remarque: l'Adaboost est un modèle qui est vraiment plus lent à faire tourner
GaussianNB
End of explanation
"""


datadis_train, datadis_test, y_train, y_test = train_test_split(datadis, y, test_size=0.33, random_state=42)
df.shape, datadis_test.shape, y_test.shape

"""
Explanation: Amélioration 1 : discrétiser les variables prédictives
Pour aller rapidement au résultat nous supprimons la représentation de la matrice de confusion et nous nous limitons à deux classifier, ceux au meilleur score: AdaBoost et DecisionTree
End of explanation
"""


from sklearn import tree
DecTree = tree.DecisionTreeClassifier(min_samples_leaf=10, min_samples_split=10)
DecTree.fit(datadis_train, y_train)
y_predDecdis = DecTree.fit(datadis_train, y_train).predict(datadis_test)


cm = confusion_matrix(y_test, y_predDecdis)
np.set_printoptions(precision=2)
print('Matrice de confusion')
print(cm)
f1disDec = f1_score(y_test, y_predDecdis)
print("f1disDec", f1disDec)

"""
Explanation: Decision Tree
End of explanation
"""


y_preddisAda= AdaBoost.fit(datadis_train, y_train).predict(datadis_test)
cm = confusion_matrix(y_test, y_preddisAda)
np.set_printoptions(precision=2)
print('Matrice de confusion')
print(cm)

f1disAda = f1_score(y_test, y_preddisAda)
print('f1disAda',f1disAda)

l = [ {"Données":"Normales", "f1Dec":f1Dec, "f1Ada":f1Ada},
      {"Données":"Discrétisées", "f1Dec":f1disDec, "f1Ada":f1disAda},]
f1=pd.DataFrame(l)
f1 = f1.set_index("Données")
f1

"""
Explanation: AdaBoost
End of explanation
"""


#visulalisation de la contraception mal prédite avec DecisionTree
mismatchesd= (y_predDec != y_test) #liste de booléen avec True si y_pred différent de y_test
mismatchesTd= mismatchesd.copy()
u = mismatchesTd[mismatchesTd==True] #je ne retiens que les y_pred différent
usort = u.index.sort_values() #je prends leurs indices
print(usort[2::])# j'exclue au hasard les deux premiers pour pouvoir avoir des arrays de même taille afin de comparer


# Visualisation de la contraception malprédite avec Adaboost
mismatchesa= (y_predAda != y_test)
mismatchesTa= mismatchesa.copy()
v=mismatchesTa[mismatchesTa==True]
vsort= v.index.sort_values()
vsort


#On compare les deux sorties
comp = (usort[2::] == vsort)
comp


"""
Explanation: Les prévisions ne sont pas excellentes : j'essaie donc de voir si ce sont les mêmes femmes qui sont toujours mal regroupées
Amélioration 2 : spécificité des individus mal prédits
End of explanation
"""


#Étude sur le nombre d'estimateurs
from sklearn.ensemble import RandomForestClassifier
from sklearn.metrics import mean_squared_error
curves = []
for n_estimators in range(1,10) :
    clf = RandomForestClassifier(n_estimators=n_estimators, min_samples_leaf=1)
    clf = clf.fit(data_train, y_train)
    erra = mean_squared_error( clf.predict(data_train), y_train)**0.5
    errb = mean_squared_error( clf.predict(data_test), y_test)**0.5
    print("mn_estimators",n_estimators, "erreur",erra,errb)
    curves.append((n_estimators, erra,errb, clf) )
plt.plot ( [c[0] for c in curves], [c[1] for c in curves], label="train")
plt.plot ( [c[0] for c in curves], [c[2] for c in curves], label="test")
plt.legend()

tic=timeit.default_timer()
# Application avec 5 estimateurs
RF= RandomForestClassifier(n_estimators=5,min_samples_leaf=1)
y_predRF= RF.fit(data_train, y_train).predict(data_test)

cm = confusion_matrix(y_test, y_predRF)
np.set_printoptions(precision=2)
print('Matrice de confusion')
print(cm)

f1RF = f1_score(y_test, y_predRF)
print('f1RF',f1RF)
toc=timeit.default_timer()
toc - tic 

"""
Explanation: Ainsi on remarque que ce n'est pas un certain groupe de femme qui a tendance à être moins bien prédit
Amélioration 3 : Utilisation d'autres modèles  et hyperparamètres
Nous utilisons à présent des modèles avec des hyperparamètres qu'il faut tester en premier pour choisir le bon nombre
Random forest
Le  paramètre à ajuster en premier dans les méthodes de RandomForest est le nombre d'estimateurs soit le nombre d'arbres dans la forêt. Ici nous n'avons pas a priori un problème de temps de résolution puisque nous n'avons que peu d'observations et de variables. Cependant l'arbitrage se fait tout de même puisque à partir d'un nombre suffisant d'estimateurs, le modèle ne fait que du surajustement. Nous cherchons donc dans un premier temps le nombre optimal d'arbres.
End of explanation
"""


#Étude sur le nombre de neighbors
from sklearn.neighbors import KNeighborsClassifier
for n_neighbors in range(1, 11):
    clf = KNeighborsClassifier(n_neighbors=n_neighbors).fit(data_train, y_train)
    y_predNKeigh = clf.predict(data_test)
    print ("KNeighbors(n_neighbors={0})".format(n_neighbors), f1_score(y_test, y_predNKeigh))

tic=timeit.default_timer()
#Application avec n_neighbors = 7
KN = KNeighborsClassifier(n_neighbors=7).fit(data_train, y_train)
y_predKNeigh = KN.predict(data_test)

cm = confusion_matrix(y_test, y_predKNeigh)
np.set_printoptions(precision=2)
print('Matrice de confusion')
print(cm)

f1KNeigh = f1_score(y_test, y_predKNeigh)
print('f1KNeigh',f1KNeigh)

toc=timeit.default_timer()
toc - tic 

#Comparaison des résultats
l = [ {"Méthode":"DécisionTree", "f1":f1Dec},
         {"Méthode":"AdaBoost", "f1":f1Ada},
         {"Méthode":"RandomForest", "f1":f1RF},
         {"Méthode":"KNeighbors", "f1":f1KNeigh}
    ]
f1=pd.DataFrame(l)
f1

"""
Explanation: KNeighbors
End of explanation
"""


# nous créons un nouveau dataframe y_nonuse en regroupant les classes short-term et long-term
y_nonuse = y.copy()
y_nonuse[y_nonuse!=1] = 0 # Utilise une contraception
y_nonuse[y_nonuse==1] = 1 # N'utilise pas de contraception

"""
Explanation: Amélioration 4 : simplifier les catégories de la variable d'intérêt
On remarque que dans les matrices de confusion le fait de ne pas avoir de contraception est à chaque fois mieux prédit. Par ailleurs pour des raisons de politiques publiques on peut penser que le décideur est plus intéressé par le fait de ne pas avoir de contraception afin de mettre des politiques en place. En effet l'enjeu sociétal doit être de maximiser l'utilisation d'un moyen de contraception vs ne pas se protéger.
On se propose donc de regarder la classification suivant deux critères uniquement : avoir une contraception ou ne pas en avoir et donc transformer un problème de classification multi-classe en un de classification binaire. Cela nous permettra donc de pouvoir introduire la notion de courbe ROC valable uniquement pour une variable de prédiction binaire.
End of explanation
"""


from sklearn.cross_validation import train_test_split
data_train, data_test, y_nutrain, y_nutest = train_test_split(data, y_nonuse, test_size=0.33, random_state=42)
df.shape, data_test.shape, y_nutest.shape

# DecisionTree
y_nupredDec = DecTree.fit(data_train, y_nutrain).predict(data_test)
cm = confusion_matrix(y_nutest, y_nupredDec)
np.set_printoptions(precision=2)
print('Matrice de confusion')
print(cm)


f1nuDec = f1_score(y_nutest, y_nupredDec)
print("f1nuDec", f1nuDec)

# AdaBoost
y_nupredAda = AdaBoost.fit(data_train, y_nutrain).predict(data_test)
cm = confusion_matrix(y_nutest, y_nupredAda)
np.set_printoptions(precision=2)
print('Matrice de confusion')
print(cm)


f1nuAda = f1_score(y_nutest, y_nupredAda)
print("f1nuAda", f1nuAda)

#KNeighbors
y_prednuKN= KN.fit(data_train, y_nutrain).predict(data_test)

cm = confusion_matrix(y_nutest, y_prednuKN)
np.set_printoptions(precision=2)
print('Matrice de confusion')
print(cm)

f1nuKN = f1_score(y_nutest, y_prednuKN)
print('f1nuKN',f1nuKN)

# RandomForest
y_nupredRF = RF.fit(data_train, y_nutrain).predict(data_test)
cm = confusion_matrix(y_nutest, y_nupredRF)
np.set_printoptions(precision=2)
print('Matrice de confusion')
print(cm)


f1nuRF = f1_score(y_nutest, y_nupredRF)
print("f1nuRF", f1nuRF)

"""
Explanation: Test sur la nouvelle variable de prédiction
Nous avons choisi de faire tourner les modèles déjà envisagés : DecisionTree, RandomForest, KNeighbors et AdaBoost.
End of explanation
"""


#Comparaison des scores
l = [ {"y":"Multiclasse", "DecTree":f1Dec, "AdaBoost":f1Ada, "Kneighbors":f1KNeigh, "RandomForest": f1RF},
      {"y":"Binaire", "DecTree":f1nuDec, "AdaBoost":f1nuAda, "Kneighbors":f1nuKN, "RandomForest":f1nuRF}
        ]
f1=pd.DataFrame(l)
f1 = f1.set_index("y")
f1

"""
Explanation: Comparaison des résultats
On met dans un tableau les scores obtenus en classant suivant la variable y (Contraception) à trois modalités (Non-Use, Short-Term, Long-Term) et celle à deux modalités (Use, Non-Use)
End of explanation
"""


#Comparaison des courbes ROC et AUC
from sklearn.metrics import roc_curve, auc
probasRF = RF.predict_proba(data_test)
probasDec = DecTree.predict_proba(data_test)
probasAda = AdaBoost.predict_proba(data_test)

fpr, tpr, thresholds = roc_curve(y_nutest, probasRF[:, 1])
fpr1, tpr1, thresholds1 = roc_curve(y_nutest, probasDec[:, 1])
fpr2, tpr2, thresholds2 = roc_curve(y_nutest, probasAda[:, 1])
roc_auc = auc(fpr, tpr)
roc_auc1 = auc(fpr1, tpr1)
roc_auc2 = auc(fpr2, tpr2)

plt.plot(fpr, tpr, label=' Courbe ROC  RF (auc = %0.2f)' % roc_auc)
plt.plot(fpr1, tpr1, label=' Courbe ROC  DecTree (auc = %0.2f)' % roc_auc1)
plt.plot(fpr2, tpr2, label=' Courbe ROC  Ada (auc = %0.2f)' % roc_auc2)

plt.plot([0, 1], [0, 1], 'k--')
plt.xlim([0.0, 1.0])
plt.ylim([0.0, 1.0])
plt.xlabel("Taux d'erreur (faux positif)")
plt.ylabel("Taux de bonne réponse (Vrai positif)")
plt.title('Comparaison de la performance des classifiers')
plt.legend(loc="lower right")

"""
Explanation: Ainsi cela corrobore notre intuition, tous les modèles sont plus performants lorsqu'on leur soumet la variable binaire d'utilisation ou non d'un moyen de contraception, quel qu'il soit. Cependant  ce qui est plus particulièrement étonnant est le changement dans le classement de la performance relative des méthodes. En effet si on avait, en prenant comme indicateur les f1_scores, les préférences suivantes : KN>Rf>DecTRee>AdaBoost, nous avons maintenant cette nette préférence : DecTree>RF=AdaBoost>KN.
Nous allons maintenant nous concentrer sur une représentation plus compète de la performance et par la même occasion éliminer le KNeighborClassifier parce qu'il y a un score nettement inférieur aux autres. Pour cela nous allons traiter les modèles avec la variable binaire mais plutôt que de considérer leur score, nous allons regarder les courbes ROC.
End of explanation
"""


#Distribution des scores du DecisionTree
from sklearn.metrics import auc, precision_recall_curve

y_testn = y_nutest.as_matrix()
y_minD = y_nupredDec.min()

y_scoreD = array( [probasDec[i,p-y_minD] for i,p in enumerate(y_nupredDec)] )
y_scoreD[:5]

positive_scores = y_scoreD[y_testn == y_nupredDec]
negative_scores = y_scoreD[y_testn != y_nupredDec]

fpr = dict()
tpr = dict()
roc_auc = dict()
nb_obs = dict()

for i in DecTree.classes_:
    fpr[i], tpr[i], _ = roc_curve(y_nutest == i, y_scoreD)
    roc_auc[i] = auc(fpr[i], tpr[i])
    nb_obs[i] = (y_nutest == i).sum()

roc_auc, nb_obs

ax = seaborn.distplot(positive_scores, rug=True, hist=True, label="+")
seaborn.distplot(negative_scores, rug=True, hist=True, ax=ax, label="-")
ax.legend()

#Distribution des scores de l'AdaBoost

y_testn = y_nutest.as_matrix()
y_minD = y_nupredAda.min()

y_scoreD = array( [probasAda[i,p-y_minD] for i,p in enumerate(y_nupredAda)] )
y_scoreD[:5]

positive_scores = y_scoreD[y_testn == y_nupredAda]
negative_scores = y_scoreD[y_testn != y_nupredAda]

fpr = dict()
tpr = dict()
roc_auc = dict()
nb_obs = dict()

for i in AdaBoost.classes_:
    fpr[i], tpr[i], _ = roc_curve(y_nutest == i, y_scoreD)
    roc_auc[i] = auc(fpr[i], tpr[i])
    nb_obs[i] = (y_nutest == i).sum()

roc_auc, nb_obs

ax = seaborn.distplot(positive_scores, rug=True, hist=True, label="+")
seaborn.distplot(negative_scores, rug=True, hist=True, ax=ax, label="-")
ax.legend()

"""
Explanation: Nous avons confirmation que le DecisionTree et l'AdaBoost ont des meilleurs prédictions. En effet la probabilité que le score d’une bonne réponse soit supérieure au score d’une mauvaise réponse est dans les deux strictement supérieure à 0.76 (cf auc)
Limite des deux modèles les plus performants : manque de discrimination
Après avoir mis en valeur par comparaison que l'AdaBoost et le DecisionTree sont les meilleurs classifieurs, nous désirons lire de plus près la distribution du score
End of explanation
"""


tic=timeit.default_timer()
from sklearn.ensemble import GradientBoostingClassifier
GB = GradientBoostingClassifier()
y_prednuGB = GB.fit(data_train, y_nutrain).predict(data_test)
toc=timeit.default_timer()
toc - tic 

cm = confusion_matrix(y_nutest, y_prednuGB)
plot_confusion_matrix(cm)
f1nuGB = f1_score(y_nutest, y_prednuGB)
print('f1nuGB',f1nuGB)
print('f1nuDec',f1nuDec)  

probasnuGB = GB.predict_proba(data_test)
fpr, tpr, thresholds = roc_curve(y_nutest, probasnuGB[:, 1])
roc_auc = auc(fpr, tpr)
print ("Area under the ROC curve GradientBoost : %f" % roc_auc)
plt.plot(fpr, tpr, label='ROC curve GradientBoost (area = %0.2f)' % roc_auc)
plt.plot([0, 1], [0, 1], 'k--')
plt.xlim([0.0, 1.0])
plt.ylim([0.0, 1.0])
fpr1, tpr1, thresholds1 = roc_curve(y_nutest, probasDec[:, 1])

# ajout pour comparer des informations obtenues avec un DecisionTree
roc_auc1 = auc(fpr1, tpr1)
print ("Area under the ROC curve DecTree: %f" % roc_auc1)
plt.plot(fpr1, tpr1, label=' Courbe ROC  DecTree (auc = %0.2f)' % roc_auc1)

plt.xlabel("Taux d'erreur (faux positif)")
plt.ylabel("Taux de bonne réponse (Vrai positif)")
plt.title('Performance comparée du GradientBoosting')
plt.legend(loc="lower right")

#Distribution des scores de l'AdaBoost

y_testn = y_nutest.as_matrix()
y_minD = y_prednuGB.min()

y_scoreD = array( [probasnuGB[i,p-y_minD] for i,p in enumerate(y_prednuGB)] )
y_scoreD[:5]

positive_scores = y_scoreD[y_testn == y_prednuGB]
negative_scores = y_scoreD[y_testn != y_prednuGB]

fpr = dict()
tpr = dict()
roc_auc = dict()
nb_obs = dict()

for i in GB.classes_:
    fpr[i], tpr[i], _ = roc_curve(y_nutest == i, y_scoreD)
    roc_auc[i] = auc(fpr[i], tpr[i])
    nb_obs[i] = (y_nutest == i).sum()

roc_auc, nb_obs

ax = seaborn.distplot(positive_scores, rug=True, hist=True, label="+")
seaborn.distplot(negative_scores, rug=True, hist=True, ax=ax, label="-")
ax.legend()

"""
Explanation: Le score des deux modèles ne sont pas des plus discriminants puisqu’il existe une aire commune entre les bonnes et les mauvaises réponses plutôt importante. 
Nous allons donc introduire un nouveau modèle, le GradientBoost, en sélectionnant toutefois pour la comparaison le DecisionTree qui a la meilleur performance de prédiction
Amélioration 5 : apport d'un nouveau modèle et analyse des features
Le GradientBoostingClassifier
End of explanation
"""


feature_name = list(data.T.index)
limit = 30
feature_importance = GB.feature_importances_[:30]
feature_importance = 100.0 * (feature_importance / feature_importance.max())
sorted_idx = np.argsort(feature_importance)
pos = np.arange(sorted_idx.shape[0]) + .5
plt.subplot(1, 2, 2)
plt.barh(pos, feature_importance, align='center')
plt.yticks(pos, feature_name)
plt.xlabel('Importance relative')
plt.title("Classement des variables selon leur importance")

"""
Explanation: Ainsi on observe une nette amélioration de la prédiction. En effet la probabilité que le score d’une bonne réponse soit supérieure au score d’une mauvaise réponse est avec le GradientBoosting quasiment de 80% (cf.AUC). Par ailleurs sa discrimination est plus importante que les deux classifieurs précédents. Cela nous encourage à utiliser cette méthode pour appréhender l'importance relative des composantes, plutôt que les RandomForest par exemple.
Importance des features
End of explanation
"""


#Discrétisation de l'âge
u = df.copy()
u["Age"]= pandas.qcut(u.Age,10)# le choix du nombre de catégories s'est effectuée selon l'allure de l'histogramme préenté au début
# Nous choisission délibérément qcut parce que proche de la distribution initiale du fait la prise en compte des quantiles
u.head()

#Application de mosaic
from statsmodels.graphics.mosaicplot import mosaic

temp1 = pandas.crosstab([u.Contraception],[u.Age])

props = lambda key: {'color':'grey' if '1' in key else 'white'}#on met en valeur plus particulièrement le fait 
#d'avoir une contraception vs ne pas en avoir puisque nous avons déjà remarqué que cela avait plus de poids dans nos données

mosaic(temp1.unstack(),properties=props, title="Répartition de la non-contraception en fonction de l'âge" )

# Défaut : 
# 1)je n'arrive pas à me débarasser de ce qui est inscrit à l'intérieur des tiles même avec une lambda function comme labellizer.
# 2)On obtient aussi l'ensemble des coordonnées certes intéressant mais qui alourdisent la lecture

"""
Explanation: Ce graphique met en évidence l'importance de quatre variables dans la prédiction : le nombre d'enfant, l'éducation du mari, l'éducation de la femme, et son âge, variable la plus prégnante.
Focus sur la relation entre l'âge et la contraception
On peut donc s'interroger sur la distribution des femmes utilisant une contraception selon leur âge. Pour se faire on  veut représenter un grpah de type "mosaic plot". Il faut donc dans un premier temps discrétiser la variable âge pour rendre le graphique lisible. Dans un deuxième temps pour pouvoir appliquer l'objet "mosaic" il faut considérer le tableau de contingence de ces variables.
End of explanation
"""


from sklearn.feature_selection import SelectKBest
from sklearn.feature_selection import chi2

data_new = SelectKBest(chi2, k=2).fit_transform(data, y)
data_new.shape

# Création de datared et des bases d'apprentissages et de test correspondantes
u= data.copy()
datared = u[["Age","Nbr Enfant", "Educ", "Educ mari", "Niveau de vie"]]
datared.head()
datared_train, datared_test, y_nutrain, y_nutest = train_test_split(datared, y_nonuse, test_size=0.33, random_state=42)
print(datared_train.shape)
print(datared_test.shape)
datared.head()

"""
Explanation: Nous pouvons voir donc qu'il y a clairement une tendance comportementale qui varie avec l'âge : ainsi après 37 ans les femmes déclareraient moins utiliser une contraception.
Amélioration 6 : réduire le nombre de variables prédictives
L'idée est d'adapter la méthode qui consiste à utiliser une ACP afin de combiner les facteurs pour augmenter la performance, à un modèle avec des variables catégorielles.
Nous allons utiliser une méthode bien plus simple qui est de réduire de le nombre de variables prédictives par la suppression des variables ayant le moins d'importance, mises en évidence par le GradientBoost.
Ainsi cinq variables ont une portée prédictive plus importante que les autres: l'âge, le nombre d'enfant, l'éducation, le niveau de vie et l'éduction du mari.
On choisit de tester ce jeu de données réduit avec à la fois de le DecisionTree et le GradientBoost afin d'avoir une idée de l'effet de ce changement sur la prédiction.
Récupération des variables à plus forte valeur prédicitive
End of explanation
"""


y_prednuredDT = DecTree.fit(datared_train, y_nutrain).predict(datared_test)
cm = confusion_matrix(y_nutest, y_prednuredDT)
plot_confusion_matrix(cm)
f1nuredDT = f1_score(y_nutest, y_prednuredDT)
print('f1nuredDec',f1nuredDT)
probasnuredDT = DecTree.predict_proba(datared_test)
fpr, tpr, thresholds = roc_curve(y_nutest, probasnuredDT[:, 1])
roc_auc = auc(fpr, tpr)
print ("Area under the ROC curve DecTReered : %f" % roc_auc)
plt.plot(fpr, tpr, label='ROC curve DecTreered(area = %0.2f)' % roc_auc)
plt.plot([0, 1], [0, 1], 'k--')
plt.xlim([0.0, 1.0])
plt.ylim([0.0, 1.0])

#Comparaison avec le modèle non réduit
print('f1nuDec',f1nuDec)
fpr1, tpr1, thresholds1 = roc_curve(y_nutest, probasDec[:, 1])
roc_auc1 = auc(fpr1, tpr1)
print ("Area under the ROC curve DecTree : %f" % roc_auc1)
plt.plot(fpr1, tpr1, label='ROC curve DecTree (area = %0.2f)' % roc_auc1)

plt.xlabel("Taux d'erreur (faux positif)")
plt.ylabel("Taux de bonne réponse (Vrai positif)")
plt.title('Performance comparée du modèle réduit avec DecisionTree')
plt.legend(loc="lower right")

"""
Explanation: DecTree sur modèle réduit
End of explanation
"""


y_prednuredGB = GB.fit(datared_train, y_nutrain).predict(datared_test)
cm = confusion_matrix(y_nutest, y_prednuredGB)
plot_confusion_matrix(cm)
f1nuredGB = f1_score(y_nutest, y_prednuredGB)
print('f1nuredGB',f1nuredGB)
print('f1nuGB',f1nuGB)
probasnuredGB = GB.predict_proba(datared_test)
fpr, tpr, thresholds = roc_curve(y_nutest, probasnuredGB[:, 1])
roc_auc2 = auc(fpr2, tpr2)
print ("Area under the ROC curve GBred : %f" % roc_auc2)
plt.plot(fpr, tpr, label='ROC curve (area = %0.2f)' % roc_auc2)
plt.plot([0, 1], [0, 1], 'k--')
plt.xlim([0.0, 1.0])
plt.ylim([0.0, 1.0])

#Comparaison avec le modèle non réduit

fpr3, tpr3, thresholds3 = roc_curve(y_nutest, probasnuGB[:, 1])
roc_auc3 = auc(fpr3, tpr3)
print ("Area under the ROC curve GB : %f" % roc_auc3)
plt.plot(fpr3, tpr3, label='ROC curve GradientBoost (area = %0.2f)' % roc_auc3)

plt.xlabel("Taux d'erreur (faux positif)")
plt.ylabel("Taux de bonne réponse (Vrai positif)")
plt.title('Performance comparée du modèle réduit avec GradientBoost')
plt.legend(loc="lower right")

"""
Explanation: GradientBoost sur modèle réduit
End of explanation
"""


