Regressione circolare: diffrazione (mod. geometrico)

Genera una rappresentazione grafica dei punti sperimentali ottenuti leggendo un file di dati CSV e, tramite la regressione circolare con modello geometrico, fornisce la circonferenza ottimale con relativo centro e raggio.

import numpy as np
from scipy import optimize
import matplotlib.pyplot as plt

La prima definizione rappresenta la funzione obiettivo nel modello geometrico, quindi la funzione che produce le coordinate di punti sulla circonferenza ottimale e infine la funzione che restituisce i punti sulla circonferenza intersezioni di questa con i segmenti radiali che collegano il centro con il punto sperimentale.

def funzione_obiettivo(argomenti, ascisse, ordinate):
    # argomenti[0] = ascissa centro,  argomenti[1] = ordinata centro, argomenti[2] = raggio,
    return np.sum((np.sqrt((ascisse-argomenti[0])**2 + (ordinate-argomenti[1])**2)-argomenti[2])**2)

def punti_cfr_ottimale(xc, yc, raggio):
    alfa = np.linspace(-np.pi, np.pi, 360)
    xp = xc + raggio*np.cos(alfa)
    yp = yc + raggio*np.sin(alfa)
    return xp, yp

def punti_cfr_corrispondenti(dati_x, dati_y, xc, yc, raggio):
    angoli = np.arctan2(dati_y-yc,dati_x-xc)
    punti = np.transpose([raggio*np.cos(angoli)+xc, raggio*np.sin(angoli)+yc])
    return punti

Letto il file di dati, si assegnano alle variabili dati_x e dati_y rispettivamente le ascisse e le ordinate.

nome_file = input("Inserire il nome del file CSV: ")
file_in = open(nome_file, "r")
coppie_dati = np.loadtxt(file_in, delimiter = ",", comments = '#', usecols = (0,1))

n = len(coppie_dati)

nparrayX_Y = np.transpose(coppie_dati)
dati_x = nparrayX_Y[0]
dati_y = nparrayX_Y[1]

file_in.close()

Si calcolano, in termini approssimativi, le coordinate del centro e il raggio per poi considerare queste stime come i valori iniziali x0 da cui partire per la ricerca del minimo della funzione obiettivo.

x_med = np.mean(dati_x)
y_med = np.mean(dati_y)
stima_iniziale_raggio = np.sqrt((dati_x[0]-x_med)**2 + (dati_y[0]-y_med)**2)

x0 = np.array([x_med, y_med, stima_iniziale_raggio])

esiti = optimize.minimize(funzione_obiettivo, x0, args = (dati_x, dati_y))

I risultati del calcolo di minimizzazione (restituiti per i tre parametri dall'array x) vengono associati alle coordinate del centro e al raggio.

xc, yc, r = esiti.x
# per disporre delle rimanenti informazioni togliere il commento all'istruzione seguente
# print(esiti)

Calcolo di 360 punti della circonferenza ottimale seguito dal calcolo dei punti di questa circonferenza allineati radialmente con le misure. Segue il calcolo dei residui e della loro deviazione standard.

xp, yp = punti_cfr_ottimale(xc, yc, r)

punti_cfr = punti_cfr_corrispondenti(dati_x, dati_y, xc, yc, r)
[puntiCfr_x, puntiCfr_y] = [punti_cfr[:,0], punti_cfr[:,1]]

residui = np.sqrt((dati_x-xc)**2+(dati_y-yc)**2) - r
somma_quadrati_residui = np.sum(residui**2)
# formula con n-3 gradi di libertà
sigma = np.sqrt(somma_quadrati_residui/(n-3))

massimo_residui = np.max(abs(residui))
x = np.linspace(0, n+1, 100)

Parte grafica: circonferenza di best-fit e dispersione delle misure.

figura = plt.figure(facecolor = 'white')
# dimensioni immagine
# plt.rcParams['figure.figsize'] = [12,7.5]
plt.axis('equal')
plt.plot([dati_x, puntiCfr_x], [dati_y, puntiCfr_y], c = 'orange')
plt.plot([np.ones(n)*xc, puntiCfr_x], [np.ones(n)*yc, puntiCfr_y], linestyle = '--', linewidth = .5, c = 'gray')
plt.plot(xp, yp, linewidth = 2, c ='r', label = 'centro = ({0:5.3f}, {1:5.3f}),\nraggio = {2:3.3f}'.format(xc, yc, r))
plt.scatter(dati_x, dati_y, c = 'b',  label = 'punti sperimentali', marker = 'o', zorder = 2)
plt.scatter(xc, yc, c = 'blue', marker = 'x')
plt.legend(loc = 'best', labelspacing = 0.5)
plt.title('Regr. circolare modello geometrico: punti sperimentali e\nbest-fit: $(x-\hat a)^2+(y-\hat b)^2=\hat r^2$')
plt.show()

Visualizzazione dei residui con un indice intero in ascissa.

ax = plt.axes()
ax.xaxis.set_major_locator(plt.MultipleLocator(1))
plt.grid(which = 'both', color = '.85', linestyle = '-', linewidth=1)
plt.xlim([0, n+1])
plt.ylim([-13/10*massimo_residui, 13/10*massimo_residui])
plt.fill_between(x, -sigma, sigma, alpha =.1, linewidth = 0, color = 'r')
plt.vlines(np.arange(1, n+1), 0, residui, linewidth = 1, color = 'orange')
plt.plot([0, n+1], [0, 0], color = 'red', linewidth = 2)
plt.scatter(np.arange(1, n+1), residui, s = 20, c = 'b', zorder = 3, label = 'residui: $e_i=\sqrt{(x_i-\hat a)^2+(y_i-\hat b)^2}-\hat r$')
plt.title("Distribuzione dei residui (modello geometrico)\n devSt $\sigma = $" + str(round(sigma,3)))
plt.xlabel('indice')
plt.ylabel('residui')
plt.legend()
plt.show()