#J.COURTIN   01/2023      Biblio 128 Probabilités et statistiques
#
#
#  "Les statistiques, c'est comme le bikini : ça donne des idées mais ça cache l'essentiel !"
#
#   Michel Colucci : Philosophe






##fonctions nécessaires à la régression
def mean(x):        #moyenne d'une variable aléatoire
    S=0.
    for i in x:
        S+=i
    return S/len(x)
    
def Sigma(x):       #ecart-type d'une variable aléatoire
    S=0.
    mu=mean(x)
    for i in x:
        S+=(i-mu)**2
    return (S/len(x))**0.5
    
def Cov(x,y):  #Covariance des variables aléatoires X et Y
    C=0.
    mux=mean(x)
    muy=mean(y)
    for i in range(len(x)):
        C+=(x[i]-mux)*(y[i]-muy)
    return C/len(x)
    
def ProdXY(x,y):  #Construction de la liste produit Xi.Yi
    XiYi=[]
    for i in range(len(x)):
        XiYi+=[x[i]*y[i]]
    return XiYi

def Bracket(x):
    S=0.             #    Opérateur de moyenne statistique
    Sum=0.           #!!! Pondération par sigma_Yi uniquement !!!
    for i in range(len(x)):
        S+=x[i]
        Sum+=1
    return S/Sum


## Calcul de la régression

def Reg_simple(x,y):
    #Calcul de a
    a= (Bracket(ProdXY(x,y)) - Bracket(x)*Bracket(y))/ \
       (Bracket(ProdXY(x,x)) - Bracket(x)**2)
    
    #calcul de b
    b=Bracket(y)-a*Bracket(x)

    #calcul de R
    RCoef=Cov(x,y)/(Sigma(x)*Sigma(y))

    
    return (a, b, RCoef)


## Mise en place de Monte Carlo
#
#  ATTENTION : "La routourne va tourner"    F.Ribéry   
#
from random import gauss


def monteCarlo(x,y,Dx,Dy,N=100):
    
    A = [];  B = [];  R = []
    
    for i in range(N):
        tab_X = [gauss(x[k],Dx[k]) for k in range(len(x))] 
        #les x aléatoires simulés
        tab_Y = [gauss(y[k],Dy[k]) for k in range(len(y))] 
        #les x aléatoires simulés
        
        (a, b, RCoef) = Reg_simple(tab_X, tab_Y) #régression
        A+=[a];   B+=[b];   R+=[RCoef]
        
    return (mean(A),mean(B),Sigma(A),Sigma(B), mean(R),Sigma(R))















###### Main  ######


##données avec incertitudes
x  = [1,2,3,4,5,6,7,8,9,10]
Dx = [0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1]

y  = [0.8,2.4,2.0,2.2,2.7,3.5,3.0,4.0,5.0,4.8]
Dy = [0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1]

N=1000


(a,b,Da,Db,R,DR) = monteCarlo(x,y,Dx,Dy,N)

print("")
print("*** Y = aX + b ***\n")       #Affichage

print("a : {:2.5f} ".format(a) + " +/- {:2.5f} ".format(Da) )
print("b : {:2.5f} ".format(b) + " +/- {:2.5f} ".format(Db) )
print("")
print("R2 : {:2.5f} ".format(R**2) + " +/- {:2.5f} ".format(2*R*DR) )
print("R : {:2.5f} ".format(R) + " +/- {:2.5f} ".format(DR) )
print("")