9b Calcolo delle linee segnalatrici di possibilità pluviometrica con R
-
Upload
riccardo-rigon -
Category
Education
-
view
4.072 -
download
9
description
Transcript of 9b Calcolo delle linee segnalatrici di possibilità pluviometrica con R
C. W
hit
e, P
art
of
“ D
rop
s of
Rai
n”,
19
67
Riccardo Rigon, Matteo Dall’Amico
Finding the LSPP with R
Tuesday, March 6, 12
Le precipitazioni Estreme
Riccardo Rigon
Obbiettivi:
2
•Calcolare le curve di possibilità pluviometrica con assegnato tempo di ritorno con R
Tuesday, March 6, 12
Le precipitazioni Estreme
Riccardo Rigon
Consideriamo le precipitazioni massime annualiQueste si trovano negli annali idrologici registrate per certe durate caratteristiche:
1h, 3h, 6h,12h 24 h e rappresentano il massimo di precipitazione cumulato sulla
prefissata durata.
3
anno 1h 3h 6h 12h 24h1 1925 50.0 NA NA NA NA2 1928 35.0 47.0 50.0 50.4 67.6
......................................
......................................
46 1979 38.6 52.8 54.8 70.2 84.247 1980 28.2 42.4 71.4 97.4 107.451 1987 32.6 40.6 64.6 77.2 81.252 1988 89.2 102.0 102.0 102.0 104.2
Tuesday, March 6, 12
Le precipitazioni Estreme
Riccardo Rigon
> read.table("PluviometriaPaperopoli_2.txt",header=TRUE) -> data1> rep(c(1,3,6,12,24),each=50) -> h> c(data1[[2]],data1[[3]],data1[[4]],data1[[5]],data1[[6]]) -> hh> plot(hh ~ h,xlab="durata",ylab="Precipitazione (mm)",main="Precipitazioni Massime a Paperopoli")
4
Tuesday, March 6, 12
Le precipitazioni Estreme
Riccardo Rigon
1 3 6 12 24
50
100
150
Precipitazioni Massime a Paperopoli
durata
Pre
cip
itazio
ne (
mm
)
Mediana
>boxplot(hh ~ h,xlab="durata",ylab="Precipitazione (mm)",main="Precipitazioni Massime a Paperopoli") 5
Tuesday, March 6, 12
Le precipitazioni Estreme
Riccardo Rigon
1 ora
Precipitazion in mm
Frequenza
20 40 60 80
05
10
15
20
25
3 ore
Precipitazion in mm
Frequenza
20 40 60 80 100
05
10
15
6 ore
Precipitazion in mm
Frequenza
40 60 80 1000
510
15
6
Tuesday, March 6, 12
Le precipitazioni Estreme
Riccardo Rigon
> par(mfrow=c(1,3))
> hist(data1[[2]], breaks=8,xlab="Precipitazion in mm", ylab="Frequenza",main="Precipitazioni massime di 1 ore a Paperopoli")
> hist(data1[[3]], breaks=8,xlab="Precipitazion in mm", ylab="Frequenza",main="Precipitazioni massime di 3 ore a Paperopoli")
> hist(data1[[4]], breaks=8,xlab="Precipitazion in mm", ylab="Frequenza",main="Precipitazioni massime di 6 ore a Paperopoli")
7
Tuesday, March 6, 12
Le precipitazioni Estreme
Riccardo Rigon
Il primo problema da risolvere con l’ausilio della teoria delle probabilità e dell’analisi statistica
E’ dunque quello di determinare, per ogni durata, la corrispondenza tra
quantili (assegnati tempi di ritorno) e altezza di precipitazione.
Per ogni durata si cercherà dunque di interpolare i dati ad una distribuzione di probabilità. La famiglia di curve candidata per questo scopo è la Curva dei valori estremi di tipo I, o curva di Gumbel
b è un parametro di forma, a un parametro di posizione (in effetti la moda)
P [H < h; a, b] = e�e�h�a
b �⇥ < h <⇥
Tuesday, March 6, 12
Le precipitazioni Estreme
Riccardo Rigon
Distribuzione di Gumbel
La media della distribuzione e data da:
E[X] = b� + a
dove:
è la costante di Eulero-Mascheroni:
� � 0.57721566490153228606
Tuesday, March 6, 12
Le precipitazioni Estreme
Riccardo Rigon
Distribuzione di Gumbel
La mediana:
La varianza :
a� b log(log(2))
V ar(X) = b2 �2
6
Tuesday, March 6, 12
Le precipitazioni Estreme
Riccardo Rigon
Metodi di adattamento dei parametrirelativi alla distribuzione di Gumbel ma con validità generale
Il metodo dei momenti applicato alla curva di Gumbel consiste allora nel
porre:
o:
�b� + a = µH
b2 �2
6 = ⇤2H
�MH [1; a, b] = µH
MH [2; a, b] = ⇥2H
Tuesday, March 6, 12
Le precipitazioni Estreme
Riccardo Rigon
Il metodo della massima verosimiglianza(maximum likelihood)
relativi alla distribuzione di Gumbel ma con validità generale
Il metodo si fonda sulla valutazione della probabilità (composta) di ottenere la
serie temporale registrata:
P [{h1, · · ·, hN}; a, b]
Nella ipotesi di indipendenza delle osservazioni, tale probabilità diviene:
P [{h1, · · ·, hN}; a, b] =N�
i=1
P [hi; a, b]
Tuesday, March 6, 12
Le precipitazioni Estreme
Riccardo Rigon
La precedente probabilità si chiama anche funzione di verosimiglianza
rappresenta ed è, evidentemente una funzione dei parametri. Per
semplificare i calcoli si definisce anche la funzione detta di log-
verosimiglianza:
13
Il metodo della massima verosimiglianza(maximum likelihood)
relativi alla distribuzione di Gumbel ma con validità generale
P [{h1, · · ·, hN}; a, b] =N�
i=1
P [hi; a, b]
log(P [{h1, · · ·, hN}; a, b]) =N�
i=1
log(P [hi; a, b])
Tuesday, March 6, 12
Le precipitazioni Estreme
Riccardo Rigon
14
Il metodo della massima verosimiglianza(maximum likelihood)
relativi alla distribuzione di Gumbel ma con validità generale
Se la serie osservata è sufficientemente lunga, si ritiene che essa debba essere
tale per cui la probabilità della sua osservazione è massima. Allora, i parametri
della curva, che ne descrive la popolazione si possono ottenere da:
�⇥ log(P [{h1,···,hN};a,b])
⇥a = 0⇥ log(P [{h1,···,hN};a,b])
⇥b = 0
Che produce un sistema di due equazioni non-lineari in due incognite.
Tuesday, March 6, 12
Le precipitazioni Estreme
Riccardo Rigon
Metodo dei minimi quadrati
Consiste nel definire lo scarto quadratico tra le misure, di ECDF, e la probabilità
di non superamento:
�2(⇥) =n�
i=1
(Fi � P [H < hi; ⇥])2
e nel minimizzarlo
15
Tuesday, March 6, 12
Le precipitazioni Estreme
Riccardo Rigon
scarto quadratico
Metodo dei minimi quadrati
Consiste nel definire lo scarto quadratico tra le misure, di ECDF, e la probabilità
di non superamento:
�2(⇥) =n�
i=1
(Fi � P [H < hi; ⇥])2
e nel minimizzarlo
15
Tuesday, March 6, 12
Le precipitazioni Estreme
Riccardo Rigon
ECDFscarto quadratico
Metodo dei minimi quadrati
Consiste nel definire lo scarto quadratico tra le misure, di ECDF, e la probabilità
di non superamento:
�2(⇥) =n�
i=1
(Fi � P [H < hi; ⇥])2
e nel minimizzarlo
15
Tuesday, March 6, 12
Le precipitazioni Estreme
Riccardo Rigon
ProbabilitàECDFscarto quadratico
Metodo dei minimi quadrati
Consiste nel definire lo scarto quadratico tra le misure, di ECDF, e la probabilità
di non superamento:
�2(⇥) =n�
i=1
(Fi � P [H < hi; ⇥])2
e nel minimizzarlo
15
Tuesday, March 6, 12
Le precipitazioni Estreme
Riccardo Rigon
⇤�2(⇥j)⇤⇥j
= 0 j = 1 · · · m
Tale minimizzazione si ottiene derivando l’espressione dello scarto rispetto
agli m parametri
Ottenendo così le m equazioni in m incognite necessarie.
16
Tuesday, March 6, 12
Le precipitazioni Estreme
Riccardo Rigon
Istogramma precipitazioni massime 1 h
Precipitazioni Orarie
Freq
uenz
a
20 40 60 80
02
46
810
1214
17
Tuesday, March 6, 12
Le precipitazioni Estreme
Riccardo Rigon
Un altro comando, che sarà riusato in seguito, è quello che fornisce la frequenza di non superamento empirica (in inglese: empirical cumulative distribution function)
> ecdf(data1h) -> x> plot(x,xlab="h[mm]",ylab="P[H<h]",main="Frequenza di non superamento")
18
Tuesday, March 6, 12
Le precipitazioni Estreme
Riccardo Rigon
20 40 60 80
0.0
0.2
0.4
0.6
0.8
1.0
Frequenza di non superamento
h[mm]
P[H<h]
●
●
●●●●●
●●●●●
●●●●●
●●●●●●●
●●●●●●
●●●●●●●
●●●●
●●
●●
19
Tuesday, March 6, 12
Le precipitazioni Estreme
Riccardo Rigon
Tra le analisi preliminari che si possono effettuare, verifichiamo se i dati orari hanno distribuzione normale. Prima bisogna standardizzare i dati, ovvero ridurli a media nulla e deviazione standard unitaria.
>h.norm=(data1h - mean(data1h))/sd(data1h)
In R esiste poi il comando qqnorm che esegue l'operazione desiderata
>qqnorm(h.norm)>abline(0,1)
il secondo comando aggiunge una linea al grafico, che, nel caso in esame, corrisponde alla distribuzione gaussiana normalizzata.
20
Tuesday, March 6, 12
Le precipitazioni Estreme
Riccardo Rigon
●
●●
●
●
●
●
●
●
●
●
●
●
●
●
●
● ●●
●
●
●
●
●
●
●
●
●
●
●●
●●
●
●
●
●●
●
●
●
●
●
●
●
●
●
●
●
●
●
−2 −1 0 1 2
−10
12
34
Normal Q−Q Plot
Theoretical Quantiles
Sam
ple
Qua
ntile
s
21
Tuesday, March 6, 12
Le precipitazioni Estreme
Riccardo Rigon
Per esercizio, il lettore potrebbe ingegnarsi a verificare l'adattamento dei
dati ad una curva di Gumbel.
Veniamo ora dunque all'argomento centrale di queste slides, la stima dei
parametri della curva di Gumbel con il metodo dei momenti. Tale metodo
consiste nel fare coincidere con i momenti della popolazione, determinati
analiticamente o numericamente generando un campione (numerico)
sufficientemente numeroso di quella popolazione, con i momenti del
campione empirico. Nel caso della curva di Gumbel sono note le
espressioni analitiche dei momenti.
22
Tuesday, March 6, 12
Le precipitazioni Estreme
Riccardo Rigon
Innanzitutto, determiniamo la media e la varianza del campione:
> m1h=mean(data1h)> m1h> v1h=var(data1h)> v1h
23
Tuesday, March 6, 12
Le precipitazioni Estreme
Riccardo Rigon
Ricordando le espressioni analitiche che legano momenti e parametri della curva di Gumbel, calcoliamo dapprima il fattore di scala:
>b1.gumbel = sqrt(6*v1h)/pi
e, il fattore di localizzazione (la moda), dopo aver definito una conveniente approssimazione del numero irrazionale di Eulero
> eulergamma= 0.577216> a1.gumbel = (m1h-b1.gumbel *eulergamma)
24
Tuesday, March 6, 12
Le precipitazioni Estreme
Riccardo Rigon
I due valori cosi' ottenuti rappresentano una prima approssimazione dei
parametri (che possiamo usare, per esempio per tracciare un qqplot() non
normalizzato^1).
In possesso dei valori dei parametri, possiamo tracciare la curva di Gumbel e
sovrapporla alla curva di non superamento empirica.
25
Tuesday, March 6, 12
Le precipitazioni Estreme
Riccardo Rigon
Per disegnare la curva di Gumbel si può procedere per due strade. Una è
quella di definire una funzione in R^2. Una seconda è quella di usare una
funzione già preparata da un esperto. Il pacchetto evd (extreme value
distributions) contiene predefinita la funzione di Gumbel. Per poterlo
caricare, si usa il comando
>load(edv)
P.S. - Il pacchetto edv non e’ distribuito con R, bisogna scaricarselo dal sito e
istallarlo. A questo proposito le interfacce di R mettono a disposizione
alcune utilities, per esempio il Package manager ...
26
Tuesday, March 6, 12
Le precipitazioni Estreme
Riccardo Rigon
Per disegnare la curva di Gumbel, si scriva allora:
> z= sort(data1h) > plot(z,pgumbel(z,loc=a1.gumbel,scale=b1.gumbel), xlab="h[mm]",ylab="P[H<h]",col="red",type="l")
27
Tuesday, March 6, 12
Le precipitazioni Estreme
Riccardo Rigon
Alla quale si vuole aggiungere, eventualmente la curva cumulata, per osservare,
per confronto visivo, se il metodo dei momenti ha lavorato correttamente.
> plot(ecdf(data1h),xlab="h [mm]",main="Frequenza di non
superamento",add=T)
cruciale la parola chiave "add=T" che permette di sovrapporre i grafici.
28
Tuesday, March 6, 12
Le precipitazioni Estreme
Riccardo Rigon
Ora, oltre al metodo dei momenti, si vuole applicare il metodo dei minimi quadrati, effettuando, alla fine una regressione lineare. Il metodo, nella versione qui applicato consiste
1-nel valutare la frequenza empirica di superamento associata a cisacuno dei dati (ad ogni hi, si associa un Fi)
2-calcolare Y_i := -\log(-\log(F_i))
3- osservare che, in un piano (h,Y), i dati, se seguissero la distribuzione di Gumbel dovrebbero disporsi su di una retta
4- interpolare linearmente questa retta
R offre numerosi strumenti per effettuare, verificare e visualizzare, queste operazioni (il già citato qqplot() è uno di questi). Una possibilità è quella di fare le operazioni ad una ad una, tuttavia le funzioni già implementate in R consentono alcune scorciatoie.
29
Tuesday, March 6, 12
Le precipitazioni Estreme
Riccardo Rigon
Si può usare, per esempio, la funzione ecdf(), finora usata solo per scopi visivi.Una sua analisi appena un po' più attenta mostra che ecdf() restituisce una sorta di funzione interpolata. Infatti, posto
> ec=ecdf(data1h)
risulta che, per esempio
> ec(25.5)
restituisce la probabilità del valore h=25.5
dunque
Fi = ec(sort(data1h))
restituisce la probabilità di non superamento empirica dei dati orari.
30
Tuesday, March 6, 12
Le precipitazioni Estreme
Riccardo Rigon
L'operazione 1 è dunque risolta con ecdf(). L'operazione 2 tradotta nel linguaggio di R è:
Y = -(log(-log(Fi)))
Questa operazione comporta un problema, poichè l'ultimo data ha, relativamente ai dati presenti Fi =1, il cui doppio logaritmo è Infinito
31
Tuesday, March 6, 12
Le precipitazioni Estreme
Riccardo Rigon
Osservate come è stato eliminato l’ultimo valore che avrebbe dato in -infinito
32
Tuesday, March 6, 12
Le precipitazioni Estreme
Riccardo Rigon
Si tratta ora di interpolare linearmente i dati con y=Yb e x = data1hb. Gli strumenti per la regressione lineare in R sono molteplici. In questo caso si userà il comando
> lsfit(X,Y) -> fts
Successivamente i vari campi prodotti da lsfit possono essere interrogati nel modo seguente:
> fts$coefficients
ovvero posponendo a fts la parola chiave $coefficients. Per gli altri "attributi" prodotti dalla regressione lineare si veda l'aiuto di lsfit() con ?lsfit()
33
Tuesday, March 6, 12
Le precipitazioni Estreme
Riccardo Rigon
Il coefficiente b2.gumbel è l'inverso del coefficiente angolare della regressione:
> b2.gumbel = fts$coefficients[[2]]^-1
Il coefficiente
> a2.gumbel = -fts$coefficients[[1]]*b2.gumbel
Analogamente a quanto fatto in precedenza, anche questi due valori permettono di graficare una curva di Gumbel e di sovrapporla alla curva cumulata di probabilità empirica.
34
Tuesday, March 6, 12
Le precipitazioni Estreme
Riccardo Rigon
Passiamo ora ad affrontare il problema della determinazione dei parametri della curva di Gumbel con il metodo della massima verosimiglianza. Sotto un esempio della funzione di massima verosimiglianza per i dati orari di precipitazione.
Out[1180]= 9-„-a H50-uL - a H50 - uL + Log@aD, -„-a H35-uL - a H35 - uL + Log@aD,-„-a H35.4-uL - a H35.4 - uL + Log@aD, -„-a H67.2-uL - a H67.2 - uL + Log@aD,-„-a H25.2-uL - a H25.2 - uL + Log@aD, -„-a H35.2-uL - a H35.2 - uL + Log@aD,-„-a H48.6-uL - a H48.6 - uL + Log@aD, -„-a H36.4-uL - a H36.4 - uL + Log@aD,-„-a H47.8-uL - a H47.8 - uL + Log@aD, -„-a H39.4-uL - a H39.4 - uL + Log@aD,-„-a H21.4-uL - a H21.4 - uL + Log@aD, -„-a H33-uL - a H33 - uL + Log@aD,-„-a H42-uL - a H42 - uL + Log@aD, -„-a H39.6-uL - a H39.6 - uL + Log@aD,-„-a H28-uL - a H28 - uL + Log@aD, -„-a H42.4-uL - a H42.4 - uL + Log@aD,-„-a H21.2-uL - a H21.2 - uL + Log@aD, -„-a H21.2-uL - a H21.2 - uL + Log@aD,-„-a H19.6-uL - a H19.6 - uL + Log@aD, -„-a H41.6-uL - a H41.6 - uL + Log@aD,-„-a H51-uL - a H51 - uL + Log@aD, -„-a H32-uL - a H32 - uL + Log@aD,-„-a H27.4-uL - a H27.4 - uL + Log@aD, -„-a H35-uL - a H35 - uL + Log@aD,-„-a H21.6-uL - a H21.6 - uL + Log@aD, -„-a H36.8-uL - a H36.8 - uL + Log@aD,-„-a H54.2-uL - a H54.2 - uL + Log@aD, -„-a H39.4-uL - a H39.4 - uL + Log@aD,-„-a H30.6-uL - a H30.6 - uL + Log@aD, -„-a H30.6-uL - a H30.6 - uL + Log@aD,-„-a H33-uL - a H33 - uL + Log@aD, -„-a H32.2-uL - a H32.2 - uL + Log@aD,-„-a H38.4-uL - a H38.4 - uL + Log@aD, -„-a H33.4-uL - a H33.4 - uL + Log@aD,-„-a H31-uL - a H31 - uL + Log@aD, -„-a H37.4-uL - a H37.4 - uL + Log@aD,-„-a H36.8-uL - a H36.8 - uL + Log@aD, -„-a H39.2-uL - a H39.2 - uL + Log@aD,-„-a H29.4-uL - a H29.4 - uL + Log@aD, -„-a H40.4-uL - a H40.4 - uL + Log@aD,-„-a H37.6-uL - a H37.6 - uL + Log@aD, -„-a H30.4-uL - a H30.4 - uL + Log@aD,-„-a H44-uL - a H44 - uL + Log@aD, -„-a H38.6-uL - a H38.6 - uL + Log@aD,-„-a H28.2-uL - a H28.2 - uL + Log@aD, -„-a H61.2-uL - a H61.2 - uL + Log@aD,-„-a H23.6-uL - a H23.6 - uL + Log@aD, -„-a H20.2-uL - a H20.2 - uL + Log@aD,-„-a H32.6-uL - a H32.6 - uL + Log@aD, -„-a H89.2-uL - a H89.2 - uL + Log@aD= 35
Tuesday, March 6, 12
Le precipitazioni Estreme
Riccardo Rigon
In R l'operazione è molto semplice. Infatti R provvede per questo il comando,
fitdistr() contenuto nel pacchetto MASS (che va caricato con load(MASS). N.B. il
pacchetto evd provvede un comando alternativo: gumbelfit())
>fitdistr(data1h,densfun=dgumbel,start=list(loc=a1.gumbel,scale=b1
.gumbel)) -> mlab
Il comando precedente dice che: i dati da interpolare sono data1h, la densità di
probabilità da usare è dgumbel (fornita dal pacchetto evd). La parola chiave
"start" è seguita dai valori iniziali dai quali parte la ricerca dei parametri
ottimali. Nel caso descritto, come valori iniziali si sono scelti quelli forniti dal
metodo dei momenti (è una pratica consolidata). 36
Tuesday, March 6, 12
Le precipitazioni Estreme
Riccardo Rigon
Come risultato abbiamo 3 coppie di parametri, tutti in un verto senso ottimi. Per distinguere quali tra questi insiemi di parametri è migliore, dobbiamo usare un criterio di confronto (un test non parametrico). Useremo test di Pearson.
37
Dopo l’applicazione dei vari metodi di adattamento ...
Tuesday, March 6, 12
Le precipitazioni Estreme
Riccardo Rigon
Il test di Pearson è un test NON parametrico e consiste:
1 - Nel suddividere il campo di probabilità in k parti, per esempio uguali
38
Il Test di Pearson
Tuesday, March 6, 12
Le precipitazioni Estreme
Riccardo Rigon
Il test di Pearson è un test NON parametrico e consiste:
2 - derivarne una suddivisione del dominio
39
Il Test di Pearson
Tuesday, March 6, 12
Le precipitazioni Estreme
Riccardo Rigon
Il test di Pearson è un test NON parametrico e consiste:
3 - Contare il numero di dati in ciascun intervallo (tra i cinque della figura)
40
Il Test di Pearson
Tuesday, March 6, 12
Le precipitazioni Estreme
Riccardo Rigon
Le operazioni richieste sono le seguenti (mostrate per il caso del metodo dei momenti):
1- sudddividere l'intervallo della probabilità in un numero sufficiente di intervalli (per esempio 5: {[0,0.2),[0.2,0.4),[0.4,0.6),[0.6,0.8),[0.8,1)})
> q=c(0.2,0.4,0.6,0.8)
2 - calcolare i valori dei quantili della distribuzione di Gumbel per i valori superiori degli intervalli (0.2,0.4,0.6,0.8)
> qgumbel(q,loc=a1.gumbel,scale=b1.gumbel) -> qi
3 - Contare il numero di elementi del campione compresi negli intervalli di h: [0,h_0.2), [h_0.2,h_0.4) e cosi' via.
> c(0,ec(qi)*51) -> no1> c(ec(qi)*51,51) -> no2> no2-no1 ->no
4 - Calcolare il numero di elementi previsti nell'intervallo prescritto dalla curva interpolante. E' facile: N * \delta p, dove, nel caso in esame, \delta p =0.2 ed N è il numero totale di elementi presenti nel campione (51)
> 0.2*length(data1h) -> deltapi> deltapi
5 - Calcolare il \Chi^2
>X1=sum((no-deltapi)^2/deltapi) 41
Tuesday, March 6, 12
Le precipitazioni Estreme
Riccardo Rigon
42
Tuesday, March 6, 12
Le precipitazioni Estreme
Riccardo Rigon
Naturalmente la procedura va ripetuta per ciascuna curva interpolante.
Per la massima verosimiglianza:
> qgumbel(q,loc=a3.gumbel,scale=b3.gumbel) -> qi
> c(0,ec(qi)*51) -> no1
> c(ec(qi)*51,51) -> no2
> no2-no1
> 0.2*length(data1h) -> deltapi
> deltapi
> X3=sum((no-deltapi)^2/deltapi)
43
Tuesday, March 6, 12
Le precipitazioni Estreme
Riccardo Rigon
0 50 100 150
0.00
0.01
0.02
0.03
0.04
Precipitazione [mm]
P[h]
1h
3h
6h
12h
24h
La procedura va anche ripetuta per ogni durata, ovvero oltre che per la durata
oraria, per le 3, 6, 12 e 24 ore.
44
Dopo aver applicato Pearson
Tuesday, March 6, 12
Le precipitazioni Estreme
Riccardo Rigon
0 50 100 150
0.0
0.2
0.4
0.6
0.8
1.0
Precipitazione [mm]
P[h]
1h
3h
6h
12h
24h
45
Dopo aver applicato Pearson
Tuesday, March 6, 12
Le precipitazioni Estreme
Riccardo Rigon
> seq(from=1, to=150,by=0.1) -> x> plot(x,pgumbel(x,loc=loc.1.gumbel,scale=scale.1.gumbel),type="l",col="red",xlab="Precipitazione [mm]",ylab="P[h]")> lines(x,pgumbel(x,loc=loc.3.gumbel,scale=scale.3.gumbel),col="blue")> lines(x,pgumbel(x,loc=loc.6.gumbel,scale=scale.6.gumbel),col="dark green")> lines(x,pgumbel(x,loc=loc.12.gumbel,scale=scale.12.gumbel),col="dark red")> lines(x,pgumbel(x,loc=loc.24.gumbel,scale=scale.24.gumbel),col="dark blue")> text(30,0.6,"1h",cex=0.8)> text(40,0.55,"3h",cex=0.8)> text(48,0.50,"6h",cex=0.8)> text(60,0.45,"12h",cex=0.8)> text(70,0.40,"24h",cex=0.8)
Ecco i comandi di R per disegnare più curve di probabilità di Gumbel sullo
stesso grafico
46
Tuesday, March 6, 12
Le precipitazioni Estreme
Riccardo Rigon
> plot(x,pgumbel(x,loc=loc.1.gumbel,scale=scale.1.gumbel),type="l",col="red",xlab="Precipitazione [mm]",ylab="P[h]")
plot( ) disegna la prima curva, con un linea continua (type=”l”) di colore
rosso (col=”red”) con gli assi, i tickmarks e le legende degli assi
> seq(from=1, to=150,by=0.1) -> x
seq( ) produce un vettore di punti del dominio
47
Tuesday, March 6, 12
Le precipitazioni Estreme
Riccardo Rigon
> text(30,0.6,"1h",cex=0.8)> text(40,0.55,"3h",cex=0.8)> text(48,0.50,"6h",cex=0.8)> text(60,0.45,"12h",cex=0.8)> text(70,0.40,"24h",cex=0.8)
text(x,y,”testo”) scrive il testo tra apice nella posizione “x” e “y”.
cex=0.8 riduce ul carattere di default di 2 decimi
lines( ) disegna le curve successive che si aggiungono alla prima.
> lines(x,pgumbel(x,loc=loc.3.gumbel,scale=scale.3.gumbel),col="blue")> lines(x,pgumbel(x,loc=loc.6.gumbel,scale=scale.6.gumbel),col="dark green")> lines(x,pgumbel(x,loc=loc.12.gumbel,scale=scale.12.gumbel),col="dark red")> lines(x,pgumbel(x,loc=loc.24.gumbel,scale=scale.24.gumbel),col="dark blue")
48
Tuesday, March 6, 12
Le precipitazioni Estreme
Riccardo Rigon
0 50 100 150
0.0
0.2
0.4
0.6
0.8
1.0
Precipitazione [mm]
P[h]
1h
3h
6h
12h
24h
Tr = 10 anni
h1 h3 h6 h12 h24
49
Dopo aver applicato Pearson
Tuesday, March 6, 12
Le precipitazioni Estreme
Riccardo Rigon
>c(qgumbel(0.9,loc=loc.1.gumbel,scale=scale.1.gumbel),qgumbel(0.9,loc=loc.3.gumbel,scale=scale.3.gumbel),qgumbel(0.9,loc=loc.6.gumbel,scale=scale.6.gumbel),qgumbel(0.9,loc=loc.12.gumbel,scale=scale.12.gumbel),qgumbel(0.9,loc=loc.24.gumbel,scale=scale.24.gumbel))-> h10
Ecco come calcolare il valore del quantile desiderato, in questo caso relativo a 10 anni di tempo di ritorno (P[H <h] =0.9.
50
Tuesday, March 6, 12
Le precipitazioni Estreme
Riccardo Rigon
Questi punti vengono successivamente interpolati in campo logaritmico per
ottenere, attraverso interpolazione lineare coefficiente ed esponente delle
curve di possibilità pluviometrica.
d = c(1,3,6,12,24)> lsfit(log(d),log(h10)) ->ft10>ft10$coefficientsIntercept X 3.8647304 0.3289624
Vista il logaritmo, per ottenere i veri coefficiente, l’intercetta va esponenziata
> exp(ft10$coefficients[[1]])[1] 47.69042
51
Tuesday, March 6, 12
Le precipitazioni Estreme
Riccardo Rigon
0 5 10 15 20 25 30 35
40
60
80
100
120
140
160
180
Linee Segnalitrici di Possibilita' Pluviometrica
h [mm]
t [o
re]
52
Si ottengono infine per interpolazione le
Tuesday, March 6, 12
Le precipitazioni Estreme
Riccardo Rigon
0.5 1.0 2.0 5.0 10.0 20.0
60
80
100
120
140
160
Linee Segnalitrici di Possibilita' Pluviometrica
t [ore]
h [
mm
]
53
Si ottengono infine per interpolazione le
Tuesday, March 6, 12
Le precipitazioni Estreme
Riccardo Rigon
>plot(dd,exp(ft20$coefficients[[1]])*dd^ft20$coefficients[[2]],type="l",xlab="t [ore]",ylab="h [mm]",main="Linee Segnalitrici di Possibilita' Pluviometrica")>lines(dd,exp(ft05$coefficients[[1]])*dd^ft05$coefficients[[2]],type="l",col="blue")>lines(dd,exp(ft10$coefficients[[1]])*dd^ft10$coefficients[[2]],type="l",col="red")> points(d,h10,pch=0)> points(d,h20,pch=1)> points(d,h05,pch=2)
Infine ecco come disegnare le curve di possibilità pluviometrica
54
Tuesday, March 6, 12
Le precipitazioni Estreme
Riccardo Rigon
Nel secondo grafico è stato creato un grafico bi-logaritmico usando l’opzione di plottaggio log=”xy” indicante che entrambi gli assi sono
in scala logaritmica.
> plot(dd,exp(ft20$coefficients[[1]])*dd^ft20$coefficients[[2]], type="l",xlab="t [ore]",ylab="h [mm]",main="Linee Segnalitrici di Possibilita' Pluviometrica",log="xy")
55
Tuesday, March 6, 12
Le precipitazioni Estreme - Addendum
Riccardo Rigon - Universita’ degli studi di Trento
R:
Il �2
dchisq(x, df, ncp=0, log = FALSE)pchisq(q, df, ncp=0, lower.tail = TRUE, log.p = FALSE)qchisq(p, df, ncp=0, lower.tail = TRUE, log.p = FALSE)rchisq(n, df, ncp=0)
x, q vector of quantiles.p vector of probabilities.n number of observations. If length(n) > 1, the length is taken to be the number required.df degrees of freedom (non-negative, but can be non-integer).ncp non-centrality parameter (non-negative).log, log.plogical; if TRUE, probabilities p are given as log(p).lower.taillogical; if TRUE (default), probabilities are P[X <= x], otherwise, P[X > x].
56
Tuesday, March 6, 12