%matplotlib inline

from jyquickhelper import add_notebook_menu
add_notebook_menu()

"""
Explanation: Evaluation Python année 2016-2017 - solution
Le répertoire data contient deux fichiers csv simulés aléatoirement dont il faudra se servir pour répondre aux 10 questions qui suivent. Chaque question vaut deux points. Le travail est à rendre pour le lundi 20 février sous la forme d'un notebook envoyé en pièce jointe d'un mail.
End of explanation
"""


import pandas
persons = pandas.read_csv("https://raw.githubusercontent.com/sdpython/actuariat_python/master/_doc/notebooks/examen/data/persons.txt",
                     sep="\t")
persons.head()                     

rend = pandas.read_csv("https://raw.githubusercontent.com/sdpython/actuariat_python/master/_doc/notebooks/examen/data/rendezvous.txt",
                     sep="\t")
rend.head()  

persons.shape, rend.shape

"""
Explanation: 1
Deux fichiers sont extraits de la base de données d'un médecin. Un fichier contient des informations sur des personnes, un autre sur les rendez-vous pris par ces personnes. Quels sont-ils ?
End of explanation
"""


gr = rend.groupby("idc").mean()
gr.head()

"""
Explanation: Le second fichier est plus volumineux et contient une variable price qui ne peut pas être reliée aux personnes. Le premier fichier est celui des personnes, le second celui des rendez-vous. La variable idc est présente dans les deux tables. C'est elle qui identifie les personnes dans les deux tables.
2
On souhaite étudier la relation entre le prix moyen payé par une personne, son âge et son genre. Calculer le prix moyen payé par une personne ?
La table des rendez-vous contient toutes l'information nécessaire. La question était un peu ambiguë. On peut déterminer le prix moyen payé par personne, ou le prix moyen des prix moyens... On va répondre à la première option car la seconde n'a pas beaucoup d'intérêt et c'est très proche du prix moyen par rendez-vous, ce que la question aurait sans doute formulé dans ce sens si telle avait été son intention. On groupe par idc et on fait la moyenne.
End of explanation
"""


join = persons.merge(gr.reset_index(), on="idc")
join.head()

"""
Explanation: 3
Faire la jointure entre les deux tables.
End of explanation
"""


join.shape

join = persons.merge(gr.reset_index(), on="idc", how="outer")
join.shape

"""
Explanation: Cette jointure est assez simple puisque la colonne partagée porte le même nom dans les deux tables. On peut néanmoins se poser la question de savoir s'il y a des personnes qui n'ont pas de rendez-vous associé et réciproquement.
End of explanation
"""


join.plot(x="age", y="price", kind="scatter")

"""
Explanation: Visiblement, ce n'est pas le cas puisqu'une jointure incluant les éléments sans correspondances dans les deux tables n'ajoute pas plus d'éléments à la jointure.
4
Tracer deux nuages de points (age, prix moyen) et (genre, prix moyen) ?
End of explanation
"""


import seaborn
g = seaborn.jointplot("age", "price", data=join, kind="reg", size=7, scatter_kws={"s": 10})

join.plot(x="gender", y="price", kind="scatter")

"""
Explanation: On peut aussi utiliser un module comme seaborn qui propose des dessins intéressants pour un statisticatien.
End of explanation
"""


import numpy
bruit = join.copy()
bruit["gx"] = bruit.gender + numpy.random.random(bruit.shape[0])/3
bruit.plot(x="gx", y="price", kind="scatter")

"""
Explanation: On ne voit pas grand chose sur ce second graphe. Une option est d'ajouter un bruit aléatoire sur le genre pour éclater le nuage.
End of explanation
"""


join[["price", "gender"]].boxplot(by="gender")

"""
Explanation: Il n'y a rien de flagrant. On peut faire un graphe moustache.
End of explanation
"""


seaborn.violinplot(x="gender", y="price", data=join, inner="quart")

"""
Explanation: C'est mieux. Un dernier. Le diagramme violon, plus complet que le précédent.
End of explanation
"""


from sklearn.linear_model import LinearRegression
lr = LinearRegression()
lr.fit(join[["age", "gender"]], join["price"])

lr.coef_

lr.intercept_

"""
Explanation: 5
Calculer les coefficients de la régression $prix_moyen \sim age + genre$.
Une régression. Le premier réflexe est scikit-learn.
End of explanation
"""


from statsmodels.formula.api import ols
lr = ols("price ~ age + gender", data=join)
res = lr.fit()
res.summary()

"""
Explanation: On utilise maintenant statsmodels qui est plus complet pour toute ce qui est un modèle linéaire.
End of explanation
"""


