## MÉTHODE MONTE-CARLO - Relation Y = f(X)
## Armel MARTIN - Septembre 2021
## distribué sous licence CC-by-nc-sa
import numpy as np
import matplotlib.pyplot as plt
import random

#### MESURE D'UNE FRÉQUENCE VIA UNE PÉRIODE
## Période mesurée
xm = 24.4e-6 # s
# Précision de la mesure de période
DeltaX = 0.1e-6*100# s
# Incertitude (distribution uniforme supposée)
uX = DeltaX / np.sqrt(3)

def f(x):
    return 1/x # Fonction inverse

## Génération aléatoire d'une série de mesures
N = 100000 # nombre de simulations à exécuter.

# En utilisant des listes, transformées ensuite en vecteurs Numpy
##liste_X=[]
##liste_Y=[]
##for i in range(0,N):
##    x = random.uniform(xm-DeltaX,xm+DeltaX) # distribution uniforme
##    liste_X.append(x)
##    liste_Y.append(f(x))
##vecX = np.array(liste_X)
##vecY = np.array(liste_Y)

## En utilisant directement des vecteurs Numpy
vecX = np.zeros(N)
for i in range(0,N):
    #vecX[i] = random.uniform(xm-DeltaX,xm+DeltaX) # spécifier l'intervalle
    vecX[i] = random.gauss(xm,uX) # spécifier moyenne et écart-type
vecY = f(vecX) # calcul vectorisé

## Bibliothèque numpy.random... ça ne donne pas la même chose (surtout normal, qui ne fonctionne pas, valeurs d'écart-type divergentes). 
##vecX = np.random.uniform(xm-DeltaX,xm+DeltaX,N) # spécifier l'intervalle, le nb N d'échantillons
##vecX = np.random.normal(xm,DeltaX,N) # spécifier l'intervalle, le nb N d'échantillons
##vecY = f(vecX)

## Statistique sur Y
Ym = np.mean(vecY)
uY = np.std(vecY,ddof=1) # ddof=1 => non biaisé
print('Ym =',Ym,' ; u(Y) =',uY)

plt.figure(1,figsize=[5,3])
plt.hist(vecX,bins='rice') #,color='grey')
plt.xlim(0,1.1*max(vecX))
plt.xlabel('X')# (10^-6 s)')
plt.ylabel('frequence')
plt.grid()
titre = "Distribution de X - Incert. relat. ="+ \
        format(100*uX/xm,'.1f')+'%'
plt.title(titre)
plt.tight_layout()
plt.savefig('hist_X_5.png')
plt.show()
#plt.clf()

plt.figure(2,figsize=[5,3])
#plt.hist(vecY,bins='rice')
plt.hist(vecY,range=(0,100000),bins=100)
plt.xlim(0,100000)#1.1*max(vecY))
plt.xlabel('Y')
plt.ylabel('frequence')
plt.grid()
titre = "Distribution de Y=f(X)=1/X"
plt.title(titre)
resultat_MC = 'MC: Moyenne ='+format(Ym,'.0f')+'; Incert. relat. ='+format(100*uY/Ym,'.1f')+'%'
resultat_diff = 'Diff: Moyenne ='+format(f(xm),'.0f')+'; Incert. relat. ='+format(100*uX/xm,'.1f')+'%'
plt.annotate(resultat_MC,(0,3200))#(40850,700))
plt.annotate(resultat_diff,(0,2200))#(40850,500))
plt.tight_layout()
plt.savefig('hist_f(X)_5.png')
plt.show()
#plt.clf()