stones.py
#
# Un petit peu de trucs rigolos sur les erreurs de mesures
#
# Vincent Legat - 2018
# Ecole Polytechnique de Louvain
#
from numpy import *
from numpy.random import randn
import matplotlib
from matplotlib import pyplot as plt
matplotlib.rcParams['toolbar'] = 'None'
# ============================= clone of histc MATLAB function ============
def histc(r,bin):
map = digitize(r,bin)
n = zeros(len(bin))
for i in map:
n[i-1] += 1
return n
# ============================= mainProgram ===============================
#
# -1- Observer qu'en répétant une expérience indépendante et
# et en effectuant la moyenne, on obtient une estimation plus précise :-)
#
# Créer un ensemble de 800 mesures de moyenne 500
# et dont l'écart type est 15
#
r = 500 + 15*randn(800)
# plt.figure()
# plt.hist(r,100,facecolor='r')
plt.figure()
n = histc(r,arange(450,550))
plt.bar(arange(len(n)),n,color='g')
plt.bar(arange(0,34),n[0:34],color='r')
plt.bar(arange(66,100),n[66:100],color='r')
plt.axis('off')
print('-1- Effectuons 800 mesures virtuelles avec un écart de 15 :-)')
print(' Moyenne : %14.7f ' % mean(r))
print(' Ecart type de la moyenne : %14.7f ' % (15/sqrt(800)) )
#
# -2- Observer qu'en prenant un ensemble de 800 moyennes de 10 mesures,
# on obtient encore un résultat encore plus précis...
#
# C'est logique, on fait 8000 mesures !
#
r = zeros(800);
for i in range(10):
r = r + 500 + 15*randn(800)
r = r/10
plt.figure()
n = histc(r,arange(450,550))
plt.bar(arange(len(n)),n,color='g')
plt.bar(arange(0,44),n[0:44],color='r')
plt.bar(arange(56,100),n[56:100],color='r')
plt.axis('off')
print('-2- Effectuons 800 moyennes de 10 mesures virtuelles avec un écart de 15 :-)')
print(' Moyenne des moyennes : %14.7f ' % mean(r))
print(' Ecart type de la moyenne : %14.7f ' % (15/sqrt(10)/sqrt(800)) )
#
# -3- En déduire une technique pour obtenir la meilleure estimation
# de la densité sur base d'une série de données expérimentales
#
# Comment calculer la meilleure estimation de n mesures
# dont la précision est differente....
# s achant que variance(aX) = a^2 * variance(X) ?
#
#
Vol = array([ 10, 8, 8, 10, 99])
ErrVol = array([ 5, 2, 5, 1, 2])
Mass = array([ 29.1,25.7,22.2,27.5,266.0])
ErrMass = array([ 0.1, 0.1, 0.1, 0.1, 0.1])
Rho = Mass / Vol
ErrRho = (ErrMass * Vol + ErrVol * Mass)/ (Vol*Vol)
VarRho = ErrRho*ErrRho
plt.figure();
plt.errorbar(arange(len(Rho)),Rho,ErrRho,
linestyle='',marker='o',
markerfacecolor='blue',markeredgecolor='blue',
color='red',capsize=5)
plt.axis('off')
coeff = 1 / (VarRho)
coeff = coeff / sum(coeff)
BestEstimate = coeff @ Rho
VarBestEstimate = (coeff*coeff) @ VarRho
EcartBestEstimate = sqrt(VarBestEstimate)
print('-3- Effectuons la meilleure combili de mesures de précisions différentes :-)');
for i in range(len(Rho)):
print(' Mesure %d et écart type : %14.7f (%9.7f)' % (i,Rho[i],sqrt(VarRho[i])))
print(' Meilleure estimation et écart type : %14.7f (%9.7f)' % (BestEstimate,EcartBestEstimate))
plt.show()
# =========================================================================