Sulla base della lettura di un file CSV, calcola i quadrati delle ascisse per linearizzare la relazione con le ordinate e quindi fornisce i coefficienti della retta di regressione, le incertezze dei parametri e il coefficiente di correlazione utilizzando la funzione linregress
del pacchetto SciPy.
Visualizza i risultati in due grafici a dispersione dove, per la piccolezza delle incertezze, appaiono solo i dati: il primo comprendente i dati linearizzati, il secondo gli stessi non linearizzati.
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import linregress
nome_file = input("Inserire il nome del file CSV: ")
# motoAcc.csv
file_in = open(nome_file, "r")
coppie_dati = np.loadtxt(file_in, delimiter = ",", comments = '#', usecols = (0,1))
Viene utilizzata la funzione linregress()
della libreria SciPy. Riportati i risultati nella tuple esiti
vengono da questa estratti gli elementi slope
e intercept
.
n = len(coppie_dati)
nparrayX_Y = coppie_dati.transpose()
dati_x = nparrayX_Y[0]
# i dati in ascissa vengono linearizzati eseguendone il quadrato
dati_x_quad = dati_x**2
dati_y = nparrayX_Y[1]
file_in.close()
esiti = linregress(dati_x_quad, dati_y)
a = esiti.slope
b = esiti.intercept
Si estraggono le stime delle varie deviazioni standard e il coefficiente di correlazione.
previsioni_y = a*dati_x_quad + b
residui = dati_y - previsioni_y
sigma = np.std(residui, ddof=2)
sigma_a = esiti.stderr
sigma_b = esiti.intercept_stderr
coeff_correlazione = esiti.rvalue
Definiscono gli estremi degli intervalli entro cui tracciare la retta e la parabola
min_x_quad = min(dati_x_quad)
max_x_quad = max(dati_x_quad)
delta_x_quad = (max_x_quad-min_x_quad)/15
min_x = min(dati_x)
max_x = max(dati_x)
delta_x = (max_x-min_x)/15
Si costruiscono le rette che individuano la fascia di incertezza (non appaiono nei grafici perché sovrapposte dalla retta e dalla parabola di regressione)
x_quad = np.linspace(min_x_quad-delta_x_quad, max_x_quad+delta_x_quad, 100)
y_quad = a*x_quad + b
y_sup_quad = (a+sigma_a)*x_quad + (b+sigma_b)
y_inf_quad = (a-sigma_a)*x_quad + (b-sigma_b)
x = np.linspace(min_x-delta_x, max_x+delta_x, 100)
y = a*x**2 + b
y_sup = (a+sigma_a)*x + (b+sigma_b)
y_inf = (a-sigma_a)*x + (b-sigma_b)
Primo grafico: dispersione dei dati, retta di regressione e fascia di incertezza (il posizionamento del testo è scelto ad hoc).
plt.grid(which = 'both', color = '.85', linestyle = '-', linewidth=1)
plt.vlines(dati_x_quad, dati_y-sigma, dati_y+sigma, linewidth = .5, color = 'b')
plt.fill_between(x_quad, y_inf_quad, y_sup_quad, alpha = .1, linewidth = 0, color='r')
plt.plot(x_quad, y_quad, color = 'red', linewidth = 2, label = 'retta di regressione')
plt.scatter(dati_x_quad, dati_y, s = 10, c = 'cornflowerblue', zorder = 3, label = 'dati rilevati')
plt.title("Moto uniformemente accelerato")
plt.xlabel('quadrato degli istanti $z=t^2$ (s$^2$)')
plt.ylabel('posizione y (m)')
plt.text(2.85, 0, 'y = ({0:6.4f})z + ({1:6.4f})\na = {0:6.4f} $\pm$ {2:6.4f}\nb = {1:6.4f} $\pm$ {3:6.4f}\n r = {4:6.3f}\naccel = {5:6.4f} m/s$^2$\n$\sigma(accel)$ = {6:6.4f} m/s$^2$'.format(a, b, sigma_a, sigma_b, coeff_correlazione, 2*a, 2*sigma_a), c = 'r')
plt.legend()
plt.show()
Parabola di regressione con sovrapposti i dati iniziali.
plt.grid(which = 'both', color = '.85', linestyle = '-', linewidth=1)
plt.vlines(dati_x, dati_y-sigma, dati_y+sigma, linewidth = .5, color = 'b')
plt.fill_between(x, y_inf, y_sup, alpha = .1, linewidth = 0, color='r')
plt.plot(x, y, color = 'red', linewidth = 2, label = 'parabola di regressione')
plt.scatter(dati_x, dati_y, s = 10, c = 'cornflowerblue', zorder = 3, label = 'dati rilevati')
plt.title("Moto uniformemente accelerato")
plt.xlabel('istanti t (s)')
plt.ylabel('posizione y (m)')
plt.text(0, .4, 'y = ({0:6.4f})$\,t^2$ + ({1:6.4f})\na = {0:6.4f} $\pm$ {2:6.4f}\nb = {1:6.4f} $\pm$ {3:6.4f}\n r = {4:6.3f}\naccel = {5:6.4f} m/s$^2$\n$\sigma(accel)$ = {6:6.4f} m/s$^2$'.format(a, b, sigma_a, sigma_b, coeff_correlazione, 2*a, 2*sigma_a), c = 'r')
plt.legend()
plt.show()