import pandas as pan #pour les données :gestion des data.frame
import numpy as np
import matplotlib.pyplot as plt
import scipy.stats.mstats as ms


df=pan.read_csv("http://math.univ-lyon1.fr/homes-www/dabrowski/nutriage.csv",sep="\t")

df=pan.read_csv("https://math.univ-lyon1.fr/~dabrowski/nutriage.csv",sep="\t")




for nom in df.keys():
    globals()[nom] = df[nom]




   ###################Exercice 1#########################
x=np.array([1,8,5,1])
y=np.concatenate(([0],1+2*np.array(range(5))))
##autre solution
y=np.array([0,1,3,5,7,9])
   
   #1 
   plt.plot(x,y)
#erreur car x et y n’ont pas même longueur.
   
   #2
 x=np.concatenate((x,[3,5]))

plt.clf()
np.cor
fig,ax=plt.subplots()
ax.plot(x,y,'.')
plt.plot(x,y,'.')
#Autre solution
plt.scatter(x,y)
#3
fig,ax=plt.subplots()
ax.plot(x,y,'.')
ax.plot(np.mean(x),np.mean(y),'r+')

#4
a,b=ms.linregress(x,y)[:2]
print("La droite de régression est y=%s x + %s" % (a,b))

help(plt.axline)

fig,ax=plt.subplots()
ax.plot(x,y,'.')
ax.plot(np.mean(x),np.mean(y),'r+')
ax.axline(xy1=(0,b),slope=a,color="green")
fig.suptitle("Droite de régression de y en fonction de x")

#Autre solution avec axline

fig,ax=plt.subplots()
ax.plot(x,y,'.')
ax.plot(np.mean(x),np.mean(y),'r+')
ax.axline(xy1=(0,b),xy2=(np.mean(x),np.mean(y)),color="green")
fig.suptitle("Droite de régression de y en fonction de x")
np.cor

#Autre solution sans axline

def abline(slope, intercept,axes,col="red"):
    x_vals = np.array(axes.get_xlim())
    y_vals = intercept + slope * x_vals
    axes.plot(x_vals, y_vals, '-',color=col)


fig,ax=plt.subplots()
ax.plot(x,y,'.')
ax.plot(np.mean(x),np.mean(y),'r+')
abline(a,b,ax,col="green")
fig.suptitle("Droite de régression de y en fonction de x")

#Exemple pour comprendre fonctionnement de la fonction
def abline(slope, intercept,axes,col="red"):
    x_vals = np.array(axes.get_xlim())
    y_vals = intercept + slope * x_vals
    axes.plot(x_vals, y_vals, '-',color=col)
    
fig,ax=plt.subplots()
ax.plot((0,1),(0,1),'.')
abline(1,0,ax,col="green")



#6
np.corrcoef(x,y)[0,1]#=-0.03873023356985952 c'est très 
#proche de 0 donc l'approximation par la droite de
# régression n'est pas pertinente

   ###################Exercice 2#########################


#1
np.mean(poids)#moyenne 66.48230088495575
np.var(poids)#variance  144.1611911661055
np.var(poids,ddof=1)#variance sans biais 144.8019075712882 (plus grand)
np.std(poids,ddof=1)#écart-type sans biais 12.033366427201003
np.median(poids)#médiane 66.0

np.mean(taille)#moyenne
np.var(taille)#variance empirique
np.var(taille,ddof=1)#variance empirique non biaisée
np.std(taille,ddof=1)#écart-type
np.median(taille)#médiane


#2
np.quantile(poids, [0.25,0.5,0.75],interpolation="lower")#python 3.9

np.quantile(poids, [0.25,0.5,0.75],method="lower")#python 3.11

#l'argument interpolation="lower" est indispensable pour obtenir la définition du cours, 
# la version par défaut est similaire à la définition de la médiane

