Fit Lineare attraverso il metodo Monte Carlo
Salve a tutti!
Come da titolo sto cercando un metodo rigoroso che consenta di fare un fit lineare (ma anche no...) utilizzando il metodo Monte Carlo. Mi sono imbattuto in questa situazione dovendo analizzare dei dati presi durante la riproduzione dell'esperienza di Thomson. Nonostante ciò le mie domande sono più prettamente statistiche ed è per questo che scrivo in questa sezione. Detto ciò vi pongo il mio problema ringraziandovi anticipatamente!
Ho delle misure $ (x_i , y_i) $ le quali dovrebbero avere un andamento lineare. Per ogni $ y_i $ ho un'incertezza, sempre la stessa, di tipo massimo: chiamiamola $ delta_y $ . In primo luogo il fatto che le incertezze siano sempre le stesse mi porta a poter assumere che anche le ipotetiche incertezze statistiche (ossia le deviazioni standard) delle singole $ y_i $ siano le stesse per ogni $ y_i $; inoltre il fatto che siano presenti su una sola grandezza (la $ y $) mi pone nelle condizioni di poter effettuare un best-fit ai minimi quadrati. Ricordo che i parametri $ A $ e $ B $ della retta $ y=A+Bx $ sono indipendenti dalle incertezze statistiche. Quanto ho detto finora è per spiegare il perché ho ritenuto opportuno usare i minimi quadrati.
PASSIAMO AL PROBLEMA VERO
L'idea è quella di generare casualmente e uniformemente le varie $ y_i $ negli opportuni intervalli di variazione $ [y_i-delta_y,y_i+delta_y] $ (infatti è questo il significato di errore massimo) e a ogni simulazione fare un fit.
Detto più semplicemente, diciamo di fare $ N $ simulazioni, avremo così $ N $ set di dati $ {x_i,y_i}_j $ con $ jin {1,...,N} $ , a cui corrisponderanno altrettanti parametri $ A_j $ e $ B_j $. Stando a quanto dice il professore tali parametri dovrebbero distribuirsi gaussianamente. Purtroppo non è quanto avviene, avendo effettuato un test del $ chi^2 $ ed avendo ottenuto un disaccordo significativo. Specifico che le distribuzioni ottenute sono state confrontate con quelle Gaussiane che minimizzano il $ chi^2 $. Sapete dirmi quali sono le distribuzioni attese? O comunque come ricavare informazioni dalle distribuzioni ottenute?
POSSIBILI SOLUZIONI (???) A CUI HO PENSATO
1. Questo metodo, che mi è stato consigliato dal prof, credo fallisca proprio nella teoria. Infatti per applicare i minimi quadrati si ipotizza che ogni $ y_i $ sia governata da una distribuzione gaussiana, mentre poi quando si applica il Monte Carlo queste vengono distribuite uniformemente.
2. Nelle formule dei minimi quadrati rientrano sostanzialmente SOMMATORIE di PRODOTTI (ad esempio $ sum_i(x_i*y_i) $ ) questo vuole dire che aumentando $ i $ ,per il teorema del limite centrale, dovrei ottenere qualcosa di "sempre più gaussiano".
Mi sono inoltre imbattuto nel problema in cui oltre a $ delta_y $ è presente anche $ delta_x $, in quel caso però non credo abbia senso utilizzare la formula dei minimi quadrati per ogni simulazione. Concluso dicendo che nel mio caso $ iin {1,...,8} $ e $ N=10000 $ .
Ringrazio nuovamente tutti per l'aiuto che spero possiate darmi!
Come da titolo sto cercando un metodo rigoroso che consenta di fare un fit lineare (ma anche no...) utilizzando il metodo Monte Carlo. Mi sono imbattuto in questa situazione dovendo analizzare dei dati presi durante la riproduzione dell'esperienza di Thomson. Nonostante ciò le mie domande sono più prettamente statistiche ed è per questo che scrivo in questa sezione. Detto ciò vi pongo il mio problema ringraziandovi anticipatamente!
Ho delle misure $ (x_i , y_i) $ le quali dovrebbero avere un andamento lineare. Per ogni $ y_i $ ho un'incertezza, sempre la stessa, di tipo massimo: chiamiamola $ delta_y $ . In primo luogo il fatto che le incertezze siano sempre le stesse mi porta a poter assumere che anche le ipotetiche incertezze statistiche (ossia le deviazioni standard) delle singole $ y_i $ siano le stesse per ogni $ y_i $; inoltre il fatto che siano presenti su una sola grandezza (la $ y $) mi pone nelle condizioni di poter effettuare un best-fit ai minimi quadrati. Ricordo che i parametri $ A $ e $ B $ della retta $ y=A+Bx $ sono indipendenti dalle incertezze statistiche. Quanto ho detto finora è per spiegare il perché ho ritenuto opportuno usare i minimi quadrati.
PASSIAMO AL PROBLEMA VERO
L'idea è quella di generare casualmente e uniformemente le varie $ y_i $ negli opportuni intervalli di variazione $ [y_i-delta_y,y_i+delta_y] $ (infatti è questo il significato di errore massimo) e a ogni simulazione fare un fit.
Detto più semplicemente, diciamo di fare $ N $ simulazioni, avremo così $ N $ set di dati $ {x_i,y_i}_j $ con $ jin {1,...,N} $ , a cui corrisponderanno altrettanti parametri $ A_j $ e $ B_j $. Stando a quanto dice il professore tali parametri dovrebbero distribuirsi gaussianamente. Purtroppo non è quanto avviene, avendo effettuato un test del $ chi^2 $ ed avendo ottenuto un disaccordo significativo. Specifico che le distribuzioni ottenute sono state confrontate con quelle Gaussiane che minimizzano il $ chi^2 $. Sapete dirmi quali sono le distribuzioni attese? O comunque come ricavare informazioni dalle distribuzioni ottenute?
POSSIBILI SOLUZIONI (???) A CUI HO PENSATO
1. Questo metodo, che mi è stato consigliato dal prof, credo fallisca proprio nella teoria. Infatti per applicare i minimi quadrati si ipotizza che ogni $ y_i $ sia governata da una distribuzione gaussiana, mentre poi quando si applica il Monte Carlo queste vengono distribuite uniformemente.
2. Nelle formule dei minimi quadrati rientrano sostanzialmente SOMMATORIE di PRODOTTI (ad esempio $ sum_i(x_i*y_i) $ ) questo vuole dire che aumentando $ i $ ,per il teorema del limite centrale, dovrei ottenere qualcosa di "sempre più gaussiano".
Mi sono inoltre imbattuto nel problema in cui oltre a $ delta_y $ è presente anche $ delta_x $, in quel caso però non credo abbia senso utilizzare la formula dei minimi quadrati per ogni simulazione. Concluso dicendo che nel mio caso $ iin {1,...,8} $ e $ N=10000 $ .
Ringrazio nuovamente tutti per l'aiuto che spero possiate darmi!
Risposte
Non sono un esperto di statistica, ma ho usato abbastanza spesso metodi Monte Carlo. Perché se a livello teorico hai deciso che le variabili sono distribuite con una distribuzione Gaussiana hai poi usato una distribuzione uniforme? Anche senza provare a fare i calcoli mi sembra plausibile aspettarsi un comportamento diverso da quello previsto dalla teoria in questo caso. Come generi i campioni?
Se non ho capito male te vuoi simulare dei dati provenienti da un modello lineare del tipo:
$ Y=alpha+betaX+epsi $
Dove $ epsi~ N(0,sigma^2) $ e quindi anche la distribuzione condizionata di Y è normale di tipo $ Y_|x~ N(alpha+betax,sigma^2) $
Quindi scelta una valore per la varianza dovresti simulare i valori di yi usando la normale con la media in funzione di xi.
$ Y=alpha+betaX+epsi $
Dove $ epsi~ N(0,sigma^2) $ e quindi anche la distribuzione condizionata di Y è normale di tipo $ Y_|x~ N(alpha+betax,sigma^2) $
Quindi scelta una valore per la varianza dovresti simulare i valori di yi usando la normale con la media in funzione di xi.
Scusate l'assenza e grazie per le risposte. Forse ho fatto un po' di confusione tra i punti reali e quelli generati. Chiamerò $ y_i $ l'i-esimo punto misurato e $ (y_i)_j $ l'i-esimo punto generato durante la j-esima simulazione.
Per apatriarca. Cerco di spiegarmi meglio, comunque sì, in sostanza è come hai detto tu, infatti questo fatto l'ho incluso tra le possibili soluzioni ma...
Ho un set di dati costituito da $ n $ punti i quali sono approssimativamente disposti su una retta. Ogni punto è stato misurato (lo conosco per via sperimentale insomma), questo significa che ogni punto è soggetto all'errore di sensibilità dello strumento. Il significato di un errore del genere è che noi non possiamo conoscere il valore vero di ciò che stiamo misurando ma siamo sicuri che la nostra misura si trova, per esempio, tra due stanghette di un righello (quello che prima ho chiamato $ [y_i-delta_y,y_i+delta_y] $) . Per quando riguarda la sua distribuzione non possiamo fare altro che dire che è uniforme nel suddetto intervallo. ALDILÀ DI OGNI FIT, dovendo simulare dei dati dovremmo generare i nostri $ y_i $ in maniera uniforme nei rispettivi intervalli. Questo è il Monte Carlo. Il "problema" arriva quando devo trovare, per ogni simulazione di dati, la migliore retta. Io ho scelto i minimi quadrati NON PESATI i quali, tra le altre ipotesi, richiedono che, per un fissato $ j $, ogni $ (y_i)_j $ generato casualmente debba essere distribuito con lo stesso parametro di larghezza $ sigma_y $ intorno a $ (y_i)_j $. ATTENZIONE. Il fatto di fare tale assunzione non vieta di generare un altro set di dati (da pensarsi come un altro esperimento) e di ripetere lo stesso procedimento.
Per Overflow94: purtroppo non riesco a capire cosa mi stai chiedendo. Se hai voglia prova a rispiegarmelo in termini più elementari. Grazie comunque!
P.S. Sapete la routine di randomizzazione di Matlab fino a che punto è efficiente?
Per apatriarca. Cerco di spiegarmi meglio, comunque sì, in sostanza è come hai detto tu, infatti questo fatto l'ho incluso tra le possibili soluzioni ma...
Ho un set di dati costituito da $ n $ punti i quali sono approssimativamente disposti su una retta. Ogni punto è stato misurato (lo conosco per via sperimentale insomma), questo significa che ogni punto è soggetto all'errore di sensibilità dello strumento. Il significato di un errore del genere è che noi non possiamo conoscere il valore vero di ciò che stiamo misurando ma siamo sicuri che la nostra misura si trova, per esempio, tra due stanghette di un righello (quello che prima ho chiamato $ [y_i-delta_y,y_i+delta_y] $) . Per quando riguarda la sua distribuzione non possiamo fare altro che dire che è uniforme nel suddetto intervallo. ALDILÀ DI OGNI FIT, dovendo simulare dei dati dovremmo generare i nostri $ y_i $ in maniera uniforme nei rispettivi intervalli. Questo è il Monte Carlo. Il "problema" arriva quando devo trovare, per ogni simulazione di dati, la migliore retta. Io ho scelto i minimi quadrati NON PESATI i quali, tra le altre ipotesi, richiedono che, per un fissato $ j $, ogni $ (y_i)_j $ generato casualmente debba essere distribuito con lo stesso parametro di larghezza $ sigma_y $ intorno a $ (y_i)_j $. ATTENZIONE. Il fatto di fare tale assunzione non vieta di generare un altro set di dati (da pensarsi come un altro esperimento) e di ripetere lo stesso procedimento.
Per Overflow94: purtroppo non riesco a capire cosa mi stai chiedendo. Se hai voglia prova a rispiegarmelo in termini più elementari. Grazie comunque!
P.S. Sapete la routine di randomizzazione di Matlab fino a che punto è efficiente?
Stiamo parlando di uno strumento digitale per cui l'errore che calcoli è dovuto alla conversione analogico/digitale? Non mi è infatti chiara la ragione per cui stai supponendo che l'errore di misura sia uniforme in un intervallo. Io ho sempre visto utilizzare distribuzioni normali per gli errori strumentali.
Venendo ora al problema vero e proprio. Se ho capito bene, tu stai generando diverse successioni \( (y_i^j \) di valori possibili in base ai tuoi valori sperimentali (scegliendo quindi ogni \(y_i^j\) uniformemente in \([ y_i - \delta, y_i + \delta ]\)). Per ogni successione campione calcoli quindi i parametri \(A_j\) e \(B_j\) con il metodo dei minimi quadrati e vorresti sapere come sono distribuiti. È esatto? Quanti campioni hai calcolato per \( A_j \) e \( B_j \)?
C'è una ragione per cui sei interessato alla loro distribuzione (che probabilmente non sarà normale) e non, ad esempio, al calcolo della loro media (che per il teorema del limite centrale è distribuita come una Gaussiana)?.
Per quanto riguarda Matlab direi che dipende dalle tue necessità. Direi che sono abbastanza ottimizzate e veloci. Nulla ti impedisce di provarle e nel caso passare ad altro se non si rivelassero abbastanza efficienti per le tue necessità.
Venendo ora al problema vero e proprio. Se ho capito bene, tu stai generando diverse successioni \( (y_i^j \) di valori possibili in base ai tuoi valori sperimentali (scegliendo quindi ogni \(y_i^j\) uniformemente in \([ y_i - \delta, y_i + \delta ]\)). Per ogni successione campione calcoli quindi i parametri \(A_j\) e \(B_j\) con il metodo dei minimi quadrati e vorresti sapere come sono distribuiti. È esatto? Quanti campioni hai calcolato per \( A_j \) e \( B_j \)?
C'è una ragione per cui sei interessato alla loro distribuzione (che probabilmente non sarà normale) e non, ad esempio, al calcolo della loro media (che per il teorema del limite centrale è distribuita come una Gaussiana)?.
Per quanto riguarda Matlab direi che dipende dalle tue necessità. Direi che sono abbastanza ottimizzate e veloci. Nulla ti impedisce di provarle e nel caso passare ad altro se non si rivelassero abbastanza efficienti per le tue necessità.
Esatto! Il problema è proprio questo. Perché sono interessato alla distribuzione dei parametri della retta? Più che altro mi è stato detto che dovrebbero essere gaussiani e quindi non trovandomi volevo capire il perché. In secondo luogo sappiamo che in una Gaussiana il valor medio è la migliore stima del valor vero di una grandezza. Nel caso di una distribuzione qualsiasi posso ovviamente calcolare il valor medio, ma non saprei che significato attribuire a questo. Ho calcolato 10.000 campioni per $ A_j $ e $ B_j $. Quello di simulare esperimenti o situazioni in cui si abbia a che fare con relazioni lineari credo sia un problema molto frequente. Per questo credevo ci fosse una teoria nella quale eventualmente fosse contemplata anche la distribuzione dei parametri che si ottengono da un fit. Probabilmente non è così, si tratterà di una qualche distribuzione della quale la gaussiana rappresenta un limite...
Ad ogni modo ho testato la seconda possibile soluzione che ho esposto nel primo post.
In effetti avendo più punti sperimentali la distribuzione ottenuta è decisamente gaussiana. (Nel caso particolare ho inventato 100 punti sperimentali, fatto tutte le simulazioni, e ho eseguito un test del $ chi^2 $ ottenendo circa 2 con 8 gradi di libertà.
Ad ogni modo ho testato la seconda possibile soluzione che ho esposto nel primo post.
"iggi":
2. Nelle formule dei minimi quadrati rientrano sostanzialmente SOMMATORIE di PRODOTTI (ad esempio ∑i(xi⋅yi) ) questo vuole dire che aumentando i ,per il teorema del limite centrale, dovrei ottenere qualcosa di "sempre più gaussiano".
In effetti avendo più punti sperimentali la distribuzione ottenuta è decisamente gaussiana. (Nel caso particolare ho inventato 100 punti sperimentali, fatto tutte le simulazioni, e ho eseguito un test del $ chi^2 $ ottenendo circa 2 con 8 gradi di libertà.
Accedi a tutti gli appunti
Tutor AI: studia meglio e in meno tempo