import numpy as np
import scipy.optimize as spo
# données expérimentales
x = np.array([1.0,2.1,3.1,3.9])
y = np.array([2.1,3.9,6.1,7.8])
# incertitudes-types sur les données expérimentales
ux = np.array([0.2,0.3,0.2,0.2])
uy = np.array([0.3,0.3,0.2,0.3])

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

# 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)))

# calcul de la valeur du "chi2 réduit" pour les paramètres ajustés
chi2r = np.sum(np.square(residual(popt,y,x)))/(x.size-popt.size)