I SEMI ANTICHI a cura di Angelo Bianchi, erborista

I SEMI ANTICHI di Angelo Bianchi

Molti cibi che arrivano sulle nostra tavola sono in gran parte prodotti dall’industria. Per soddisfare ed anche per imporre certi criteri esclusivamente commerciali, viene ritenuto necessario pianificare le colture in modo che si abbiano piante adatte alla raccolta, alla conservazione ed alla trasformazione delle derrate e che tengano conto delle esigenze dei processi di lavorazione e commercializzazione delle grosse catene commerciali.

        Questa prassi assai diffusa non prevede, se non in minima parte, la salvaguardia dell’ambiente, la biodiversità genetica.  Il risultato di questa tendenza impone, tra l’altro, che i piccoli coltivatori siano dipendenti dalle sementi ( in larga parte ibridi e OGM) controllate dai grandi complessi agroindustriali ( per es. Monsanto ed altri) e costretti poi di conseguenza ad utilizzare i prodotti di sintesi di cui le stesse industrie sono produttrici. Infatti tali sementi avendo in qualche modo perso la naturale robustezza e rusticità, a causa delle manipolazioni atte a privilegiarne la produttività, necessitano per essere coltivate, trattamenti massicci di pesticidi per combatterne le avversità e necessitano anche apporti notevoli di concimi chimici per sostenerne la crescita e lo sviluppo.
         Tra i piccoli agricoltori si è diffusa la volontà di controbattere questa tendenza dell’industria agroalimentare e si è cercato di recuperare varietà di sementi che possano garantire la biodiversità, il sapore e la salubrità del cibo, e , non ultima, la possibilità di autoprodursi le sementi necessarie, cercando di svincolarsi, per quanto possibile, dalla morsa dei condizionamenti imposti dall’agroindustria.
          Questi agricoltori, insieme ad alcuni ricercatori universitari, hanno cercato di recuperare semi di orticole e di cereali accantonati, nei tempi recenti, perché meno produttivi e meno adatti alle colture intensive anche se notevolmente resistenti alle avversità ed in  grado di produrre piante particolarmente ricche di valori nutritivi.
           Questa ricerca riguarda soprattutto le varietà locali adatte per uno specifico terreno e clima, assolutamente diverse delle sementi costruite per la monocoltura che invece sono le stesse per differenti climi, terreni e latitudini.
            Nei casi di particolari avversità climatiche oppure di avversità causate da parassiti e crittogame, le monocolture vengono praticamente distrutte in modo assai grave, in quanto tutte le piante sono sottoposte allo stesso rischio, mentre nel caso di una promisquità di diverse varietà ci saranno alcune piante che soccomberanno, ma altre invece che riusciranno a resistere e quindi a garantire in ogni caso il raccolto. Alcuni agricoltori si sono attrezzati anche per fornire alla popolazione locale prodotti provenienti da queste varietà “antiche” ed a garantire quindi alimenti assai più ricchi di valori nutritivi, esenti da residui tossici e rispettosi dell’ambiente e della biodiversità.
              Questa piccola rete “commerciale” inoltre permette di limitare i trasporti che creano non pochi danni all’ambiente e notevoli spese che vanno ad incidere in maniera sostanziale sul prezzo del prodotto.
               Attraverso questi semplici meccanismi si può cercare in qualche modo di modificare le abitudini e gli acquisti delle persone, in gran parte condizionati dalla pubblicità e dalle mire di chi ha come scopo solo l’interesse economico.
Angelo Bianchi
Erborista
Annunci

FLORA E VEGETAZIONE NEL VOLTERRANO del dott. Juri David Bettini

Il seguente articolo del dott. Juri Bettini è stato trasferito da “Il Sillabario” 4-1996  e  1-1997, inserti della Comunità di Pomarance.
juri10001

Per vedere la tab. N.1 cliccare su:

ARTICOLO1

juri0002

CONTROLLA IL CONTO! ESEMPIO DI ANALISI STATISTICA CON EXCEL SU PIOGGE E TEMPERATURE MENSILI NEL VOLTERRANO a cura del dott. Piero Pistoia

POST IN VIA DI COSTRUZIONE: LA VIA SI FA NELL’ANDARE! e, con Foerster e, per certi versi, Bruner, la cultura non si comunica ma si costruisce insieme.

Per questo lavoro si ringrazia il dott. David Bettini per aver fornito nel 1996, circa venti anni fa, i  dati mensili di piogge e temperature del Volterrano per trenta anni a partire dal 1956.

RIFLESSIONI E CONTI  SU PIOGGE E TEMPERATURE MENSILI DAL 1956 AL 1986, CALCOLO E VARIAZIONE DEGLI  INDICI CLIMATICI DI BAGNOULS E GAUSSEN COL TEMPO,  ATTRAVERSO MEDIE MOBILI, FATTORI STAGIONALI ED ALTRO SULLE SERIE STORICHE a cura del dott. Piero Pistoia

Ecco i ‘conti’ con Excel che sviluppai appunto una ventina d’anni fa:

PIOGGE_VOLTERRA 1956-1986

Sulla cartella di lavoro di Excel appaiono 7 grafici (N°: 1-2; 3; 4; 7; 8; 9; 11) e 7 fogli di lavoro (N°: 18; 1; 2; 3; 5; 6; 7), alla rinfusa a zibaldone, su cui discuteremo nel seguito. Intanto il lettore potrà trovare una guida alla ‘lettura’ dei precedenti processi (individuando le formule utilizzate riportate) su altri esempi di statistica, anche  con l’uso di Excel, dello stesso autore, per es. nel post “UN PARZIALE PERCORSO DI BASE SULL’ANALISI DI UNA SERIE STORICA REALE…”, uno zibaldone di statistica e linguaggi informatici (la statistica ‘raccontata’ con Excel, con il Basic, il Mathematica di Wolfram e l’R), ecc..

DA CORREGGERE!

Intanto cerchiamo di ricavare dal link i dati di pioggia e temperatura in Excel, per memorizzarle in un file .csv, leggibile da R. Iniziamo con i dati di pioggia. Si clicca sul link: PIOGGE_VOLTERRA_1956_1986, riportato sopra.

Appare una console di Excel con riportate varie colonne, nominate sopra; si va al Foglio 1 (7 colonne e un grafico degli Effetti Stagionali) e si isola intanto la seconda colonna  che contiene 360 dati mensili della pioggia da 1956 al 1986, eliminando le altre, memorizzandola poi col suffisso .CSV, nel file VOLTERRA_PIOGGE_1956_1986_A (ricordarsi in quale Memoria di Massa), con le seguenti opzioni: Si mantenga il formato corrente – Tipo di carattere occidentale – Separatore di campo la “virgola”.

Si entra nella console di R e si guarda se il processo funziona.

setwd(“X:/”)

in X va la lettera della memoria di massa; nel nostro caso setwd(“I:/”) o ultimamente G-

Si fanno tentativi per far leggere da  R  i dati con suffisso .xls

Altri scritti e argomentazioni seguiranno successivamente

ECCO UN PRIMO ‘ASSAGGIO’ ORGANIZZATO DEL RACCONTARE  L’ARTICOLARSI DEL CONTENUTO DEL LINK IN XLS CON IL LINGUAGGIO DI R, INIZIANDO IN PARTICOLARE DALLA COLONNA DUE CHE CONTIENE I DATI .

In definitiva tenterò di raccontare una storia statistico-informatica usando un linguaggio informatico alternativo; i diversi trucchi informatici con i loro risvolti logico-razionali necessari aumenteranno la possibilità di tentare veloci prove diversificate, attivando una maggiore concentrazione sui concetti statistici, visti da punti di vista diversi, insieme alla loro memorizzazione e alla loro assimilazione, oltre ad un intenso ammaestramento informatico.

Come in un gioco di puzzle a più vie, si attivano costruzioni con i ‘mattoni’  statistici a più possibilità e si costruiscono gli stessi concetti statistici con ‘mattoni’ informatici diversi logicamente, individuando percorsi diversi  e il percorso spesso diventa una facility  per l’apprendimento.

SCRIPTS IN R E COMMENTI

setwd(“I:/”)
dataset=read.csv(“PIOGVOLTR0.CSV”, header=T, dec=”.”, sep=”;”). La lettera ‘0’  è uno zero!

#Da tener d’occhio l’attributo dec=”.” dell’argomento !?

#COME SI COSTRUISCE IL FILE PIOGVOLTR0.CSV a partire da foglio 1 di Excel:

# 1 – Dal post “CONTROLLA IL CONTO…” si carica, cliccando sul link
#”PIOGGE_VOLTERRA 1956_1986″, l’analisi in .xls su questi dati di cui parla il post
# e questo al fine di costruire dal foglio 1 un data frame leggibile da R.

# 2 – E’ un foglio di Excel, ma siamo in Open Office; mi pongo sul foglio 1
# (siglato anche col nome del link)

# 3 . Cambio le intestazioni originali e le sostituisco col nome di variabili neutre
#X.1 X.2 …X.7 per le sette colonne del foglio; ne chiariremo nel proseguio
#il contenuto.

