www.

ingam

bunit

n.wor

dpres

s.com

Calcolo delle curve di possibilitàpluviometrica con R

www.ingambunitn.wordpress.com

9 maggio 2013

Si riportano i comandi per creare le curve di possibilità pluviometrica conR. Per qualche strano motivo lo script che ne esce da risultati uguali a quellodel professore solo per la durata di un ora.

Per prima cosa occorre pulire lo workspace, onde evitare problemi com-putazionali derivanti da precedenti elaborazioni. Lo facciamo con il comando

rm(list=ls(all=TRUE))

Si procede poi con l’impostazione della workdirectory, ovvero della car-tella dove inseriamo tutti i dati innerenti il progetto

setwd("/cartella/dove/ho/messo/i/files")

Nella workdirectory si vanno ad inserire i dati relativi alla pluviometria.Ora occorre ordinare al programma di caricarli

read.table("pluviometria.txt",header=TRUE) ->data1

La stringa header=TRUE sta ad indicare che il file condiente un header eche, per evitare problemi computazionali, esso non deve essere letto.

Vado ora a creare due vettori contenenti la colonna anni e la colonnarelativa alle precipitazioni di durata 1 ora

data1h <- data1[,2]anni <- data1[,1]

sto leggendo, rispettivamente, la prima e la seconda colonna del file pluviometria.txt.Potrebbe essere che vi siano delle righe di dati mancanti, ovvero dati

di tipologia NA. Occorre eliminarli per evitare problemi computazionali inseguito.

1

www.

ingam

bunit

n.wor

dpres

s.com

anni <-anni[!is.na(data1h)]data1h <-data1h[!is.na(data1h)]

Posso ora creare dei grafici relativi alle precipitazioni massime di durata1 ora. Per prima cosa creiamo un grafico semplice ed un istogramma

plot(anni,data1h,xlab="anno",ylab="h [mm]",main="Precipitazioni maxannuali di durata 1h",col="purple ",lwd="5",type="l")hist(data1h)

Andiamo ora a costruire la frequenza empirica di non superamento (ecdf)

ecdf(data1h) -> ecdf1hx <- seq(from=0,to=100,by=1)ecdf1h(x)plot(ecdf1h,xlab="h[mm]",ylab="Fr[H<h]",main="ecdf dati 1h")

Tale grafico può essere ben rappresentato dalle curve di Gumbel. Perprima cosa carico il pacchetto evd che mi permette di usarle

library("evd")

Per poter disegnare la curva di Gumbel bisogna prima stimarne i para-metri caratteristici a e b. Per farlo esistono tre modi:

• Metodo dei momenti

• Metodo di massima verosimiglianza

• Metodo dei minimi quadrati

Cominciamo con il metodo dei momenti. Occorre determinare media evarianza della distribuzione per poi andare a determinare i parametri dellacurva di Gumbel

mean(data1h)->m1hvar(data1h)->var1heulergamma <- 0.5772156649b1hm <- sqrt(6*var1h)/pia1hm <- m1h-b1hm*eulergamma

Ora che conosco i parametri della curva di Gumbel, posso plottarla

plot(ecdf1h,xlab="h[mm]",ylab="Fr[H<h]",main="ecdf dati 1h")lines(x,pgumbel(x,loc=a1hm,scale=b1hm),lwd=2, col="red")

2

www.

ingam

bunit

n.wor

dpres

s.com

Procedo ora al calcolo dei parametri con il metodo di massima verosimi-glianza. Per procedere devo prima caricare il pacchetto MASS

library("MASS")

Utilizzando le funzioni di questo pacchetto posso stimare i parametri

fitdistr(data1h,densfun=dgumbel,start=list(loc=a1h,scale=b1h)) ->ml1hml1h$estimate["loc"] <- a1hvml1h$estimate["scale"] <- b1hv

Ora posso andare a disegnare la curva di Gumbel

lines(x,pgumbel(x,loc=a1hv,scale=b1hv),lwd=2, col="green")

Calcoliamo infine i parametri con il metodo dei minimi quadrati

sort(data1h) -> sorteddata1hecdf1h(sorteddata1h)Y <- -(log(-log(ecdf1h(sorteddata1h))))Y <- Y[-50]X <- sorteddata[-50]lm(Y ~ X) ->l01hb1hq <- 1/l01h$coefficients[2]a1hq <- - l01h$coefficients[1]*b1hq

Il comando Y <- Y[-50] sta ad indicare che tolgo l’ultimo elemento delvettore Y (in questo caso perché è un dato che vale infinito). Se l’ultimoelemento non è in cinquantesimo occorre, ovviamente, cambiare il parametrotra parentesi quadre.

Vado a disegnare anche questa curva di Gumbel

lines(x,pgumbel(x,loc=a1hl,scale=b1hl),lwd=2, col="blue")

Ora devo scegliere quale delle tre curve di Gumbel calcolate approssimameglio la mia distribuzione: uso il test di Pearson. Per prima cosa definiscoi quantili

q <- c(0.2,0.4,0.6,0.8)

Poi comincio a stimare il valore del χ2 dei tre metodi.Metodo dei momenti

3

www.

ingam

bunit

n.wor

dpres

s.com

