report-temperatura/algoritmo.md

Algoritmo di fitting: Minimi Quadrati Non Lineari con Trust Region Reflective

1. Introduzione

Il raffreddamento newtoniano produce curve della forma:

T(t) = T_\infty + A \cdot e^{-(t - t_0)/\tau}

I parametri da stimare sono A, \tau e t_0. Nonostante la parte esponenziale sia linearizzabile con un logaritmo, la presenza di T_\infty e di t_0 rende il modello intrinsecamente non lineare nei parametri: non esiste una trasformazione algebrica che lo riduca a una regressione lineare. È quindi necessario un algoritmo iterativo di ottimizzazione non lineare.

---

2. Minimi Quadrati Non Lineari

L'obiettivo è trovare il vettore di parametri \boldsymbol{\theta} che minimizza la somma dei quadrati dei residui pesati:

\min_{\boldsymbol{\theta}} \; S(\boldsymbol{\theta}) = \sum_{i=1}^{N} \left(\frac{y_i - f(x_i,\, \boldsymbol{\theta})}{\sigma_i}\right)^2

dove:

• y_i sono le temperature misurate,
• f(x_i, \boldsymbol{\theta}) è il modello (esponenziale semplice o doppio),
• \sigma_i è l'incertezza associata al campione i (usata per pesare i punti).

A differenza dei minimi quadrati lineari, qui S(\boldsymbol{\theta}) non ha soluzione analitica chiusa. Il problema viene risolto iterativamente sfruttando il Jacobiano J_{ij} = \partial f(x_i, \boldsymbol{\theta}) / \partial \theta_j, che descrive come ogni parametro influenza i residui localmente.

---

3. Trust Region Reflective (TRF)

scipy.optimize.curve_fit usa per default il metodo Trust Region Reflective, progettato per problemi di minimi quadrati non lineari — in particolare con vincoli di bound (parametri limitati inferiormente/superiormente).

Idea di base: regione di fiducia

Ad ogni iterazione k, invece di cercare il minimo globale di S, l'algoritmo risolve un sottoproblema locale entro una regione di fiducia di raggio \Delta_k:

\min_{\mathbf{p}} \; m_k(\mathbf{p}) = S(\boldsymbol{\theta}_k) + \mathbf{g}_k^T \mathbf{p} + \tfrac{1}{2}\mathbf{p}^T B_k \mathbf{p} \text{soggetto a } \|\mathbf{p}\| \leq \Delta_k

dove \mathbf{g}_k = J^T \mathbf{r} è il gradiente e B_k \approx J^T J è l'approssimazione hessiana (formula di Gauss–Newton).

Aggiornamento adattivo del raggio

Dopo ogni passo proposto \mathbf{p}_k, si calcola il rapporto di riduzione:

\rho_k = \frac{S(\boldsymbol{\theta}_k) - S(\boldsymbol{\theta}_k + \mathbf{p}_k)}{m_k(\mathbf{0}) - m_k(\mathbf{p}_k)}

• \rho_k \approx 1: il modello quadratico è accurato → il passo viene accettato e \Delta_k aumenta.
• \rho_k \ll 1 o negativo: il modello è impreciso → il passo viene rifiutato e \Delta_k si riduce.

Questo meccanismo garantisce convergenza anche quando il punto iniziale è lontano dal minimo.

Gestione dei bound ("Reflective")

Il termine Reflective indica la strategia per rispettare i vincoli di box (l_j \leq \theta_j \leq u_j): quando il passo proposto supera un bound, viene riflesso lungo quella direzione, mantenendo il sottoproblema risolvibile in modo efficiente senza proiezioni esplicite.

---

4. Pesi (Weighted Least Squares)

Il parametro sigma di curve_fit permette di assegnare incertezze diverse ai punti. Campioni con \sigma_i grande contribuiscono poco a S(\boldsymbol{\theta}) e vengono quasi ignorati dal fit.

Questa tecnica è usata per escludere zone di dati corrotti senza eliminarle fisicamente: ai campioni nella zona problematica si assegna \sigma_i = 10^{10}, riducendo il loro peso a zero in modo pratico. I campioni affidabili mantengono \sigma_i = 1.

---

5. Stima iniziale dei parametri

I minimi quadrati non lineari sono sensibili al punto di partenza: un'inizializzazione lontana dal minimo globale può far convergere l'algoritmo verso un minimo locale. Per questo motivo i parametri iniziali vanno scelti con cura, usando:

• Conoscenza fisica del processo (es. costante di tempo dell'ordine di decine di secondi per un oggetto metallico in aria),
• Fit preliminari su sotto-finestre dei dati (es. il risultato del fit sull'intero intervallo come punto di partenza per il doppio esponenziale),
• Ispezione visiva della curva.

---

6. Metriche di qualità

Coefficiente di determinazione R²

R^2 = 1 - \frac{\sum_i (y_i - \hat{y}_i)^2}{\sum_i (y_i - \bar{y})^2}

Misura la frazione di varianza dei dati spiegata dal modello. Valori vicini a 1 indicano un fit eccellente.

Incertezze sui parametri

curve_fit restituisce la matrice di covarianza C = (J^T J)^{-1} \sigma^2_{\text{res}}. La deviazione standard del parametro \theta_j è:

\sigma_{\theta_j} = \sqrt{C_{jj}}

Queste incertezze quantificano la sensibilità del fit alla variabilità dei dati e dipendono sia dalla qualità del modello sia dalla struttura dei residui.