R, l’analisi delle serie storiche partendo da Copenaghen

Eccoci di nuovo a parlare di R, linguaggio ed un ambiente open source per il calcolo statistico e la rappresentazione grafica in ambito statistico.

Anche questa volta faccio frutto del “dovere” per scrivere questo articolo per “piacere”.

Nel titolo leggete Copenaghen, mi spiego meglio.

Questa volta parliamo di come effettuare analisi delle serie storiche con il software R.

Per analizzare una serie storica dobbiamo averne una. Io che ogni tanto sbarco sul sito del settimanale internazionale ho memoria di una serie storica il cui grafico è apparso in un post sul sito di internazionale.

Ho cercato l’articolo e l’ho ritrovato: http://www.internazionale.it/home/?p=12879 .

(Potremmo anche riportare l’intero articolo dato che Internazionale fa uso di Creative Commons Attribution-Non-Commercial-Share Alike 2.5.)

Bene, ho pensato, perchè mai cercare i dati delle vendite di una gelateria alla ricerca di trend e pattern stagionali quando possiamo cercare serie storiche sulla concentrazione della co2 visto che siamo nel dopo Copenaghen ed il problema del riscaldamento terrestre e più che mai di interesse?

Ho cercato quindi tramite i motori di ricerca un documento che contenesse i dati sui rilevamenti di co2 di diversi anni campionati mensilmente ed ho trovato questo documento: http://cdiac.ornl.gov/ftp/trends/co2/maunaloa.co2 .

Proviamo quindi a studiare con R i dati relativi alle concentrazioni mensili di Co2 dal 2000 al 2008.

Prima di poter lavorare sui dati dobbiamo importarli in R e per farlo copiamo i valori mensili in un’unica riga in un foglio di calcolo; io ho usato Calc.

Una volta trascritti i dati in un’unica riga, salviamo il foglio in formato csv e scegliamo, al momento di salvare, la Tabulazione come delimitatore di campo.

Adesso possiamo importare i dati in R e iniziarne lo studio.

Con il seguente comando importiamo i dati in R e assegniamo ad un oggetto data frame i valori della serie storica.

>A=read.table(file=”/home/andrea/Documents/Flandoli/homework2/SeriStoricaCo2.csv”)

Per visualizzare una serie storica abbiamo bisogno di avere un vettore da passare come argomento al comando ts.plot().

Quindi digitiamo:

>X=t(A[1,])

cosi da trasporre il vettore.

X ora conterrà i 108 valori corrispondenti ai valori mensili di co2 dal 2000 al 2008.

Ora usiamo il comando ts.plot.

>ts.plot(X)

A questo punto ci apparirà il grafico della serie storica.

ts.plot(X)

ts.plot(X) - X contiene la serie storica delle concentrazioni di co2 dal 2000 al 2008

Già da una analisi visiva si desume un trend per cui i valori crescono quasi linearmente e una periodicità stagionale con un periodo vicino a 12 mesi.

Con R possiamo produrre il grafico della funzione di autocorrelazione della serie storica con il comando acf(X), ottenendo il grafico della figura seguente.

>acf(X)

Grafico della funzione di autoccrelazione

acf(X) - grafico della funzione di autocorrelazione con R

Come è giusto che sia, in zero ovvero con una traslazione di 0, la correlazione è pari a 1.

Il grafico della funzione di autocorrelazione mostra che vi è una alta correlazione tra la serie storica e la sua traslata di 11 e 12 unità e questo ci è utile per impostare poi il modello di previsione lineare. I picchi sono oltre l’intervallo di confidenza quindi possiamo prenderli in considerazione.

L’impressione sulla periodicità della serie storica non era errata ma tutto sommato era evidente, ed in questo modo ne abbiamo avuto conferma tramite la funzione di autocorrelazione calcolata da R.

Se assegniamo i risultati del comando acf(X) ad una variabile Y con il comando

>Y<-acf(X)

oltre ad avere il grafico della funzione di autocorrelazione avremo in Y i valori della funzione di autocorrelazione che appunto risultano corrispondenti al grafico e valere 0.711 per un periodo di 11 mesi e 0.706 per un periodo di 12 mesi.

Grafico della funzione di autocorrelazione in R

Grafico della funzione di autocorrelazione in R

Avendo 108 valori nella serie storica e volendo visualizzare la presenza di un trend lineare, possiamo costruire un vettore con gli indici temporali ed effettuando la regressione lineare visualizzare intercetta e coefficiente così da poter costruire una retta da sovrapporre al grafico della serie storica che evidenzi il trend.

