Fit simultanei
Supponiamo di volere fittare una retta $f(x) = a x + b$ a dei dati $(x_i,y_i)$.
Tipicamente si stabilisce che le migliori stime per $a$ e $b$ sono quelle per cui si minimizza la somma dei quadrati degli scarti dalla media:
$S = \sum_{i=1}^{n} (f(x_i) - y_i)^2$
Ciò è equivalente a minimizzare la varianza:
$\sigma^2 = \frac{S}{n-1}$
Supponiamo ora di volere avere un ulteriore set di $n$ dati $(x_j^,y_j^,)$, e di volere eseguire due fit simultanei, che trovino le migliori stime per $a$ e $b$ dati i due set di dati, ma senza volerli unire in unico set.
Siano $\sigma_1^2$ e $\sigma_2^2$ le varianze dei residui per $a$ e $b$ dati.
Qual è la quantità da minimizzare questa volta?
Io rifacendomi alla teoria su come si combinano i risultati di due esperimenti diversi, pensavo alla media armonica delle varianze:
$\sigma^2 = \frac{2}{ \frac{1}{ \sigma_1^2 } + \frac{1}{ \sigma_2^2 }}$.
Ma se mischiassi i due set di dati a formarne uno solo, il risultato da minimizzare sembra essere la media aritmetica delle due varianze (per $n$ grande).
Qual è la soluzione corretta?
A me sembra che i due metodi siano entrambi legittimi, tuttavia se una delle due misure è molto più accurata dell'altra, mischiandole semplicemente si perde molto, invece la media armonica tiene maggiormente in considerazione la misura più accurata, senza buttare del tutto la misura meno accurata.
E' corretto?
Tipicamente si stabilisce che le migliori stime per $a$ e $b$ sono quelle per cui si minimizza la somma dei quadrati degli scarti dalla media:
$S = \sum_{i=1}^{n} (f(x_i) - y_i)^2$
Ciò è equivalente a minimizzare la varianza:
$\sigma^2 = \frac{S}{n-1}$
Supponiamo ora di volere avere un ulteriore set di $n$ dati $(x_j^,y_j^,)$, e di volere eseguire due fit simultanei, che trovino le migliori stime per $a$ e $b$ dati i due set di dati, ma senza volerli unire in unico set.
Siano $\sigma_1^2$ e $\sigma_2^2$ le varianze dei residui per $a$ e $b$ dati.
Qual è la quantità da minimizzare questa volta?
Io rifacendomi alla teoria su come si combinano i risultati di due esperimenti diversi, pensavo alla media armonica delle varianze:
$\sigma^2 = \frac{2}{ \frac{1}{ \sigma_1^2 } + \frac{1}{ \sigma_2^2 }}$.
Ma se mischiassi i due set di dati a formarne uno solo, il risultato da minimizzare sembra essere la media aritmetica delle due varianze (per $n$ grande).
Qual è la soluzione corretta?
A me sembra che i due metodi siano entrambi legittimi, tuttavia se una delle due misure è molto più accurata dell'altra, mischiandole semplicemente si perde molto, invece la media armonica tiene maggiormente in considerazione la misura più accurata, senza buttare del tutto la misura meno accurata.
E' corretto?
Risposte
Grazie per la risposta.
Ok.
Supponiamo che una misura ti dia $5.10 \pm 0.12$ e l'altra $5.15 \pm 0.03$.
La seconda misura è più accurata, in quanto ha un errore più piccolo, il che sicuramente vuol dire che la sua varianza era minore. Stando ai libri di Analisi degli Errori (per esempio il Taylor), se voglio stabilire la migliore stima di questa grandezza, avendo le due misure, non scarto la prima perchè meno accurata, ma:
- per il valore best faccio una media pesata tra le due misure;
- per la varianza prendo la metà della media armonica delle varianze (se le misure fossero 3 prenderei 1/3, e così via).
In questo modo ottengo una barra d'errore che è più piccola della più accurata delle misure, grazie all'informazione congiunta delle due.
Supponiamo adesso che io sia lo sperimentatore che ha svolto entrambi gli esperimenti, e anzichè analizzarli separatamente li voglia in qualche modo unire. Se mischiassi i dati mi sembra che perderei la maggiore accuratezza data dalla seconda misura.
Un'alternativa mi pare quella dei fit simultanei, minimizzando la media armonica delle due varianze (mi sono accorto ora che in realtà dovrei minimizzare la sua metà), il che mi pare concettualmente quasi equivalente a combinare le due misure a posteriori.
Qui ho ridotto il problema per potermi concentrare sull'aspetto statistico del problema. Nella realtà la cosa è complicata dal fatto che i due set di dati hanno dei parametri sicuramente comuni, mentre altri possono essere diversi, da cui il non poterli mischiare. Quindi vorrei fare un fit simultaneo, imponendo che i parametri comuni siano uguali nei due set e gli altri siano liberi.
"Sergio":
Perché dividi per \(n-1\)?
Nel modello su cui stai ragionando, \(f(x)=ax+b\), i parametri sono due, \(a\) e \(b\), quindi si divide per \(n-2\).
Ok.
"Sergio":
Ora come ora posso solo dire che la minimizzazione è ovviamente separata (due set di dati, due fit, ciascun fit scaturisce da una minimizzazione) e se una misura è più "accurata" dell'altra ("accurata" in che senso?) tale rimane.
Supponiamo che una misura ti dia $5.10 \pm 0.12$ e l'altra $5.15 \pm 0.03$.
La seconda misura è più accurata, in quanto ha un errore più piccolo, il che sicuramente vuol dire che la sua varianza era minore. Stando ai libri di Analisi degli Errori (per esempio il Taylor), se voglio stabilire la migliore stima di questa grandezza, avendo le due misure, non scarto la prima perchè meno accurata, ma:
- per il valore best faccio una media pesata tra le due misure;
- per la varianza prendo la metà della media armonica delle varianze (se le misure fossero 3 prenderei 1/3, e così via).
In questo modo ottengo una barra d'errore che è più piccola della più accurata delle misure, grazie all'informazione congiunta delle due.
Supponiamo adesso che io sia lo sperimentatore che ha svolto entrambi gli esperimenti, e anzichè analizzarli separatamente li voglia in qualche modo unire. Se mischiassi i dati mi sembra che perderei la maggiore accuratezza data dalla seconda misura.
Un'alternativa mi pare quella dei fit simultanei, minimizzando la media armonica delle due varianze (mi sono accorto ora che in realtà dovrei minimizzare la sua metà), il che mi pare concettualmente quasi equivalente a combinare le due misure a posteriori.
"Sergio":
Francamente non capisco dove vuoi andare a parare. Se hai due set di dati relativi a uno stesso fenomeno, perché non vorresti unirli? Se me ne fosse chiaro il motivo mi sarebbe più facile rispondere.
Qui ho ridotto il problema per potermi concentrare sull'aspetto statistico del problema. Nella realtà la cosa è complicata dal fatto che i due set di dati hanno dei parametri sicuramente comuni, mentre altri possono essere diversi, da cui il non poterli mischiare. Quindi vorrei fare un fit simultaneo, imponendo che i parametri comuni siano uguali nei due set e gli altri siano liberi.
So che ci sono varie assunzioni, per cominciare le volevo dare tutte per scontate.
Non è che non voglio dire esattamente qual è il mio caso, ma per mettermi a descrivere esattamente la situazione dovrei fare un trattato di 2 pagine, dopo il quale si potrebbe pure discutere per mesi sulla validità delle assunzioni fatte, sia dal punto di vista statistico che dal punto di vista fisico. Per questo sto proponendo un caso semplice, che abbia tutte le buone proprietà assunte dal Taylor, quindi parlavo di combinare due misure assolutamente identiche, tranne che per lo scatter dei punti.
A posteriori ho detto come si dovrebbero combinare. Il dubbio è come fare se invece voglio fare un fit simultaneo come descritto.
Se proprio volessimo complicare la cosa (ma non credo ce ne sia bisogno), si può pensare che i due set di misure siano relativi allo stesso fenomeno, debbano avere la stessa pendenza, ma l'intercetta è affetta da un offset strumentale, che varia di misura in misura in maniera incognita. In questo modo imporrei un fit simultaneo sui due set, con coefficiente angolare uguale per entrambi e intercette totalmente libere.
Non è che non voglio dire esattamente qual è il mio caso, ma per mettermi a descrivere esattamente la situazione dovrei fare un trattato di 2 pagine, dopo il quale si potrebbe pure discutere per mesi sulla validità delle assunzioni fatte, sia dal punto di vista statistico che dal punto di vista fisico. Per questo sto proponendo un caso semplice, che abbia tutte le buone proprietà assunte dal Taylor, quindi parlavo di combinare due misure assolutamente identiche, tranne che per lo scatter dei punti.
A posteriori ho detto come si dovrebbero combinare. Il dubbio è come fare se invece voglio fare un fit simultaneo come descritto.
Se proprio volessimo complicare la cosa (ma non credo ce ne sia bisogno), si può pensare che i due set di misure siano relativi allo stesso fenomeno, debbano avere la stessa pendenza, ma l'intercetta è affetta da un offset strumentale, che varia di misura in misura in maniera incognita. In questo modo imporrei un fit simultaneo sui due set, con coefficiente angolare uguale per entrambi e intercette totalmente libere.
Ringrazio per la pazienza e la pretesa di rigore matematico, ma penso che la mia richiesta sia molto più semplice.
Se l'offset strumentale è costante lungo ciascuno dei singoli set di misure, non influenza la pendenza delle rette, per cui penso di potere assumere che questa sia condizionata da errori casuali, con errori sistematici trascurabili.
Non ho chiesto una ricetta universalmente valida, ma ho fatto due esempi specifici:
1. doppio fit lineare con parametri comuni [nota]Ci si può chiedere perchè non mischiare i due set di dati, ma in linea di principio deve essere possibile fare questa scelta; comunque una motivazione è data dal fatto che se li mischio prima, o dopo, ottengo risultati diversi, almeno le barre d'errore[/nota];
2. doppio fit lineare con pendenza comune, e intercette indipendenti.
"Sergio":
Se non sbaglio, quando facevi riferimento al Taylor intendevi il suo capitolo 7, in cui dice «in tutto questo capitolo assumeremo che tutti gli errori sistematici sono stati identificati e ridotti ad un livello trascurabile».
Sbaglierò, ma non mi pare che l'ipotesi di "un offset strumentale che varia di misura in misura in maniera incognita" sia coerente con quell'assunto.
Se l'offset strumentale è costante lungo ciascuno dei singoli set di misure, non influenza la pendenza delle rette, per cui penso di potere assumere che questa sia condizionata da errori casuali, con errori sistematici trascurabili.
"Sergio":
Se esiste una ricettina facile e valida per qualsiasi situazione, senza "complicare le cose", non la conosco.
Non ho chiesto una ricetta universalmente valida, ma ho fatto due esempi specifici:
1. doppio fit lineare con parametri comuni [nota]Ci si può chiedere perchè non mischiare i due set di dati, ma in linea di principio deve essere possibile fare questa scelta; comunque una motivazione è data dal fatto che se li mischio prima, o dopo, ottengo risultati diversi, almeno le barre d'errore[/nota];
2. doppio fit lineare con pendenza comune, e intercette indipendenti.
Perfetto, penso che sia proprio questa la risposta che cercavo.
Dando un'occhiata a wikipedia (http://it.wikipedia.org/wiki/Analisi_della_varianza),
noto che vengono considerati due contributi alla varianza totale:
- la varianza between, che tiene conto delle differenze tra i vari set di dati;
- la varianza within, che è la varianza media.
Per quanto riguarda la varianza within, mi pare che sia la media aritmetica delle varianze (se i set di dati hanno tutti la stessa lunghezza, come nel mio caso).
Mi chiedevo se, sotto certe ipotesi, ha senso usare una media armonica anzichè aritmetica.
Dando un'occhiata a wikipedia (http://it.wikipedia.org/wiki/Analisi_della_varianza),
noto che vengono considerati due contributi alla varianza totale:
- la varianza between, che tiene conto delle differenze tra i vari set di dati;
- la varianza within, che è la varianza media.
Per quanto riguarda la varianza within, mi pare che sia la media aritmetica delle varianze (se i set di dati hanno tutti la stessa lunghezza, come nel mio caso).
Mi chiedevo se, sotto certe ipotesi, ha senso usare una media armonica anzichè aritmetica.
Vediamo cosa ho capito finora.
Sia il primo set di dati $(x_i,y_i)$, $1<=i<=n$, e il secondo set di dati $(x'_j,y'_j)$, $1<=j<=n$.
Supponiamo che sia riuscito ad individuare le seguenti rette da una regressione con coefficiente angolare comune:
$r_1$: $ y = 2x$
$r_2$: $ y = 2x + 1$
La varianza complessiva è data da:
$\sigma = \frac{1}{2n-3} ( \sum_{i=1}^{n} (y_i - 2x_i )^2 + \sum_{j=1}^{n} (y'_j - 2x'_j - 1)^2 )$
Questa può essere riscritta come la somma di una componente between ed una within.
Fino a qua ci sono?
Sia il primo set di dati $(x_i,y_i)$, $1<=i<=n$, e il secondo set di dati $(x'_j,y'_j)$, $1<=j<=n$.
Supponiamo che sia riuscito ad individuare le seguenti rette da una regressione con coefficiente angolare comune:
$r_1$: $ y = 2x$
$r_2$: $ y = 2x + 1$
La varianza complessiva è data da:
$\sigma = \frac{1}{2n-3} ( \sum_{i=1}^{n} (y_i - 2x_i )^2 + \sum_{j=1}^{n} (y'_j - 2x'_j - 1)^2 )$
Questa può essere riscritta come la somma di una componente between ed una within.
Fino a qua ci sono?
"Sergio":
Conviene partire dalle devianze (poi le varianze si ottengono dividendo per i gradi di libertà) e la devianza totale è \(\sum_{i=1}^{2n} (y_i-\bar{y})^2\): somma dei quadrati degi scarti dalla media.
La devianza totale si scompone in \(\sum_{i=1}^{2n}(\hat{y}_i-\bar{y})^2+\sum_{i=1}^{2n}(y_i-\hat{y}_i)^2\), dove gli \(\hat{y}_i\) sono i valori "teorici" o "predetti", cioè \(2x\) o \(2x+1\), secondo il valore di \(s\).
La prima è la devianza spiegata (between), la seconda è quella residua (within, l'errore).
Per ottenere le varianze si divide per i gradi di libertà, che in questo caso sono \(2n-1\) per la devianza totale, \(3-1=2\) per quella spiegata (numero dei parametri meno uno), \(2n-3\) per quella residua.
Mi sa che c'è un equivoco di fondo: quando parlo di varianza, sottointendo "dei residui rispetto al modello".
Ad esempio nel caso di un semplice fit lineare:
[tex]\sigma^2 = \frac{1}{n-2} \sum_{i=1}^{n} (y_i - m x_i - q )^2[/tex]
dove $m$ e $q$ sono i parametri determinati dalla regressione.
Uso questa quantità come misura dello scostamento dei dati rispetto al modello.
Per cui la domanda originale era come mediare queste "varianze" nel caso di fit simultanei su set multipli.
Ora mi sembra di capire che quella che io chiamo varianza, sarebbe la varianza residua.
Da quello che hai scritto, la varianza residua totale sarebbe la media aritmetica delle singole varianze residue.
Ora è giusto?
"Sergio":
Indicando le devianze residue con \(D_1,D_2\), è ovvio che \(\displaystyle \frac{D_1+D_2}{2n-3}\ne\frac{D_1}{n-3}+\frac{D_2}{n-3} \).
Ok, è chiaro, mi sono espresso in modo impreciso perchè dò per scontato [tex]n \gg 1[/tex].
Quello che si sta facendo così però è di "mettere tutto insieme". Io pensavo di potere dare maggiore individualità ai singoli set di dati, calcolando le singole varianze residue, e minimizzando una varianza residua complessiva che non sia la pura media aritmetica (approssimativamente), ma qualcosa di diverso, per esempio suggerivo la media armonica. Questo mi sembra avere particolarmente senso proprio nel caso in cui diversi set hanno varianze residue diverse, ma non sono uno statistico quindi se dici che non esiste nulla di simile ci credo.
"Sergio":
Interviene inoltre un problema: non dovresti ottenere varianze residue diverse per i due sottoinsiemi di dati.
Se le ottieni, un assunto dell'analisi della covarianza - ma più in generale della stima col metodo dei minimi quadrati - viene violato e si deve ricorrere ad altri metodi.
Il motivo è che se diversi sottoinsiemi di dati presentano varianze residue diverse l'incertezza nella stima dei coefficienti è maggiore di quella che puoi calcolare, ottieni cioè stime dei coefficienti meno precise di quello che sembra, e quindi si richiedono azioni correttive.
Intuitivamente posso capire che se voglio determinare una retta è i punti a sinistra hanno molta meno incertezza rispetto ai punti a destra, non sarebbe equo contarli tutti alla stessa maniera, e per questo si può ricorrere ai minimi quadrati pesati.
Tuttavia, se invece ho due set che ricoprono lo stesso dominio lungo $x$ e hanno lo stesso numero di punti, il fatto che uno abbia una varianza residua minore dell'altro non mi disturberebbe affatto, anzi mi inviterebbe a combinare i due set in modo non banale, "fidandomi di più" del set con varianza residua minore (si potrebbe usare il metodo della media armonica che ho già proposto più volte, oppure usare una media pesata delle varianze residue, ...; certo, bisognerebbe però capire il significato preciso dell'operazione che si sceglie di fare)
"Sergio":
Se hai già dei risultati, hai varianze residue stimate su cui effettuare calcoli; se però devi "mettere insieme" i dati elementari e decidere come e quali risultati ottenere, secondo l'obiettivo dell'analisi, non hai le stime delle varianze residue, quindi .... la media armonica di cosa potresti calcolare? Di qualcosa che non hai?
Il problema è tutt'altro: come procedere per ottenere quelle stime, come verificare che siano attendibili ecc.
Questo potrebbe essere fatto iterativamente, cioè parto con delle stime dei modelli, calcolo le varianze residue individuali, le combino come voglio, e poi correggo il modello cercando di minimizzare la "varianza residua combinata".
Entrando più in dettaglio nel mio caso, le curve non sono rette, bensì funzioni piuttosto sofisticate che dipendono da vari parametri [nota]
[/nota], e il metodo che uso è quello delle Catene MonteCarlo. Nel caso in cui faccio un fit su un solo set, la varianza residua la utilizzo come varianza della likelihood function (gaussiana). Le ampiezze delle distribuzioni a posteriori dei parametri ovviamente dipendono dalla varianza della likelihood function.
Nel caso di fit a più curve, posso imporre che i parametri siano gli stessi per tutte le curve, o alcuni siano uguali e altri varino (dipende dalle assunzioni fisiche che faccio, che poi possono risultare più o meno supportate dai risultati dell'analisi). Pensavo che però scegliere come varianza della likelihood function la varianza residua "semplice" potesse non essere la scelta ottimale, specialmente se alcuni set di dati fossero migliori di altri. Comunque, come dicevo, la discussione si fa molto sofisticata a questo punto.
"Sergio":
Le varianze residue svolgono il ruolo di iperparametri e, se vengono stimate separatamente, si rientra nell'approccio detto "bayesiano empirico". Che vanta una storia gloriosa, ma ormai si preferisce un approccio bayesiano "pieno".
Se per bayesiano "pieno" si intende che la stessa varianza della likelihood function viene stimata iterativamente attraverso una catena (al momento non ricordo i dettagli), allora applico anche questo metodo. Tuttavia immagino che faccia qualche differenza la struttura con cui decido di valutare la sigma.
"Sergio":
Uno dei (tanti) vantaggi dell'inferenza bayesiana è la libertà assoluta nella modellazione del fenomeno.
D'altra parte, modellare una struttura complessa della varianza non è uno scherzo.
Potrei consigliarti qualche recente, e potente, strumento per l'inferenza bayesiana. Se ti interessa.
Certo che mi interessa, grazie. (Quasi tutte le analisi di dati che faccio ricorrono all'inferenza bayesiana, ma direi di averne una conoscenza solo essenziale)
Grazie, proverò presto ad usarlo.