GIULIANO GHILLI, UN APPASSIONATO COLLEZIONISTA, dell’insegnante Maria Chiara Bianchi Burgassi; note del coordinatore

 

NOTE DEL COORDINATORE
Articolo, al tempo Pubblicato sulla “Comunità di Pomarance” n.1 2006, che riproponiamo su questo blog, con note del coordinatore (NDC) Piero Pistoia, come ricordo  dello studioso della Natura e di altro, Giuliano Ghilli
Foto di una ‘magica’ settaria, una pietra di origine e interpretazione complesse e piuttosto criptica, scattata da piero pistoia ad una mostra di Ghilli a Pomarance nel 2008, raccolta dallo stesso Giuliano sulle biancane e crete a calanchi alle pendici della Verna (Gr).

A ‘CACCIA’ DI SEPTARIE IN UN  POSTO DI RACCOLTA IN VAL DI CECINA (Buriano, Le Fogliare), 2009

A destra Giuliano Ghilli, con accanto l’accademico dott. P. Orlandi dell’Università di Pisa, prof. di Mineralogia e Gemmologia; all’estrema sinistra lo studioso e autore della foto Massimo Magni, di cui abbiamo pubblicato su questo blog recentemente l’articolata ricerca sulle Septarie della Val di Cecina; fra parentesi abbiamo pubblicato anche un altro punto di vista sulle Septarie scritto da una studiosa di cristalloterapia Eva Saroglia, alla ricerca di un ‘invariante culturale’ regolativo, come da nostro progetto (leggere premessa al blog); l’ultimo è un personaggio di rilievo del Gruppo Mineralogico di Cecina (Livorno) ed altro.

Ma veniamo all’articolo; in prima istanza, per leggere le tre pagine dello scritto dell’insegnante Maria Chiara Burgassi cliccare in successione sui tre links seguenti:

 

Mail0911

Mail0913

Mail0915

 

Successivamente, per facilitarne la lettura, con calma, previa ristampa dei tre files, cercheremo di inserirli, sempre in pdf, in uno solo, con una opzione dello scanner, ovvero di riportare per esteso l’articolo sul blog, scannerizzando di nuovo le tre pagine stampate, questa volta in JPG, da inserire in successione. La via della conversione dei testi da PDF a ODT modificabile, sembra poco affidabile nonostante le garanzie degli svariati convertitori e certamente non percorribile se esistono nel documento intermezzi diversi (es. foto) . Faremo dei tentativi e vedremo…

In base alla precedente nota, intanto cliccando sul link sotto, si legge l’intero articolo.

GHILLI0001

Proviamo ora a riportare i tre files in jpg, inserendoli poi in successione.

 

Annunci

PIOGGE: prove su analisi dei dati in mm/mese, LARDERELLO 51-97 a cura del dott. piero pistoia

#Questo articolo permette ad ogni lettore di accedere ai dati #sperimentali su cui rifare tentativi e prove di analisi statistica, #come facciamo noi o contemporaneamente a noi, di divertirsi #a scegliere percorsi, interpretazioni e  integrazioni alternative #su essi e correggere quei passaggi che non convincono, #cercando di risolvere i problemi da noi posti ed altri che #sorgeranno nel processo, usando il potente linguaggio R #orientato al calcolo statistico, dotato di una moltitudine di #strumenti per lo stesso obbiettivo.

#Il linguaggio R fa molto bene ciò che riesce a fare!

#Chi volesse ASSAGGI SUCCESSIVI in divenire di questo post, #clicchi su questo file in odt e con copia-incolla, lo trasferisca #sulla console di R, dopo aver cambiato l’argomento del setwd.

PIOGGE_A_LARDERELLO_1951_1957_MARZO_PARTE_PRIMA

Il seguente link invece attiva un file di premessa che, trasferito su R,  gira senza dare errore sul mio computer con setwd(I:/)

PIOGGIA1_18_2_LARDEREL_1951-97 in pdf

#PROVE INIZIALI: attenzione il post è in costruzione; ci possono essere errori!!!!

#PREMESSA

#Problema da risolvere: come far leggere il file dati da qualsiasi #lettore? I files devono essere col suffisso xls.

Larderello_pioggia_csv1 file xls

 

#1- Si deve caricare sulla console di Microsoft Excel il file #precedente (si esplicita e si carica in Excel), cliccandoci sopra.

#2 – Si deve memorizzare poi da Excel (file-salva con nome)  #con l’opzione CSV di excel (delimitato da separatore di #elenco), sul proprio Hd o su qualsiasi altro supporto di massa #di proprietà del lettore interessato, col nome “Larderello_pioggia_csv1”; verrà memorizzato col suffisso CSV.

3 -Si proceda a seguire il post.

 

____________________________________________

getwd()
setwd(“I:/”) #ogni lettore deve sostituire ad I il suo dispositivo #di massa in cui vuole memorizzare il file CSV o txt dei dati
#getwd()

