#!/usr/bin/python3
# -*- coding:Utf-8 -*

"""
regression.py - Version 1.4 22/05/2022

Regression linéaire avec la méthode des moindres carrés et estimation de
l'incertitude type sur a et b à l'aide de la méthode de Monte-Carlo

Basé sur le fichier regression.csv (codage utf-8 avec un séparateur ; et ,
comme séparateur décimal) contenant les points (X,Y) et sY, l'incertitude type
sur chaque valeur de Y

Licence MIT, Damien Riou
"""

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

plt.close('all')

# Nombre d'échantillon
N = 20000

# Import et traitement des données du fichier regression.csv
X, Y, sY = pd.read_csv('regression.csv', delimiter=';', decimal=',',
                       dtype='float', encoding='utf-8').transpose().to_numpy()
n = len(X)
print('Données avec {} points :\nX : {}\nY : {}'.format(n, X, Y))

# Calcul des sommes d'intérêt pour la méthode des moindes carrés
print('\nRégression linéaire par la méthode des moindres carrés...')
s_x = np.sum(X)
s_y = np.sum(Y)
s_xy = np.sum(X*Y)
s_x2 = np.sum(X**2)

# Régression linéaire
a = (n*s_xy-s_x*s_y)/(n*s_x2-s_x**2)
b = (s_y-a*s_x)/n
Ymod = a*X+b
print('Modèle : y = a x + b')
print('Coefficient a : ' + np.format_float_scientific(a, precision=2))
print('Coefficient b : ' + np.format_float_scientific(b, precision=2))

# Représentation graphique des points et du modèle
fig = plt.figure(1)
fig.suptitle('Modélisation affine des données $y=a\cdot x+b$')
(ax1, ax2) = fig.subplots(2, 1, sharex=True)
ax1.plot(X, Ymod, label='Modèle')
ax1.errorbar(X, Y, sY, fmt='.', label='Points exp.')
ax1.set_ylabel('$y$')
ax1.legend()
ax2.plot(X, Y-Ymod, 'og', label='Résidus')
ax2.set_xlabel('$x$')
ax2.set_ylabel('Résidu')

# N régressions linéaires pour avoir les incertitudes types sur a et b
print('\nEstimation des incertitudes types sur les coefficients avec {} \
droites...'.format(N))
an = np.zeros(N)  # Coefficient a, N fois
bn = np.zeros(N)  # Coefficient b, N fois
for j in range(N):
    ybrui = Y + sY*np.random.randn(n)  # Données bruitées avec Monte-Carlo
    s_y = np.sum(ybrui)  # Actualisation des données pour les moindres carrés
    s_xy = np.sum(X*ybrui)
    an[j] = (n*s_xy-s_x*s_y)/(n*s_x2-s_x**2)  # Moindres carrés, N fois
    bn[j] = (s_y-a*s_x)/n
anmoy = np.mean(an)  # Moyennes et écart-types
anstd = np.std(an)
bnmoy = np.mean(bn)
bnstd = np.std(bn)
print('Modèle : y = a x + b')
print('Coefficient a : ' + np.format_float_scientific(anmoy, precision=2))
print('Incertitude sur a : ' + np.format_float_scientific(anstd, precision=1))
print('Coefficient b : ' + np.format_float_scientific(bnmoy, precision=2))
print('Incertitude sur b : ' + np.format_float_scientific(bnstd, precision=1))

# Histogrammes de an et bn
fig = plt.figure(2)
fig.suptitle('Distributions des coefficients $a$ et de $b$')
(ax1, ax2) = fig.subplots(1, 2, sharey=True)
ax1.hist(an, bins=40)
ax1.set_xlabel('$a$')
ax1.set_ylabel('Occurence des valeurs de $a$ et de $b$')
ax2.hist(bn, bins=40)
ax2.set_xlabel('$b$')
plt.show()