#!/usr/bin/python3
# -*- coding:Utf-8 -*

"""
Régression linéaire avec montecarlo

PSI* Lycée Pierre Corneille

"""

import numpy as np
import matplotlib.pyplot as plt

## Situation 1 : Y=aX (exemple pendule simple l= g*(T/2Pi)^2 )
X=np.array([1,2,3,4,5])
Y=np.array([2.04,4.01,5.95,7.97,10.02])
uX=np.array([0.01,0.01,.01,0.1,0.1])
uY=np.array(Y/100)

# On calcule Y/X => valeur de a
a = Y/X
# On visualise le graphe (permet de vérifier à l'oeil que la régression sera acceptable) :
plt.close()    # ferme une eventuelle fenetre graphique precedemment ouverte
plt.plot(X ,Y , 'k+')  # trace des points ( y en ordonnee , x en abscisse ) , croix noires
plt.xlabel('nom  axe des abscisses')   # pour legender l'axe des abscisses
plt.ylabel('nom axe des ordonnees')   # pour legender l'axe des ordonnees
plt.title('titre')   # pour indiquer le titre du graphe
plt.errorbar (X , Y , xerr = uX, yerr = uY , color = 'r' , fmt ='none' , capsize = 5)  # trace des barres d'incertitudes , en rouge
plt.show()  # ouvre la fenetre graphique

# N régressions linéaires pour avoir les incertitudes types sur a et b
N = 10000 # Nombre d'échantillons

listea=[]
for i in range(N):
    Xbruit =np.random.uniform(X-uX,X+uX)  # Données bruitées avec Monte-Carlo
    Ybruit =np.random.uniform(Y-uY,Y+uY)
    listea.append(np.mean(Ybruit/Xbruit))  # On fait la liste des pentes trouvées

anmoy = np.mean(listea)
anstd=np.std(listea)
print('Modèle : y = a x avec a =',anmoy,' et incertitude = ',anstd)


## Situation 2 : Y=aX+b (exemple l=l0 + m * (g/k))
X=np.array([1,2,3,4,5])
Y=np.array([2.54,4.51,6.45,8.47,10.52])
uX=np.array([0.01,0.01,.01,0.1,0.1])
uY=np.array(Y/100)

a,b=np.polyfit(X,Y,1) #tuple qui donne a et b
# On visualise le graphe (permet de vérifier à l'oeil que la régression sera acceptable) :
plt.close()    # ferme une eventuelle fenetre graphique precedemment ouverte
plt.plot(X ,Y , 'k+')  # trace des points ( y en ordonnee , x en abscisse ) , croix noires
plt.xlabel('nom  axe des abscisses')   # pour legender l'axe des abscisses
plt.ylabel('nom axe des ordonnees')   # pour legender l'axe des ordonnees
plt.title('titre')   # pour indiquer le titre du graphe
plt.errorbar (X , Y , xerr = uX, yerr = uY , color = 'r' , fmt ='none' , capsize = 5)  # trace des barres d'incertitudes , en rouge
plt.show()  # ouvre la fenetre graphique

# N régressions linéaires pour avoir les incertitudes types sur a et b
N = 10000 # Nombre d'échantillons
listea=[]
listeb=[]
for i in range(N):
    Xbruit = np.random.uniform(X-uX,X+uX)  # Données bruitées avec Monte-Carlo
    Ybruit = np.random.uniform(Y-uY,Y+uY)  # ici loi uniforme...
    a,b=np.polyfit(Xbruit,Ybruit,1)
    listea.append(a)  # On fait la liste des pentes trouvées
    listeb.append(b)

amoy = np.mean(listea)
astd=np.std(listea)
bmoy = np.mean(listeb)
bstd=np.std(listeb)
residus=(Y-amoy*X-bmoy)/uY
print('Modèle : y = ax+b avec a =',amoy,' b= ',bmoy)
print('Incertitudes = ',astd,bstd)
print('residus = ',residus)