Alla scoperta di bboot, il package R che implementa il block stationary bootstrap

Alla scoperta di bboot, il package R che implementa il block stationary bootstrap

Una delle tecniche di ricampionamento con reimmissione più utilizzate per fare inferenza statistica su un parametro è il metodo bootstrap. Un importante limite di questa metodologia è che non può essere utilizzata ogniqualvolta si ha un campione y = [y(1), y(2), …, y(n)] di numerosità n i cui dati risultano essere correlati tra di loro. Infatti, il campionamento singolo da y per ottenere l’i-esimo campione bootstrap  y(i)* = [y(i,1)*, y(i,2)*, …, y(i,N)*] di lunghezza N, comporta la rottura del legame esistente tra osservazioni correlate quando i dati non sono indipendenti. Per ovviare a questo problema, in letteratura è stata proposta la tecnica block stationary bootstrap (BSB), un’estensione della metodologia base che campiona da y a blocchi piuttosto che lavorando sulle singole unità. 

Il block stationary bootstrap risulta essere particolarmente indicato per la simulazione di serie storiche stazionarie, ove le loro statistiche non dipendendo dal tempo t, ma dalla distanza temporale tra le osservazioni che costituiscono y. Quindi, se L viene scelto in modo tale che ogni singolo blocco sia sufficientemente esaustivo per descrivere l’autocorrelazione esistente nei dati, allora i blocchi possono essere liberamente accostati per ottenere una serie storica y(i)* che rispecchi le caratteristiche statistiche di quella originale.

Per MATLAB è già disponibile una funzione che implementa il BSB (stationaryBB, scaricabile tramite il MATLAB file transfer), mentre per R ancora non esiste, pertanto abbiamo deciso di realizzare bboot, un package avente il compito di fornire ai programmatori R una soluzione pronta all’uso del block stationary bootstrap tramite l’omonima funzione blockboot(). La necessità di questo package nasce dal progetto AgrImOnIA. Data la serie storica di un inquinante, volevamo simulare i dati mancanti tramite il BSB, così da misurare le prestazioni in validazione dei modelli ARIMA utilizzati per la ricostruzione dei missing values. I dati mancanti sono dovuti ai malfunzionamenti delle centraline adibite alla misurazione delle concentrazioni di inquinanti nell’aria. Prima che un guasto venga risolto possono passare alcuni giorni, di conseguenza i missing values occorrono nelle serie storiche non in maniera isolata, ma a blocchi; da ciò deriva la necessità di utilizzare il block stationary bootstrap piuttosto che la tecnica base di ricampionamento.

La funzione MATLAB che implementa il BSB è implementata come ricampionamento con ripetizione, ovvero i blocchi non sono disgiunti tra di loro, un comportamento indesiderato per quello che è il fine di utilizzo del block stationary bootstrap per questo caso specifico. Pertanto, blockboot() è stata implementata come ricampionamento senza ripetizione, ovvero la selezione dei blocchi è tale che all’iterazione j-esima non possa essere scelto un blocco contenente indici di osservazioni che sono stati precedentemente scelti. Ovviamente, questo comportamento si verifica solo quando N, ossia la lunghezza del singolo campione bootstrap y(i)*, è minore di n, cioè la lunghezza del campione originale; infatti, se N è maggiore di n, allora il vincolo di non sovrapposizione deve essere rilassato affinché y* raggiunga la lunghezza desiderata avendo a disposizione solo n dati (è necessario scegliere delle osservazioni che sono state precedentemente selezionate).

In questo momento bboot non è ancora disponibile su CRAN, l’archivio ufficiale di R, tuttavia è possibile consultare il seguente repository GitHub:

https://github.com/lamferzon/bboot

Per l’installazione è sufficiente eseguire in R i seguenti comandi:

>> devtools::install_github(“lamferzon/block-bootstrap-for-R”)

>> library(bboot)

Il package è ben documentato per assicurare un più vasto utilizzo, pertanto si rimanda all’help di blockboot() per avere ulteriori informazioni in merito ai parametri della funzione e al suo funzionamento. Degna di nota la possibilità di ricevere in output non sono gli indici delle osservazioni che non sono state campionate, ma anche un summary contenente, per ognuna delle simulazioni, il valore di alcuni parametri prestazionali:

 – num. iterations: numero di iterazioni dell’algoritmo necessarie per ottenere una serie storica simulata di lunghezza N;

 – num. blocks: numero di blocchi creati;

