Bonjour, j'essaie de comprendre comment faire une simulation Monte Carlo sur python. J'ai quelques bases sur ce langage de programmation mais j'ai vraiment du mal à comprendre le programme à partir de

N=10000

#Import des bibliothèques
import numpy as np
import matplotlib.pyplot as plt

# Données expérimentales - création de listes : t et A
t = np.array([1.466, 2.383, 2.717, 3, 3.366, 3.717, 4.180, 4.5, 4.850, 6, 7 , 8, 9, 10, 11, 12, 13.02, 14, 15, 16, 17, 18, 19, 20])   #temps en min
A = np.array([1.023, 0.960, 0.941, 0.925, 0.908, 0.883, 0.862, 0.848, 0.830, 0.777, 0.733, 0.690, 0.650, 0.612, 0.578, 0.546, 0.514, 0.487, 0.459, 0.434, 0.410, 0.386, 0.364, 0.343])# absorbance correspondante

# tracé de la courbe donnant A en fonction de t : 
plt.figure(1)
plt.plot(t,A,'x')         #nuage de points avec t en abscisse, A en ordonnée en forme de "." bleu
plt.xlabel("temps")       # légende sur l'axe des abscisses
plt.ylabel("Absorbance")       # légende sur l'axe des ordonnées
plt.grid()         		#fait apparaître un quadrillage
plt.title("Ordre 0 : A = f(t)")        #pour donner un titre au graphique
plt.show

# de même, ajouter les courbes permettant de vérifier l'ordre 1 et l'ordre 2




# Régression linéaire :
p0 = np.polyfit(t, A, 1)		# régression linéaire de A en fonction de t
p1 = np.polyfit(t,np.log(A),1 )		# régression linéaire pour l'ordre 1
p2 = np.polyfit(t,1/A,1 )		# régression linéaire pour l'ordre 2


# ajout des résidus pour l'ordre 0 : 
z = (A-np.polyval(p0,t))
plt.figure(4)
plt.plot(t, z,'bo')               
plt.xlabel('t')
plt.ylabel(r'Résidus normalisés de $A$' )
plt.title('Résidus - ordre 0')
plt.grid()
plt.show()

# de même, ajouter les résidus pour l'ordre 1 et l'ordre 2




#Détermination de l'incertitude-type sur les paramètres (a, b) de la droite de régression à l'ordre 1 :
N = 1000          #Nombre de simulations souhaitées
liste_a = []        #On crée une liste vide qui sera remplie avec les diverses valeurs de la pente a de la droite de régression
liste_b = []	    #On crée une liste vide qui sera remplie avec les diverses valeurs de l'ordonnée à l'origine b

deltaA = 0.01       # Demi-étendue associée à chaque observation de A

for i in range(0,N):   #Création d'une boucle for permettant de générer les N expériences virtuelles
    A_sim = np.random.uniform(A - deltaA,A + deltaA) #Génère des valeurs aléatoires de A selon une loi uniforme, avec une demi-étendue deltaA
    lnA_sim = np.log(A_sim)
    p_sim = np.polyfit(t, lnA_sim, 1)  #Régression linéaire de lnA_sim en fonction de t ; 
    liste_a.append(p_sim[0])  	   #ajout de la pente de p_sim à la fin de la liste
    liste_b.append(p_sim[1])  #ajout de l'ordonnée à l'origine de p_sim à la fin de la liste
    
#Calcul de la valeur moyenne et écart type :
a_moy = np.mean(liste_a)      # Moyenne arithmétique des valeurs calculées de a
u_a = np.std(liste_a, ddof=1) # Ecart-type des valeurs calculées de a
b_moy = np.mean(liste_b)     # Moyenne arithmétique des valeurs calculées de b
u_b = np.std(liste_b, ddof=1)	 # Ecart-type des valeurs calculées de b

# affichage des résultats
print('meilleur estimateur de la pente (L/mmol):',a_moy)
print('incertitude-type sur a (L/mmol):',u_a)

print('meilleur estimateur de l ordonnée à l origine (L/mmol):',b_moy)	# afficher le meilleur estimateur de b
print('incertitude-type sur b (L/mmol):',u_b)	# afficher l'incertitude type sur b