# 4 – Salvo questo foglio in xls con l’opzione di Open Office “Testo CSV(.csv)”
# nel file, in questo caso, PIOGVOLTR0 che ricaricheremo con R; ognuno
#può chiamarlo col nome che vuole e memorizzarlo sul disco che vuole, basta
#ricordasi il nome; noi abbiamo chiamato il file come detto e lo memoriziamo in
# un disco rimovibile indicato con I o G. Durante la memorizzazione rispondiamo
#alle domande, separatore di campo “;” , separatore di testo Apice ‘ ecc.

# 5 – Si copi il presente testo sulla console di R e controlla che giri; il nostro gira! Per un po’!! Vedere nel proseguo dell’argomentazione.

Oppure, per  preparare i racconti successivi, si memorizzano tutte le sette colonne (cioè l’intera pagina di excel), ponendo al posto delle intestazioni che figuravano nella pagina, le seguenti sette una per colonna: X.1,X.2….X.7, alla testa di esse.

 

PREMESSA E DESCRIZIONE DEL PROCESSO

ECCO IL NOSTRO PROGETTO IPOTETICO CON OUTPUTS

attach(dataset)
X.1

Dovrebbe apparire  la prima colonna

Dovrebbero apparire i 6 valori iniziali del dataset

head(dataset)

X.2[1:6]

In effetti appaiono tre NA iniziali prima dei valori del 1956; per il resto Ok; con

X.2[4:375], si pensa di correggere; appaiono i 372 valori effettivi delle piogge misurate! La corroborazione dell’ipotesi sembrava abbastanza scontata. Il dataframe col processo descritto portava, nominando la seconda variabile di colonna, ai dati da analizzare! Il processo aveva funzionato.

Invece da qui la sorpresa non prevista. Se cerco di memorizzare  X.2[4-375] in una variabile (es.,dataset2 ), cioè

dataset2=X.2[4:375] e guardo i dati in essa (battendo dataset2 sulla console di R) i contenuti sono stranamente completamente cambiati! ts.plot(dataset2) dà un grafico che sembra non riguardare più la nostra prova, comunque continuiamo la loro analisi per cercare anche di capire. Il processo comunque è un modello abbastanza trasferibile ad altri insieme di dati! Poi con calma cercheremo di entrare in possesso di questi dati reali in qualche modo, al limite copiandoli direttamente in un vettore  di R. Insomma nel trasferimento del contenuto conosciuto di una variabile ad un’altra….entriamo in una zona caotica imprevista, almeno per ora. VEDREMO poi!

_________________________________________________________

ECCO IL RIQUADRO DELLO  SBAGLIO NON VOLUTO!

Da una revisione sui processi nella sintassi del  richiamo del file in .CSV, fra gli altri attributi, c’era dec (separatore dei decimali) =; qui fra virgolette avevo inserito il punto. Dovevamo dire che invece c’era la virgola: dec=”,”. Ora il programma gira. Col punto il programma pensava che le scritte dei numeri decimali fossero testo,  attivando le funzioni factor e level.

ECCO IL RIQUADRO DELLA CAUSA

I am sorry. Troppe direttrici culturali, aperte e diverse fra loro, da ordinare, creano una situazione entropica densa e dispersiva che, coniugata al tempo che ora scorre rapido ed ai problemi sempre più numerosi in quest’ultimo scorcio della vita, aumentano una richiesta di concentrazione poco spiegabile per un lavoro semplicemente hobbistico, in totale assenza di contributi, che vengono invece elargiti, per la cultura, da organi sociali; per passare il tempo insomma! Vedrò.

Comunque , ‘si parva licet componere magnis’ (Virgilio ‘Georgiche, IV), ritengo che un apporto culturale che non si esaurisca in ‘racconti e descrizioni’, ma proceda con punti interrogativi, caso proprio di questo blog, non sia inferiore a quello dei molti centri attivi benificiati da contributi, ed anche mi sembra che le stesse ‘lectiones magistrales’ dei gruppi di eccellenza, che si susseguono nei palchi, spesso, nel migliore dei casi, si ‘spengano’ in se stesse, comunicando poco, e, negli altri casi si riducano ad escamotages pubblicitari per ottenere finanziamenti.

In attesa …per i lettori curiosi, se ci sono, ho intanto ricopiato direttamente in un vettore di R i 372 dati delle piogge mensili, su cui sarà possibile accedere, nell’immediato, alla loro analisi, ‘divertendosi’ nel controllare conti e processi, cioè nel fare cultura

dataset=c(89.2,32,68,106,39.8,69.4,33.4,25.2,64.4,41.4,123.8,34.4,
66.8,100.2,17,109.8,159.8,23.6,27,5.6,2.4,74.2,100.8,79,
47.4,27,133,108.8,36,30.6,22,21,12.8,163.8,67,134.4,
61.2,35.6,133.2,69.2,120.8,25.4,4.6,46,46.2,74.6,73.4,149.2,
57.8,112.8,113.4,62,9.6,91.4,64.6,5.2,179.6,214.6,132.2,210.2,
107.8,30,0.6,152.8,37.6,73,18.4,1.5,87.6,160.4,135.8,86.4,
50.4,46.6,139.6,35.6,35.6,21.6,37.6,50,64.8,161.6,196.8,35.4,
147.6,89.4,75.3,74.2,45.8,89.2,28.8,38.8,94.6,95.2,100.4,88.8,
4.6,77.6,157.6,58.8,49.6,32.2,84,69.5,34,193,78,118.4,
127.8,16.8,91,107.4,50.4,36,80,24.8,119.8,2.6,240.2,104.8,
97.8,88.4,28.4,28.4,36.2,28.6,67.2,52.2,122.8,173.4,259,59.6,
49.2,42.2,56.4,11.4,84.2,83.6,19.4,75.8,99.2,42.8,99.8,62.6,
57,149.8,29.8,49,156,91.4,28,65.8,35.6,75,148.2,68,
89,156.2,85.8,40,56.6,43.4,38.4,83.8,128.6,26.4,157.2,106.4,
104.8,98.2,140.6,79.6,66.6,63.8,18.6,22.4,75,28.2,106.4,95,
81,51.4,55,53,122.4,67.2,15.8,2.2,54.8,17.4,149,8.6,
85.4,91.8,43,106.8,51.4,26,71,49,136.6,76,46.8,67.2,
70.8,33,13,48.2,17.2,56.2,38.8,22.6,171,80,73.2,30.6,
35.4,105.2,66,80.6,52.8,27.2,12.6,61.4,88.8,87.6,64.2,16.6,
16.2,39.6,101.6,92.2,47.6,74.4,11.4,92.8,51.2,125.2,103.8,102.2,
12.4,55,92.8,75.8,28.6,32.6,86,94.4,110.8,154.6,69.4,175.4,
61.2,108,54,15,101,12,29.4,100.2,33.2,31.4,65.4,60.4,
102.2,96.6,68,172.2,66.8,67,87.6,23.2,18.6,88,38.8,103,
127.4,75.8,50.8,84.8,1.4,80,9.4,131.4,84.8,110.4,126.6,108,
50.2,20.4,125.8,34,102.6,66.6,22.8,13,1,142.4,182.2,79.8,
64.2,25.6,85.2,74.4,31.4,13.4,45,15.8,99.8,172.2,1.6,152.2,
33,40.8,64.2,23.2,55.6,34.4,3.4,59.4,70.4,142.6,210.4,100,
23,199.6,128,65.6,22.8,34.2,8,192,12.4,73.2,17.6,37,
72.4,65.6,83.8,91.4,161,87.2,3.3,131.4,91,139.6,142.8,86,
93.8,67.6,147.2,9.2,94.4,25,15,72.2,0.4,58,117.8,32.6,
102.4,113.4,125.6,122,1.4,121.8,70,7.4,26.8,13,63.6,53.2)

Comunque lo sviluppo a partire da questo vettore verrà inserito in links  successivi.  Al termine aggiungeremo anche l’outputs senza errori come garanzia del lavoro

DATASET_PIOGGE_VOLT_24_5_ stag1_stag2_OUT_OK1

DATASET_PIOGGE_VOLT_24_5_ stag1_stag2_OUT

 Seguono i 5 grafici del link precedente a quello sopra, costruiti nell’ OUTPUT relativo. Dovremo aggiungere anche altri tre grafici  relativi all’output del nuovo ultimo links.

Nella serie dataset  ipotizziamo che sia assente un trend, come suggerito dal GRAF.1

GRAF_PIOG_VOLT0001

GRAF_PIOG_VOLT0002

GRAF_PIOG_VOLT0003

GRAF_PIOG_VOLT0004

GRAF_PIOG_VOLT0005

GRAF.4 analisi

I links successivi riguarderanno la Stagionalità, Fattori stagionali, Effetti stagionali.

piog_volt_5_8_graf0001

piog_volt_5_8_graf0002

piog_volt_5_8_graf0003

piog_volt_5_8_graf0004

piog_volt_5_8_graf0005

piog_volt_5_8_graf0006

 

_______________________________________________

Prove statistiche  e tentativi su ipotesi al fine di costruire un modello trasferibile.

#Ipotizzo che con una media mobile a tre termini pesata, #possa eliminare in buona parte i randoms da questa strana #serie dataset2 (già senza il trend come dal grafico).

