Markov Chain MonteCarlo errore
Studiando la teoria delle Markov Chain MonteCarlo, ho trovato un espressione per l'errore, nel caso di catena stazionaria, che non mi convince:
[tex]\sigma^2 = var \{ x_i \} + 2 \sum_{k=0}^{ \infty} cov \{ x_i, x_{i+k} \} = \gamma_{0} + 2 \sum_{k=0}^{ \infty} \gamma_k[/tex]
dove
[tex]\gamma_k = \frac{1}{n} \sum_{i=1}^{n-k} (x_i - \mu)(x_{i+k}- \mu)[/tex]
essendo $n$ il numero di punti del vettore [tex]\textbf{x}[/tex], $\mu$ è la media di [tex]\textbf{x}[/tex].
Così facendo mi risulta sempre che $\sigma=0$.
Infatti è facile notare che risulta:
[tex]\sigma^2 = \frac{1}{n} \left [ \sum_{i=1}^{n} (x_i - \mu ) \right ]^2[/tex]
ed è chiaro che l'ultima sommatoria è nulla per definizione.
Dove sbaglio?
[tex]\sigma^2 = var \{ x_i \} + 2 \sum_{k=0}^{ \infty} cov \{ x_i, x_{i+k} \} = \gamma_{0} + 2 \sum_{k=0}^{ \infty} \gamma_k[/tex]
dove
[tex]\gamma_k = \frac{1}{n} \sum_{i=1}^{n-k} (x_i - \mu)(x_{i+k}- \mu)[/tex]
essendo $n$ il numero di punti del vettore [tex]\textbf{x}[/tex], $\mu$ è la media di [tex]\textbf{x}[/tex].
Così facendo mi risulta sempre che $\sigma=0$.
Infatti è facile notare che risulta:
[tex]\sigma^2 = \frac{1}{n} \left [ \sum_{i=1}^{n} (x_i - \mu ) \right ]^2[/tex]
ed è chiaro che l'ultima sommatoria è nulla per definizione.
Dove sbaglio?
Risposte
In termini più avanzati, che forse risultano più chiari a chi si occupa di statistica, quello che non riesco tuttora a capire è il così detto teorema del limite centrale per catene di Markov:
[tex]\sqrt{n} \frac{ \bar{t}_n - E(t( \theta))}{ \tau} \to N(0,1)[/tex]
cioè (se interpreto bene) la differenza tra la media campionaria di un parametro e la sua distribuzione limite tende ad una distribuzione normale con ampiezza [tex]\frac{ \tau}{ \sqrt{n}}[/tex], dove [tex]\tau^2 = \sigma^2 + \sum_{k=1}^{ + \infty} \rho_k[/tex].
Il problema sta nella definizione di [tex]\tau^2[/tex]: esso non è la varianza, bensì la varianza più la somma di tutte le autocovarianze shiftate. In base al ragionamento riportato nel post precedente, supportato da controlli pratici su matlab, mi risulta che la somma di tutte le autocovarianze shiftate di un segnale è pari a meno la metà della sua varianza, così che mi risulta sempre [tex]\tau = 0[/tex].
Ovviamente è assurdo, ci deve essere qualcosa che fraintendo, ma non capisco bene cosa.
[tex]\sqrt{n} \frac{ \bar{t}_n - E(t( \theta))}{ \tau} \to N(0,1)[/tex]
cioè (se interpreto bene) la differenza tra la media campionaria di un parametro e la sua distribuzione limite tende ad una distribuzione normale con ampiezza [tex]\frac{ \tau}{ \sqrt{n}}[/tex], dove [tex]\tau^2 = \sigma^2 + \sum_{k=1}^{ + \infty} \rho_k[/tex].
Il problema sta nella definizione di [tex]\tau^2[/tex]: esso non è la varianza, bensì la varianza più la somma di tutte le autocovarianze shiftate. In base al ragionamento riportato nel post precedente, supportato da controlli pratici su matlab, mi risulta che la somma di tutte le autocovarianze shiftate di un segnale è pari a meno la metà della sua varianza, così che mi risulta sempre [tex]\tau = 0[/tex].
Ovviamente è assurdo, ci deve essere qualcosa che fraintendo, ma non capisco bene cosa.
Il quadrato va alla parentesi tonda, non all'intera sommatoria
Grazie per la risposta, ma quella è la varianza. Invece $\sigma^2$ è definita come la varianza, cioè la somma di tutti i quadrati, più una serie di doppi prodotti, che sono le autocovarianze shiftate. Per cui mi risulta che la $\sigma^2$ è proprio il quadrato di quella sommatoria, che è zero.
Forse in effetti la scelta di chiamare $\sigma$ questo termine fa confondere, ma il libro da cui ho tratto la formula la chiamava così, e non ho pensato a cambiare notazione.
Forse in effetti la scelta di chiamare $\sigma$ questo termine fa confondere, ma il libro da cui ho tratto la formula la chiamava così, e non ho pensato a cambiare notazione.
come l'hai scritta tu è per definizione pari a zero.
Il fatto è che la tua $\sigma^2$ deve essere la somma delle varianze $\gamma_0$ più tutte le autocovarianza a tutti i lag $\gamma_k$. Anche se la sommatoria dovrebbe partire da 1 anzichè da 0, perchè $\gamma_0$ altri non è che la varianza.
Per me hai sbagliato a guardare la formula. Non vedo altre spiegazioni. Dove l'hai presa?
Il fatto è che la tua $\sigma^2$ deve essere la somma delle varianze $\gamma_0$ più tutte le autocovarianza a tutti i lag $\gamma_k$. Anche se la sommatoria dovrebbe partire da 1 anzichè da 0, perchè $\gamma_0$ altri non è che la varianza.
Per me hai sbagliato a guardare la formula. Non vedo altre spiegazioni. Dove l'hai presa?
Grazie per il link.
In pratica però dice le stesse cose che ho trovato sul libro relativamente a questa parte.
In particolare nella prima diapositiva sui windows estimators usa la stessa formula, solo che al posto dei doppi prodotti mette un termine in $+t$ e uno in $-t$ che sono uguali tra loro, e al posto di $x_i$ usa $g(x_i)$, quindi non cambia nulla. Continuo a non capire perchè la somma di tutte le autocovarianze (inclusa la varianza) non è 0.
P.s.: riguardo al fatto che la sommatoria partiva da 1 sono d'accordo.
In pratica però dice le stesse cose che ho trovato sul libro relativamente a questa parte.
In particolare nella prima diapositiva sui windows estimators usa la stessa formula, solo che al posto dei doppi prodotti mette un termine in $+t$ e uno in $-t$ che sono uguali tra loro, e al posto di $x_i$ usa $g(x_i)$, quindi non cambia nulla. Continuo a non capire perchè la somma di tutte le autocovarianze (inclusa la varianza) non è 0.
P.s.: riguardo al fatto che la sommatoria partiva da 1 sono d'accordo.
Leggi bene.
quando fa la sommatori tra -t a +t, il termine t=0 non è nient'altro che la classica varianza!!!!!
quando fa la sommatori tra -t a +t, il termine t=0 non è nient'altro che la classica varianza!!!!!
E che ho detto io?
Non capisco proprio deve sia il problema.
Per capire. tu come trovi che $\sigma^2=1/n [\sum (x_i-\mu)]^2$? fammi vedere i calcoli. Perchè non è corretto
Per capire. tu come trovi che $\sigma^2=1/n [\sum (x_i-\mu)]^2$? fammi vedere i calcoli. Perchè non è corretto
Facciamo per semplicità che la catena contenga solo 3 termini:
$(x_1, x_2, x_3)$
Calcolo le varie autocovarianze:
$\gamma_0 = ((x_1 - \mu)^2 + (x_2 - \mu)^2 + (x_3 - \mu)^2)/3$
$\gamma_1 = ((x_1 - \mu)(x_2- \mu) + (x_2 - \mu)(x_3 - \mu))/3$
$\gamma_2 = ((x_1 - \mu)(x_3 - \mu))/3$
Dunque:
$\sigma^2 = \gamma_0 + 2 \gamma_1 + 2 \gamma_2 =$
$ = 1/3 [(x_1 - \mu)^2 + (x_2 - \mu)^2 + (x_3 - \mu)^2 + 2(x_1 - \mu)(x_2- \mu) + 2(x_2 - \mu)(x_3 - \mu) + 2(x_1 - \mu)(x_3 - \mu)] = $
$ = 1/3 [(x_1 - \mu) + (x_2 - \mu) + (x_3 - \mu)]^2$
Dov'è che sbaglio?
$(x_1, x_2, x_3)$
Calcolo le varie autocovarianze:
$\gamma_0 = ((x_1 - \mu)^2 + (x_2 - \mu)^2 + (x_3 - \mu)^2)/3$
$\gamma_1 = ((x_1 - \mu)(x_2- \mu) + (x_2 - \mu)(x_3 - \mu))/3$
$\gamma_2 = ((x_1 - \mu)(x_3 - \mu))/3$
Dunque:
$\sigma^2 = \gamma_0 + 2 \gamma_1 + 2 \gamma_2 =$
$ = 1/3 [(x_1 - \mu)^2 + (x_2 - \mu)^2 + (x_3 - \mu)^2 + 2(x_1 - \mu)(x_2- \mu) + 2(x_2 - \mu)(x_3 - \mu) + 2(x_1 - \mu)(x_3 - \mu)] = $
$ = 1/3 [(x_1 - \mu) + (x_2 - \mu) + (x_3 - \mu)]^2$
Dov'è che sbaglio?
Niente?
Scusate se risollevo la questione, ma sto cercando ovunque, su libri, internet, non riesco proprio a risolvere la questione.
A me risulta che la sommatoria di tutte le covarianze shiftate di un segnale, inclusa la classica varianza è 0.
E Matlab me lo conferma pure.
Le soluzioni a cui ho pensato sono:
- si devono considerare i quadrati delle covarianze, ma questo non risulta dalla formula data sul libro di statistica e su internet;
- la media va stimata in un altro modo, ma come? In ogni caso penso che qualsiasi stima della media mi porterà ad avere covarianze sia positive che negative, mentre ho la sensazione che debbano essere tutte positive per avere senso il teorema, quindi ritorno all'opzione dei quadrati.
A me risulta che la sommatoria di tutte le covarianze shiftate di un segnale, inclusa la classica varianza è 0.
E Matlab me lo conferma pure.
Le soluzioni a cui ho pensato sono:
- si devono considerare i quadrati delle covarianze, ma questo non risulta dalla formula data sul libro di statistica e su internet;
- la media va stimata in un altro modo, ma come? In ogni caso penso che qualsiasi stima della media mi porterà ad avere covarianze sia positive che negative, mentre ho la sensazione che debbano essere tutte positive per avere senso il teorema, quindi ritorno all'opzione dei quadrati.
A me viene in mente che $\gamma_1$ va diviso per 2 e non per 3, mentre $\gamma_3$ va diviso per uno
volevo dire $\gamma_2$ e non $\gamma_3$
"niandra82":
A me viene in mente che $\gamma_1$ va diviso per 2 e non per 3, mentre $\gamma_2$ va diviso per uno
Ci avevo pensato in un primo momento, eppure il libro non solo scrive la formula con $n$ al denominatore, che potrebbe essere un errore di stampa, ma subito dopo specifica:
It might be thought that one should divide by $n-k$ instead of $n$, but the large $k$ terms are already very noisy so dividing by $n-k$ only makes a bad situation worse.
Non è certo il massimo della chiarezza questa frase, però sembra che il problema non sia cosa metto a denominatore, anche se mettendo $n-k$ almeno evito che faccia automaticamente $0$.
Le autocorrelazioni non si considerano mai tutte quante, ma ci si ferma quando perdono di significatività. Per essere più chiari la sommatoria della prima formula del post, non dovrebbe essere tra meno e più infinito, ma da più o meno c, dove c è un valore "ragionevole". Non ha senso considerare autocorrelazioni stimate con 1 o 2 valori.
Il libro ti dice che per diminuire il rumore conviene usare sempre n, ma e vuoi essere preciso devi cambiare il denominatore.
Penso che il problema possa essere questo
Il libro ti dice che per diminuire il rumore conviene usare sempre n, ma e vuoi essere preciso devi cambiare il denominatore.
Penso che il problema possa essere questo
Se non ricordo male, nelle dispense che ti avevo indicato, questa cosa è spiegata
"niandra82":
Le autocorrelazioni non si considerano mai tutte quante, ma ci si ferma quando perdono di significatività. Per essere più chiari la sommatoria della prima formula del post, non dovrebbe essere tra meno e più infinito, ma da più o meno c, dove c è un valore "ragionevole". Non ha senso considerare autocorrelazioni stimate con 1 o 2 valori.
Probabilmente hai ragione, ma mi stupisco che un libro apposito per le MCMC, dove l'autore sembra preciso in maniera fin troppo stressante commetta un errore così grossolano. Soprattutto ho trovato decine di fonti che commetterebbero lo stesso errore (magari copia e incolla?).
"niandra82":
Se non ricordo male, nelle dispense che ti avevo indicato, questa cosa è spiegata
Nelle dispense sostiene che considerare tutti lag fino ad infinito crea dei problemi, ma da nessuna parte viene evidenziata l'incosistenza che ho trovato io, ed essendo una cosa piuttosto semplice, mi sarei aspettato di trovarla almeno in una fonte.
Oltretutto devo dire che comunque considerando tutti i lag, l'aspetto della funzione di autocorrelazione ha una forma dipendente dal tipo di segnale che sto studiando, anche ad alti lag, per quanto stimati con 1 o 2 punti.
Comunque non so che dire per ora, ho bisogno di indagare meglio, anche se si dovesse decidere di tagliare ad esempio quando tagliare, ecc.. Visto quanto ci sto impazzendo su, aggiornerò quando avrò il tutto più chiaro, chissà possa servire ad altri. In ogni caso grazie per l'aiuto.
Il criterio per tagliare è la significatività. Se un'autocorrelazione non è significativa, allora la puoi assumere pari a zero. Arrivati a un lag c, troverai sicuramente che lag del tipo c+k, cn k positivo e intero saranno non significativi e quindi la sommatoria non è più infinita.
Visto che stiamo parlando di Markov chain Monte carlo, assumo che tu le stia studiando per stimare modelli bayesiani e quindi, per un probema semplice deve essere almeno di 2000 simulazioni. Se la catena è fatta bene, deve essere poco correlata, e di solito l'autocorrelazioni si annullano per lag > 5 più o meno, se invece è fatta male, puoi arrivare ad autocorrelazioni anche di 30 o 40, ma non di più. Se hai autocorrelazioni con lag maggiori di 50 allora devi rivedere l'algoritmo, perchè non funziona benissimo.
Quindi nella tua sommatoria il valore limite c è di 50. e l'ultima autocorrelazione la stimi con 1950 osservazioni. Adesso se la sommatoria la dividi per 2000 o 1950, non è che cambia poi molto.
Visto che stiamo parlando di Markov chain Monte carlo, assumo che tu le stia studiando per stimare modelli bayesiani e quindi, per un probema semplice deve essere almeno di 2000 simulazioni. Se la catena è fatta bene, deve essere poco correlata, e di solito l'autocorrelazioni si annullano per lag > 5 più o meno, se invece è fatta male, puoi arrivare ad autocorrelazioni anche di 30 o 40, ma non di più. Se hai autocorrelazioni con lag maggiori di 50 allora devi rivedere l'algoritmo, perchè non funziona benissimo.
Quindi nella tua sommatoria il valore limite c è di 50. e l'ultima autocorrelazione la stimi con 1950 osservazioni. Adesso se la sommatoria la dividi per 2000 o 1950, non è che cambia poi molto.
In realtà ho posto il problema in questi termini, cioè di MCMC, perchè pensavo fosse più noto in questo campo, ma non è qui che lo sto applicando.
Quello che cerco di fare è di stimare l'errore di un modello che si adatta a una serie temporale, assumendo che i residui presentino del "rumore colorato". Le autocorrelazioni non vanno a zero dopo un certo lag, ma si verificano due casi:
- se il residuo è piatto, ho che il grafico delle autocorrelazioni è molto caotico, con alcune leggermente più significative a certi lag, probabilmente corrispondenti a disturbi periodici del modello, dovuti a effetti strumentali;
- se il residuo presenta un trend, il grafico delle autocorrelazioni è una funzione convessa, con valori alti positivi a piccoli lag e grandi lag, valori alti negativi per i lag centrali.
Quello che cerco di fare è di stimare l'errore di un modello che si adatta a una serie temporale, assumendo che i residui presentino del "rumore colorato". Le autocorrelazioni non vanno a zero dopo un certo lag, ma si verificano due casi:
- se il residuo è piatto, ho che il grafico delle autocorrelazioni è molto caotico, con alcune leggermente più significative a certi lag, probabilmente corrispondenti a disturbi periodici del modello, dovuti a effetti strumentali;
- se il residuo presenta un trend, il grafico delle autocorrelazioni è una funzione convessa, con valori alti positivi a piccoli lag e grandi lag, valori alti negativi per i lag centrali.