from statsmodels.api import OLS
join2 = join.copy()
join2["cst"] = 1
lr = OLS(join2["price"], join2[["age", "gender", "cst"]])
res = lr.fit()
res.summary()

"""
Explanation: Ou encore (après avoir ajouté une constante).
End of explanation
"""


from patsy import dmatrices
y, X = dmatrices("price ~ age + gender" , join, return_type="matrix")
y = numpy.ravel(y)
lr = LinearRegression(fit_intercept=False)
lr.fit(X, y)

lr.coef_, lr.intercept_

X[:2]

"""
Explanation: On peut aussi définir la régression sous la forme de formule avec le module patsy.
End of explanation
"""


rend["date2"] = pandas.to_datetime(rend.date)
rend.dtypes

"""
Explanation: 6
On souhaite étudier le prix d'une consultation en fonction du jour de la semaine. Ajouter une colonne dans la table de votre choix avec le jour de la semaine.
On convertit d'abord la colonne date (chaîne de caractères au format date) avec la fonction to_datetime.
End of explanation
"""


rend["weekday"] = rend.date2.dt.weekday
rend.head()

"""
Explanation: Et on récupère le jour de la semaine avec weekday.
End of explanation
"""


rend[["price", "weekday"]].boxplot(by="weekday")

"""
Explanation: 7
Créer un graphe moustache qui permet de vérifier cette hypothèse.
On réutilise le code d'une question précédente.
End of explanation
"""


rend["constant"] = rend["date2"].min()
rend["jour"] = rend["date2"] - rend["constant"]
rend.head(n=2)

"""
Explanation: C'est clairement plus cher le dimanche.
8
Ajouter une colonne dans la table de votre choix qui contient 365 si c'est le premier rendez-vous, le nombre de jour écoulés depuis le précédent rendez-vous. On appelle cette colonne $delay$. On ajoute également la colonne $1/delay$.
Pour commencer, on convertit la date en nombre de jours depuis la première date.
End of explanation
"""


rend["jouri"] = rend.jour.apply(lambda d: d.days)
rend.head(n=2)

"""
Explanation: On convertit en entier.
End of explanation
"""


diff = rend.sort_values(["idc", "jouri"])["jouri"].diff()
rend["diff"] = diff
rend.head(n=2)

"""
Explanation: On trie par patient et jour puis on effectue la différence.
End of explanation
"""


first = rend[["idc", "date"]].groupby("idc", as_index=False).min()
first["j365"] = 365
first.head(n=2)

"""
Explanation: Il reste à traiter le premier jour ou plutôt le premier rendez-vous. On le récupère pour chaque patient.
End of explanation
"""


tout = rend.merge(first, on=["idc", "date"], how="outer")
tout[["idc", "jouri", "date", "j365"]].head(n=5)

"""
Explanation: Puis on fait une jointure.
End of explanation
"""


tout["delay"] = tout.j365.fillna(tout.jouri)
tout[["idc", "jouri", "date", "j365", "delay"]].head(n=8)

"""
Explanation: Il ne reste plus qu'à remplacer les NaN par jouri.
End of explanation
"""


tout["delay1"] = 1/ (tout["delay"] + 1)
tout[["delay", "delay1"]].head()

"""
Explanation: Finalement, il faut ajouter une colonne $1/delay$. Comme des patients ont parfois deux rendez-vous le même jour, pour éviter les valeurs nulles, on ajoute la colonne $1/(1+delay)$. On aurait pu également pour éviter les valeurs nulles considérer le temps en secondes et non en jour entre deux rendez-vous.
End of explanation
"""


rend2 = rend.sort_values(["idc", "jouri"]).reset_index(drop=True).copy()
rend2["diff"] = rend2["jouri"].diff()
rend2.loc[rend2.idc != rend2.idc.shift(1), "diff"] = 365
rend2.head()

"""
Explanation: 8 - réponse plus courte
On garde l'idée de la fonction diff et on ajoute la fonction shift.
End of explanation
"""


mat = tout.merge(persons, on="idc")

"""
Explanation: 9
Calculer les coefficients de la régression $prix \sim age + genre + delay + 1/delay + jour_semaine$.
L'âge ne fait pas partie de la table tout. Il faut faire une jointure pour le récupérer.
End of explanation
"""


mat[["age", "gender", "delay", "delay1", "weekday", "price"]].corr()

seaborn.clustermap(mat[["age", "gender", "delay", "delay1", "weekday", "price"]].corr(), figsize=(5,5))

