#J.COURTIN   01/2023      Biblio 128 Probabilités et statistiques






##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):



    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]   #A modifier
#
# 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("")