#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Tue May 18 18:00:07 2021

@author: tex
"""

import pylab as py # import fonctions graphiques
import numpy as np
import matplotlib.pyplot as plt # pour le tracé de l'histogramme

#Fixer les valeurs de CB0, VB0, VA0, sCB, sVB et sVA
CB0=0.100
VB0=5.0
VA0=20.0
sCB=0.003
sVB=0.2
sVA=0.1

# Générer 10000 valeurs de CB par tirage aléatoire dans une gaussienne 
#de moyenne CB0 et d’écart-type sCB. Idem pour VB et VA
N=1000
CB=CB0+sCB*np.random.randn(N)
VB=VB0+sVB*np.random.randn(N)
VA=VA0+sVA*np.random.randn(N)

CA=CB*VB/VA

#Calculer la moyenne et l’écart-type des valeurs de CA
def moy(l):
    N=len(l)
    s=0
    for i in range(N):
        s=s+l[i]
    return s/N

def ecart(l):
    m=moy(l)
    s=0
    N=len(l)
    for i in range(N):
        s=s+(l[i]-m)**2
    return np.sqrt(s/(N-1)) 

print(CB0*VB0/VA0)
print(moy(CA))
print(ecart(CA))


#Retrouve-t-on l’estim la règle de calcul d’incertitude de type B ?
CA0=VB0*CB0/VA0
sCA=CA0*np.sqrt((sVB/VB0)**2+(sVA/VA0)**2+(sCB/CB0)**2)
print(CA0)
print(sCA)

#Tracer l’histogramme de distribution des CA
plt.figure(4,figsize= [10, 2])

plt.subplot(1,4,1)
plt.hist(CB, bins=20)
plt.xlabel("$C_B$")

plt.subplot(1,4,2)
plt.hist(VB, bins=20)
plt.xlabel("$V_B$")

plt.subplot(1,4,3)
plt.hist(VA, bins=20)
plt.xlabel("$V_A$")

plt.subplot(1,4,4)
plt.hist(CA, bins=20)
plt.xlabel("$C_A$")


#Créer un vecteur x d’abscisses [0.1, 0.2, 0.3, … 0.9, 1]
N=10
vecx=np.linspace(0.1,1,N)

#Fixer les paramètres théoriques: atheo = 0 et btheo = 1
# Générer le vecteur ytheo
atheo=0
btheo=1
ytheo=atheo+btheo*vecx
sexp=0.05

# Générer le vecteur y des « mesures » bruitées
y=atheo+btheo*vecx+sexp*np.random.randn(N)
                                    
#Tracer la droite théorique (noir) et les « mesures »
plt.plot(vecx,ytheo)
plt.plot(vecx,y,'*')
plt.show()



#Calculer les coefficients 𝑎 et 𝑏 de la régression 
def sum1(l):
    N=len(l)
    s=0
    for i in range(N):
        s=s+l[i]
    return s

def sum2(l1,l2):
    N=len(l1)
    s=0
    for i in range(N):
        s=s+l1[i]*l2[i]
    return s

# Tracer la droite de régression (vert) sur la figure
    # Calculer la prédiction ypred par la droite pour un xpred=0,57
# Tracer en bleu le point (xpred; ypred) sur la figure
b=(N*sum2(vecx,y)-sum1(vecx)*sum1(y))/(N*sum1(vecx**2)-(sum1(vecx))**2)
a=moy(y)-b*moy(vecx)
xpred=0.57
ypred=a+b*xpred
plt.plot(vecx,a+b*vecx,color='g')
plt.plot(vecx,ytheo)
plt.plot(vecx,y,'*')
plt.plot([xpred],[ypred],'o')
plt.show()

#Générer 10000 valeurs de a,b et ypred grâce à une boucle
Np=10000
s=np.sqrt((sum1(y**2)-a*sum1(y)-b*sum2(vecx,y))/(N-2))
t=2
t1=0
for i in range(N):
    t1=t1+(vecx[i]-moy(vecx))**2
sb=t*s/np.sqrt(t1)
sa=t*s*np.sqrt((1/float(N))+(moy(vecx)**2/t1))
amc=a+sa*np.random.randn(Np)
bmc=b+sb*np.random.randn(Np)
vecypred=np.zeros(Np)
for i in range(Np):
    vecypred[i]=amc[i]+bmc[i]*xpred

plt.figure(3,figsize= [10, 2])

plt.subplot(1,3,1)
plt.hist(amc, bins=20)
plt.xlabel("a")

plt.subplot(1,3,2)
plt.hist(bmc, bins=20)
plt.xlabel("b")

plt.subplot(1,3,3)
plt.hist(vecypred, bins=20)
plt.xlabel("ypred")