yt2=dataset2; n2=length(yt2); mbt2=c()
for(t in 2: n2-1){mbt2[t]=(yt2[t-1]+2*yt2[t]+yt2[t+1])/4}
ts.plot(mbt2)

#Ho smussato dall’originale (dataset2 )i randoms per cui avrei ottenuto
#una serie senza i randoms (con plausibili stagionalità+ciclo; assente il trend
#iniziale come si vedeva ad occhio dal grafico di dataset2.
#Così se applico una media mobile di ordine dodici a mbt2 (primo caso), mi aspetto di trovare una serie con solo il ciclo.
#Potevo applicare la media mobile 12 direttamente su dataset2, smussandola
#di stagionalità e forse anche di randoms, e poi, togliendo
#la stagionalità+ randoms (questa nuova serie) da dataset2 avrei ottenuto solo
# il ciclo, perché all’inizio non possedeva trend .
#Protocollo sperimentale per  il modello: da confrontare questi due processi che dovrebbero condurre ambedue al ciclo.

#Questa serie con solo il ciclo la chiamo mbt2_12.

#PROVIAMO IL CASO 1:
yt=mbt2; n1=length(mbt2)-2; mbt2_12=c()
for(t in 7: n1-6){mbt2_12[t]=(yt[t-6]/2+yt[t-5]+yt[t-4]+yt[t-3]+yt[t-2]+
yt[t-1]+yt[t]+yt[t+1]+yt[t+2]+yt[t+3]+yt[t+4]+yt[t+5]+yt[t+6]/2)/12}

ts.plot(mbt2)
#stagionalità + randoms
lines(mbt2_12)

#CICLO; è la serie mbt2 smussata della stagionalità
#se tolgo stagionalità + ciclo – ciclo, mbt2-mbt2_12, ottengo la stagionalità.

#Stiamo pensando anche di scrivere direttamente i 372 valori

#della di X.2[4:375], nella variabile dataset2 con

#dataset2=c(………)

 

LA VERSIONE CORRETTA + I RELATIVI  CINQUE GRAFICsetwd(“F:/”)

dataset=read.csv("PIOGVOLTR0.CSV", header=T, dec=",", sep=";")

#Da osservare l'attributo dec=","! nell'argomento della funzione "read.csv".
 
#COME SI COSTRUISCE IL FILE PIOGVONTR.CSV a partire da foglio 1 di Excel.
 
par(ask=T) 
attach(dataset) 
X.1 
#Stampa i 372 dati della 1a colonna delle 7 colonne del DATA.FRAME 
#contenuto nel file PIOGVOLTR0.CSV (indicazione degli anni a partire dal 1956 
#intervallati da NA head(dataset) # X.1 X.2 X.3 X.4 X.5 X.6 X.7 
#1 1956 89,2 -3,3 70,62580645 #2 NA 32 2,0 73,94193548 #3 NA 68 10,2 83,02258065 
#4 NA 106 -2,0 72,30322581 #5 NA 39,8 -9,6 62,67741935 #6 NA 69,4 -21,0 52,52903226 
#Stampa i valori delle prime sei righe del data.frame, costituito da 
#sette colonne. 
X.2[4:375] 
#Stampa i 372 dati della 2a colonna in mm di pioggia mensili 
#partendo dall'anno 1956 con gennaio  
X.2[1:6] 
#Stampa i primi sei valori della seconda colonna 
#[1] 89,2 32 68 106 39,8 69,4
dataset=ts(dataset) 
#considera il data.frame dataset come una serie storica dataset2=X.2[1:372] 
#prende 1 valori da 1 a 372 del vettore X.2 e li mette 
#nella variabile dataset2 
dataset2=ts(dataset2) 
#dataset2 è una serie storica ts.plot(dataset2) 
#GRAF.1 #Stampa la serie storica dataset2; sembra assente il trend. PIOGGE_Vo_GRAF0001yt2=dataset2; n2=length(yt2); mbt2=c() 
for(t in 2: n2-1){mbt2[t]=(yt2[t-1]+2*yt2[t]+yt2[t+1])/4} 
ts.plot(mbt2)
# GRAF.2 
#disegno il grafico di mbt2, GRAF.2, cioè i dati originali (senza trend) 
#privati anche dei random s.l.

GRAF.2
 
PIOGGE_Vo_GRAF0002#Penso di smussare cioè dataset2 dai randoms s.l.; nel vettore mbt2 è plausibile 
#siano contenuti dati relativi a stagionalità e ciclo. 
#Ho smussato dall'originale (dataset2 )i randoms s.l. per cui avrei ottenuto 
#una serie senza i randoms (con plausibili stagionalità+ciclo); assente il trend 
#iniziale come si vedeva ad occhio dal grafico di dataset2. 
#Così se applico una media mobile di ordine dodici a mbt2 mi aspetto di trovare 
#una serie con solo il ciclo (PRIMO CASO). 
#Potevo applicare la media mobile 12 direttamente su dataset2 (SECONDO CASO), 
#smussandola dalla stagionalità e forse anche dai randoms (?), e poi, togliendo 
#la stagionalità + randoms (questa nuova serie) da dataset2 avrei ottenuto solo 
# il ciclo. 
#Da confrontare questi due processi che dovrebbero condurre ambedue al ciclo. 
#Questa seconda serie con solo il ciclo la chiamo mbt2_12. 
#PROVIAMO IL PRIMO CASO (applico una media modile 12 su mbt2: 
yt=mbt2; n1=length(mbt2)-2; mbt2_12=c() 
for(t in 7: n1-6){mbt2_12[t]=(yt[t-6]/2+yt[t-5]+yt[t-4]+yt[t-3]+yt[t-2]+ yt[t-1]+
yt[t]+yt[t+1]+yt[t+2]+yt[t+3]+yt[t+4]+yt[t+5]+yt[t+6]/2)/12} 
ts.plot(mbt2_12)
#GRAF.3 Disegno il grafico del ciclo

#GRAF.3
 
PIOGGE_Vo_GRAF0003
ts.plot(mbt2) # + CICLO1 sovrapposto. 
#Disegno il grafico di mbt2 (stagionalità + ciclo) e sovrappongo mbt2_12 (ciclo): 
lines(mbt2_12) 
#è la serie mbt2 smussata della stagionalità 
#Insieme al grafico mbt2 (stagionalità+ciclo) sovrappongo il ciclo: GRAF.4 PIOGGE_Vo_GRAF0005#Se tolgo il ciclo da stagionalità + ciclo (mbt2), ottengo mbt2-mbt2_12, 
#cioè la stagionalità. 
#Sorge il problema che mbt2 e mbt2_12 debbono avere la stessa lunghezza

#per poterli sottrarre 
#FACCIO DELLE PROVE PER RENDERE I VETTORI LUNGHI UGUALE 
#Controllo mbt2 
length(mbt2) # 371 = 12 
head(mbt2) # NA 1.00 70.75 177.50 202.75 159.00 
mbt2=mbt2[2:(length(mbt2)-1)] 
#Controllo mbt2_12 
length(mbt2_12) # 363 
head(mbt2_12) 
# NA NA NA NA NA NA 
#Impongo che mbt2 e mbt2_12 abbiano la stessa lunghezza per sottrarli 
mbt2_12=mbt2_12[7: (length(mbt2_12)-6)] 
length(mbt2_12) 

#Proviamo il SECONDO CASO (applico direttamente la Mb12 su dataset2) 
#per il calcolo del ciclo 
yt=dataset2; n1=length(dataset2); mbt2_12_0=c() 
for(t in 7: n1-6){mbt2_12_0[t]=(yt[t-6]/2+yt[t-5]+yt[t-4]+yt[t-3]+yt[t-2]+ yt[t-1]+
yt[t]+yt[t+1]+yt[t+2]+yt[t+3]+yt[t+4]+yt[t+5]+yt[t+6]/2)/12} 
ts.plot(mbt2_12_0) 
#Disegno il grafico del ciclo nel secondo modo: GRAF.5 PIOGGE_Vo_GRAF0006 PIOGGE_Vo_GRAF0007#Da confrontare i due grafici del ciclo mbt2_12_0 e mbt2_12. 
#I due cicli praticamente coincidono. Ma nel secondo sono rimasti 
#più errori randoms (?), nel senso che la Mb12 sugli originali 
#praticamente elimina da essi solo la stagionalità. Per cui
#nella serie nuova rimarrà ciclo+randoms e se sottraggo 
#Ciclo+Randoms dall'originale (dataset2) otterrei 
#direttamente la Stagionalità. #Potremmo fare un test statistico per controllo. 
#CALCOLO DELLA STAGIONALITA' nel proseguo (SECONDA PARTE) 