par(ask=T)

#Non era così immediato far caricare ad R il nostro file dati #originale relative alle piogge cadute a Larderello dal 1951 al #1997, fornito dal CNR, gentilmente concesso dal dott. #Alessandro Bettini che gestisce la Stazione Meteo di #Castenuovo V.C. (da precisare). Da un primo controllo sembra #che un solo dato non sia disponibile (sostituito da -999). Da #prima si carica il file in Excel in 7 colonne di cui quella dei #mm/giorno di pioggia è la quarta. Imponiamo che la quarta #colonna sia numerica con “,” (?) il separatore di decimali è “.” #quello delle migliaia.
#Con procedimento riportato all’inizio salviamo 4  colonne in #.CSV #di Excel con dati delimitati da separatore di elenco.

# INTANTO: carichiamo e prepariamo il file dati da studiare:
dataset_lard_51_97=read.csv(“Larderello_pioggia_csv1.CSV”, sep=”,”)
dataset_lard_51_97=dataset_lard_51_97[1:4]
dataset_lard_51_97
attach(dataset_lard_51_97)
dataset=dataset_lard_51_97[,4]
dataset=ts(dataset)
head(dataset)
tail(dataset)

_____________________________

INTERMEZZO

#MODO ALTERNATIVO DI CARICARE,
#ATTRAVERSO EXCEL, FILES CON SUFFISSO
#.TXT CARICABILE DA R
#Partendo da un file con suffisso .xls è possibile memorizzare in #un file testo “tab delimited“(opzione di Excel 2000), invece #che in CSV di Excel; file testo che è pure visto da R. Perché #funzioni nel nostro post, va caricato col nome “Larderel_pioggia_csv1.txt

dataset_lard_51_97 = read.csv(“Larderel_pioggia_csv1.txt”, sep=”\t”)

#o togliendo il sep e sostituendo a read.csv, read.csv2, funziona #ancora

dataset_lard_51_97 = read.csv2(“Larderel_pioggia_csv1.txt)

#>dataset_lard_51_97=read.csv(“Larderel_pioggia_csv1.txt”, #sep=”\t”)
#>dataset_lard_51_97=read.csv2(“Larderel_pioggia_csv1.txt”)

FINE INTERMEZZO

______________________________

#dataset
#Time Series:
#Start = 1
#End = 17154; poichè non è divisibile per 30 – si considerano
#tutti i mesi di trenta giorni – si aggiungono 6 valori finali,
#suggeriti dai valori terminali del mese di dicembre precedente.
#I dati complessivi giornalieri della pioggia caduta in mm/giorno
#ammonta così a 17160; 1760/30 = 572 mesi di 30 giorni;
#572/12 = 47 anni circa.
#Frequency = 1

piog_mensili=matrix(dataset,ncol=30,byrow=T)
medie_piog_mensili1=rowMeans(piog_mensili)
sum_piog_mensili2=medie_piog_mensili1*30
vmesi=sum_piog_mensili2

ts.plot(vmesi) #appare un dato anomalo in vmesi[301]

fix(vmesi) #serve a controllare sul video i grandi files anche
#per correggerli.

vmesi[301]

# Si corregge facendo direttamente a mano in itinere la media #dei due valori vicini e poi su richiesta si memorizza il vettore #corretto

ts.plot(vmesi) #dovrebbe essere corretto il dato mancante

#A partire da qui si deve ancora  correggere!

____________________________________________

#E’ possibile anche costruire un vettore di 572 medie mensili #della pioggia a Larderello richiamabile (da pensare)

#Il programma R salta le righe con “#”, per cui il lavoro #presentato può essere. direttamente riportato sulla console di #R con copia incolla; il programma girerà autonomamente!

getwd()
#setwd(“X:/”) #X è il nome del dispositivo di massa su cui #memorizzare il file “Larderello_pioggia_csv1” che dovrà #avere suffisso .CSV
getwd()

#CERCANDO DI RIASSUMERE
dataset_lard_51_97=read.csv(“Larderello_pioggia_csv1.CSV”, sep=”,”) #OK, ma da controllare questo sep!
dataset_lard_51_97=dataset_lard_51_97[1:4]
#dataset_lard_51_97
attach(dataset_lard_51_97)
dataset=dataset_lard_51_97[,4]
dataset=ts(dataset)
head(dataset)
tail(dataset)

#dataset

piog_mensili=matrix(dataset,ncol=30,byrow=T)
medie_piog_mensili1=rowMeans(piog_mensili)
sum_piog_mensili2=medie_piog_mensili1*30
vmesi=sum_piog_mensili2

ts.plot(vmesi) #appare un dato anomalo in vmesi[301]

fix(vmesi) #serve a controllare sul video i grandi files anche
#per correggerli.

vmesi[301]

# Si corregga in divenire (in fieri, in corso d’opera) facendo #direttamente a man0 la media dei due vicini.

