# -*- coding: utf-8 -*-
"""
Regression lineaire avec :
affichage graphe avec barres d'erreur,
affichage residus,
calcul parametres ajustement lineaire,
calcul incertitude sur les parametres.

a changer : valeurs de x, y , u_x, u_y, N iterations
"""

## Bibliotheques
import numpy as np                 # Pour sqrt, log, cos...
import matplotlib.pyplot as plt    # graphes



## Valeurs et incertitudes de y
y = np.array([2.2,1.65,0.82,0.57,0.35])
u_y = 0.1       # exemple d'incertitude identique pour chaque mesure
# Valeurs et incertitudes de x
x=np.array([1.18, 1.01, 0.82, 0.74, 0.69])
u_x=np.array([0.0023, 0.0017, 0.0011, 0.0009, 0.0008]) # exemple d'incertitudes differentes



## Graphe avec barres d'erreur
plt.figure(figsize=(8,6))
plt.xlabel('x')
plt.ylabel('y')

plt.errorbar(x, y,xerr=u_x ,yerr=u_y,fmt='b+',zorder=2,label='Données')
plt.legend()
plt.show()



## Ajustement des donnees par une fonction affine
p = np.polyfit(x,y,1)         # calcul de l'ajustement
a = p[0]  # pente de la droite d'ajustement
b = p[1]  # ordonnee a l'origine de la droite d'ajustement

# Graphe avec droite d'ajustement
plt.figure(figsize=(8,6))
plt.xlabel('x')
plt.ylabel('y')

xfit = np.linspace(min(x),max(x),2)
plt.plot(xfit, b + a*xfit, 'r', label='Régression linéaire')
plt.errorbar(x, y,xerr=u_x ,yerr=u_y,fmt='b+',zorder=2,label='Données')
plt.legend()
plt.show()



## Verification validite de la regression lineaire avec les residus
plt.subplot(211)
res = (y-np.polyval(p,x))      # Residus
plt.plot(x,res,'ro')           # Graphe des residus
plt.errorbar(x, res,yerr=u_y,fmt='r+',zorder=2,label='Mesures')
plt.grid(True)
plt.ylabel('résidus')
plt.xlabel('x')
plt.subplot(212)
z=(y-np.polyval(p,x))/u_y      # Residus normalisés
plt.plot(x,z,'ro')             # Graphe des residus normalises
plt.grid(True)
plt.ylabel('résidus normalisés')
plt.xlabel('x')
plt.show()

#1. sur les residus : la valeur 0 est bien dans les barres d'erreurs -> ajustement valide
#2. sur les residus normalises : z<2 pour tous les points -> ajustement valide



## Determination des incertitudes sur les parametres de modelisation par methode Monte-Carlo
N = 1000
L_a = [] # liste qui contiendra les valeurs de pente
L_b = [] # liste qui contiendra les valeurs d'ordonnee a l'origine

for i in range(0,N):
    l = len(x) 
    my = y + u_y*np.sqrt(3)*np.random.uniform(-1,1,l) #le sqrt(3) vient du fait qu'on prend une distribution uniforme sur la graduation de y / tirage de y
    mx = x + u_x*np.sqrt(3)*np.random.uniform(-1,1,l) #le sqrt(3) vient du fait qu'on prend une distribution uniforme sur la graduation de x / tirage de x
    p_MC=np.polyfit(mx,my,1) # Ajustement apres tirage aleatoire en chaque point
    L_a.append(p_MC[0]) 
    L_b.append(p_MC[1])

u_a = np.std(L_a,ddof=1) #incertitude sur la pente
u_b = np.std(L_b,ddof=1) #incertitude sur l'ordonnee a l'origine

print('a= {:.3e} +/- {:.1e} unité'.format(a,u_a),'et b={:.3e} +/- {:.1e} unité'.format(b,u_b))