Approximation et extrapolation par moindres carrés

import numpy as np
import numpy.linalg as linalg
import matplotlib.pyplot as plt

# Chargement des données
data = np.loadtxt("data/temperatures.csv", delimiter=",")
year = data[:,0]
temp = data[:,1]
N = len(year)

Approximation par moindres carrés¶

\theta = (H^TH)^{-1} H^Ty

# Normalisation des abscisses entre 0 et 1
debut = year[0]
fin   = year[-1]
x = (year-debut) / (fin-debut)

# Données sous forme de vecteur colonne
y = np.transpose(temp)

# Matrice H
H = np.transpose(np.array([ x**0, x**1, x**2 ]))

# Estimation des paramètres θ
HT     = np.transpose(H)
HTH    = np.dot(HT, H)
HTHinv = linalg.pinv(HTH)
HTy    = np.dot(HT, y)
theta  = np.dot(HTHinv, HTy)

print(theta)

# Calcul du polynôme z
z = np.dot(H,theta)

[-0.19123016 -0.780101    1.98398694]

# Affichage
plt.figure(figsize=(10,6))
plt.plot(year, temp, label="Mesures")
plt.plot(year, z, label="Approximation")
plt.xlabel("Année")
plt.ylabel("Température")
plt.legend()
plt.grid()
plt.show()

Extrapolation¶

Pour extrapoler la courbe (c’est-à-dire calculer ses valeurs au délà des données), il suffit de calculer le polynôme sur une autre plage des abscisses, donc en redéfinissant $H$ .

year2 = np.arange(debut,2070)
x = (year2-debut)/(fin-debut)
H = np.transpose(np.array([ x**0, x**1, x**2 ]))
z = np.dot(H,theta)

# Affichage
plt.figure(figsize=(10,6))
plt.plot(year, temp, label="Mesures")
plt.plot(year2, z, label="Extrapolation")
plt.xlabel("Année")
plt.ylabel("Température")
plt.legend()
plt.grid()
plt.show()

Si vous avez le moral accroché, vous pouvez connaître les conséquences issues de cette extrapolation sur l’article Wikipedia Changement climatique... 😱

Données :