>Z<-c(1:108)

>reg<-lm(X~Z) // con questo comando R applica il modello di regressione lineare tra la variabile Z che contiene i valori temporali e X che contiene la serie storica ed in più associamo il risultato prodotto da R ad un oggetto reg.

>cor(X,Z) // con il comando cor è possibile calcolare la correlazione tra due variabili
>summary(reg) // con questo comando visualizziamo una serie di informazioni risultati dal calcolo della regressione lineare effettuato da R
F<-(0:1080)/10
b=reg$coefficients[1]
//associamo il valore dell’intercetta ad un oggetto b
a=reg$coefficients[2]
G=a*F+b
//costruiamo la retta con coefficiente angolare a
ts.plot(X) //
lines(F,G) // il comando lines permette di sovrapporre un grafico ad un altro grafico precedentemente plottato difatti se abbiamo chiuso la finestra del grafico precedentemente plottato il comano restituirà un errore.

Con la sequenza di comandi appena scritti otteniamo il seguente risultato:

Serie di comandi per sovrapporre un retta ottenuta tramite la regressione lineare al grafico della serie storica

Ora possiamo costruire un modello che descriva i dati della serie storica.

Ho scelto di costruire e confrontare tre modelli lineari che descrivano i dati visto anche il grafico della funzione di autocorrelazione il quale mostra una periodicità alta tra un periodo di 11 e 12.

1)xk=a1x(k-1)+a12x(k-12)+b+ek

2)xk=a1x(k-1)+a11x(k-11)+b+ek

3)xk=a1x(k-1)+a11x(k-11)+a12x(k-12)+b+ek

modello2 dei dati

modello 3 dei dati

La varianza spiegata nel terzo modello, che prende in considerazione una periodicità 11 e 12 insieme è maggiore.

Dai p-value molto vicini al valore zero si deduce che le ipotesi che le variabili non siano utili come predittori sono quasi nulle.

Possiamo usare quindi i valori a1, a2, a3 e b presi dalla regressione lineare del terzo modello per effettuare le previsioni poiché il terzo modello spiega maggiore varianza dei dati x0 da 13 a n e “fidandomi” anche del grafico della funzione di autocorrelazione.

Previsione dei valori futuri

Associando ad un vettore P i valori della serie storica presenti nel vettore X con un ciclo for possiamo calcolare i valori dei mesi del 2009 utilizzando come modello di previsione lineare quello costruito con i coefficienti e l’intercetta del terzo modello di descrizione dei dati della serie storica che ho scelto.


>x01<-X[13:n]
>x11<-X[12:(n-1)]
>x21<-X[1:(n-12)]
>REG1<-lm(x01~x11+x21)
>summary(REG1)
>x02<-X[13:n]
>x12<-X[12:(n-1)]
>x22<-X[2:(n-11)]
>REG2<-lm(x02~x12+x22)
>summary(REG2)

>x03<-X[13:n]
>x13<-X[12:(n-1)]
>x23<-X[2:(n-11)]

>x33<-X[1:(n-12)]
>REG3<-lm(x03~x13+x23+x33)
>summary(REG3)

>P<-X
>length(P)
>a11<-REG3$coefficients[2]
>a21<-REG3$coefficients[3]
>a31<-REG3$coefficients[4]
>b<-REG3$coefficients[1]
>for (k in (n+1) : (n+12) ) {
P[k]=a11*P[k-1]+a21*P[k-12]+a31*P[k-11]+b
}
>length(P) // ora p dovrebbe risultarvi composto da 120 valori, i primi 108 della serie + dodici dalla previsione per

// il 2009 basata sul modello lieneare

ts.plot(P,col=”red”)
lines(X,col=”blue”)

Dal grafico seguente si distingue la previsione per l’anno 2009 di valori di co2 che è di colore rosso.

Previsione con modello lineare

Previsione con modello lineare

Tra le  mlle potenzialità di R c’è la possibilità di calcolare le previsioni con il metodo di Holt-Winters.

Possiamo usare i comandi di R per predire l’anno 2009 e generare un grafico che sovrapponga la previsione con il metodo di Holt-Winters (in giallo) con quella calcolata con il modello lineare ottenendo il risultato visibile nella figura seguente.

