#!/usr/bin/env python # coding: utf-8 # # **1. Incertitude-type composée avec Python** # # Imaginons que l'on mesure un angle $\theta$ et que l'on détermine sa loi de probabilité par une méthode de type B (on se restreint cette année au cas d'une loi uniforme continue ou bien normale). # # On souhaiterait déterminer l'incertitude-type sur la valeur de $\sin(\theta)$ mais malheureusement cela ne fait pas partie des deux cas de figure vus en cours (forme linéaire ou bien loi de puissance). La méthode de Monte Carlo nous offre, grâce à Python, un moyen simple de calculer $u(\sin(\theta))$. # # Cette méthode consiste à tirer aléatoirement un grand nombre de valeurs de $\theta$ (un million de fois par exemple) en suivant la loi de probabilité attendue, puis à calculer à chaque tirage la valeur de $\sin(\theta)$. # # La dispersion des valeurs de $\sin(\theta)$ peut-être représentée avec un histogramme et l'incertitude-type $u(\sin(\theta))$ évaluée en calculant l'écart-type de la liste des $\sin(\theta)$. # # Pour mettre en oeuvre cette méthode, il faut importer les bibliothèques ```numpy``` (pour les calculs mathématiques et les tirages aléatoires) et ```matplotlib.pyplot``` (pour le tracé de l'histogramme). On utilisera le module ```random``` de ```numpy``` qui permet de tirer aléatoirement un nombre avec une distribution uniforme continue ```np.random.uniform()``` ou bien une loi normale ```np.random.normal()```. # In[ ]: import numpy as np import matplotlib.pyplot as plt N = 1000000 # Nombre de tirages aléatoires theta = 50 # valeur de l'angle en degrés precision = 5 # précision en degrés (cas d'une distribution uniforme continue) u_theta = 5 # incertitude-type en degrés (cas d'une distribution gaussienne) # Création d'un tableau de valeurs tirées aléatoirement. D'abord avec une distribution uniforme continue theta_rand = np.random.uniform(theta - precision, theta + precision, N) # on tire N valeurs aléatoires suivant une loi uniforme continue dans l'intervalle [theta - precision,theta + precision] # theta_rand est un tableau contenant N valeurs (type np.ndarray) # puis avec une distribution gaussienne #theta_rand = np.random.normal(theta, u_theta, N) # on tire N valeurs aléatoires suivant une loi normale de moyenne theta et d'écart-type u_theta #on calcule le sinus de chaque angle tiré aléatoirement sin_theta_rand = np.sin(np.radians(theta_rand)) # np.radians() convertit les degrés en radians # sin_theta_rand est un tableau de N valeurs contenant le sinus de chaque angle tiré aléatoirement N_c = 100 # calcul de l'histogramme des valeurs de theta plt.hist(theta_rand, N_c, alpha = 0.5) # plt.hist calcule l'histogramme des valeurs de theta_rand avec N_c classes réparties uniformément entre la valeur min et la valeur max # le nombre alpha, compris entre 0 et 1, détermine la transparence de l'histogramme (1 = opaque, 0 = transparent) # affichage de l'histogramme des valeurs de theta plt.xlabel('theta') # légende en l'abscisse plt.ylabel('Frequence') # légende en ordonnée plt.title('Dispersion des valeurs de theta') # titre du graphe plt.grid(True) # une grille se superpose à l'histogramme plt.show() # affichage de l'histogramme # calcul de l'histogramme des valeurs de sin_theta plt.hist(sin_theta_rand, N_c, alpha = 0.5) # affichage de l'histogramme des valeurs de sin_theta plt.xlabel('sin_theta') plt.ylabel('Frequence') plt.title('Dispersion des valeurs de sin(theta)') plt.grid(True) plt.show() # calcul de l'incertitude-type u(sin_theta) print(np.std(sin_theta_rand,ddof=1)) # Vous remarquez que même si $\theta$ est distribué de manière continue, ce n'est pas le cas de $\sin(\theta)$. Certaines valeurs sont plus probables que d'autres. # # Grâce à la fonction ```np.std()```, on détermine rapidement l'incertitude-type $u(\sin(\theta))$. La méthode de Monte-Carlo suppose de connaître à l'avance la loi de probabilité de la grandeur mesurée. Elle a l'avantage de donner un résultat aussi précis que l'on veut, à condition de choisir un nombre de tirages $N$ suffisamment élevé. La puissance de calcul des processeurs actuels autorise sans problème à aller jusqu'au million de tirages. # # Le calcul de $u(\sin(\theta))$ sera particulièrement utile pour le TP d'aujourd'hui. Afin d'éviter d'avoir à répéter toutes ces instructions pour chaque angle, nous pouvons définir une fonction Python qui fera automatiquement le calcul d'incertitude-type composée pour un ensemble d'angles $\theta$ rassemblés dans un tableau numpy. # # In[ ]: import numpy as np # ITC_sin est une fonction qui prend en entrée un tableau contenant des valeurs d'angle en degrés (theta est un np.ndarray) et un tableau de même taille contenant l'amplitude de la plage d'incertitude # associée à chacune de ces valeurs (cas d'une distribution uniforme continue). Ex : theta = ([0, 20, 40, 60, 80]) et delta = ([0.5, 0.5, 0.7, 1, 1.5]) # ITC_sin renvoie un tableau de même taille que theta contenant l'incertitude-type sur chaque valeur de sin(theta). Cette fonction généralise le code précédent en l'appliquant à un ensemble de valeurs d'angles # rentrées dans un tableau plutôt qu'à un angle unique def ITC_sin(theta, precision): N = 100000 # Nombre de tirages aléatoires u_sin = np.zeros(len(theta)) # on prédéfinit un tableau de même taille que theta, remplit initialement de zéros for k in range(len(theta)): theta_rand = np.random.uniform(theta[k] - precision[k], theta[k] + precision[k], N) # pour chaque rang k, on tire aléatoirement N valeurs dans l'intervalle [theta[k] - precision[k],theta[k] + precision[k]] sin_theta_rand = np.sin(np.radians(theta_rand)) # on calcule la distribution des valeurs de sinus u_sin[k] = np.std(sin_theta_rand) # on calcule l'incertitude-type u(sin(theta[k])), que l'on affecte au rang k du tableau u_sin return u_sin # une fois la boucle terminée, u_sin est remplie avec les incertitudes-type de chacun des sin(theta[k]) # Exemple r = np.array([0, 7.2, 13.4, 19.5, 25.2, 30.6, 35.25, 38.6, 41]) # tableau des angles en degrés precision = np.array([0.15, 0.15, 0.15, 0.15, 0.15, 0.15, 0.15, 0.25, 0.5]) # tableau des précisions en degrés ITC_sin(r, precision) # tableau des u(sin(theta[k])) # # **2. Regression linéaire avec Python** # # Nous cherchons à déterminer le modèle affine $y=ax+b$ qui approche le mieux possible les différents points ($x_i,y_i$) obtenus expérimentalement. # # Pour cela, nous pouvons utiliser la fonction ```np.polyfit(x, y, 1)``` qui détermine le polynôme de degré 1 qui s'approche le mieux possible des points expérimentaux. ```np.polyfit(x, y, 1)``` renvoie dans un tableau le coefficient directeur $a$ et l'ordonnée à l'origine $b$ en utilisant une méthode de calcul appelée **méthode des moindres carrés ordinaire** (MCO). Cette méthode consiste à calculer pour chaque point ($x_i,y_i$) l'écart par rapport au modèle, aussi appelé le **résidu** : # # $$r_i=y_i-(ax_i+b)$$ # # et à chercher les valeurs de $a$ et $b$ qui minimisent la somme des carrés des résidus : # # $$\sum_{i}r_{i}^{2}=\sum_{i}\left[y_{i}-(ax_{i}+b)\right]^{2}$$ # # C'est également la méthode utilisée par votre calculatrice pour effectuer une regression linéaire avec le mode "STAT". On illustre la méthode sur l'exemple ci-dessous. # In[ ]: # Importation des bibliothèques import numpy as np import matplotlib.pyplot as plt get_ipython().run_line_magic('config', "InlineBackend.figure_format='retina' # optionnel : pour améliorer la résolution des figures") i = np.array([]) # tableau des valeurs de i en degrés, à remplir avec vos valeurs r = np.array([]) # tableau des valeurs de r en degrés, à remplir avec vos valeurs precision_r = np.array([]) # tableau des précisions pour chaque valeur de r, à remplir avec vos valeurs sini = np.sin(np.radians(i)) # calcul du sinus de chaque angle i sinr = np.sin(np.radians(r)) # calcul du sinus de chaque angle r u_sin_r = ITC_sin(r, precision_r) # tableau des u(sinr) p = np.polyfit(sini, sinr, 1) # p est un tableau contenant les coefficients de l'ajustement # p[0] : pente ; p[1] : ordonnée à l'origine #print(p) sinr_mod = np.polyval(p, sini) # calcul des ordonnées des points "modèle" residus = sinr - sinr_mod # calcul des résidus RN = residus / u_sin_r # calcul des résidus normalisés plt.figure(figsize = (17,5)) # création d'une fenêtre graphique plt.subplot(1, 3, 1) # 1er graphe plt.errorbar(sini, sinr, yerr = u_sin_r, fmt = 'bo',label = "Points expérimentaux") # points expérimentaux plt.plot(sini, sinr_mod, 'c--',label = "Modèle affine") # points issus de l'ajustement plt.xlabel(r"$\sin(i)$",fontsize=15) plt.ylabel(r"$\sin(r)$",fontsize=15) plt.xticks(fontsize=15) plt.yticks(fontsize=15) plt.grid(), plt.legend(loc = 'best',fontsize=15) plt.subplot(1, 3, 2) # 2ème graphe plt.errorbar(sini, residus, yerr = u_sin_r, fmt = 'bo') # résidus avec barres d'incertitude-type plt.plot([0, 1], [0, 0], 'c--') # pour mieux visualiser la droite correspondant à un résidu nul plt.xlabel(r"$\sin(i)$") plt.ylabel(r"résidus"), plt.ticklabel_format(axis = 'y', style = 'sci', scilimits = (0,0)) plt.grid() plt.subplot(1, 3, 3) # 3ème graphe plt.plot(sini, RN, 'bo') # résidus normalisés plt.fill_between([0, 1], y1 = -2, y2 = 2, color = 'cyan', alpha = .1) # pour mieux visualiser le domaine des En acceptables plt.xlabel(r"$\sin(i)$",fontsize=15) plt.ylabel(r"résidus normalisés",fontsize=15), plt.ylim(-3,3) plt.xticks(fontsize=15) plt.yticks(fontsize=15) plt.grid() plt.show() # Pour déterminer si le modèle est compatible avec les points expérimentaux, on se posera les questions suivantes : # # * les écarts normalisés sont-ils tous inférieurs à $2$ ? Si ce n'est pas le cas, soit le modèle est incorrect, soit les incertitudes ont été sous-évaluées, soit il y a une/des valeur(s) aberrante(s). # * les résidus sont-ils distribués aléatoirement autour de $0$ ? Si ce n'est pas le cas, cela signifie peut-être qu'un modèle autre qu'affine pourrait mieux convenir. # * les incertitudes relatives sont-elles suffisamment faibles ? Si les barres d'erreur sont trop larges, cela suppose que de nombreux modèles autre qu'affine pourraient convenir. Dans ce cas, on ne peut pas conclure. # # # **3. Déterminer l'incertitude sur les paramètres de la regression linéaire** # # Une fois le modèle affine obtenu, on peut se demander comment les incertitudes sur les valeurs des $y_i$ affectent l'incertitude sur les valeurs de $a$ et $b$. Pour le savoir, on peut à nouveau utiliser une méthode de Monte-Carlo. # # Dans ce cas, connaissant la loi de probabilité de chacun des $y_i$, on tire aléatoirement un ensemble de points $(x_i,y_i)$ et on détermine les valeurs de $a$ et $b$ correspondantes avec ```np.polyfit()```. On répète ensuite l'opération un très grand nombre de fois avec des tirages aléatoires de points expérimentaux. On obtient alors deux tableaux avec des valeurs de $a$ et des valeurs de $b$. Il ne reste plus qu'à déterminer les incertitudes-type $u(a)$ et $u(b)$ avec ```np.std()```. # In[ ]: N = 100000 # nombre de tirages aléatoires a = np.zeros(N) # on prédéfinit deux tableaux pour les valeurs de a et b qu'on remplit initialement avec des zéros b = np.zeros(N) i = np.array([]) # tableau des valeurs de i en degrés, à remplir avec vos valeurs r = np.array([]) # tableau des valeurs de r en degrés, à remplir avec vos valeurs sini = np.sin(np.radians(i)) # calcul du sinus de chaque angle i sinr = np.sin(np.radians(r)) # calcul du sinus de chaque angle r precision_r = np.array([]) # tableau des precisions pour chaque valeur de r, à remplir avec vos valeurs for k in range(N): r_rand = r + np.random.uniform(-precision_r, precision_r) # on tire aléatoirement des valeurs de r en prenant pour chaque valeur la plage d'incertitude # delta_r qui lui correspond sinr_rand = np.sin(np.radians(r_rand)) # on passe au sinus p = np.polyfit(sini, sinr_rand, 1) # on calcule les paramètres de regression en prenant les mêmes sin(i) à chaque fois car on néglige l'incertitude sur i a[k] = p[0] # on remplit le k-ième rang du tableau a avec le coefficient directeur (premier élément du tableau renvoyé par polyfit) b[k] = p[1] # on remplit le k-ième rang du tableau b avec l'ordonnée à l'origine (deuxième élément du tableau renvoyé par polyfit) #Calcul de la moyenne a_moy = np.mean(a) #Calcul de l'écart-type u_a = np.sqrt(N / (N - 1)) * np.std(a) #Affichage du résultat print("a = " + str(a_moy) + r" ; u(a) = " + str(u_a)) # Proposez enfin une méthode pour déterminer la valeur de l'indice de réfration du prisme avec son incertitude-type associée.