#3
help(pan.DataFrame)
tp=pan.DataFrame({"taille" : taille, "poids": poids})
#ne marche pas : pan.DataFrame([taille,poids])
#Autre solution:
pan.DataFrame(np.transpose(np.array([taille,poids])))
#4
tp.describe()#résumé des variables stats: moyenne,
# écart-type, quantile
#           taille       poids
#count  226.000000  226.000000
#mean   163.960177   66.482301
#std      9.003368   12.033366
#min    140.000000   38.000000
#25%    157.000000   57.250000
#50%    163.000000   66.000000
#75%    170.000000   75.000000
#max    188.000000   96.000000

#Remarquer que ce ne sont pas le même 1er quartile que pour le cours. 


#5
x=taille
y=poids
fig,ax=plt.subplots()
ax.plot(x,y,'b.')
ax.plot(np.mean(x),np.mean(y),'r+')
ax.set_xlim((135,195))
ax.set_ylim((35,100))

#6

a,b=ms.linregress(x,y)[:2]
print("la droite de Régression est y=",a,"x+",b)
ax.axline((0,b), slope=a,color='g')
fig.suptitle("Droite de régression de y en fonction de x")

#autre solution
ax.axline((np.mean(x),np.mean(y)), slope=a,color='g')

#Dernière solution
ax.axline(((0,b), (np.mean(x),np.mean(y)),color='g')



#Un point au dessus de la droite de régression est une personne dont 
#le poids est au dessus de celui attendu pour sa taille.
#En langage courant, on dirait au "dessus de la moyenne", même si ce n'est pas exactement une valeur moyenne
#que donne la droite de régression.

 ###################Exercice 3#########################
 
#1
np.mean(age[sexe=='F'])#74.34042553191489

#2
y1=poids[sexe=='F']
x1=age[sexe=='F']

y2=poids[sexe=='H']
x2=age[sexe=='H']
fig,ax=plt.subplots()
ax.plot(x1,y1,'r.')
ax.plot(x2,y2,'b.')
ax.set_xlim((60,95))
ax.set_ylim((35,100))

#3
a1,b1=ms.linregress(x1,y1)[:2]
print("la droite de Régression pour les femmes est y=",a1,"x+",b1)

a2,b2=ms.linregress(x2,y2)[:2]
print("la droite de Régression pour les hommes est y=",a2,"x+",b2)
ax.axline((0,b1), slope=a1,color='pink')
ax.axline((0,b2), slope=a2,color='cyan')
   ###################Exercice 4#########################

import statsmodels.api as sm

#1
X=sm.add_constant(taille)#besoin de donner vecteur constant pour avoir b
lmdf=sm.OLS(poids,X);resdf=lmdf.fit()
print(resdf.summary())

#On trouve a comme coeff en face de taille soit 0.8429 et b comme coeff en face de const: -71.7196
# Le résumé donne aussi l'incertitude sur a et b sous forme d'écart type, 
# de p-valeur et d'intervalle de confiance (cf. fin du cours)
# L'indication P>|t| 0 indique que les coefficients sont significativement non nulle
# Pour a, cela indique la pertinence de la régression.

#2
y= a* taille + b
fig,ax=plt.subplots()
fig,ax=plt.subplots()
ax.plot(taille,poids,'b.')
ax.plot(taille,resdf.fittedvalues,'g')

#3 

y= a* taille + b
np.var(y)
np.var(poids)
np.var(y)/np.var(poids)#0.3977289595467793

np.var(resdf.fittedvalues)/np.var(poids)
#On retrouve R-squared:                       0.398
#Plus c'est proche de 1, plus le nuage de point est proche de la droite de régression, et meilleur 
#est l'utilisation de la droite de régression seule pour approcher les points
#Plus il est petit et plus il y a d'aléas autour des données dans la direction perpendiculaire à la droite de régression

#Commentaire pour exo 3
X1=sm.add_constant(x1)#besoin de donner vecteur constant pour avoir b
lmdf1=sm.OLS(y1,X1);resdf1=lmdf1.fit()
print(resdf1.summary())
#0.5% de la variance capturée par la régression, P>|t| = 0.397 élévé > 0.05 pour a, donc signe d'une régression linéaire peu significative  

X2=sm.add_constant(x2)#besoin de donner vecteur constant pour avoir b
lmdf2=sm.OLS(y2,X2);resdf2=lmdf2.fit()
print(resdf2.summary())
#0.1% de la variance capturée par la régression, P>|t| = 0.816 élévé pour a, donc signe d'une régression linéaire peu significative  

#4
prediction=resdf.get_prediction().summary_frame(alpha=0.05)#1-alpha=0.95 niveau de confiance pour laprediction
fig, ax = plt.subplots(figsize=(8, 6))
ax.plot(taille, poids, "o", label="data")
ax.plot(taille,  prediction["mean"], label="OLS",color="orange")#droite de régression
ax.plot(taille, prediction["obs_ci_lower"], color="red")#borne inf de l'intervalle
ax.plot(taille, prediction["obs_ci_upper"], color="red")#borne sup de l'intervalle
ax.legend(loc="best")
fig.suptitle("Régression de poids par rapport à taille avec intervalle de prédiction (rouge)")
#facultatif
ax.plot(taille, prediction["mean_ci_lower"], color="green")#borne inf de l'incertitude sur la droite
ax.plot(taille, prediction["mean_ci_upper"], color="green")#borne sup de l'incertitude sur la droite

#voir les résultats calculés
prediction
prediction.keys()

   ###################Exercice 5#########################

Air=pan.read_csv("/home/Enseignement/Cours/StatInfo/TP2024/pollution.csv",sep="\t",na_values="-")

Air=pan.read_csv("https://math.univ-lyon1.fr/~dabrowski/pollution.csv",sep="\t",na_values="-")

#1
PO=pan.DataFrame({'PM10': Air['PM10'], 'NO2' : Air['NO2']    })
PO2=PO.dropna()

#2
y=np.log(PO2['PM10'])
x=np.log(PO2['NO2'])

X=sm.add_constant(x)#besoin de donner vecteur constant pour avoir b
lmdf=sm.OLS(y,X);resdf=lmdf.fit()
print(resdf.summary())
#Les coefficients sont significativement non nuls, donc la régression est pertinente
#26.5% de la variance des données est capturée par la régression.

prediction=resdf.get_prediction().summary_frame(alpha=0.05)#1-alpha=0.95 niveau de confiance pour laprediction
fig, ax = plt.subplots(figsize=(8, 6))
ax.plot(x, y, "o", label="data")
ax.plot(x,  prediction["mean"], label="OLS",color="orange")#droite de régression
ax.plot(x, prediction["obs_ci_lower"], color="red")#borne inf de l'intervalle
ax.plot(x, prediction["obs_ci_upper"], color="red")#borne sup de l'intervalle


#3
#1er essai
fig, ax = plt.subplots(figsize=(8, 6))
ax.plot(np.exp(x),np.exp( y), "o", label="data")
ax.plot(np.exp(x), np.exp(prediction["mean"]), label="OLS",color="orange")#droite de régression
ax.plot(np.exp(x), np.exp(prediction["obs_ci_lower"]), color="red")#borne inf de l'intervalle
ax.plot(np.exp(x), np.exp(prediction["obs_ci_upper"]), color="red")#borne sup de l'intervalle
#c'est insatisfaisant car les courbes ne sont pas parcourues dans l'ordre

#Solution
fig, ax = plt.subplots(figsize=(8, 6))
ax.plot(np.exp(x),np.exp( y), "o", label="data")
ax.plot(np.exp(np.sort(x)), np.exp(np.sort(prediction["mean"])), label="OLS",color="orange")#droite de régression
ax.plot(np.exp(np.sort(x)), np.exp(np.sort(prediction["obs_ci_lower"])), color="red")#borne inf de l'intervalle
ax.plot(np.exp(np.sort(x)), np.exp(np.sort(prediction["obs_ci_upper"])), color="red")#borne sup de l'intervalle



#Les deux polluants croissent simultanément (avec la pollution automobile) 
#mais l'intervalle d'incertitude est assez large (correspondant au 26.5% de la variance capturée, valeur assez faible dans le résumé),
# il est donc pertinent de mesurer les deux, qui ne s'accumulent pas dans l'air de la même manière.