>x.new<-ts(data=X,frequency=12) // il comando HoltWinters ha bisogno che all’oggetto contenente la serie storica sia  //associato il parametro che ne indica la periodicità quindi in realtà stiamo confrontando il modello lineare con quello di //Holt-Winters considerando una frequenza pari a 12
>HW<-HoltWinters(x.new)
>predHW<-predict(HW,12)
>C<-X
>for (k in 1 : 12 ) {
C[n+k]=predHW[k]
}
>ts.plot(P,col=”red”)
>lines(C,col=”yellow”)
>lines(X,col=”blue”)

Si nota bene già graficamente che le due previsioni per quanto simili in grandi linee non sono le stesse.

Training Set e Set Test per il confronto tra il modello lineare e quello di Holt-Winters

Avendo a disposizione dati per molto anni possiamo prendere i valori della serie storica fino al 2007 e calcolare le previsioni per il 2008 con il terzo modello lineare di cui prima e con il metodo Holt-Winters per metterle a confronto con i dati veri della serie storica del 2008.

Quindi partendo da n=n1+n2 dati della serie ho usiamo n1=96 dati per training set e n2=12 come test set.

>X1<-X[1:96]
>n1<-length(X1)
>P1<-X1
>for (k in (n1+1) : (n1+12) ) {
P1[k]=a11*P1[k-1]+a21*P1[k-12]+a31*P1[k-11]+b
}

>x1.new<-ts(data=X1,frequency=12)
>HW1<-HoltWinters(x1.new)
>predHW1<-predict(HW1,12)
>C1<-X1
>for (k in 1 : 12 ) {
C1[n1+k]=predHW1[k]
}

length(C1)

P12<-P1[97:108]
X12<-X[97:108]
C12<-C1[97:108]

>ts.plot(P12,col=”red”)
>lines(C12,col=”yellow”)
>lines(X12,col=”blue”)

Plottando soltanto i valori del 2008, quelli veri e quelli previsti con i due modelli, otteniamo il seguente risultato mostrato in figura dove con il colore blu sono rappresentati i dati della vera serie storica con il rosso i dati previsti con il terzo modello lineare e in giallo i dati previsti con il modello Holt-Winters.

Confronto previsioni

Confronto previsioni

Di prima impressione i dati previsti con il modello lineare sembrano essere più aderenti a quelli della serie storica vera.

In tutti i casi conviene calcolare la deviazione standard sigma usandola come metodo per verificare l’errore e per vedere quale dei due modelli di previsione si è comportato meglio in questo caso presentando una deviazione minore.

>sigma1<-sqrt(mean((P12-X12)^2))
>sigmaHW<-sqrt(mean((C12-X12)^2))
>sigma1
>sigmaHW

Sigma calcolato con R vale circa 0.49 per il modello lineare e 0.58 con il modello HW quindi l’intuizione grafica era corretta.

Intervallo di confidenza

Il valore di sigma può essere utilizzato per stabilire un intervallo di confidenza entro il quale accettare future previsioni.

Ipotizzando errori gaussiani e un largo margine di confidenza possiamo ipotizzare che il valore vero futuro ricada in una previsione pk + o – un valore 3*sigma scelto tra quelli precedentemente calcolato.

Scegliendo il sigma calcolato con il modello lineare si può decidere che i prossimi valori della serie storica molto probabilmente rientreranno in un valore p compreso tra pk+1.47 e pk-1.47.

Abbiamo visto soltanto qualche altra piccola funzionalità offerta da R.

Cercando dati a proposito della serie storica sulla co2 ho trovato un articolo molto interessante che riguarda lo strano crescere del trend delle rilevazioni noaa http://globalwarming.blog.meteogiornale.it/2009/06/17/anomalie-oceaniche-lo-strano-caso-dei-dati-noaa/ . (sempre in tema Copenaghen)

Vi ringrazio dell’attenzione e vi lascio esplorare R sperando, che al di là della propria maggiore o minore esperienza con il calcolo statistico, riusciate a sfruttare questo ambiente open source, a presto

Andrea Stani

R, open source per il calcolo statistico.

Ho pensato di unire con quest’articolo l’utile e il dilettevole.

Devo, e con piacere perchè lo sto trovando stimolante, studiare un pò di probabilità insieme alla soluzione software R.

Dal sito del progetto leggiamo che R è un linguaggio ed un ambiente per il calcolo statistico e la rappresentazione grafica in ambito statistico.

Sul sito del progetto troviamo tutto quanto può servire per approfondire le conoscenze di R: screenshot, wiki, manuali, faq e tanto altro.

Prima di tutto è importante dire che R è un progetto open source rilasciato con licenza GPL v.2 ed è multipiattaforma quindi disponibile per linux, mac os x e windows.