ts.plot(vmesi) #dovrebbe essere corretto il dato mancante

 

UN VELOCE EXCURSUM DA RIVEDERE ED AGGIORNARE  ATTRAVERSO I LINKS INTERNI CONTINUAMENTE AMPLIATI.

#Potrei riclassificare in una matrice vmesi per righe, a partire #da gennaio 1951, di 12 colonne, dove le somme dei valori di #ogni riga rappresenterebbero la pioggia caduta nei diversi #anni a partire da 1951 e in ogni colonna andrebbero i valori #dei mesi per ogni anno; avremo così  – in quanto abbiamo #corretto il numero dati per la presenza dei mesi diversi da 30 #giorni – circa 47 righe complete e 12 colonne corrispondenti

#ai valori dei dodici mesi di ogni anno. La media di ogni riga #rappresenterebbe la pioggia media di ogni anno (per 12 #otterrei le piogge in ogni anno) e la media di ogni colonna #sarebbe la media di pioggia caduta in tutti i 47 i gennaio, i #febbraio ecc.

#Le piogge cadute nei 47 anni, ottenute col comando #rowMeans, costituirebbero una serie storica annuale da #studiare a parte, magari inserendola nel post costruito per #questo studio già elaborato.

#Potremmo farlo ancora con la funzione

#matrix (piog_anni=matrix(vmesi, ncol=12, byrow=T), o #anche con la funzione ts (che ha come #argomenti: file,start,   #e frequency), la quale  raggruppa i dati con i valori di ogni #mese nella stessa colonna. Nella tabella apparirebbero il

#nome dei mesi ad ogni colonna e degli anni ad ogni riga; #siamo così in grado di prendere i 47 valori di ogni mese (una #volta ogni dodici) per farne la media #piog_anni.ts1=ts(vmesi, #start=1951, frequency=12).

#Questo per un controllo incrociato.
#Poi con subas =vmesi.ts1[seq(1, length(vmesi), by=12)] si #raccolgono i dati di gennaio per i 47 anni (o quanti sono) e ne #fa la media; per gli altri mesi la media si calcola con un for #integrando il comando precedente con la sostituzione di un i #ad 1.

#media_mesi=c()
#for(i in 1:12) {media_mesi[i]=mean(vmesi.ts1[seq(i, #length(vmesi), by=12)])}

#Da confrontare poi media-mesi con colMeans comando usato #con matrix.
#Seguiranno le istruzioni descritte. Si passerà poi ai filtri!

 

#Da confrontare poi media-mesi con mediacol comando usato #con matrix.
#In tal modo si viene a costruire l ‘EFFETTO STAGIONALE.

#Le dodici medie di colonna relative ai 47 mesi dello stesso #nome (gennaio, febbraio…dicembre), è conseguenza #dell’ipotesi plausibile che, all’interno di ogni anno #debba oscillare un’onda di pioggia legata alle stagioni

#(componente stagionale della serie storica).
#Osservando il grafico del vettore vmesi invece non si notano #trends e ad occhio anche le eventuali oscillazioni sembrano #mascherate, nonostante l’acf del vettore (autocorrelazione) #risulti significativo. Se esiste un’oscillazione di ordine dodici la #riassumiamo dalle medie di ogni mese nel corso dei 47 anni. #Potremmo anche spezzare i dati in 5 parti (4 da dieci anni e #una da 7) per studiarle separatamente onde controllarne le #variazioni nel tempo (come cambia la pioggia nelle varie #stagioni col tempo?)

#Si passerà poi ai filtri!

#Un modo alternativo più ortodosso di studiare una #serie storica mensile detrendizzata  infatti è quello, #suggerito sempre dall’ipotesi accennata suddetta, di calcolare #una media mobile di ordine 12 su tutti dati o sui singoli pezzi #per estrarre dai dati originali un’onda di periodo opportuno     #(stagionalità) costruendo così una seconda serie che, poi #sottratta dall’originale, ne lascia una residuale contenente #ancora trend (se esisteva prima), e qualche ciclo. Alla fine #rimarrà così solo da testare i residui.

#Seguono le istruzioni descritte.

piog_anni=matrix(vmesi, ncol=12, byrow=T)

m_anni=rowMeans(vmesi)

pioganni=m_anni*12

piog_anni.ts1=ts(vmesi, start=1951, frequency=12)

vmesi.ts1 = piog_anni.ts1

subas =vmesi.ts1[seq(1, length(vmesi), by=12)]

 

media_mesi=c()
for(i in 1:12) {media_mesi[i]=mean(vmesi.ts1[seq(i, length(vmesi), by=12)])}

ts.plot(media-mesi)

DA CONTROLLARE, CORREGGERE ESPANDERE, CONTINUARE….con calma, avendo altri post avviati da gestire.

SEGUIRE I LINKS INTERNI CONTINUAMENTE AGGIORNATI

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 CURRICULUMOK

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