##  simulation Monte Carlo  pour evaluer l'incertitude sur les paramètres du modèle  ici régression linéaire  :  pente et ordonnée à l'origine


## Importation des bibliothèques utiles
import numpy as np                 # pour le traitement vectoriel des données
import matplotlib.pyplot as plt    # pour les tracés
import numpy.random as rd

## Saisie des valeurs expérimentales
v = np.array([* ,* ,*,*,*,*])      #  !!!!à compléter !!!!6 vitesses en mm/s
deltaf = np.array([*,*,*,*,*,*]) # !!!!! à compléter !!! 6 écarts fréquentiels  en Hz
u_v = np.array([0.1]*6)          # !!!!  6 incertitudes-type sur les vitesses  en mm/s
u_deltaf = np.array([*]*6)     # !!!!!  à compléter !!!! 6 incertitudes-type sur les écarts fréquentiels  en Hz




## Ajustement affine (polynome de degré 1)
X , Y = deltaf, v
p = np.polyfit(X, Y, 1) # p est un tableau contenant les coefficients de l'ajustement
                        # p[0] : pente ; p[1] : ordonnée à l'origine
Yreg = p[0]*deltaf + p[1]     # calcul des ordonnées des points "modèle"

## Affichage du résultat:
print("Pente = ", p[0])
print("Ordonnée à l'origine = ", p[1])


## Tracé
plt.errorbar(X, Y, xerr=u_deltaf, yerr=u_v, fmt='bo', label="points expérimentaux")
plt.plot(X, Yreg, 'c--', label="Modèle")
plt.xlabel("deltaf [Hz]")
plt.ylabel("v [mm/s]")
plt.legend()
plt.show()


## Simulation Monte-Carlo !

NMC = 10000                           # nombre d'expériences simulées soit 10000
Amc, Bmc = [], []                     # initialisation des listes dans lesquelles on va stocker
                                      # les valeurs de A (pente) et B (ordonnée à l'origine) correspondant à chaque expérience simulée

for i in range(NMC):                  # pour chaque expérience
  # on simule la mesure des 6 vitesses
  vmc = v + rd.normal(0, u_v, size=6)
  ##on simule la mesure des 6 écarts de fréquences
  deltafmc = deltaf + rd.normal(0,u_deltaf,size=6)
  # on reprend l'ajustement affine avec cette série de valeurs pour v et deltaf
  p = np.polyfit(deltafmc, vmc, 1)
  # on stocke les valeurs des paramètres d'ajustement dans les listes Amc et Bmc
  Amc.append(p[0])
  Bmc.append(p[1])


## Analyse statistique de l'ensemble des valeurs de A et de B obtenues
Amoy = np.mean(Amc)
u_A = np.std(Amc, ddof = 1)
print("Estimation de la célérité :")
print(f" Valeur mesurée : c/80 = {Amoy}")
print(f" Incertitude-type : u(c)/80 = {u_A}\n")

Bmoy = np.mean(Bmc)
u_B = np.std(Bmc, ddof = 1)
print("Estimation de l'ordonnée à l'origine :")
print(f" Valeur mesurée : B = {Bmoy}")
print(f" Incertitude-type : u(B) = {u_B}")