"""
Explanation: Ensuite retour à scikit-learn et plutôt le second statsmodels pour effectuer des tests sur les coefficients du modèle. On regarde d'abord les corrélations.
End of explanation
"""


seaborn.pairplot(mat[["age", "gender", "delay", "delay1", "weekday", "price"]],
                 plot_kws={"s": 10}, size=1)

"""
Explanation: Si le jeu de données n'est pas trop volumineux.
End of explanation
"""


feat = mat[["age", "gender", "delay", "delay1", "weekday", "price"]]
g = seaborn.PairGrid(feat.sort_values("price", ascending=False), x_vars=feat.columns[:-1], 
                     y_vars=["price"], size=5, aspect=.25)
g.map(seaborn.stripplot, size=3, orient="h", palette="Reds_r", edgecolor="gray")

"""
Explanation: Un dernier pour la route.
End of explanation
"""


lr = LinearRegression()
lr.fit(mat[["age", "gender", "delay", "delay1", "weekday"]], mat["price"])
lr.coef_

from statsmodels.formula.api import ols
lr = ols("price ~ age + gender + delay + delay1 + weekday", data=mat)
res = lr.fit()
res.summary()

"""
Explanation: Régression
End of explanation
"""


lr_moy = LinearRegression()
lr_moy.fit(join[["age", "gender"]], join["price"])
lr_moy.coef_, lr_moy.intercept_

pred_moy = lr_moy.predict(join[["age", "gender"]])
join["pred_moy"] = pred_moy
join.head()

"""
Explanation: 10
Comment comparer ce modèle avec le précédent ? Implémentez le calcul qui vous permet de répondre à cette question.
Nous pourrions comparer les coefficients $R^2$ (0.57, 0.61) des régressions pour savoir quelle est la meilleur excepté que celle-ci ne sont pas calculées sur les mêmes données. La comparaison n'a pas de sens et il serait dangeraux d'en tirer des conclusions. Les valeurs sont de plus très proches. Il faut comparer les prédictions. Dans le premier cas, on prédit le prix moyen. Dans le second, on prédit le prix d'une consultation. Il est alors possible de calculer une prédiction moyenne par patient et de comparer les erreurs de prédiction du prix moyen. D'un côté, la prédiction du prix moyen, de l'autre la prédiction du prix d'une consultation agrégé par patient.
End of explanation
"""


err1 = ((join.pred_moy - join.price)**2).sum() / join.shape[0]
err1

join.plot(x="price", y="pred_moy", kind="scatter")

lrc = LinearRegression()
feat = mat[["age", "gender", "delay", "delay1", "weekday", "price", "idc"]].copy()
lrc.fit(feat[["age", "gender", "delay", "delay1", "weekday"]], feat["price"])
lrc.coef_, lrc.intercept_

predc = lrc.predict(feat[["age", "gender", "delay", "delay1", "weekday"]])
feat["predc"] = predc
feat.head()

"""
Explanation: On calcule l'erreur.
End of explanation
"""


agg = feat[["idc","predc", "price"]].groupby("idc").mean()
agg.head()

err2 = ((agg.predc - agg.price)**2).sum() / agg.shape[0]
err2

"""
Explanation: On agrège les secondes prédictions.
End of explanation
"""


agg.plot(x="price", y="predc", kind="scatter")

"""
Explanation: Le second modèle est clairement meilleur.
End of explanation
"""


temp = join.sort_values("price").reset_index(drop=True).reset_index(drop=False)
temp.head(n=1)

temp.plot(x="index", y=["price", "pred_moy"])

temp2 = agg.sort_values("price").reset_index(drop=True).reset_index(drop=False)
temp2.head(n=1)

ax = temp.plot(x="index", y="price", figsize=(14,4), ylim=[60,120])
temp.plot(x="index", y="pred_moy", linewidth=1, ax=ax, ylim=[60,120])
temp2.plot(x="index", y="predc", ax=ax, linewidth=0.6, ylim=[60,120])

"""
Explanation: La seconde régression utilise une information dont on ne dispose pas au niveau agrégé : le jour de la semaine et un précédent graphe a clairement montré que c'était une variable importante. Un dernier graphe pour comparer les deux prédictions en montrant les prédictions triées par prix à prédire.
End of explanation
"""


temp["erreur"] = temp.pred_moy - temp.price
temp2["erreur2"] = temp2.predc - temp2.price
ax = temp[["erreur"]].sort_values("erreur").reset_index(drop=True).plot()
temp2[["erreur2"]].sort_values("erreur2").reset_index(drop=True).plot(ax=ax)
ax.plot([0,1000], [0,0], "r--")