___________________________________________ OUTPUTS di R 
> rm(list=ls(all=TRUE)) 
> setwd("F:/") 
> dataset=read.csv("PIOGVOLTR0.CSV", header=T, dec=",", sep=";") 
>#COME SI COSTRUISCE IL FILE PIOGVONTR.CSV a partire da foglio 1 di Excel 
> par(ask=T) 
> attach(dataset) 
>#The following objects are masked from dataset (pos = 3): 
>X.1, X.2, X.3, X.4, X.5, X.6, X.7  
> X.1 
[1] NA NA NA 1956 NA NA NA NA NA NA NA NA NA NA NA 
[16] 1957 NA NA NA NA NA NA NA NA NA NA NA 1958 NA NA 
[31] NA NA NA NA NA NA NA NA NA 1959 NA NA NA NA NA 
[46] NA NA NA NA NA NA 1960 NA NA NA NA NA NA NA NA 
[61] NA NA NA 1961 NA NA NA NA NA NA NA NA NA NA NA 
[76] 1962 NA NA NA NA NA NA NA NA NA NA NA 1963 NA NA 
[91] NA NA NA NA NA NA NA NA NA 1964 NA NA NA NA NA 
[106] NA NA NA NA NA NA 1965 NA NA NA NA NA NA NA NA 
[121] NA NA NA 1966 NA NA NA NA NA NA NA NA NA NA NA 
[136] 1967 NA NA NA NA NA NA NA NA NA NA NA 1968 NA NA 
[151] NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA 
[166] NA NA NA NA NA NA 1970 NA NA NA NA NA NA NA NA 
[181] NA NA NA 1971 NA NA NA NA NA NA NA NA NA NA NA 
[196] 1972 NA NA NA NA NA NA NA NA NA NA NA 1973 NA NA 
[211] NA NA NA NA NA NA NA NA NA 1974 NA NA NA NA NA 
[226] NA NA NA NA NA NA 1975 NA NA NA NA NA NA NA NA 
[241] NA NA NA 1976 NA NA NA NA NA NA NA NA NA NA NA 
[256] 1977 NA NA NA NA NA NA NA NA NA NA NA 1978 NA NA 
[271] NA NA NA NA NA NA NA NA NA 1979 NA NA NA NA NA 
[286] NA NA NA NA NA NA 1980 NA NA NA NA NA NA NA NA 
[301] NA NA NA 1981 NA NA NA NA NA NA NA NA NA NA NA 
[316] 1982 NA NA NA NA NA NA NA NA NA NA NA 1983 NA NA 
[331] NA NA NA NA NA NA NA NA NA 1984 NA NA NA NA NA 
[346] NA NA NA NA NA NA 1985 NA NA NA NA NA NA NA NA 
[361] NA NA NA 1986 NA NA NA NA NA NA NA NA NA NA NA 
> #Stampa i 372 dati della 1a colonna delle 7 colonne del DATA.FRAME 
> #contenuto nel file PIOGVOLTR0.CSV (indicazione degli anni a partire dal 1956 
> #intervallati da NA 
>  head(dataset) X.1 X.2 X.3 X.4 X.5 X.6 X.7 
1 NA NA NA NA NA NA NA 
2 NA NA NA NA NA NA NA 
3 NA NA NA NA NA NA NA 
4 1956 89.2 NA NA NA -3.3 70.62581 
5 NA 32.0 NA NA NA 2.0 73.94194 
6 NA 68.0 NA NA NA 10.2 83.02258 
> # X.1 X.2 X.3 X.4 X.5 X.6 X.7 
> #1 1956 89,2 -3,3 70,62580645 
> #2 NA 32 2,0 73,94193548 > 
>#3 NA 68 10,2 83,02258065 
> #4 NA 106 -2,0 72,30322581 
> #5 NA 39,8 -9,6 62,67741935 
> #6 NA 69,4 -21,0 52,52903226 
> #Stampa i valori delle prime sei righe del data.frame, costituito da 
> #sette colonne. 
> X.2[4:375] 
[1] 89.2 32.0 68.0 106.0 39.8 69.4 33.4 25.2 64.4 41.4 123.8 34.4 
[13] 66.8 100.2 17.0 109.8 159.8 23.6 27.0 5.6 2.4 74.2 100.8 79.0 
[25] 47.4 27.0 133.0 108.8 36.0 30.6 22.0 21.0 12.8 163.8 67.0 134.4 
[37] 61.2 35.6 133.2 69.2 120.8 25.4 4.6 46.0 46.2 74.6 73.4 149.2 
[49] 57.8 112.8 113.4 62.0 9.6 91.4 64.6 5.2 179.6 214.6 132.2 210.2 
[61] 107.8 30.0 0.6 152.8 37.6 73.0 18.4 1.5 87.6 160.4 135.8 86.4 
[73] 50.4 46.6 139.6 35.6 35.6 21.6 37.6 50.0 64.8 161.6 196.8 35.4 
[85] 147.6 89.4 75.3 74.2 45.8 89.2 28.8 38.8 94.6 95.2 100.4 88.8 
[97] 4.6 77.6 157.6 58.8 49.6 32.2 84.0 69.5 34.0 193.0 78.0 118.4 
[109] 127.8 16.8 91.0 107.4 50.4 36.0 80.0 24.8 119.8 2.6 240.2 104.8 
[121] 97.8 88.4 28.4 28.4 36.2 28.6 67.2 52.2 122.8 173.4 259.0 59.6 
[133] 49.2 42.2 56.4 11.4 84.2 83.6 19.4 75.8 99.2 42.8 99.8 62.6 
[145] 57.0 149.8 29.8 49.0 156.0 91.4 28.0 65.8 35.6 75.0 148.2 68.0 
[157] 89.0 156.2 85.8 40.0 56.6 43.4 38.4 83.8 128.6 26.4 157.2 106.4 
[169] 104.8 98.2 140.6 79.6 66.6 63.8 18.6 22.4 75.0 28.2 106.4 95.0 
[181] 81.0 51.4 55.0 53.0 122.4 67.2 15.8 2.2 54.8 17.4 149.0 8.6 
[193] 85.4 91.8 43.0 106.8 51.4 26.0 71.0 49.0 136.6 76.0 46.8 67.2 
[205] 70.8 33.0 13.0 48.2 17.2 56.2 38.8 22.6 171.0 80.0 73.2 30.6 
[217] 35.4 105.2 66.0 80.6 52.8 27.2 12.6 61.4 88.8 87.6 64.2 16.6 
[229] 16.2 39.6 101.6 92.2 47.6 74.4 11.4 92.8 51.2 125.2 103.8 102.2 
[241] 12.4 55.0 92.8 75.8 28.6 32.6 86.0 94.4 110.8 154.6 69.4 175.4 
[253] 61.2 108.0 54.0 15.0 101.0 12.0 29.4 100.2 33.2 31.4 65.4 60.4 
[265] 102.2 96.6 68.0 172.2 66.8 67.0 87.6 23.2 18.6 88.0 38.8 103.0 
[277] 127.4 75.8 50.8 84.8 1.4 80.0 9.4 131.4 84.8 110.4 126.6 108.0 
[289] 50.2 20.4 125.8 34.0 102.6 66.6 22.8 13.0 1.0 142.4 182.2 79.8 
[301] 64.2 25.6 85.2 74.4 31.4 13.4 45.0 15.8 99.8 172.2 1.6 152.2 
[313] 33.0 40.8 64.2 23.2 55.6 34.4 3.4 59.4 70.4 142.6 210.4 100.0 
[325] 23.0 199.6 128.0 65.6 22.8 34.2 8.0 192.0 12.4 73.2 17.6 37.0 
[337] 72.4 65.6 83.8 91.4 161.0 87.2 3.3 131.4 91.0 139.6 142.8 86.0 
[349] 93.8 67.6 147.2 9.2 94.4 25.0 15.0 72.2 0.4 58.0 117.8 32.6 
[361] 102.4 113.4 125.6 122.0 1.4 121.8 70.0 7.4 26.8 13.0 63.6 53.2 
> #Stampa i 372 dati della 2a colonna in mm di pioggia mensili > #partendo dall'anno 1956 con gennaio (i dati del data.frame > #partivano da tre mesi prima (da ottobre 1955) > > X.2[1:6] [1] NA NA NA 89.2 32.0 68.0 > #Stampa i primi sei valori della seconda colonna > #[1] 89,2 32 68 106 39,8 69,4 > X.2=as.number(X.2) Errore: non trovo la funzione "as.number" > > dataset=ts(dataset) > #considera il data.frame dataset come una serie storica > > dataset2=X.2[1:372] > #prende 1 valori da 1 a 372 del vettore X.2 e li mette > #nella variabile dataset2 > > dataset2=ts(dataset2) > #dataset2 è una serie storica > > ts.plot(dataset2) #GRAF.1 Aspetto per confermare cambio pagina... > #Stampa la serie storica dataset2; sembra assente il trend. PIOGGE_Vo_GRAF0001> yt2=dataset2; n2=length(yt2); mbt2=c() 
> for(t in 2: n2-1){mbt2[t]=(yt2[t-1]+2*yt2[t]+yt2[t+1])/4} > > ts.plot(mbt2)# GRAF.2 Aspetto per confermare cambio pagina... > #disegno il grafico di mbt2, cioè i dati originali (senza trend) > #privati anche dei random s.l. > #Penso di smussare cioè dataset2 dai randoms s.l.; nel vettore mbt2 è plausibile > #siano contenuti dati relativi a stagionalità e ciclo. PIOGGE_Vo_GRAF0002 > #Ho smussato dall'originale (dataset2 )i randoms s.l. per cui avrei ottenuto > #una serie senza i randoms (con plausibili stagionalità+ciclo; assente il trend > #iniziale come si vedeva ad occhio dal grafico di dataset2. > #Così se applico una media mobile di ordine dodici a mbt2 mi aspetto di trovare > #una serie con solo il ciclo (PRIMO CASO). > > #Potevo applicare la media mobile 12 direttamente su dataset2 (SECONDO CASO), > #smussandola dalla stagionalità e forse anche dai randoms, e poi, togliendo > #la stagionalità + randoms (questa nuova serie) da dataset2 avrei ottenuto solo > # il ciclo. > #Da confrontare questi due processi che dovrebbero condurre ambedue al ciclo. > > #Questa serie con solo il ciclo la chiamo mbt2_12. > > #PROVIAMO IL PRIMO CASO (applico una media modile 12 su mbt2: > yt=mbt2; n1=length(mbt2)-2; mbt2_12=c() > for(t in 7: n1-6){mbt2_12[t]=(yt[t-6]/2+yt[t-5]+yt[t-4]+yt[t-3]+yt[t-2]+ + yt[t-1]+yt[t]+yt[t+1]+yt[t+2]+yt[t+3]+yt[t+4]+yt[t+5]+yt[t+6]/2)/12} > > ts.plot(mbt2_12)#GRAF.3 Aspetto per confermare cambio pagina... > Disegno il grafico del ciclo Errore: unexpected symbol in "Disegno il" PIOGGE_Vo_GRAF0003 > ts.plot(mbt2) Aspetto per confermare cambio pagina... > #Disegno il grafico di mbt2 (stagionalità + ciclo) e sovrappongo mbt2_12 (ciclo): > lines(mbt2_12) #è la serie mbt2 smussata della stagionalità > #Insieme al grafico mbt2 (stagianalità+ciclo) sovrappongo il ciclo: GRAF.4 PIOGGE_Vo_GRAF0005 > #Se tolgo il ciclo da stagionalità + ciclo (mbt2), ottengo mbt2-mbt2_12, > #cioè la stagionalità. > #Sorge il problema che mbt2 e mbt2_12 debbono avere la stessa lunghezza > #per poterli sottrarre > > #FACCIO DELLE PROVE PER RENDERE I VETTORI LUNGHI UGUALE > > #Controllo mbt2 > length(mbt2) [1] 371 > # 371 = 12 > head(mbt2) [1] NA NA NA NA 55.3 68.5 > # NA 1.00 70.75 177.50 202.75 159.00 > mbt2=mbt2[2:(length(mbt2)-1)] > > #Controllo mbt2_12 > length(mbt2_12) [1] 363 > # 363 > head(mbt2_12) [1] NA NA NA NA NA NA > # NA NA NA NA NA NA > #Impongo che mbt2 e mbt2_12 abbiano la stessa lunghezza per sottrarli > mbt2_12=mbt2_12[7: (length(mbt2_12)-6)] > length(mbt2_12) [1] 351 > > #Proviamo il SECONDO CASO (applico direttamente la Mb12 su dataset2) > #per il calcolo del ciclo > > yt=dataset2; n1=length(dataset2); mbt2_12_0=c() > for(t in 7: n1-6){mbt2_12_0[t]=(yt[t-6]/2+yt[t-5]+yt[t-4]+yt[t-3]+yt[t-2]+ + yt[t-1]+yt[t]+yt[t+1]+yt[t+2]+yt[t+3]+yt[t+4]+yt[t+5]+yt[t+6]/2)/12} > > ts.plot(mbt2_12_0) Aspetto per confermare cambio pagina... > #Disegno il grafico del ciclo nel secondo modo. PIOGGE_Vo_GRAF0006> #Da confrontare i due grafici del ciclo mbt2_12_0 e mbt2_12. 
> #I due cicli praticamentte coincidono. nel secondo sembra siano rimasti più 
> #errori randoms. Potremmo fare un test statistico per controllo.