Quindi se conoscete solo windows e siete vissuti con esso, sappiate che non siete costretti a cambiare il vostro sistema operativo per utilizzare R.

Mi sono accorto che uno strumento per il calcolo statistico è di interesse per molti laureati che pur non avendo studiato nell’ambito del ramo dell’informazione conoscono la statistica più che me ma non sanno di strumenti che gli facilitino il calcolo e la rappresentazione.

Con questo articolo voglio provare ad introdurre ad R coloro che potrebbero trovare utile questo strumento.

Possiamo quindi scaricare R dalla pagina dei download del progetto scegliendo la locazione da cui scaricare il software più vicina a noi; naturalmente essendo distribuito con licenza gpl troviamo disponibili anche i sorgenti e nel momento in cui scrivo, l’ultima versione di R disponibile è la 2.10.0 rilasciata il 26 ottobre 2009.

Io personalmente ho installato R su Debian Lenny tramite il pacchetto “r-base” presente nel repository di lenny e mi ritrovo quindi installata la versione 2.7.1.

In tutti i casi possiamo scaricare da quella che viene indicata come rete CRAN, i cui server sono a vostra disposizione tramite i collegamenti della sezione download, l’utima versione binaria già compilata per alcune delle più diffuse distribuzioni linux rese disponibili di recente tra cui troviamo Debian Lenny, Ubuntu Karmic, Fedora 11 e Suse 11.1.

Mi perdonino gli utenti windows qualora si sentano trascurati ma immagino che l’installazione su lo stesso windows non sia differente da quella di altri software; se avete problemi provate a comunicarmelo e cercherò d’essere d’aiuto per quanto mi è possibile; con molta probabilità per windows bastera scaricare l’installer e cliccare due volte su di esso per avviare l’installazione di R che sarà poi disponibile nel menu di avvio di windows.

Se stiamo usando linux lanciamo R da terminale con il comando “R”.

Per questa volta cerchiamo di arrivare a tracciare una semplice densità Gaussiana.

Ci accorgiamo che R funziona benissimo senza alcuna interfaccia grafica ma questo non deve spaventare perchè non è così ostico da usare come può sembrare a chi non è abituato a software con interfaccia testuale.

Appena avviato R (nella versione che ho installato su Debian, ma non differisce molto per le altre), ci compare un messaggio che ci comunica la versione che stiamo utilizzando, ci dice che il software è un software libero e che è rilasciato senza alcuna garanzia, e ci propone di utilizzare dei primi comandi quali:

  • license()
  • help()
  • help.start()
  • q()

(Installando ora R su Ubuntu 9.10, a qualche giorno di distanza dalla prima bozza di questo articolo, il messaggio di benvenuto che appare una volta avviato R non “consiglia” da subito i comandi di help ma poco importa, questi ci sono sempre; interessante dovrebbe essere il suggerimento ad installare il pacchetto revolution-r il quale è una versione ottimizzata per Ubuntu del pacchetto r-revolution-revobase che dovrebbe contenre funzioni per il controllo della gestione del multithreading. Ora come ora non so niente su come funzionino queste funzioni ed in tutti casi se usate windows penso che questa parentesi non vi riguardi per niente)

Bene se ci siamo spaventati digitiamo “q()” senza le virgolette e premiamo il tasto invio, se ci viene chiesto di salvare qualcosa, digitiamo “n” per ora e saremo fuori del tutto.

Siamo usciti da R in tal caso, la paura è passata e possiamo riprovarci sapendo che non stiamo facendo nulla di male.

Se usiamo il comando license() vedremo il tipo di licenza con cui è distribuita la versione di R che stiamo utilizzando. Nel mio caso per la versione 2.7.1 e mi viene detto che questa versione è rilasciata con GPL 2.

(nei repository di Karmik è disponibile la versione 2.9.2)

Per avere informazioni sui comandi disponibili abbiamo a disposizione due tipi di help.

Se digitiamo solo help() ci viene descritto il comando che ci permette di cercare informazioni senza appogiarci ad un browser  mentre con help.start() verrà aperta con il nostro browser predefinito una pagina html all’interno della quale troviamo un link “Search Engine & Keywords” con il quale possiamo facilmente ricercare ciò che ci interessa.

Sulla stessa pagina html che ci compare con il comando help.start() troviamo con il link packages la documentazione a tutti i pacchetti che compongono R e le funzioni/comandi da questi fornite.

R è composto da diversi pacchetti che forniscono i vari comandi/funzioni ed è possibile estendere le sue funzionalità installando tanti altri pacchetti provenienti da persone che contribuiscono allo sviluppo del software; con il comando contributors() possiamo scorrere i nomi e i contatti mail di chi ha contribuito allo sviluppo.