"""
Explanation: C'est finalement un peu plus visible sur le graphe précédent (nuage de points) mais aussi un peu trompeur du fait de la superposition des points. Une dernière remarque. En machine learning, nous avons l'habitude d'apprendre un modèle sur une base d'apprentissage et de tester les prédictions sur une autre. Dans notre cas, nous avons appris et prédit sur la même base. Ce type de tester est évidemment plus fiable. Mais nous avons comparé ici deux erreurs d'apprentissage moyennes et c'est exactement ce que l'on fait lorsqu'on compare deux coefficients $R^2$.
Un dernier graphe pour la route obtenu en triant les erreurs par ordre croissant.
End of explanation
"""


dummies = pandas.get_dummies(mat.weekday)
dummies.columns=["lu", "ma", "me", "je", "ve", "sa", "di"]
dummies.head()

"""
Explanation: Le second modèle fait des erreurs moins importantes surtout côté négatif. Il sous-estime moins la bonne valeur.
Traitement spécifique de la variable catégorielle weekday
Le second modèle ne prend pas en compte le dimanche comme jour de la semaine. weekday est une variable catégorielle. Contrairment au genre, elle possède plus de deux modalités. Il serait intéressant de la traiter comme un ensemble de variable binaire et non une colonne de type entier.
End of explanation
"""


mat2 = pandas.concat([mat, dummies.drop("lu", axis=1)], axis=1)
mat2.head(n=1)

lr = ols("price ~ age + gender + delay + delay1 + ma + me + je + ve + sa + di", data=mat2)
res = lr.fit()
res.summary()

"""
Explanation: On supprime une modalité pour éviter d'avoir une matrice corrélée avec la constante et on ajoute des variables au modèle de régression.
End of explanation
"""


lrc = LinearRegression()
feat2 = mat2[["age", "gender", "delay", "delay1", "price", "idc", "di"]].copy()
lrc.fit(feat2[["age", "gender", "delay", "delay1", "di"]], feat2["price"])
lrc.coef_, lrc.intercept_

predc = lrc.predict(feat2[["age", "gender", "delay", "delay1", "di"]])
feat2["predc"] = predc
feat2.head()

agg2 = feat2[["idc","predc", "price"]].groupby("idc").mean()
agg2.head()

err2d = ((agg2.predc - agg2.price)**2).sum() / agg2.shape[0]
err2d

"""
Explanation: On vérifie que le coefficient pour dimanche n'est clairement significatif contrairement aux autres dont la probabilité d'être nul est élevée. Le médecin appliquerait une majoration de 20 euros le dimanche. Le coefficient $R^2$ est aussi nettement plus élevé. On construit les prédictions d'un second modèle en ne tenant compte que de la variable dimanche.
End of explanation
"""


agg2.plot(x="price", y="predc", kind="scatter")

temp2 = agg2.sort_values("price").reset_index(drop=True).reset_index(drop=False)
ax = temp2.plot(x="index", y="price", figsize=(14,4), ylim=[60,120])
temp2.plot(x="index", y="predc", linewidth=1, ax=ax, ylim=[60,120])

"""
Explanation: Nettement, nettement mieux.
End of explanation
"""


from category_encoders import PolynomialEncoder
encoder = PolynomialEncoder(cols=["weekday"])
encoder.fit(rend[["weekday"]], rend["price"])
pred = encoder.transform(rend[["weekday"]])
conc = pandas.concat([rend[["weekday"]], pred], axis=1)
conc.head()

"""
Explanation: Une variable catégorielle en une seule colonne ?
Un jeu de données peut rapidement croître s'il est étendu pour chaque variable catégorielle. On peut utiliser le module category_encoders ou statsmodels.
End of explanation
"""


copy = rend[["weekday", "price"]].copy()
gr = copy.groupby("weekday", as_index=False).mean()
gr

feat3 = mat[["age", "gender", "delay", "delay1", "price", "idc", "weekday"]]
feat3 = feat3.merge(gr, on="weekday", suffixes=("", "_label"))
feat3.head()

"""
Explanation: Ou associer la valeur de la cible à prédire pour chaque jour de la semaine.
End of explanation
"""


lr = ols("price ~ age + gender + delay + delay1 + weekday", data=feat3)
res = lr.fit()
res.summary()

"""
Explanation: On apprend avec la variable weekday :
End of explanation
"""


lr = ols("price ~ age + gender + delay + delay1 + price_label", data=feat3)
res = lr.fit()
res.summary()

"""
Explanation: Et on compare avec la variable price_label :
End of explanation
"""