Curriculum di Piero Pistoia: PIERO PISTOIA CURRICULUM2




DA CONTINUARE

RIFLESSIONI SU UN ESEMPIO DI ANALISI STATISTICA CON R: MM DI PIOGGIA CADUTI IN CIASCUN ANNO DAL 1907 AL 1990 NELLA ZONA DI LARDERELLO (POMARANCE, Pisa) ED ALTRO; precisazioni su Regressioni lineari, Intervalli di confidenza ed altro, oltre al modo più o meno ortodosso di procedere; del dott. Piero Pistoia

PREMESSA

Post in via di sviluppo; fase iniziale di un processo che si costruisce nell’andare.

CURRICULUM DI PIERO PISTOIA

PIERO PISTOIA CURRICULUM2

Questi ‘racconti’ di Statistica applicata di fatto sono modelli che automatizzano percorsi nell’analisi dei dati.

R è un programma potente e gratuito, continuamente aggiornato in tempo reale nelle principali Università del mondo. Sono disponibili centinaia di manuali gratuiti e migliaia di packages per tutte le esigenze.

Per vedere gli outputs anche grafici,  le istruzioni (quando pronte girano!) possono essere riportate, con copia incolla, per es, prima sul Blocco Note o direttamente sulla console di R. Nel procedere potremmo anche decidere di alternare le istruzioni ai grafici (vedremo). Nota bene: prima di incollare in R, è necessario ripulire il piano di lavoro di R, premendo dal MENU MODIFICA ‘Pulisci console’ e poi dal MENU VARIE ‘Rimuovi tutti gli oggetti’

Questi primi dati di Larderello sono stati ripresi dalla “COMUNITA’  DI POMARANCE 1-1991” consegnati alla Redazione da Mauro Fanfani, allora Tecnico presso il laboratorio ENEL.

Quelli di Volterra furono forniti dal dott. Juri Bettini.

—————————————-

RIFLESSIONI SU UN ESEMPIO DI ANALISI STATISTICA CON R: MM DI PIOGGIA CADUTI IN OGNI ANNO DAL 1907 AL 1998 NELLA ZONA DI LARDERELLO (POMARANCE, Pisa). POSSIBILITA’ DI TRASFERIRE, CON POCHE MODIFICHE, IL PROCESSO SU ALTRE ANALOGHE SERIE STORICHE (Es., Volterra 1956-1986)

PRECISAZIONI SU REGRESSIONI LINEARI, MEDIE MOBILI, TESTS STATISTICI, INTERVALLI DI CONFIDENZA, FORECASTS ED ALTRO, OLTRE AL MODO PIU’ O MENO ORTODOSSO DI PROCEDERE.
dott. Piero Pistoia

 Versioni del lavoro

PIOGGE_ANN_FORECAST_20_1_7.45 pdf

PIOGGE_ANN_20_1_LAPREVISIONE_ore17.45 in odt

PIOGGE_ANN_OUTPUTS_19.10

Si evidenziano i  files richiamati e si copiano sulla console di R: iniziano automaticamente a girare; se incontrano un grafico si fermano in attesa di un tasto premuto o un click del mouse, permettendo di stampare o memorizzare il grafico stesso. Volendo possiamo anche con copia-incolla trasferire pezzi dell’articolo, per es., da un grafico all’altro. Infine è possibile ricopiare le istruzione direttamente sulla console di R.

——————————-

PARTE PRIMA CON INSERIMENTO DI OTTO PAGINE DI GRAFICI (FIG. 1-8) RIPRESI DAGLI OUTPUTS

Si ‘guarda’ dentro i dati ‘provocandoli’ con le istruzioni di R!

#ESEMPIO DI ANALISI STATISTICA CON IL LINGUAGGIO R: MILLIMETRI
#DI PIOGGIA CADUTI IN CIASCUN ANNO DAL 1907 Al 1990
#A LARDERELLO E DAL 1956 AL 1986 A VOLTERRA.
#PRECISAZIONI SU INTERVALLI DI CONFIDENZA ED ALTRA STATISTICA
#dott. Piero Pistoia
#IN QUESTO PRIMO INTERVENTO SI PROPONGONO UNA SERIE DI AZIONI
#MESSE IN ATTO PER PORRE PUNTI INTERROGATIVI AI DATI
#DA ANALIZZARE. NEI GRAFICI SI ‘LEGGONO’ LE RISPOSTE RESE DAI DATI.

#I dati originali erano a scansione mensile, ma noi inizieremo
#ad analizzare i dati annuali ottenuti sommando i 12
#valori mensili delle piogge per ogni anno.

#Se volessimo esprime un’ipotesi sulla pioggia caduta nell’anno
#in un certo intervallo di anni, quella più immediata
#è che la quantità di essa costruisca nel tempo una serie
#stazionaria, senza uniformità interne, (assenza di trend,
#fluttuazioni di circa uguale ampiezza) perché sosterremmo
#che sommare le entità mensili ridurrebbe le oscillazioni
#stagionali di tutti i livelli, mentre non è così immediato
#pensare a fenomeni atmosferici e/o astronomici che
#possano attivare oscillazioni annuali periodi o non periodiche
#(come, per es. i cicli undecennali delle macchie solari…)

library(“UsingR”); library(“TTR”);library(“tseries”);
library(“graphics”);library(“forecast”)
#Le librerie, se non esistono già nella biblioteca attuale della
#versione, devono essere prima caricate dal Menù.

par(ask=T)
par(mfrow=c(2,2))