In particolare i pacchetti (d’ora in poi package) più utilizzati sono base, graphics e stats.

Proviamo a utilizzare per ora la ricerca tramite pagina html ricordando che R parla inglese e quindi se cerchiamo quali sono i comandi che riguardano la distribuzione di Gauss o Normale dobbiamo cercare “normal“.

Cerchiamo la parola chiave “normal” e troviamo i comandi che hanno a che vedere con la distribuzione normale o di Gauss.

Avviata la ricerca ci compariranno vari risultati con abbinata una descrizione e tra questi troviamo: Normal con descrizione “The Normal Distribution“.

Cliccando sul link Normal ci appare una pagina descrittiva nella quale viene detto che i comandi inerenti la distribuzione normale sono forniti dal package stats e ne viene spiegata la sintassi.

Leggiamo ad esempio che il comando dnorm() ha quattro parametri che sono: una variabile di tipo vettore dei quantili per cui vogliamo che la densità sia definita, la moda, la varianza e un variabile logica chiamata log che può valere TRUE o FALSE e permette di ottenere, se posta a TRUE, come risultati del comando il logaritmo delle probabilità al posto delle probabilità stesse.

Se l’argomento log viene omesso il comando funzionerà parimenti e considererà log=FALSE.

Ad es. una invocazione del comando potrebbe essere dnorm(x,0,1) ed in tal caso se abbiamo definito un vettore x di quantili “abbastanza denso” composto magari da 2000 punti tra i valori -4 e 4 otterremo come output del comando una serie di valori che compongono il vettore delle densità di probabilità corrispettive per i vari quantili forniti al comando e che rappresentati su di un grafico ci darebbero una bella distribuzione normale o di gauss.

Prima di utilizzare il comando però abbiamo bisogno di definire quantomeno un vettore di punti x; vediamo quindi come fare.

Se vogliamo definire un vettore di 8001 punti compreso tra -4 e 4. dobbiamo assegnarlo alla variabile x ad esempio con il seguente comando:

x<-(-4*1000:4*1000)/1000

se avessimo voluto 801 punti avremmo dovuto effettuare l’assegnamento cosi:

x<-(-4*100:4*100)/100

o anche

x=(-4*100:4*100)/100

se ora vogliamo vedere il vettore di valori che è stato assegnato alla variabile x digitiamo x in R e premiamo invio.

Cia appariranno tante righe di valori e ogni riga inzierà con un valore tra parentesi crescente del tipo [1] per la prima riga che non è un valore del vettore delle x ma l’indice che identifica il valore nel vettore; se infatti voglia sapere il valore del cinquecentosessantesimo quantile del vettore x possiamo digitare x[560] e premere il tasto invio.

A questo punto digitanto i comandi a seguire avremo un vettore di quantili x e un vettore di densità y accoppiati; se non facciamo attenzione con i valori che assegnamo R ci dirà che la dimensione dei due vettori non corrisponde:

x<-(-4*1000:4*1000)/1000

assegnamo ad y i valori della densità con il prossimo comando

y<-dnorm(x,0,1)

e poi se vogliamo visualizzare questi punti su di un piano cartesiano digitiamo:

plot(x,y)

Dovrebbe ora apparire il grafico della nostra normale standard con moda e media uguali a 1 e deviazione standard pari a 1.

Quando conosciamo già il nome di un comando, ad esempio dnorm per dnorm(), e non ricordiamo alcuni parametri/argomenti possiamo cercare subito la sintassi di questo comando digitando help(dnorm) e premendo invio, senza dover di nuovo cercare la voce normal.

Se proviamo infatti a digitare help(plot) troveremo informazioni su come dire al comando plot di inserire dei titoli per gli assi o anche colorare il grafico se il nero non ci piace.

Potremmo ad esempio decidere di usare il comando cosi:

plot(x,y,col=”red”,xlab=”quantili”,ylab=”densità”)

e otterremmo in tal caso una bella curva di gauss rossa con gli assi x e y intitolati rispettivamente quantili e densità.

Schermata di R con digitati i comandi per disegnare una gaussiana standard

Plot di una curva di gauss disegnata con R

Magari cercherò di scrivere qualcos’altro a proposito di R prossimamente, in tutti i casi chi conosce la statistica più di me a questo punto dovrebbe potersi divertire parecchio.

Sperando al solito di non aver commesso troppi errori vi auguro buon divertimento con R,

Andrea Stani