#-------------------------------------------------------------------------------
# Calcul d'incertitudes, loi normale, méthode Monte Carlo et régression linéaire
#-------------------------------------------------------------------------------
# Bibliothèques utilisées
import numpy as np
import scipy.optimize
from scipy.optimize import curve_fit
import matplotlib
import matplotlib.pyplot as plt

## Paramètres modifiables
N_exp = 1000 # nombre de tirages aléatoires (à modifier éventuellement)
n_pnt = 5 # nombre de points d'expérience (à modifier)

## Valeurs mesurées (attention, la taille des listes doit correspondre à n_pnt !)
valY = np.array([2,4,6,8,10]) # a remplacer par vos valeurs
valX = np.array([1,2,3,4,5]) # a remplacer par vos valeurs
sY = np.array([0.1,0.1,0.1,0.1,0.1]) # a remplacer par vos valeurs
sX = np.array([0.1,0.1,0.1,0.1,0.1]) # a remplacer par vos valeurs

## Définitions des fonctions utiles
def make_data(n) :
    tabY = np.zeros(n)
    tabX = np.zeros(n)

    for k in range(n) :
        tabY[k] = float(np.random.normal(valY[k],sY[k],1))
        tabX[k] = float(np.random.normal(valX[k],sX[k],1))

    return tabY,tabX

def func(x, a): # regression linéaire de la forme y = a*x (+b ?)
    return a*x #+b

## Tirages aléatoires et regressions linéaires
valA = np.zeros(N_exp)
# valB = np.zeros(N_exp)

for k in range(N_exp) :
    Y,X = make_data(n_pnt)
    popt, pcov = curve_fit(func, X, Y)
    valA[k] = float(popt[0])
    # valB[k] = float(popt[1]) # exemple de récupération d'autres paramètres de l'ajustement
    plt.plot(X,Y,'x') # tracé des points aléatoires
    plt.plot(X,func(X,*popt),'-') # tracé de chaque régression linéaire

## Calcul des indicateurs statistiques
Am = np.mean(valA) # moyenne de la liste
sigma = np.std(valA,ddof=1) # écart-type de la liste
sA = sigma/np.sqrt(N_exp) # incertitude-type prenant en compte le nombre de points expérimentaux
# Bm = np.mean(valA) # idem pour l'ordonée à l'origine éventuellement
# sigma = np.std(valB,ddof=1)
# sB = sigma/np.sqrt(N_exp)
print("Valeur moyenne pente : %2.4f"  %Am)
print("Incertitude-type pente : %2.4f"  %sA)
# print("Valeur moyenne oao : %2.4f"  %Bm)
# print("Incertitude type oao : %2.4f"  %sB)

# Tracés pour vérification (à commenter si execution trop longue)
plt.plot(valX,valY,'o',color='black',markersize=10) # tracé des points expérimentaux
plt.plot(valX,func(valX,Am),'--',color='black',linewidth=2) # ajustement "moyen"

plt.show()