– mean blocks length: lunghezza media dei blocchi; è bene notare che di default la lunghezza di ogni singolo blocco viene scelta casualmente da una distribuzione poissoniana di parametro λ (da specificare in input);

– num. out-of-indexes: numero di iterazioni (perse) in cui è stato selezionato un blocco contenente osservazioni i cui indici superano quello dell’ultimo dato nel campione y;

– num. rejected indexes: numero di interazioni (perse) in cui è stato selezionato un blocco contenente osservazioni i cui indici sono già stati selezionati;

– execution time: tempo necessario per generare y*.

Di blockboot() è stata svolta anche un’analisi prestazionale. I grafici in Figura 2 e 3 mostrano l’andamento della lunghezza media dei blocchi e del tempo di esecuzione in funzione del rapporto N/n in due casi:

1. blocks rejection disattivata per N maggiore di n (Figura 2): accettazione di blocchi il cui indice iniziale è già stato selezionato come primo indice per blocchi che già sono stati creati.

2. blocks rejection attivata per N maggiore di n (Figura 3): rifiuto di blocchi il cui indice iniziale è già stato selezionato come primo indice per blocchi che già sono stati creati.

L’analisi è stata svolta con K = 5, L = 3 e distribuzione poissoniana delle lunghezze L.

Nel primo caso, è possibile notare che il tempo di esecuzione (punti in arancione nel grafico) aumenta esponenzialmente all’avvicinarsi a 1 del rapporto N/n, ossia quando si vuole generare una simulazione la cui lunghezza coincide esattamente con quella del campione originale; i blocchi non possono sovrapporsi, pertanto ogni singolo dato di y dev’essere scelto una e una sola volta affinché N sia uguale a n. Il numero di indici che può essere scelto come idx(0) alla j-esima iterazione diventa sempre più piccolo all’aumentare di j e la lunghezza L (scelta casualmente) deve essere tale da non provocare una sovrapposizione con i blocchi che già sono stati creati, quindi il numero di iterazioni perse diventa sempre più grande, provocando pertanto un aumento del tempo di esecuzione; nonostante ciò, l’algoritmo arriva a convergenza in un tempo finito. Per N maggiore di n, invece, non c’è alcun vincolo sulla scelta del primo indice idx(0) perché il blocks rejection è disattivato. A ogni iterazione le osservazioni sono sempre tutte disponibili, pertanto non possono verificarsi rejected indexes; ciò si traduce in un tempo di esecuzione prossimo a zero. Per quanto riguarda la lunghezza media dei blocchi (punti in verde nel grafico), essa risulta essere sempre vicina alla scelta L tranne per N prossimo a n, situazione in cui la lunghezza dei blocchi diventa sempre più piccola all’aumentare di j: la probabilità che molti dei dati nelle vicinanze di idx(0) alla j-esima iterazione siano già stati inseriti in blocchi che già sono stati creati è elevata, quindi è necessario che L sia piccolo (o addirittura uguale a 1) affinché il blocco corrente non venga rifiutato.

Infine, per il secondo caso valgono in generale le considerazioni fatte precedentemente, ma c’è una sostanziale differenza: per N maggiore di n il tempo di esecuzione aumenta all’aumentare del rapporto N/n. Questo comportamento è causato dall’impossibilità di scegliere come primo indice idx₀ alla j-esima iterazione un’osservazione che è già primo indice di un blocco che è già stato generato in una delle iterazioni precedenti dell’algoritmo (il blocks rejection è attivato), di conseguenza il numero di osservazioni disponibili diminuisce progressivamente all’aumentare di j, provocando un aumento del numero di iterazioni perse e quindi del tempo necessario per generare una delle k simulazioni y(i)* di lunghezza N.

Dott. Lorenzo Leoni
Studente magistrale di Ingegneria Informatica, indirizzo DSDE
Università degli Studi di Bergamo
l.leoni2@studenti.unibg.it