import numpy as np
from math import *
import scipy.optimize as spo
import matplotlib.pyplot as plt

# données expérimentales
x = np.array([sin(),sin(),sin(),sin(),sin(),sin(),sin(),sin(),sin()]) # à modifier avec vos données
y= np.array([sin(0*pi/180),sin(10*pi/180),sin(20*pi/180),sin(30*pi/180),sin(40*pi/180),sin(50*pi/180),sin(60*pi/180),sin(70*pi/180),sin(80*pi/180)]) # à modifier avec vos données
# incertitudes-types sur les données expérimentales
ux = np.array([0,0,0,0,0,0,0,0,0]) # à modifier avec vos données
uy = np.array([0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1]) # à modifier avec vos données

# fonction f décrivant la courbe à ajuster aux données
def f(x,p):
    n,b = p
    return n*x+b
# dérivée de la fonction f par rapport à la variable de contrôle x
def Dx_f(x,p):
    n,b = p
    return n
def affine (x, n, b ):
    return n * x + b

# fonction d'écart pondérée par les erreurs
def residual(p, y, x):
    return (y-f(x,p))/np.sqrt(uy**2 + (Dx_f(x,p)*ux)**2)

# estimation initiale des paramètres
# elle ne joue généralement aucun rôle
# néanmoins, le résultat de l'ajustement est parfois aberrant
# il faut alors choisir une meilleure estimation initiale
p0 = np.array([0,0])

# on utilise l'algorithme des moindres carrés non-linéaires
# disponible dans la biliothèque scipy (et indirectement la
# bibliothèque Fortran MINPACK qui implémente l'algorithme
# de Levenberg-Marquardt) pour déterminer le minimum voulu
result = spo.leastsq(residual, p0, args=(y, x), full_output=True)

# on obtient :
# les paramètres d'ajustement optimaux
popt = result[0];
# la matrice de variance-covariance estimée des paramètres
pcov = result[1];
# les incertitudes-types sur ces paramètres
upopt = np.sqrt(np.abs(np.diagonal(pcov)))
ymod=[]
for i in range(len(x)) :
    ymod.append(affine(x[i],popt[0],popt[1]))
plt.plot(x,y,marker= '.', linestyle='')
plt.errorbar(x, y, xerr = ux, yerr = uy,
  fmt = 'none', capsize = 5, ecolor = 'black', zorder = 1)
plt.plot(x,ymod)
plt.title('y = '+str(popt[0])+'x +' +str(popt[1]))
plt.xlabel("axe des abscisses")
plt.ylabel("axe des ordonnées")
print("l'incertitude-type sur la pente de la droite est u(n) = " + str(upopt[0]) )
print( "l'incertitude-type sur l'ordonnée à l'origine est u(b) = " + str(upopt[1]) )
plt.show()