#!/usr/bin/env python3
# -*- coding: utf-8 -*-

import numpy as np
import matplotlib.pyplot as plt


############### Données expérimentales   ###############

# Données mesurées



x=np.array([0,0.9,1.8,2.6,3.3,4.4,5.2,6.1,6.5,7.4]) 
y=np.array([5.9,5.4,4.4,4.6,3.5,3.7,2.8,2.8,2.4,1.5]) 


# Incertitudes type mesurées (elles peuvent éventuellement être nulles)
Delta_x = np.array([.03162277660,.03162277660,.04472135954,.03535533906,.07071067814,.1118033989,.1290994449,.2236067977,.7453559922,1.000000000])
Delta_y=np.array([1.000000000,.7453559922,.5000000000,.3535533906,.2236067977,.2236067977,.1195228609,.1195228609,.1000000000,.04472135954])

# Possibilité de changer les labels des axes
plt.figure(figsize=(8,6))
plt.xlabel('x')
plt.ylabel('y')



############### Code la régression ###################


N=10000

ta,tap=[],[]
tb,tbp=[],[]

for i in range(0,N):
    l = len(x) 
    mx=np.random.normal(x,Delta_x,l)
    my=np.random.normal(y,Delta_y,l)
    p=np.polyfit(mx,my,1)
    pp=np.polyfit(my,mx,1)
    ta.append(p[0])
    tb.append(p[1])
    tap.append(pp[0])
    tbp.append(pp[1])
    
a = np.mean(ta)
b = np.mean(tb)
ap = np.mean(tap)
bp = np.mean(tbp)

u_a = np.std(ta)
u_b = np.std(tb)
u_ap = np.std(tap)
u_bp = np.std(tbp)

cov=0
for i in range(N):
    cov=cov+(ta[i-1]-a)*(tb[i-1]-b)
cov=cov/N

covp=0
for i in range(N):
    covp=covp+(tap[i-1]-ap)*(tbp[i-1]-bp)
covp=covp/N



yfit = a*x + b

res = yfit-y 

xfit=ap*y+bp
resp=xfit-x

n=len(x)
pas=abs(max(x)-min(x))/n
x1=min(x)-pas
x2=max(x)+pas
y1=min(y)-abs(a*pas)
y2=max(y)+abs(a*pas)
xfit=np.linspace(x1,x2,n)
yfit=np.linspace(y1,y2,n)
ybestfit=-0.4823
plt.axis([x1,x2,y1,y2])
plt.plot(xfit, b + a*xfit, 'r', label='Régression linéaire y=ax+b')

plt.plot(ap*yfit+bp, yfit, 'y', label='Régression linéaire x=ap.y+bp')
plt.errorbar(x,y,xerr=Delta_x,yerr=Delta_y,fmt='b+',zorder=2,label='Mesures')
plt.legend()
plt.show()

print('pour y=ax+b:')
print ("Pente :", a ,", Incertitude type :", u_a)
print ("Ordonnée à l'origine :", b ,", Incertitude type :", u_b)
print("Covariance:",cov)

print('pour x=ap.x+bp:')
print ("Pente :", ap ,", Incertitude type :", u_ap)
print ("Ordonnée à l'origine :", bp ,", Incertitude type :", u_bp)
print("Covariance:",covp)

print('Comparaison: a=',a, '1/ap=',1/ap)
print('b=',b,' -bp/ap=',-bp/ap)