#Da prove preliminari si evidenzia che nella serie originale di
#Larderello esistono almeno tre autliers: nel 1913
#(636 sostituito da 736.2), nel 1915 (1505.5 sostituito
#da 1145) e 1972 (539.2 sostituito da 739.2), pensando ad una
#serie più armonica internamente (gli autliers, essendo ‘accideni’
#isolati pesano meno nel forecast).

#Con pochi cambiamenti possiamo analizzare meccanicamente diverse
#serie storiche fornendo il vettore dati:
#proponiamo la serie di Larderello di 84 dato a partire dal 1907
#e quella di Volterra di 30 dati a partire dal 1956

#::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::

Piogge_annuali_1907_1990=

c(1174,915.7,1204.1,1203.6,891.5,950.6,836.2, 960.6,1145.5,
1142.2,1156.3,876.4,1045.9,968.0,750.8, 878.3,
719.5, 720.7, 1045.7,1198.1,904.2,1142.0,958.0,
914.4,952.1,1141.8,1023.1,966.6,1013.0,823.1,1138.1,
763.1,1102,1083.0,949.3,1054.2,741.2,804.9,814.8,964.6,
1236.2,829.6,1086.2,782.9,1153.8,877.8,780.4,799.4,800.7,
863.8,760.7,821.6,753,1175.6,873,857.4,1192.6,1276.4,
1117.3,1251.2,742.4,1023.8,1107.4,845.6, 769.4,739.2,
782.7,799.2,985.6,1163.4,826.0,880.4,1067.7,923.6,903.6,
856.6,855.6,1010.4,842.1,891.7,1061.8,783,739.3,1009.6)
#si chiude o si apre il vettore dati per analizzare un’altra serie

piogann907990=ts(Piogge_annuali_1907_1990) #si chiude per altro
#vettore dati dello stesso nome yt0
yt0=piogann907990 #si chiude o si apre
t=c(907:990)#si chiude o si apre
#::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
#Di seguito analizziamo la serie di Volterra di trenta dati dal 1956.
#Da prove precedenti consideriamo 2 autliers 1960 (1253)e 1984 (1156)
#interpolando con 1000 e 950 rispettivamente

#yt0=c(727,766,804,839,1000,892,866,968,957,1002,1042,727,954,1012,869,678,851,655,698,858,988,
#671,932,991,841,781,837,813,950,733,821) piogge annuali Volterra 1956-1986
#t=c(956-986) #si apre o si chiude
#yt0=ts(yt0) #si apre o si chiude

#Fig. 1 -> quattro grafici sulla stessa pagina: 1)(yt0,Time); 2)acf(yt0);
#3)hist(yt0); 4)punti dati e retta di regressione

#3)hist(yt0); 4)punti dati e retta di regressione

piogge_lard_FIG0001

ts.plot(yt0)

#Osservando il grafico dei dati ad occhio ci turba notare un certa tendenza
#delle piogge a diminuire nel tempo,
#falsificando le aspettative se il fenomeno fosse corroborato.

acf(yt0)

hist(yt0)

#Osservando il grafico dei coefficienti di auto-correlazione h,
#sembra chiaramente che all’interno della serie
#non ci siano legami di causalità avvallando le aspettative.

#Applichiamo comunque un modello di regressione lineare semplice
#ai dati per un controllo.
# INIZIO DELL’ELABORAZIONE DEI DATI

# Si prova a cercare il package Using-R
#library(“UsingR”)
#library(“TTR”)

# Fitta i dati delle piogge col tempo in anni e fornisce gli autputs

t=c(1:length(yt0))
result=simple.lm(t,yt0)
summary(result)

#FIG. 2 -> #plot (result) fitta 4 grafici per l’analisi dei residui:

piogge_lard_FIG0002

#1)Nel primo (Residual,fitted) si osserva la diffusione dei dati intorno
#alla linea y=0 e il trend che non è ovvio. 2)Nel secondo (Normal qqplot),
#i residui sono gaussiani se questo grafico segue da vicino la linea
#tratteggiata.

plot(result)
par(mfrow=c(1,1)) #da ora si plotta grafico per grafico.
fit=lm(yt0~t)
attributes(fit)
summary(fit)
fit$coefficients
fit$resid

#plot(lm(yt0~t))
pre=predict(fit)# valori sulla retta
pre # valori predetti

#FIG. 3 -> acf(pre)

piogge_lard_FIG0003

acf(pre)
# Plotta i dati e la linea di regressione
ts.plot(yt0)
abline(fit)
#Opppure:plot(t,yt0); abline(result)
yt1=yt0-pre

#FIG. 4 -> nel piano cartesiano “yt0,time” viene disegnato l’intervallo
#di confidenza al 90% per la media. Ci sono due tipi di intervallo,
quello di ‘confidenza per la media’ e quello per la ‘predizione individuale’,
#naturalmente più grande. Consideriamo questi intervalli al 90%. Nel secondo
#vengono riportati ambedue gli intervalli

piogge_lard_FIG0004

 

# INTERVALLI DI CONFIDENZA
#1
summary(fit)
#2
predict(fit,data.frame(t=sort(t)),level=.90, interval=”confidence”)
#3
ci.lwr=predict(fit,data.frame(t=sort(t)),level=.90,interval=”confidence”)[,2]
#4
points(sort(t),ci.lwr,type=”l”)
#5
ci.upr=predict(fit,data.frame(t=t),level=.90, interval=”confidence”,add=T)[,3]
#6
points(sort(t),ci.upr,type=”l”)

#FIG. 5 -> punti-dati con i due intervalli, di confidenza e di previsone, al 90%

piogge_lard_FIG0005

# Si possono usare comandi di più alto livello del package UsingR di VENABLE
fit=simple.lm(t,yt0)
summary(fit)
simple.lm(t,yt0,show.ci=T,conf.level=0.90,pred=T) #fa un po’ casino!

#FACCIO DELLE PROVE DI ‘SMUSSAMENTO’

#Elimino una parte del contenuto dei files piogann908990
# (yt0) con una media mobile di ordine 5, pesata con 1,2,3,2,1 (Media Mobile 3*3)
#con comandi di basso livello e osservo quello che accade
#(confrontare con media  mobile semplice 5)
#—————————

#SMUSSO5 pesato 1,2,3,2,1 (3*3)
m5=c()
for(i in 2:length(yt0)){
m5[i] =(yt0[i-2] +2* yt0[i-1] + 3*yt0[i] +2* yt0[i+1] + yt0[i+2])/9}
#prova a plottare m5
m5=ts(m5)
m5 # smusso5 la yt0 3*3; elimino dall’originale ciò che non è random!?

L’idea era solo sperimentare en passant, smussando con m5!

Da notare è che d’improvviso ho trovato cancellati i pesi dell’m5,                                                                           sbilanciando risultati e grafici!!

#——————————
#Secondo ‘pezzo’ per cambiare vettore dati da sottoporre al prog.
#::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::

#Attenzione! devo fare modifiche per scegliere fra i due vettori dati
#da analizzare
m5=m5[3:82] #per Volterra (31 dati) va da 3 a 29; PER LARDERELLO (84 DATI)
#DA 3 A 82
n5=length(m5)
yt=yt0[3:82] #in yt ci sono da 3 a 82 dati yt<yt0; per Volterra 3-29
#::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::

par=(mfrow=c(2,2))

yt
nyt=length(yt)
nyt
par(mfrow=c(2,2))

ts.plot(m5)
acf(m5)

yt_m5=yt-m5 #m5 e yt hanno la stessa lunghezza per cui si possono sottrarre
yt_m5 #originale meno m5_smussata

ts.plot(yt_m5)
lines(m5)

acf(yt-m5) # auto-correlazione di yt-smussamento 3*3

#Analizzo il vettore dati yt-m5 cercando la retta di regressione
t1=c(1:length(m5))
fit1=lm(yt_m5~t1)
summary(fit1)

La retta di regressione non esiste!

#FACCIO DELLE PROVE DI ‘SMUSSAMENTO’

——————————————————

# FIG. 6 -> 4 grafici sulla stesso piano: 1)graf. m5 (smusso 5 pesato 12321,3*3);
# 2)acf(m5); 3)graf. yt1=yt-m5 (yt è yt0 di lunghezza pari a m5 per sottrarre);
# 4)acf(yt1).

piogge_correzione_18_10001

OSSERVAZIONI FUORI NORMA

Da osservare l’acf della serie yt-m5 : il picco fuori range può essere considerato dovuto a un caso accidentale come 1 su 20? e le oscillazioni smorzate dell’ acf di m5 non suggeriscono niente? Le armoniche di Fourier potrebbero ‘suonare’ qualcosa? Ritengo che non siano fuori luogo (anche se azzardate e deboli, ma cariche di significati scientifici, K. Popper,  ipotesi ‘tentative’ o Tentative Theory) su oscillazioni super-annuali  (di periodo circa 5, 11..?) relative ad armoniche rilevanti (si controllino anche i grafici di Fig. 7)

_________________________________

#TENTATIVI DI SMUSSAMENTO con SMA per raccogliere informazioni su                                                        #eventuali cicli super-annuali spesso periodici da analizzare con Fourier

piog3=SMA(yt0,n=3) #opero un semplice smussamento di ordine 3
piog5=SMA(yt0,n=5) # opero ora un semplice smussamento di ordine 5
piog8=SMA(yt0,n=8)
#piog11=SMA(yt0,n=11)
piog12=SMA(yt0,n=12) #opero un semplice smussamento di ordine 12
#piog15=SMA(yt0,n=15)

