Sulla base di un certo numero di dati casuali $n=30$, calcola i coefficienti della retta di regressione e visualizza i risultati in un grafico a dispersione.
Importa le funzioni dei moduli NumPy e, per la parte grafica, Matplotlib.
import numpy as np
import matplotlib.pyplot as plt
Codice del tutto simile al notebook esempio2.ipynb
n = 30
narray = []
for i in range(n):
x_rnd = np.random.uniform(0,50)
coppia_dati = np.array([x_rnd, (2+np.random.uniform(-20,20)+(-3.5)*x_rnd)])
narray.append(coppia_dati)
nparray = np.array(narray)
nparrayX_Y = nparray.transpose()
dati_x = nparrayX_Y[0]
dati_y = nparrayX_Y[1]
XY = dati_x * dati_y
quadrato_x = dati_x**2
a = (n*sum(XY) - sum(dati_x)*sum(dati_y)) / (n*sum(quadrato_x) - sum(dati_x)**2)
b = (sum(dati_y) - a*sum(dati_x))/n
Calcolo delle dimensioni degli assi: utili per definire la finestra di visualizzazione dei risultati.
min_x = min(dati_x)
max_x = max(dati_x)
min_y = min(dati_y)
max_y = max(dati_y)
delta_x = (max_x-min_x)/15
delta_y = (max_y-min_y)/15
Definisce 100 punti equispaziati in ascissa e calcola le corrispondenti ordinate sulla retta di regressione $y=ax+b$.
x = np.linspace(min_x-delta_x, max_x+delta_x, 100)
y = a*x + b
In corrispondenza dei valori "misurati" delle ascisse, calcola l'ordinata prevista dal modello
previsioni_y = a*dati_x + b
A seconda del segno della pendenza definisce la posizione del riquadro dei risultati
if a>=0:
posizione_txt = max_x - 4*delta_x
elif a<0:
posizione_txt = min_x + delta_x
La parte grafica riporta la dispersione dei dati, la linea di tendenza o regressione e i segmenti che rappresentano i residui. In riquadro i risultati.
plt.grid(which='both', color = '.85', linestyle='-', linewidth=1)
plt.vlines(dati_x, dati_y,previsioni_y, linewidth = 1, color = 'orange')
plt.plot(x,y, color = 'red', linewidth = 2, label = 'retta di regressione')
plt.scatter(dati_x, previsioni_y, s = 2, c = 'b', zorder = 2)
plt.scatter(dati_x, dati_y, s = 10, c = 'cornflowerblue', zorder = 3, label = 'dati rilevati')
plt.title("Dati e regressione lineare")
plt.xlabel('grandezza X')
plt.ylabel('grandezza Y')
plt.text(posizione_txt, min_y+delta_y/2, 'y = ('+str(round(a,2))+')x + ('+str(round(b,2))+')',c='r')
plt.legend()
plt.show()