qgumbel(q,loc=a1hm,scale=b1hm) -> qi1hmc(0,ecdf1h(qi1hm))*length(data1h) -> no11hmc(ecdf1h(qi1hm),1)*length(data1h) -> no21hmno21hm-no11hm->no1hm0.2*length(data1h) -> deltapi1hmX1hm <- sum((no1hm-deltapi1hm)^2/deltapi1hm)

Metodo di massima verosimiglianza

qgumbel(q,loc=a1hv,scale=b1hv) -> qi1hvc(0,ecdf1h(qi1hv))*length(data1h) -> no11hvc(ecdf1h(qi1hv),1)*length(data1h) -> no21hvno21hv-no11hv->no1hv0.2*length(data1h) -> deltapi1hvX1hv <- sum((no1hv-deltapi1hv)^2/deltapi1hv)

Metodo dei minimi quadrati

qgumbel(q,loc=a1hq,scale=b1hq) -> qi1hqc(0,ecdf1h(qi1hq))*length(data1h) -> no11hqc(ecdf1h(qi1hq),1)*length(data1h) -> no21hqno21hq-no11hq->no1hq0.2*length(data1h) -> deltapi1hqX1hq <- sum((no1hq-deltapi1hq)^2/deltapi1hq)

Abbiamo, finalmente, i tre outputs che ci servono: X1hm, X1hv, X1hq. Ilminore dei tre ci dice quale dei tre metodi associato stima meglio i parametridella curva di Gumbel.

Si deve ripetere l’intera operazione per i vari intervalli di tempo perdeterminare i parametri delle curve di Gumbel relativi ai vari intervalli.

Calcolo delle curve di possibilità pluviometrica

Supponiamo di aver fatto tutti gli script per le varie durate e di aver trovatoi valori buoni dei parametri a e b per le varie ore.

a1h <- 31.59047b1h <-9.073165a3h<-39.18705b3h<-10.69436a6h<-48.10011b6h<-12.39947a12h<-60.5621

4

www.

ingam

bunit

n.wor

dpres

s.com

b12h<-18.00545a24h<-75.84853b24h<-22.43485

Tracciamo le curve di Gumbel per le varie durate

x <- seq(from=1, to=200, by=1)plot(x,pgumbel(x,a1h,b1h),type="l",col="black")lines(x,pgumbel(x,a3h,b3h),type="l",col="red")lines(x,pgumbel(x,a6h,b6h),type="l",col="blue")lines(x,pgumbel(x,a12h,b12h),type="l",col="darkgreen")lines(x,pgumbel(x,a24h,b24h),type="l",col="skyblue")

Possiamo anche disegnare la densità di probabilità legata alle curve diGumbel

plot(x,dgumbel(x,a1h,b1h),type="l",col="black")lines(x,dgumbel(x,a3h,b3h),type="l",col="red")lines(x,dgumbel(x,a6h,b6h),type="l",col="blue")lines(x,dgumbel(x,a12h,b12h),type="l",col="darkgreen")lines(x,dgumbel(x,a24h,b24h),type="l",col="skyblue")

Supponiamo di voler tracciare le curve di possibilità pluviometrica contempo di ritorno 10 anni, 50 anni e 100 anni.

durations <- c(1,3,6,12,24)quantiles <- c(1-1/10,1-1/50,1-1/100)h1 <- c(qgumbel(quantiles,a1h,b1h))h3 <- c(qgumbel(quantiles,a3h,b3h))h6 <- c(qgumbel(quantiles,a6h,b6h))h12 <- c(qgumbel(quantiles,a12h,b12h))h24 <- c(qgumbel(quantiles,a24h,b24h))

Procediamo con la curva a 10 anni

ten <- c(h1[1],h3[1],h6[1],h12[1],h24[1])plot(durations,ten,log="xy",add=T)lm(log(ten) ~ log(durations)) -> lr10lr10$coefficients[1]lr10$coefficients[2]z <- seq(from=0.99,to=50,by=0.25)plot(z,exp(lr10$coefficients[1])*z^lr10$coefficients[2],type="l",log="xy",col="dark red",ylim=c(40,300))points(durations,ten)

5

www.

ingam

bunit

n.wor

dpres

s.com

La curva a 50 anni

fifty <- c(h1[2],h3[2],h6[2],h12[2],h24[2])#lines(durations,fifty,log="xy",add=T)lm(log(fifty) ~ log(durations)) -> lr50lr50$coefficients[1]lr50$coefficients[2]z50 <- seq(from=0.99,to=50,by=0.25)lines(z,exp(lr50$coefficients[1])*z50^lr50$coefficients[2],type="l",log="xy",col="dark green",ylim=c(40,300))points(durations,fifty)

La curva a 100 anni

hundred <- c(h1[3],h3[3],h6[3],h12[3],h24[3])#plot(durations,hundred,log="xy",add=T)lm(log(hundred) ~ log(durations)) -> lr100lr100$coefficients[1]lr100$coefficients[2]z100 <- seq(from=0.99,to=50,by=0.25)lines(z100,exp(lr100$coefficients[1])*z^lr100$coefficients[2],type="l",log="xy",col="blue",ylim=c(40,300))points(durations,hundred)

6