#ts.plot(piog3)
#ts.plot(piog5)
#ts.plot(piog8)
#ts.plot(piog12)

#smussa l’originale yt0 con 4 ordini (5 pesato 12312; 5;8;12)
par(mfrow=c(2,2))

#4 risultati dello smussamento di yt0: m5, 5, 8, 12

#FIG. 7 -> 4 grafici di smussamento sulla stessa pagina: 1)yt0 smussato m5;
#2)yt0 smussato 5; 3)yt0 smussato 8; 4)yt0 smussato 12.

piogg_18_10005

ts.plot(yt0)
lines(m5,col=”orange”,lwd=2)
ts.plot(yt0)
lines(piog5,col=”orange”,lwd=2)
ts.plot(yt0)
lines(piog8,col=”blue”,lwd=2)
ts.plot(yt0)
lines(piog12,col=”red”,lwd=4)

#Osservando le ‘lisciature’ della curva_dati originale, sarebbe

#interessante sottoporre i vettori dati degli #smussamenti alla

#funzione del PERIODOGRAMMA (scritta dal presente autore

#e testata e riportata nei diversi posts relativi ad esempi statistici a

#cura dello scrivente: in particolare vedere “UN PARZIALE

#PERCORSO DI BASE SULL’ANALISI DI UNA SERIE STORICA

#REALE” dove tale funzione in .jpg è trasferibile e funzionante);

# per la ricerca di armoniche rilevanti.

#Non sembra così strano pensare ad influenze, per lo più

#periodiche,  del Cosmo esterno sulla nostra atmosfera, con periodi

#maggiori di un anno terrestre. Anche l’Analisi di Fourier  scritta

#col Mathematica di Wolfram, in più occasioni utilizzata dallo

#stesso autore, potrebbe essere uno strumento di controllo efficace.

#Si potrebbe anche riportare la Funzione _Periodogramma

#direttamente su questo post, per ogni possibile, anche futuro,

#utilizzo su qualsiasi vettore_dati.

par(mfrow=c(1,1))

#FIG. 8 -> 4 grafici di smussamento sullo stesso piano cartesiano (yt0,Time)

 

piogge_correzione_18_10006

#Confronto 6 curve in un unico grafico (yt0,pre,m5,piog5,piog8,piog12)
ts.plot(yt0)
lines(pre)
lines(m5,col=”orange”,lwd=2)
lines(piog5,col=4,lwd=2)
lines(piog8,col=”blue”,lwd=2)
lines(piog12,col=”red”,lwd=4)
#Riflettiamo sullo smussamento interessante 12: sembra che esso tolga cicli super-annuali
#da investigarne in qualche modo le armoniche con fourier
# LA STATISTICA DI DURBIN WATSON

library(tseries)
y=yt0
n=length(y)
result=c()
result1=c()
for(t in 2:n){
result[t]=(y[t]-y[t-1])^2
}
result=result[2:n]
a=sum(result)
for(t in 1:n)
result1[t]=y[t]
b=sum(y)
dw=a/b
dw

#Possiamo usare anche funzioni di più alto livello caricando qualche #libreria nuova

#(‘lmtest’; attenzione alla versione di R!)

#PARTE SECONDA E IL PACKAGE ‘forecast’; si aggiungono
#6 pagine di grafici ripresi dagli outputs relativi

#Ci permettiamo di prendere un momento di riflessione in questo
#scorcio iniziale della PARTE SECONDA del nostro lavoro di
#’lettura’ dei dati. La retta di regressione apparsa ad occhio sui dati,
#dopo la regressione con i minimi quadrati di fatto ‘spiega’ ben poco
#del nostro campione (R-quadro trascurabile), ma se consideriamo
#le ‘statistiche campionarie T’ relative ai suoi coefficienti e le
#inseriamo sulle ascisse delle rispettive distribuzioni di Student,
#conosciute sotto l’ipotesi nulla, ci accorgiamo che probabilmente
#quella retta ‘esiste anche nell’Universo di tutti i campioni riferiti
#alla nostra ricerca. Così,nonostante questa retta non così ovvia,
#riteniamo ipoteticamente che: “se eliminiamo dai dati originali
#yt0, i valori della retta (pre), yt1=yt0-pre, già e ‘ricondizioniamo’
#il campione aggiungendo la media dei valori (cioè la media di yt0
#o di pre), otteniamo un valore migliore dell’originale almeno per
#l’applicazione del nostro modello di forecast (Simple Exponential
#Smoothing), cioè si rispetterebbero meglio i criteri di applicazione
#(omogeneità delle fluttuazioni col tempo, assenza di correlazioni
#interne…insomma migliore Stazionarietà)”. Tale ipotesi verrebbe
#controllata a posteriori con l’analisi dei residui. Basterà attribuire i
#valori della variabile ‘campione-condizionato’ alla variabile
#originale yt0 e procedere a far girare il programma successivo di
#previsione( yt0=campione-condizionato). Tale prova la qualche                                                                                 #lettore, se vorrà.

#yt1 era il vettore dei dati originali (yt0) meno le ordinate della retta
#(pre)

Myt0=mean(yt0); Mpre=mean(pre); Myt1=mean(yt1)

Myt0

Mpre

Myt1

#Dal plot di yt1

#Possiamo notare dal plot di yt1 che esiste un livello costante
#situato in corrispondenza ad valore di zero per la media; le
#fluttuazioni sembrano rimanere grossolanamente costanti nel
#tempo. E’ possibile probabilmente considerare la serie stazionaria,
#per cui si può tentare di usare un modello addittivo e uno
#’Smoothing Esponenziale Semplice’ per il forecast, forse con
#maggiore ‘matching’ di prima. Naturalmente bisogna
#’ricondizionare’ yt1 (centrata intorno alla media zero) alla serie
#originale. Riteniamo di poterlo fare aggiungendo ai valori di yt1, la
#media dei valori predetti ovvero dei valori originali; procederemo
#poi al solito a controllare questa scelta dal risultato (rilevanza del
#forecast)

Mpre=mean(pre)

yt0_condizionato=yt1+Mpre

#plotto yt0_condizionato insieme a yt0

FIG. 9′

ts.plot(yt0)
lines(yt0_condizionato,col=”blue”,lwd=2)

#Pongo:

#yt0=yt0_condizionato #per immettere nel processo_forecast i dati
#da cui si è tolto  il trend rettilineo di dubbia rilevanza

#Noi però intanto riguardiamo il grafico originale yt0 da sottoporre,

#con A. Coghian, (A LITTLE BOOK OF R FOR TIME SERIES; release

#0.2), ad un processo di previsione.

ts.plot(yt0)

piogge_forecast0001

——————————————-

 

 

#PARTE SECONDA: LA PREVISIONE
#PREMESSA SUL ‘SIMPLE EXPONENTIAL SMOOTHING (SES)’
par(ask=T)
par(mfrow=c(1,1))

#Piogge_annuali_1907_1990=

c(1174,915.7,1204.1,1203.6,891.5,950.6,836.2, 960.6,1145.5,
1142.2,1156.3,876.4,1045.9,968.0,750.8, 878.3,
719.5, 720.7, 1045.7,1198.1,904.2,1142.0,958.0,
914.4,952.1,1141.8,1023.1,966.6,1013.0,823.1,1138.1,
763.1,1102,1083.0,949.3,1054.2,741.2,804.9,814.8,964.6,
1236.2,829.6,1086.2,782.9,1153.8,877.8,780.4,799.4,800.7,
863.8,760.7,821.6,753,1175.6,873,857.4,1192.6,1276.4,
1117.3,1251.2,742.4,1023.8,1107.4,845.6, 769.4,739.2,
782.7,799.2,985.6,1163.4,826.0,880.4,1067.7,923.6,903.6,
856.6,855.6,1010.4,842.1,891.7,1061.8,783,739.3,1009.6)
#Ridefiniamo yt0
#yt0= Piogge_annuali_1907_1990

#FIG. 9
#ts.plot(yt0)

#Le righe dei commenti associati al programma che segue devono essere
#scorciate (con cancelletti all’inizio) se vogliamo far girare il tutto dalla console di R!
#come accade nella prima parte di questo post; in caso contrario vanno riscritte le linee di programma
#in successione sulla console di R

#Il metodo del Semplice ‘Lisciare’ Esponenziale è raccomandato quando 1) i dati non
#presentano un trend ‘esplicito’ e, 2) non è coglibile una stagionalità.

#Da precisare che non si ha ‘trend esplicito’ anche quando si ha, sì, un cambiamento di livello
#o media nel corso del tempo fuori del random, ma tale evento si presenta poco regolare e
#intuibile, scarsamente razionalizzabile in un trend ‘esplicito’

#Sono almeno tre i possibili interventi nei processi di forecast del metodo SES , due
#praticamente immediati, il metodo ‘ingenuo’ , o ‘naive methods’, ‘average methods’
#(per prevedere è importante solo l’ultima osservazione che rimane costante nel tempo
#di previsione), l’altra, il metodo della ‘Media’ (i valori predetti sono tutti uguali alla media
#dei valori che servono per prevedere) e il terzo, più mediato, usa medie pesate, che tengono conto
#dei valori pesandoli in maniera diversa; quelli più lontani sempre meno di quelli più vicini
#secondo una curva esponenziale.

#Fra i modelli del terzo tipo utilizzeremo due modelli di HoltWinters con le funzioni di R HoltWinters()
#e forecast.HoltWinters(); per es., per Larderello: valori annuali 1907-1990; per Volterra valori
#annuali 1956-1986.

#RIASSUMIAMO IL PROCESSO

#Poiché riteniamo di avere una serie storica a cui si può applicare un modello addittivo con fluttuzioni
#pressochè costanti pressochè costanti e che possiede un livello che si mantiene, pensiamo, uguale
#nel tempo (valore medio) e nessuna stagionalità, possiamo usare uno dei tanti modelli SES (Simple
#Exponential Smoothing), per brevi previsioni future, cioè un mezzo per stimare i livelli previsti nei
#punti correnti col tempo.

#Se usiamo però un SES predittivo che utilizza la funzione HoltWinters() di R, sono necessari,
#per questo modello due parametri da immettere nell’argomento, beta e gamma che nella fattispecie
#dobbiamo uguagliare a FALSE (si cambia il loro valore, per es., per casi con trend e/o stagionalità),
#mentre calcola da sé il terzo parametro alfa, che controlla lo smoothing (alfa vicino a zero significa
#che diamo più peso ai dati recenti). Come outputs abbiamo il valore di alfa e tutti i valori previsti
#corrispondenti ai dati delle misurazioni originali (nella fattispecie 84) con cui è possibile costruire
#la curva di previsione all’interno dei dati.

#La funzione di HoltWinters() memorizza gli outputs in una lista di ‘oggetti’ da richiamare contenuta
#nella variabile ‘yt0previsti’ da noi scelta. Così la funzione HW() scrive i valori previsti (ed altro) in
#uno dela lista degli ‘oggetti’ (fitted, SSE…) contenuti nella variabile yt0previsti calcolata, per cui
#essi si troveranno nella variabile yt0previsti$fitted. Gli ‘errori della previsione’ si calcolano
#sottraendo dai valori predetti quelli osservati e l’oggetto SSE (Sum of Squared Errors),
#rappresenta una misura dell’accuratezza del forecast, fornito sempre dalla HW( ), che
# invece verrà richiamato da yt0previsti$SSE.

#E’ interessante notare che è necessario scegliere un valore iniziale di partenza per il calcolo
#dei livelli previsti; se questo valore non viene esplicitato il programma sceglie il primo valore
#dei dati.

#Ci sono anche dei processi per aiutarci a scegliere questo valore.

#Possiamo immettere però anche questo primo valore scelto da noi, aggiungendo il
#parametro ‘1.start’ nell’argomento della funzione HoltWinters(),

#HoltWinters(yt0, beta=FALSE, gamma=False, 1.start=X)

yt0previsti=HoltWinters(yt0,beta=FALSE, gamma= FALSE)
yt0previsti

#HoltWinters exponential smoothing without trend and without seasonal
#component.

yt0previsti$fitted

plot(yt0previsti)

FIG: 10

piogge_forecast0002

#La funzione HoltWinters() calcola i forecasts solo per i dati originali (nella
#fattispecie quelli da 1907 a 1990). Se vogliamo prevedere alcuni dati futuri
#(es., 8) è necessario scaricare il package di R ‘forecast’ , tramite library(‘forecast’),
#il quale contiene la funzione nuova ‘forecast.HoltWinters()’, che memorizzerà l’output
#nella variabile da noi scelta yt0previsti2 e aggiungeremo nel suo argomento
#’yt0previsti, h=8′
yt0previsti2 = forecast.HoltWinters(yt0previsti, h = 8)

#Con la funzione plot.forecast messa disposizione dalla libreria scaricata possiamo stampare
#le predizioni ottenute con la funzione forecast.HoltWinters():

plot.forecast (yt0previsti2)

#La funzione forecast.HoltWinters() dà la previsione per ogni anno con un intervallo di
#predizione all’80% (fascia arancionee) e al 95% (fascia gialla).

library(“forecast”)

FIG. 11

piogge_forecast0003

#La funzione HoltWinters() calcola i forecasts solo per i dati originali (nella
#fattispecie quelli da 1907 a 1990). Se vogliamo prevedere alcuni dati futuri
#(es., 8) è necessario scaricare il package di R ‘forecast’ , tramite library(‘forecast’),
#il quale contiene la funzione nuova ‘forecast.HoltWinters()’, che memorizzerà l’output
#nella variabile da noi scelta yt0previsti2 e aggiungeremo nel suo argomento
#’yt0previsti, h=8′
yt0previsti2 = forecast.HoltWinters(yt0previsti, h = 8)

#Con la funzione plot.forecast messa disposizione dalla libreria scaricata possiamo stampare
#le predizioni ottenute con la funzione forecast.HoltWinters():

plot.forecast (yt0previsti2)

#La funzione forecast.HoltWinters() dà la previsione per ogni anno con un intervallo di
#predizione all’80% (fascia arancionee) e al 95% (fascia gialla).

acf(yt0previsti2$residuals, lag.max=20)

FIG. 12

piogge_forecast0004
#Osservando l’acf dei residui può accadere che qualche ordinata superi la zona permessa
#(dove c’è assenza di autocorrelazione). Può essere che uno o due segmenti escano dai
#’significance bounds’. restando incerti se, fra 1-20 lags, ci sia qualche correlazione significativa.
#Per una conferma attiviamo il test in R di Ljung-Box-test che usa la funzione “ Box-test()”

#della libreria(stats):

#Box-test(yt0previsti2$residuals, lag=20, type=”Ljung-Box”)
#IL TEST DI LJUNG-BOX

#Se il p-value supera 0.05 di poco ci sarà una piccola certezza di qualche correlazione interna,
#che decideremo se sia dovuta al caso.

#E’ buona idea controllare anche la costanza della varianza e la normalità della distribuzione
#dei residui. Per la costanza della varianza facciamo il plot dei residui con il tempo.

plot.ts(yt0previsti2$residuals)

#Appare un grafico che presenta una successive di fluttuazioni col tempo. Può accadere a vista
#che in alcune parti del tempo l’ampiezza di esse risulti leggermente divers da quella di altre parti.
#Sarà da decidere se riteniamo che sia circa costante sostenendo l’ipotesi della uguaglianza
#della varianza col tempo

FIG: 13

piogge_forecast0005

#Per controllare infine se gli errori di previsione siano a distribuzione gaussiana con media zero,
#faremo in un istogramma (area coperta=1) degli errori con sovrapposta una curva gaussiana
#con media zero e deviazione standard pari a quella della distribuzione degli errori di previsione.
#Per far questo, di seguito viene riportata la routine di Avril Coghian di cui l’autore si sente debitore
#di questo scritto.

#Basta, dalla console di R, richiamare la funzione, scritta ad hoc, plotForecastErrors() con l’argomento
#’yt0previsti2$residuals’:

#INIZIO FUNZIONE PERSONALE
#plotForecastErrors=function(forecasterrors)

plotForecastErrors=function(forecasterrors)
{
#faccio un istogramma degli errori di previsioni

mybinsize=IQR(forecasterrors)/4
mysd=sd(forecasterrors)
mymin=min(forecasterrors)-mysd*5
mymax=max(forecasterrors)+mysd*3

#Genero una distribuzione di dati gaussiana con media zero e standard deviation
#pari a mysd

mynorm=rnorm(10000,mean=0, sd=mysd)
mymin2=min(mynorm)
mymax2=max(mynorm)

if (mymin2<mymin){mymin=mymin2}
if (mymax2>mymax){mymax=mymax2}

#Faccio un istogramma rosso degli errori di previsione con un gaussiana
#sovrapposta.

mybins=seq(mymin,mymax,mybinsize)
hist(forecasterrors, col=”red”, freq=FALSE,breaks=mybins)

#freq=F assicura che l’area sotto l’istogramma è = 1.
#genera una gaussian con media zero e standard deviation mysd

myhist = hist(mynorm,plot=FALSE, breaks=mybins)

#plotto una gaussiana blu sopra l’istogramma ddegli errori di previsione:

points(myhist$mids,myhist$density, type=”l”,col=”blue”,lwd=2)
}
#FINE FUNZIONE
plotForecastErrors(yt0previsti2$residuals )

#Possiamo usare questa funzione per costruire un istogramma con gaussiana
#degli errori di previsione per qualsiasi predizione delle
#piogge!
FIG. 14

 

 

piogge_forecast0006

e cambiando scala verticale:

piogge_correzione_18_10007

#Da controllare se la distribuzione degli errori è grossolanamente centrata

#sullo zero, se è più o meno distribuita come una gaussiana, se è ‘skewed’
#da una parte rispetto alla curva normale, ecc.

#Comunque, a mio parere, è da sottolineare che questi nostri racconti di
#statistica nel loro articolarsi, non di rado si imbattono in scelte volontarie
#ed intuitive su grandezze in gioco e questo suggerisce che non è necessario
#poi in generale un atteggiamento troppo intransigente e restrittivo nelle
#valutazioni sulla rilevanza delle ipotesi.

——————————————-

40 916.8477 916.8477