this is for holding javascript data
Maria Katherina Dal Barco deleted folder R
almost 9 years ago
Commit id: 4027493588333ba37db34f11c5b6fd6c2a5c6933
deletions | additions
diff --git a/R/.placeholder.md b/R/.placeholder.md
deleted file mode 100644
index 3cb958d..0000000
--- a/R/.placeholder.md
+++ /dev/null
...
this is a placeholder for the directory
diff --git a/R/R.R b/R/R.R
deleted file mode 100644
index be52d64..0000000
--- a/R/R.R
+++ /dev/null
...
#ANALISI PLUVIOETRICA STAZIONE METEOROLOGICA DI SPORMAGGIORE
library("zoo")
library("evd")
seq(from=as.Date("1921-01-01"),to=as.Date("2002-12-31"),by="years")->ymd
ymd
xrange<-range(ymd)
xrange
colors <-rainbow(6)
list.files ()
setwd("~/Desktop")
mio <- read.table("massimi di precipitazione.txt",header=TRUE,fill=TRUE)
j<-0
plot(xrange,c(3,140),xlab="anno",ylab="h [mm]",type="n")
for(i in 1:7){
dtq<-mio[,i]
time(dtq)
zoo(dtq) -> ptr
time(ptr) <- ymd
lines(ptr, col=colors[j])
j=j+1
legend(-19000,145,c("15 min","1 ora","3 ore","6 ore","12 ore","24 ore"),lty=c(1,1,1,1,1,1), lwd=c(2.5,2.5,2.5,2.5,2.5,2.5), col=c("violet","red","yellow","green","cyan","blue"))
}
list.files ()
setwd("~/Desktop")
data1 <- read.table("massimi di precipitazione.txt",header=TRUE)
rep(c(1,3,6,12,24,0.25),each=82) -> h
c(data1[[2]],data1[[3]],data1[[4]],data1[[5]],data1[[6]],data1[[7]]) -> hh
boxplot(hh ~ h,xlab="durata",ylab="precipitazione[mm]")
#Osservazione 15 minuti
list.files ()
setwd("~/Desktop")
read.table("massimi di precipitazione.txt",header=TRUE) -> data
data[,7] -> data15m
data[,1] -> anni
#rimuovo i valori NA
anni <- anni[!is.na(data15m)]
data15m <- data15m[!is.na(data15m)]
data15m
hist(data15m, breaks=40)
plot(anni,data15m,xlab="anno",ylab="h [mm]",main="Precipitazioni max annuali di durata 15m",col="black",lwd="1",type="l")
ecdf(data15m) -> ecdf15m
ecdf15m
plot(ecdf15m)
mean(data15m)
#metodo dei momenti
mean(data15m) -> m15m
var(data15m) -> v15m
eulergamma <- 0.5772156649
b1.gumbel.15m <- sqrt(6*v15m)/pi
a1.gumbel.15m <- (m15m - b1.gumbel.15m*eulergamma)
z15 <- sort(data15m)
plot(z15,pgumbel(z15,loc=a1.gumbel.15m, scale=b1.gumbel.15m), col="blue", type="l", xlab="h[mm]", ylab="P[H
plot(ecdf15m, add=T)
z15
#metodo dei minimi quadrati
ec15 <- ecdf(sort(data15m))
ec15
z15 <- z15[-38]
Y15 <- -log(-log(ec15(z15)))
Y15
lsfit(z15,Y15) -> fts15
fts15
b2.gumbel.15m <- fts15$coefficients[[2]]^-1
a2.gumbel.15m <- -fts15$coefficients[[1]]*b2.gumbel.15m
#metodo di massima verosimiglianza
library("MASS")
fitdistr(data15m, densfun=dgumbel, start=list(loc=a1.gumbel.15m, scale=b1.gumbel.15m)) -> mlab.15m
mlab.15m
a1.gumbel.15m
b1.gumbel.15m
a2.gumbel.15m
b2.gumbel.15m
a3.gumbel.15m <- 9.7837073
b3.gumbel.15m <- 3.0869798
#test di pearson
plot(z15, pgumbel(z15, loc=a1.gumbel.15m, scale=b1.gumbel.15m), type="l", col="blue")
plot(z15, pgumbel(z15, loc=a2.gumbel.15m, scale=b2.gumbel.15m), type="l", col="purple")
plot(z15, pgumbel(z15, loc=a3.gumbel.15m, scale=b3.gumbel.15m), type="l", col="red")
q <- c(0.2, 0.4, 0.6, 0.8)
qgumbel(q, loc=a1.gumbel.15m, scale=b1.gumbel.15m) -> qi_1.15m
qi_1.15m
data15m
c(0,ec15(qi_1.15m))*38 -> no1_1.15m
c(ec15(qi_1.15m)*38,38) -> no2_1.15m
no2_1.15m - no1_1.15m -> no_1.15m
0.2*length(data15m) -> deltapi_1.15m
X1.15m <- sum((no_1.15m - deltapi_1.15m)^2/deltapi_1.15m)
X1.15m
qgumbel(q, loc=a2.gumbel.15m, scale=b2.gumbel.15m) -> qi_2.15m
qi_2.15m
c(0,ec15(qi_2.15m)*38) -> no1_2.15m
c(ec15(qi_2.15m)*38,38) -> no2_2.15m
no2_2.15m - no1_2.15m -> no_2.15m
0.2*length(data15m) -> deltapi_2.15m
X2.15m <- sum((no_2.15m - deltapi_2.15m)^2/deltapi_2.15m)
X2.15m
qgumbel(q, loc=a3.gumbel.15m, scale=b3.gumbel.15m) -> qi_3.15m
qi_3.15m
c(0,ec15(qi_3.15m))*38 -> no1_3.15m
c(ec15(qi_3.15m)*38,38) -> no2_3.15m
no2_3.15m - no1_3.15m -> no_3.15m
0.2*length(data15m) -> deltapi_3.15m
X3.15m <- sum((no_3.15m - deltapi_3.15m)^2/deltapi_3.15m)
X3.15m
#Osservazione 1 ora
read.table("massimi di precipitazione.txt",header=TRUE) -> data
data[,2] -> data1h
data[,1] -> anni
#rimuovo i valori NA
anni <- anni[!is.na(data1h)]
data1h <- data1h[!is.na(data1h)]
data1h
hist(data1h, breaks=40)
plot(anni,data1h,xlab="anno",ylab="h [mm]",main="Precipitazioni max annuali di durata 1h",col="black",lwd="1",type="l")
ecdf(data1h) -> ecdf1h
ecdf1h
plot(ecdf1h)
mean(data1h)
#metodo dei momenti
library("evd")
mean(data1h) -> m1h
var(data1h) -> v1h
eulergamma <- 0.5772156649
b1.gumbel.1h <- sqrt(6*v1h)/pi
a1.gumbel.1h <- (m1h-b1.gumbel.1h*eulergamma)
z1 <-sort(data1h)
plot(z1,pgumbel(z1,loc=a1.gumbel.1h, scale=b1.gumbel.1h), col="blue", type="l", xlab="h[mm]", ylab="P[H
plot(ecdf1h, add=T)
z1
#metodo dei minimi quadrati
ec1 <- ecdf(sort(data1h))
ec1
z1 <- z1[-55]
Y1 <- -log(-log(ec1(z1)))
Y1
lsfit(z1,Y1) -> fts1
fts1
b2.gumbel.1h <- fts1$coefficients[[2]]^-1
a2.gumbel.1h <- -fts1$coefficients[[1]]*b2.gumbel.1h
#metodo di massima verosimiglianza
library("MASS")
fitdistr(data1h, densfun=dgumbel, start=list(loc=a1.gumbel.1h, scale=b1.gumbel.1h)) -> mlab.1h
mlab.1h
a1.gumbel.1h
b1.gumbel.1h
a2.gumbel.1h
b2.gumbel.1h
a3.gumbel.1h <- 16.6102029
b3.gumbel.1h <- 4.5843401
#test di pearson
plot(z1, pgumbel(z1, loc=a1.gumbel.1h, scale=b1.gumbel.1h), type="l", col="blue")
plot(z1, pgumbel(z1, loc=a2.gumbel.1h, scale=b2.gumbel.1h), type="l", col="purple")
plot(z1, pgumbel(z1, loc=a3.gumbel.1h, scale=b3.gumbel.1h), type="l", col="red")
q <- c(0.2, 0.4, 0.6, 0.8)
qgumbel(q, loc=a1.gumbel.1h, scale=b1.gumbel.1h) -> qi_1.1h
qi_1.1h
c(0,ec1(qi_1.1h))*55 -> no1_1.1h
c(ec1(qi_1.1h)*55,55) -> no2_1.1h
no2_1.1h - no1_1.1h -> no_1.1h
0.2*length(data1h) -> deltapi_1.1h
X1.1h <- sum((no_1.1h - deltapi_1.1h)^2/deltapi_1.1h)
X1.1h
qgumbel(q, loc=a2.gumbel.1h, scale=b2.gumbel.1h) -> qi_2.1h
qi_2.1h
c(0,ec1(qi_2.1h))*55 -> no1_2.1h
c(ec1(qi_2.1h)*55,55) -> no2_2.1h
no2_2.1h - no1_2.1h -> no_2.1h
0.2*length(data1h) -> deltapi_2.1h
X2.1h <- sum((no_2.1h - deltapi_2.1h)^2/deltapi_2.1h)
X2.1h
qgumbel(q, loc=a3.gumbel.1h, scale=b3.gumbel.1h) -> qi_3.1h
qi_3.1h
c(0,ec1(qi_3.1h))*55 -> no1_3.1h
c(ec1(qi_3.1h)*55,55) -> no2_3.1h
no2_3.1h - no1_3.1h -> no_3.1h
0.2*length(data1h) -> deltapi_3.1h
X3.1h <- sum((no_3.1h - deltapi_3.1h)^2/deltapi_3.1h)
X3.1h
#Osservazione 3 ore
read.table("massimi di precipitazione.txt",header=TRUE) -> data
data[,3] -> data3h
data[,1] -> anni
#rimuovo i valori NA
anni <- anni[!is.na(data3h)]
data3h <- data3h[!is.na(data3h)]
data3h
hist(data3h, breaks=40)
plot(anni,data3h,xlab="anno",ylab="h [mm]",main="Precipitazioni max annuali di durata 3h",col="black",lwd="1",type="l")
ecdf(data3h) -> ecdf3h
ecdf3h
plot(ecdf3h)
mean(data3h)
#metodo dei momenti
library("evd")
mean(data3h) -> m3h
var(data3h) -> v3h
eulergamma <- 0.5772156649
b1.gumbel.3h <- sqrt(6*v3h)/pi
a1.gumbel.3h <- (m3h - b1.gumbel.3h*eulergamma)
z3 <- sort(data3h)
plot(z3,pgumbel(z3,loc=a1.gumbel.3h, scale=b1.gumbel.3h), col="blue", type="l", xlab="h[mm]", ylab="P[H
plot(ecdf3h, add=T)
z3
#metodo dei minimi quadrati
ec3 <- ecdf(sort(data3h))
ec3
z3 <- z3[-55]
Y3 <- -log(-log(ec3(z3)))
Y3
lsfit(z3,Y3) -> fts3
fts3
b2.gumbel.3h <- fts3$coefficients[[2]]^-1
a2.gumbel.3h <- -fts3$coefficients[[1]]*b2.gumbel.3h
#metodo di massima verosimiglianza
library("MASS")
fitdistr(data3h, densfun=dgumbel, start=list(loc=a1.gumbel.3h, scale=b1.gumbel.3h)) -> mlab.3h
mlab.3h
a1.gumbel.3h
b1.gumbel.3h
a2.gumbel.3h
b2.gumbel.3h
a3.gumbel.3h <- 24.7465611
b3.gumbel.3h <- 6.1060564
#test di pearson
plot(z3, pgumbel(z3, loc=a1.gumbel.3h, scale=b1.gumbel.3h), type="l", col="blue")
plot(z3, pgumbel(z3, loc=a2.gumbel.3h, scale=b2.gumbel.3h), type="l", col="purple")
plot(z3, pgumbel(z3, loc=a3.gumbel.3h, scale=b3.gumbel.3h), type="l", col="red")
q <- c(0.2, 0.4, 0.6, 0.8)
qgumbel(q, loc=a1.gumbel.3h, scale=b1.gumbel.3h) -> qi_1.3h
qi_1.3h
data3h
c(0,ec3(qi_1.3h))*55 -> no1_1.3h
c(ec3(qi_1.3h)*55,55) -> no2_1.3h
no2_1.3h - no1_1.3h -> no_1.3h
0.2*length(data3h) -> deltapi_1.3h
X1.3h <- sum((no_1.3h - deltapi_1.3h)^2/deltapi_1.3h)
X1.3h
qgumbel(q, loc=a2.gumbel.3h, scale=b2.gumbel.3h) -> qi_2.3h
qi_2.3h
c(0,ec3(qi_2.3h)*55) -> no1_2.3h
c(ec3(qi_2.3h)*55,55) -> no2_2.3h
no2_2.3h - no1_2.3h -> no_2.3h
0.2*length(data3h) -> deltapi_2.3h
X2.3h <- sum((no_2.3h - deltapi_2.3h)^2/deltapi_2.3h)
X2.3h
qgumbel(q, loc=a3.gumbel.3h, scale=b3.gumbel.3h) -> qi_3.3h
qi_3.3h
c(0,ec3(qi_3.3h))*55 -> no1_3.3h
c(ec3(qi_3.3h)*55,55) -> no2_3.3h
no2_3.3h - no1_3.3h -> no_3.3h
0.2*length(data3h) -> deltapi_3.3h
X3.3h <- sum((no_3.3h - deltapi_3.3h)^2/deltapi_3.3h)
X3.3h
#Osservazione 6 ore
read.table("massimi di precipitazione.txt",header=TRUE) -> data
data[,4] -> data6h
data[,1] -> anni
#rimuovo i valori NA
anni <- anni[!is.na(data6h)]
data6h <- data6h[!is.na(data6h)]
data6h
hist(data6h, breaks=40)
plot(anni,data6h,xlab="anno",ylab="h [mm]",main="Precipitazioni max annuali di durata 6h",col="black",lwd="1",type="l")
ecdf(data6h) -> ecdf6h
ecdf6h
plot(ecdf6h)
mean(data6h)
#metodo dei momenti
library("evd")
mean(data6h) -> m6h
var(data6h) -> v6h
eulergamma <- 0.5772156649
b1.gumbel.6h <- sqrt(6*v6h)/pi
a1.gumbel.6h <- (m6h - b1.gumbel.6h*eulergamma)
z6 <- sort(data6h)
plot(z6,pgumbel(z6,loc=a1.gumbel.6h, scale=b1.gumbel.6h), col="blue", type="l", xlab="h[mm]", ylab="P[H
plot(ecdf6h, add=T)
z6
#metodo dei minimi quadrati
ec6 <- ecdf(sort(data6h))
ec6
z6 <- z6[-56]
z6 <- z6[-55]
Y6 <- -log(-log(ec6(z6)))
Y6
lsfit(z6,Y6) -> fts6
fts6
b2.gumbel.6h <- fts6$coefficients[[2]]^-1
a2.gumbel.6h <- -fts6$coefficients[[1]]*b2.gumbel.6h
#metodo di massima verosimiglianza
library("MASS")
fitdistr(data6h, densfun=dgumbel, start=list(loc=a1.gumbel.6h, scale=b1.gumbel.6h)) -> mlab.6h
mlab.6h
a1.gumbel.6h
b1.gumbel.6h
a2.gumbel.6h
b2.gumbel.6h
a3.gumbel.6h <- 34.5098162
b3.gumbel.6h <- 6.8655925
#test di pearson
plot(z6, pgumbel(z6, loc=a1.gumbel.6h, scale=b1.gumbel.6h), type="l", col="blue")
plot(z6, pgumbel(z6, loc=a2.gumbel.6h, scale=b2.gumbel.6h), type="l", col="purple")
plot(z6, pgumbel(z6, loc=a3.gumbel.6h, scale=b3.gumbel.6h), type="l", col="red")
q <- c(0.2, 0.4, 0.6, 0.8)
qgumbel(q, loc=a1.gumbel.6h, scale=b1.gumbel.6h) -> qi_1.6h
qi_1.6h
data6h
c(0,ec6(qi_1.6h))*56 -> no1_1.6h
c(ec6(qi_1.6h)*56,56) -> no2_1.6h
no2_1.6h - no1_1.6h -> no_1.6h
0.2*length(data6h) -> deltapi_1.6h
X1.6h <- sum((no_1.6h - deltapi_1.6h)^2/deltapi_1.6h)
X1.6h
qgumbel(q, loc=a2.gumbel.6h, scale=b2.gumbel.6h) -> qi_2.6h
qi_2.6h
c(0,ec6(qi_2.6h)*56) -> no1_2.6h
c(ec6(qi_2.6h)*56,56) -> no2_2.6h
no2_2.6h - no1_2.6h -> no_2.6h
0.2*length(data6h) -> deltapi_2.6h
X2.6h <- sum((no_2.6h - deltapi_2.6h)^2/deltapi_2.6h)
X2.6h
qgumbel(q, loc=a3.gumbel.6h, scale=b3.gumbel.6h) -> qi_3.6h
qi_3.6h
c(0,ec6(qi_3.6h))*56 -> no1_3.6h
c(ec6(qi_3.6h)*56,56) -> no2_3.6h
no2_3.6h - no1_3.6h -> no_3.6h
0.2*length(data6h) -> deltapi_3.6h
X3.6h <- sum((no_3.6h - deltapi_3.6h)^2/deltapi_3.6h)
X3.6h
#Osservazione 12 ore
read.table("massimi di precipitazione.txt",header=TRUE) -> data
data[,5] -> data12h
data[,1] -> anni
#rimuovo i valori NA
anni <- anni[!is.na(data12h)]
data12h <- data12h[!is.na(data12h)]
data12h
hist(data12h, breaks=40)
plot(anni,data12h,xlab="anno",ylab="h [mm]",main="Precipitazioni max annuali di durata 12h",col="black",lwd="1",type="l")
ecdf(data12h) -> ecdf12h
ecdf12h
plot(ecdf12h)
mean(data12h)
#metodo dei momenti
library("evd")
mean(data12h) -> m12h
var(data12h) -> v12h
eulergamma <- 0.5772156649
b1.gumbel.12h <- sqrt(6*v12h)/pi
a1.gumbel.12h <- (m12h - b1.gumbel.12h*eulergamma)
z12 <- sort(data12h)
plot(z12,pgumbel(z12,loc=a1.gumbel.12h, scale=b1.gumbel.12h), col="blue", type="l", xlab="h[mm]", ylab="P[H
plot(ecdf12h, add=T)
z12
#metodo dei minimi quadrati
ec12 <- ecdf(sort(data12h))
ec12
z12 <- z12[-56]
Y12 <- -log(-log(ec12(z12)))
Y12
lsfit(z12,Y12) -> fts12
fts12
b2.gumbel.12h <- fts12$coefficients[[2]]^-1
a2.gumbel.12h <- -fts12$coefficients[[1]]*b2.gumbel.12h
#metodo di massima verosimiglianza
library("MASS")
fitdistr(data12h, densfun=dgumbel, start=list(loc=a1.gumbel.12h, scale=b1.gumbel.12h)) -> mlab.12h
mlab.12h
a1.gumbel.12h
b1.gumbel.12h
a2.gumbel.12h
b2.gumbel.12h
a3.gumbel.12h <- 48.039083
b3.gumbel.12h <- 11.613895
#test di pearson
plot(z12, pgumbel(z12, loc=a1.gumbel.12h, scale=b1.gumbel.12h), type="l", col="blue")
plot(z12, pgumbel(z12, loc=a2.gumbel.12h, scale=b2.gumbel.12h), type="l", col="purple")
plot(z12, pgumbel(z12, loc=a3.gumbel.12h, scale=b3.gumbel.12h), type="l", col="red")
q <- c(0.2, 0.4, 0.6, 0.8)
qgumbel(q, loc=a1.gumbel.12h, scale=b1.gumbel.12h) -> qi_1.12h
qi_1.12h
data12h
c(0,ec12(qi_1.12h))*56 -> no1_1.12h
c(ec12(qi_1.12h)*56,56) -> no2_1.12h
no2_1.12h - no1_1.12h -> no_1.12h
0.2*length(data12h) -> deltapi_1.12h
X1.12h <- sum((no_1.12h - deltapi_1.12h)^2/deltapi_1.12h)
X1.12h
qgumbel(q, loc=a2.gumbel.12h, scale=b2.gumbel.12h) -> qi_2.12h
qi_2.12h
c(0,ec12(qi_2.12h)*56) -> no1_2.12h
c(ec12(qi_2.12h)*56,56) -> no2_2.12h
no2_2.12h - no1_2.12h -> no_2.12h
0.2*length(data12h) -> deltapi_2.12h
X2.12h <- sum((no_2.12h - deltapi_2.12h)^2/deltapi_2.12h)
X2.12h
qgumbel(q, loc=a3.gumbel.12h, scale=b3.gumbel.12h) -> qi_3.12h
qi_3.12h
c(0,ec12(qi_3.12h))*56 -> no1_3.12h
c(ec12(qi_3.12h)*56,56) -> no2_3.12h
no2_3.12h - no1_3.12h -> no_3.12h
0.2*length(data12h) -> deltapi_3.12h
X3.12h <- sum((no_3.12h - deltapi_3.12h)^2/deltapi_3.12h)
X3.12h
#Osservazione 24 ore
read.table("massimi di precipitazione.txt",header=TRUE) -> data
data[,6] -> data24h
data[,1] -> anni
#rimuovo i valori NA
anni <- anni[!is.na(data24h)]
data24h <- data24h[!is.na(data24h)]
data24h
hist(data24h, breaks=40)
plot(anni,data24h,xlab="anno",ylab="h [mm]",main="Precipitazioni max annuali di durata 24h",col="black",lwd="1",type="l")
ecdf(data24h) -> ecdf24h
ecdf24h
plot(ecdf24h)
mean(data24h)
#metodo dei momenti
library("evd")
mean(data24h) -> m24h
var(data24h) -> v24h
eulergamma <- 0.5772156649
b1.gumbel.24h <- sqrt(6*v24h)/pi
a1.gumbel.24h <- (m24h - b1.gumbel.24h*eulergamma)
z24 <- sort(data24h)
plot(z24,pgumbel(z24,loc=a1.gumbel.24h, scale=b1.gumbel.24h), col="blue", type="l", xlab="h[mm]", ylab="P[H
plot(ecdf24h, add=T)
z24
#metodo dei minimi quadrati
ec24 <- ecdf(sort(data24h))
ec24
z24 <- z24[-56]
Y24 <- -log(-log(ec24(z24)))
Y24
lsfit(z24,Y24) -> fts24
fts24
b2.gumbel.24h <- fts24$coefficients[[2]]^-1
a2.gumbel.24h <- -fts24$coefficients[[1]]*b2.gumbel.24h
#metodo di massima verosimiglianza
library("MASS")
fitdistr(data24h, densfun=dgumbel, start=list(loc=a1.gumbel.24h, scale=b1.gumbel.24h)) -> mlab.24h
mlab.24h
a1.gumbel.24h
b1.gumbel.24h
a2.gumbel.24h
b2.gumbel.24h
a3.gumbel.24h <- 65.346178
b3.gumbel.24h <- 18.501232
#test di pearson
plot(z24, pgumbel(z24, loc=a1.gumbel.24h, scale=b1.gumbel.24h), type="l", col="blue")
plot(z24, pgumbel(z24, loc=a2.gumbel.24h, scale=b2.gumbel.24h), type="l", col="purple")
plot(z24, pgumbel(z24, loc=a3.gumbel.24h, scale=b3.gumbel.24h), type="l", col="red")
q <- c(0.2, 0.4, 0.6, 0.8)
qgumbel(q, loc=a1.gumbel.24h, scale=b1.gumbel.24h) -> qi_1.24h
qi_1.24h
data24h
c(0,ec24(qi_1.24h))*56 -> no1_1.24h
c(ec24(qi_1.24h)*56,56) -> no2_1.24h
no2_1.24h - no1_1.24h -> no_1.24h
0.2*length(data24h) -> deltapi_1.24h
X1.24h <- sum((no_1.24h - deltapi_1.24h)^2/deltapi_1.24h)
X1.24h
qgumbel(q, loc=a2.gumbel.24h, scale=b2.gumbel.24h) -> qi_2.24h
qi_2.24h
c(0,ec24(qi_2.24h)*56) -> no1_2.24h
c(ec24(qi_2.24h)*56,56) -> no2_2.24h
no2_2.24h - no1_2.24h -> no_2.24h
0.2*length(data24h) -> deltapi_2.24h
X2.24h <- sum((no_2.24h - deltapi_2.24h)^2/deltapi_2.24h)
X2.24h
qgumbel(q, loc=a3.gumbel.24h, scale=b3.gumbel.24h) -> qi_3.24h
qi_3.24h
c(0,ec24(qi_3.24h))*56 -> no1_3.24h
c(ec24(qi_3.24h)*56,56) -> no2_3.24h
no2_3.24h - no1_3.24h -> no_3.24h
0.2*length(data24h) -> deltapi_3.24h
X3.24h <- sum((no_3.24h - deltapi_3.24h)^2/deltapi_3.24h)
X3.24h
#15minuti
plot(z15, pgumbel(z15, loc=a1.gumbel.15m, scale=b1.gumbel.15m), col ="blue", type="l" , xlab="h[mm]", ylab="P[H
lines(z15, pgumbel(z15, loc=a2.gumbel.15m, scale=b2.gumbel.15m), col ="red", type="l" , xlab="h[mm]", ylab="P[H
lines(z15, pgumbel(z15, loc=a3.gumbel.15m, scale=b3.gumbel.15m), col ="green", type="l" , xlab="h[mm]", ylab="P[H
plot(ecdf15m, add=T)
legend(10.25,0.25,c("momenti","minimi quadrati","massima verosimiglianza "),lty=c(1,1,1), lwd=c(2.5,2.5,2.5), col=c("blue","red","green"))
#1h
plot(z1, pgumbel(z1, loc=a1.gumbel.1h, scale=b1.gumbel.1h), col ="blue", type="l" , xlab="h[mm]", ylab="P[H
lines(z1, pgumbel(z1, loc=a2.gumbel.1h, scale=b2.gumbel.1h), col ="red", type="l" , xlab="h[mm]", ylab="P[H
lines(z1, pgumbel(z1, loc=a3.gumbel.1h, scale=b3.gumbel.1h), col ="green", type="l" , xlab="h[mm]", ylab="P[H
plot(ecdf1h, add=T)
legend(19,0.25,c("momenti","minimi quadrati","massima verosimiglianza "),lty=c(1,1,1), lwd=c(2.5,2.5,2.5), col=c("blue","red","green"))
#3h
plot(z3, pgumbel(z3, loc=a1.gumbel.3h, scale=b1.gumbel.3h), col ="blue", type="l" , xlab="h[mm]", ylab="P[H
lines(z3, pgumbel(z3, loc=a2.gumbel.3h, scale=b2.gumbel.3h), col ="red", type="l" , xlab="h[mm]", ylab="P[H
lines(z3, pgumbel(z3, loc=a3.gumbel.3h, scale=b3.gumbel.3h), col ="green", type="l" , xlab="h[mm]", ylab="P[H
plot(ecdf3h, add=T)
legend(26,0.25,c("momenti","minimi quadrati","massima verosimiglianza "),lty=c(1,1,1), lwd=c(2.5,2.5,2.5), col=c("blue","red","green"))
#6h
plot(z6, pgumbel(z6, loc=a1.gumbel.6h, scale=b1.gumbel.6h), col ="blue", type="l" , xlab="h[mm]", ylab="P[H
lines(z6, pgumbel(z6, loc=a2.gumbel.6h, scale=b2.gumbel.6h), col ="red", type="l" , xlab="h[mm]", ylab="P[H
lines(z6, pgumbel(z6, loc=a3.gumbel.6h, scale=b3.gumbel.6h), col ="green", type="l" , xlab="h[mm]", ylab="P[H
plot(ecdf6h, add=T)
legend(35,0.25,c("momenti","minimi quadrati","massima verosimiglianza "),lty=c(1,1,1), lwd=c(2.5,2.5,2.5), col=c("blue","red","green"))
#12h
plot(z12, pgumbel(z12, loc=a1.gumbel.12h, scale=b1.gumbel.12h), col ="blue", type="l" , xlab="h[mm]", ylab="P[H
lines(z12, pgumbel(z12, loc=a2.gumbel.12h, scale=b2.gumbel.12h), col ="red", type="l" , xlab="h[mm]", ylab="P[H
lines(z12, pgumbel(z12, loc=a3.gumbel.12h, scale=b3.gumbel.12h), col ="green", type="l" , xlab="h[mm]", ylab="P[H
plot(ecdf12h, add=T)
legend(54,0.25,c("momenti","minimi quadrati","massima verosimiglianza "),lty=c(1,1,1), lwd=c(2.5,2.5,2.5), col=c("blue","red","green"))
#24h
plot(z24, pgumbel(z24, loc=a1.gumbel.24h, scale=b1.gumbel.24h), col ="blue", type="l" , xlab="h[mm]", ylab="P[H
lines(z24, pgumbel(z24, loc=a2.gumbel.24h, scale=b2.gumbel.24h), col ="red", type="l" , xlab="h[mm]", ylab="P[H
lines(z24, pgumbel(z24, loc=a3.gumbel.24h, scale=b3.gumbel.24h), col ="green", type="l" , xlab="h[mm]", ylab="P[H
plot(ecdf24h, add=T)
legend(76,0.25,c("momenti","minimi quadrati","massima verosimiglianza "),lty=c(1,1,1), lwd=c(2.5,2.5,2.5), col=c("blue","red","green"))
#Analisi complessiva
#Si scelgono le coppie di valori (a,b) corrispondenti a x_i.nh minore
a.15m <- 9.84838642862797
b.15m <- 2.86134571842947
a.1h <- 16.6102029
b.1h <- 4.5843401
a.3h <- 24.7465611
b.3h <- 6.1060564
a.6h <- 34.5756040136424
b.6h <- 6.58846882748384
a.12h <- 47.5046985231325
b.12h <- 11.3393627812782
a.24h <- 64.8047813576029
b.24h <- 18.3000502827714
#Funzione di probabilità
seq(from=1, to=200, by=0.1) -> x
plot(x, pgumbel(x, loc=a.15m, scale=b.15m), type="l", xlab="h[mm]", ylab="P[h]", col="purple")
lines(x, pgumbel(x, loc=a.1h, scale=b.1h), col= "blue")
lines(x, pgumbel(x, loc=a.3h, scale=b.3h), col="green")
lines(x, pgumbel(x, loc=a.6h, scale=b.6h), col="yellow")
lines(x, pgumbel(x, loc=a.12h, scale=b.12h), col="orange")
lines(x, pgumbel(x, loc=a.24h, scale=b.24h), col="red")
text(10,0.97,"15min", cex=0.8)
text(20,0.81,"1h",cex=0.8)
text(30,0.7,"3h",cex=0.8)
text(40,0.55,"6h",cex=0.8)
text(55,0.43,"12h",cex=0.8)
text(70,0.35,"24h",cex=0.8)
legend(150,0.37,c("15 minuti ","1 ora","3 ore","6 ore", "12 ore", "24 ore"),lty=c(1,1,1,1,1,1), lwd=c(2.5,2.5,2.5,2.5,2.5,2.5), col=c("purple","blue","green","yellow","orange","red"))
#Densità di probabilità delle curve di Gumbel
plot(x,dgumbel(x,a.15m,b.15m),type="l",col="purple")
lines(x,dgumbel(x,a.1h,b.1h), type="l", col="blue")
lines(x,dgumbel(x,a.3h,b.3h),type="l",col="green")
lines(x,dgumbel(x,a.6h,b.6h),type="l",col="yellow")
lines(x,dgumbel(x,a.12h,b.12h),type="l",col="orange")
lines(x,dgumbel(x,a.24h,b.24h),type="l",col="red")
text(20,0.12,"15min", cex=0.8)
text(20,0.083,"1h",cex=0.8)
text(30,0.062,"3h",cex=0.8)
text(42,0.055,"6h",cex=0.8)
text(50,0.035,"12h",cex=0.8)
text(70,0.023,"24h",cex=0.8)
legend(140,0.13,c("15 minuti ","1 ora","3 ore","6 ore", "12 ore", "24 ore"),lty=c(1,1,1,1,1,1), lwd=c(2.5,2.5,2.5,2.5,2.5,2.5), col=c("purple","blue","green","yellow","orange","red"))
#Curve di possibilità pluviometrica con tempo di ritorno 5,10,25,50,100 anni
p <- function(x){1-1/x}
p(5)
length(x)
rep(0.8,1991) -> y
length(y)
lines(x,y, col="black")
#5anni
c(qgumbel(0.8,loc=a.15m,scale=b.15m),
qgumbel(0.8,loc=a.1h,scale=b.1h),
qgumbel(0.8,loc=a.3h,scale=b.3h),
qgumbel(0.8,loc=a.6h,scale=b.6h),
qgumbel(0.8,loc=a.12h,scale=b.12h),
qgumbel(0.8,loc=a.24h,scale=b.24h)) -> h5
d <- c(0.25,1,3,6,12,24)
lsfit(log(d),log(h5)) -> ft5
ft5$coefficients[[2]]->n5
n5
exp(ft5$coefficients[[1]])->atp5
atp5
p <- function(x){1-1/x}
p(10)
length(x)
#10anni
c(qgumbel(0.9,loc=a.15m,scale=b.15m),
qgumbel(0.9,loc=a.1h,scale=b.1h),
qgumbel(0.9,loc=a.3h,scale=b.3h),
qgumbel(0.9,loc=a.6h,scale=b.6h),
qgumbel(0.9,loc=a.12h,scale=b.12h),
qgumbel(0.9,loc=a.24h,scale=b.24h)) -> h10
d <- c(0.25,1,3,6,12,24)
lsfit(log(d),log(h10)) -> ft10
ft10$coefficients[[2]]->n10
n10
exp(ft10$coefficients[[1]])->atp10
atp10
p <- function(x){1-1/x}
p(25)
length(x)
#25anni
c(qgumbel(0.96,loc=a.15m,scale=b.15m),
qgumbel(0.96,loc=a.1h,scale=b.1h),
qgumbel(0.96,loc=a.3h,scale=b.3h),
qgumbel(0.96,loc=a.6h,scale=b.6h),
qgumbel(0.96,loc=a.12h,scale=b.12h),
qgumbel(0.96,loc=a.24h,scale=b.24h)) -> h25
d <- c(0.25,1,3,6,12,24)
lsfit(log(d),log(h25)) -> ft25
ft25$coefficients[[2]]->n25
n25
exp(ft25$coefficients[[1]])->atp25
atp25
p <- function(x){1-1/x}
p(50)
length(x)
#50anni
c(qgumbel(0.98,loc=a.15m,scale=b.15m),
qgumbel(0.98,loc=a.1h,scale=b.1h),
qgumbel(0.98,loc=a.3h,scale=b.3h),
qgumbel(0.98,loc=a.6h,scale=b.6h),
qgumbel(0.98,loc=a.12h,scale=b.12h),
qgumbel(0.98,loc=a.24h,scale=b.24h)) -> h50
d <- c(0.25,1,3,6,12,24)
lsfit(log(d),log(h50)) -> ft50
ft50$coefficients[[2]]->n50
n50
exp(ft50$coefficients[[1]])->atp50
atp50
p <- function(x){1-1/x}
p(100)
length(x)
#100anni
c(qgumbel(0.99,loc=a.15m,scale=b.15m),
qgumbel(0.99,loc=a.1h,scale=b.1h),
qgumbel(0.99,loc=a.3h,scale=b.3h),
qgumbel(0.99,loc=a.6h,scale=b.6h),
qgumbel(0.99,loc=a.12h,scale=b.12h),
qgumbel(0.99,loc=a.24h,scale=b.24h)) -> h100
d <- c(0.25,1,3,6,12,24)
lsfit(log(d),log(h100)) -> ft100
ft100$coefficients[[2]]->n100
n100
exp(ft100$coefficients[[1]])->atp100
atp100
#INSIEME piano logaritmico
plot(d,exp(ft5$coefficients[[1]])*d^ft5$coefficients[[2]],type="l", col="purple", xlab="durata [h]",ylab="h [mm]",main="LineeSegnalitrici di Possibilità Pluviometrica")
lines(d,exp(ft10$coefficients[[1]])*d^ft10$coefficients[[2]],type="l",col="blue")
lines(d,exp(ft25$coefficients[[1]])*d^ft25$coefficients[[2]],type="l",col="green")
lines(d,exp(ft50$coefficients[[1]])*d^ft50$coefficients[[2]],type="l",col="orange")
lines(d,exp(ft100$coefficients[[1]])*d^ft100$coefficients[[2]],type="l",col="red")
legend(17,40,c("5 anni","10 anni","25 anni","50 anni","100 anni "),lty=c(1,1,1,1,1), lwd=c(2.5,2.5,2.5,2.5,2.5), col=c("purple","blue","green","orange","red"))
points(d,h5,pch=0)
points(d,h10,pch=1)
points(d,h25,pch=2)
points(d,h50,pch=3)
points(d,h100,pch=4)
#INSIEME piano normale
plot(d,exp(ft5$coefficients[[1]])*d^ft5$coefficients[[2]], type="l", col="purple", xlab="durata [h]", ylab="h [mm]",main="Linee Segnalatrici di Possibilita' Pluviometrica", log="xy")
lines(d,exp(ft10$coefficients[[1]])*d^ft10$coefficients[[2]],type="l",col="blue", log="xy")
lines(d,exp(ft25$coefficients[[1]])*d^ft25$coefficients[[2]],type="l",col="green", log="xy")
lines(d,exp(ft50$coefficients[[1]])*d^ft50$coefficients[[2]],type="l",col="orange", log="xy")
lines(d,exp(ft100$coefficients[[1]])*d^ft100$coefficients[[2]],type="l",col="red", log="xy")
legend(6,27,c("5 anni","10 anni","25 anni","50 anni","100 anni "),lty=c(1,1,1,1,1), lwd=c(2.5,2.5,2.5,2.5,2.5), col=c("purple","blue","green","orange","red"))
points(d,h5,pch=0)
points(d,h10,pch=1)
points(d,h25,pch=2)
points(d,h50,pch=3)
points(d,h100,pch=4)
diff --git a/R/R.txt b/R/R.txt
deleted file mode 100644
index be52d64..0000000
--- a/R/R.txt
+++ /dev/null
...
#ANALISI PLUVIOETRICA STAZIONE METEOROLOGICA DI SPORMAGGIORE
library("zoo")
library("evd")
seq(from=as.Date("1921-01-01"),to=as.Date("2002-12-31"),by="years")->ymd
ymd
xrange<-range(ymd)
xrange
colors <-rainbow(6)
list.files ()
setwd("~/Desktop")
mio <- read.table("massimi di precipitazione.txt",header=TRUE,fill=TRUE)
j<-0
plot(xrange,c(3,140),xlab="anno",ylab="h [mm]",type="n")
for(i in 1:7){
dtq<-mio[,i]
time(dtq)
zoo(dtq) -> ptr
time(ptr) <- ymd
lines(ptr, col=colors[j])
j=j+1
legend(-19000,145,c("15 min","1 ora","3 ore","6 ore","12 ore","24 ore"),lty=c(1,1,1,1,1,1), lwd=c(2.5,2.5,2.5,2.5,2.5,2.5), col=c("violet","red","yellow","green","cyan","blue"))
}
list.files ()
setwd("~/Desktop")
data1 <- read.table("massimi di precipitazione.txt",header=TRUE)
rep(c(1,3,6,12,24,0.25),each=82) -> h
c(data1[[2]],data1[[3]],data1[[4]],data1[[5]],data1[[6]],data1[[7]]) -> hh
boxplot(hh ~ h,xlab="durata",ylab="precipitazione[mm]")
#Osservazione 15 minuti
list.files ()
setwd("~/Desktop")
read.table("massimi di precipitazione.txt",header=TRUE) -> data
data[,7] -> data15m
data[,1] -> anni
#rimuovo i valori NA
anni <- anni[!is.na(data15m)]
data15m <- data15m[!is.na(data15m)]
data15m
hist(data15m, breaks=40)
plot(anni,data15m,xlab="anno",ylab="h [mm]",main="Precipitazioni max annuali di durata 15m",col="black",lwd="1",type="l")
ecdf(data15m) -> ecdf15m
ecdf15m
plot(ecdf15m)
mean(data15m)
#metodo dei momenti
mean(data15m) -> m15m
var(data15m) -> v15m
eulergamma <- 0.5772156649
b1.gumbel.15m <- sqrt(6*v15m)/pi
a1.gumbel.15m <- (m15m - b1.gumbel.15m*eulergamma)
z15 <- sort(data15m)
plot(z15,pgumbel(z15,loc=a1.gumbel.15m, scale=b1.gumbel.15m), col="blue", type="l", xlab="h[mm]", ylab="P[H
plot(ecdf15m, add=T)
z15
#metodo dei minimi quadrati
ec15 <- ecdf(sort(data15m))
ec15
z15 <- z15[-38]
Y15 <- -log(-log(ec15(z15)))
Y15
lsfit(z15,Y15) -> fts15
fts15
b2.gumbel.15m <- fts15$coefficients[[2]]^-1
a2.gumbel.15m <- -fts15$coefficients[[1]]*b2.gumbel.15m
#metodo di massima verosimiglianza
library("MASS")
fitdistr(data15m, densfun=dgumbel, start=list(loc=a1.gumbel.15m, scale=b1.gumbel.15m)) -> mlab.15m
mlab.15m
a1.gumbel.15m
b1.gumbel.15m
a2.gumbel.15m
b2.gumbel.15m
a3.gumbel.15m <- 9.7837073
b3.gumbel.15m <- 3.0869798
#test di pearson
plot(z15, pgumbel(z15, loc=a1.gumbel.15m, scale=b1.gumbel.15m), type="l", col="blue")
plot(z15, pgumbel(z15, loc=a2.gumbel.15m, scale=b2.gumbel.15m), type="l", col="purple")
plot(z15, pgumbel(z15, loc=a3.gumbel.15m, scale=b3.gumbel.15m), type="l", col="red")
q <- c(0.2, 0.4, 0.6, 0.8)
qgumbel(q, loc=a1.gumbel.15m, scale=b1.gumbel.15m) -> qi_1.15m
qi_1.15m
data15m
c(0,ec15(qi_1.15m))*38 -> no1_1.15m
c(ec15(qi_1.15m)*38,38) -> no2_1.15m
no2_1.15m - no1_1.15m -> no_1.15m
0.2*length(data15m) -> deltapi_1.15m
X1.15m <- sum((no_1.15m - deltapi_1.15m)^2/deltapi_1.15m)
X1.15m
qgumbel(q, loc=a2.gumbel.15m, scale=b2.gumbel.15m) -> qi_2.15m
qi_2.15m
c(0,ec15(qi_2.15m)*38) -> no1_2.15m
c(ec15(qi_2.15m)*38,38) -> no2_2.15m
no2_2.15m - no1_2.15m -> no_2.15m
0.2*length(data15m) -> deltapi_2.15m
X2.15m <- sum((no_2.15m - deltapi_2.15m)^2/deltapi_2.15m)
X2.15m
qgumbel(q, loc=a3.gumbel.15m, scale=b3.gumbel.15m) -> qi_3.15m
qi_3.15m
c(0,ec15(qi_3.15m))*38 -> no1_3.15m
c(ec15(qi_3.15m)*38,38) -> no2_3.15m
no2_3.15m - no1_3.15m -> no_3.15m
0.2*length(data15m) -> deltapi_3.15m
X3.15m <- sum((no_3.15m - deltapi_3.15m)^2/deltapi_3.15m)
X3.15m
#Osservazione 1 ora
read.table("massimi di precipitazione.txt",header=TRUE) -> data
data[,2] -> data1h
data[,1] -> anni
#rimuovo i valori NA
anni <- anni[!is.na(data1h)]
data1h <- data1h[!is.na(data1h)]
data1h
hist(data1h, breaks=40)
plot(anni,data1h,xlab="anno",ylab="h [mm]",main="Precipitazioni max annuali di durata 1h",col="black",lwd="1",type="l")
ecdf(data1h) -> ecdf1h
ecdf1h
plot(ecdf1h)
mean(data1h)
#metodo dei momenti
library("evd")
mean(data1h) -> m1h
var(data1h) -> v1h
eulergamma <- 0.5772156649
b1.gumbel.1h <- sqrt(6*v1h)/pi
a1.gumbel.1h <- (m1h-b1.gumbel.1h*eulergamma)
z1 <-sort(data1h)
plot(z1,pgumbel(z1,loc=a1.gumbel.1h, scale=b1.gumbel.1h), col="blue", type="l", xlab="h[mm]", ylab="P[H
plot(ecdf1h, add=T)
z1
#metodo dei minimi quadrati
ec1 <- ecdf(sort(data1h))
ec1
z1 <- z1[-55]
Y1 <- -log(-log(ec1(z1)))
Y1
lsfit(z1,Y1) -> fts1
fts1
b2.gumbel.1h <- fts1$coefficients[[2]]^-1
a2.gumbel.1h <- -fts1$coefficients[[1]]*b2.gumbel.1h
#metodo di massima verosimiglianza
library("MASS")
fitdistr(data1h, densfun=dgumbel, start=list(loc=a1.gumbel.1h, scale=b1.gumbel.1h)) -> mlab.1h
mlab.1h
a1.gumbel.1h
b1.gumbel.1h
a2.gumbel.1h
b2.gumbel.1h
a3.gumbel.1h <- 16.6102029
b3.gumbel.1h <- 4.5843401
#test di pearson
plot(z1, pgumbel(z1, loc=a1.gumbel.1h, scale=b1.gumbel.1h), type="l", col="blue")
plot(z1, pgumbel(z1, loc=a2.gumbel.1h, scale=b2.gumbel.1h), type="l", col="purple")
plot(z1, pgumbel(z1, loc=a3.gumbel.1h, scale=b3.gumbel.1h), type="l", col="red")
q <- c(0.2, 0.4, 0.6, 0.8)
qgumbel(q, loc=a1.gumbel.1h, scale=b1.gumbel.1h) -> qi_1.1h
qi_1.1h
c(0,ec1(qi_1.1h))*55 -> no1_1.1h
c(ec1(qi_1.1h)*55,55) -> no2_1.1h
no2_1.1h - no1_1.1h -> no_1.1h
0.2*length(data1h) -> deltapi_1.1h
X1.1h <- sum((no_1.1h - deltapi_1.1h)^2/deltapi_1.1h)
X1.1h
qgumbel(q, loc=a2.gumbel.1h, scale=b2.gumbel.1h) -> qi_2.1h
qi_2.1h
c(0,ec1(qi_2.1h))*55 -> no1_2.1h
c(ec1(qi_2.1h)*55,55) -> no2_2.1h
no2_2.1h - no1_2.1h -> no_2.1h
0.2*length(data1h) -> deltapi_2.1h
X2.1h <- sum((no_2.1h - deltapi_2.1h)^2/deltapi_2.1h)
X2.1h
qgumbel(q, loc=a3.gumbel.1h, scale=b3.gumbel.1h) -> qi_3.1h
qi_3.1h
c(0,ec1(qi_3.1h))*55 -> no1_3.1h
c(ec1(qi_3.1h)*55,55) -> no2_3.1h
no2_3.1h - no1_3.1h -> no_3.1h
0.2*length(data1h) -> deltapi_3.1h
X3.1h <- sum((no_3.1h - deltapi_3.1h)^2/deltapi_3.1h)
X3.1h
#Osservazione 3 ore
read.table("massimi di precipitazione.txt",header=TRUE) -> data
data[,3] -> data3h
data[,1] -> anni
#rimuovo i valori NA
anni <- anni[!is.na(data3h)]
data3h <- data3h[!is.na(data3h)]
data3h
hist(data3h, breaks=40)
plot(anni,data3h,xlab="anno",ylab="h [mm]",main="Precipitazioni max annuali di durata 3h",col="black",lwd="1",type="l")
ecdf(data3h) -> ecdf3h
ecdf3h
plot(ecdf3h)
mean(data3h)
#metodo dei momenti
library("evd")
mean(data3h) -> m3h
var(data3h) -> v3h
eulergamma <- 0.5772156649
b1.gumbel.3h <- sqrt(6*v3h)/pi
a1.gumbel.3h <- (m3h - b1.gumbel.3h*eulergamma)
z3 <- sort(data3h)
plot(z3,pgumbel(z3,loc=a1.gumbel.3h, scale=b1.gumbel.3h), col="blue", type="l", xlab="h[mm]", ylab="P[H
plot(ecdf3h, add=T)
z3
#metodo dei minimi quadrati
ec3 <- ecdf(sort(data3h))
ec3
z3 <- z3[-55]
Y3 <- -log(-log(ec3(z3)))
Y3
lsfit(z3,Y3) -> fts3
fts3
b2.gumbel.3h <- fts3$coefficients[[2]]^-1
a2.gumbel.3h <- -fts3$coefficients[[1]]*b2.gumbel.3h
#metodo di massima verosimiglianza
library("MASS")
fitdistr(data3h, densfun=dgumbel, start=list(loc=a1.gumbel.3h, scale=b1.gumbel.3h)) -> mlab.3h
mlab.3h
a1.gumbel.3h
b1.gumbel.3h
a2.gumbel.3h
b2.gumbel.3h
a3.gumbel.3h <- 24.7465611
b3.gumbel.3h <- 6.1060564
#test di pearson
plot(z3, pgumbel(z3, loc=a1.gumbel.3h, scale=b1.gumbel.3h), type="l", col="blue")
plot(z3, pgumbel(z3, loc=a2.gumbel.3h, scale=b2.gumbel.3h), type="l", col="purple")
plot(z3, pgumbel(z3, loc=a3.gumbel.3h, scale=b3.gumbel.3h), type="l", col="red")
q <- c(0.2, 0.4, 0.6, 0.8)
qgumbel(q, loc=a1.gumbel.3h, scale=b1.gumbel.3h) -> qi_1.3h
qi_1.3h
data3h
c(0,ec3(qi_1.3h))*55 -> no1_1.3h
c(ec3(qi_1.3h)*55,55) -> no2_1.3h
no2_1.3h - no1_1.3h -> no_1.3h
0.2*length(data3h) -> deltapi_1.3h
X1.3h <- sum((no_1.3h - deltapi_1.3h)^2/deltapi_1.3h)
X1.3h
qgumbel(q, loc=a2.gumbel.3h, scale=b2.gumbel.3h) -> qi_2.3h
qi_2.3h
c(0,ec3(qi_2.3h)*55) -> no1_2.3h
c(ec3(qi_2.3h)*55,55) -> no2_2.3h
no2_2.3h - no1_2.3h -> no_2.3h
0.2*length(data3h) -> deltapi_2.3h
X2.3h <- sum((no_2.3h - deltapi_2.3h)^2/deltapi_2.3h)
X2.3h
qgumbel(q, loc=a3.gumbel.3h, scale=b3.gumbel.3h) -> qi_3.3h
qi_3.3h
c(0,ec3(qi_3.3h))*55 -> no1_3.3h
c(ec3(qi_3.3h)*55,55) -> no2_3.3h
no2_3.3h - no1_3.3h -> no_3.3h
0.2*length(data3h) -> deltapi_3.3h
X3.3h <- sum((no_3.3h - deltapi_3.3h)^2/deltapi_3.3h)
X3.3h
#Osservazione 6 ore
read.table("massimi di precipitazione.txt",header=TRUE) -> data
data[,4] -> data6h
data[,1] -> anni
#rimuovo i valori NA
anni <- anni[!is.na(data6h)]
data6h <- data6h[!is.na(data6h)]
data6h
hist(data6h, breaks=40)
plot(anni,data6h,xlab="anno",ylab="h [mm]",main="Precipitazioni max annuali di durata 6h",col="black",lwd="1",type="l")
ecdf(data6h) -> ecdf6h
ecdf6h
plot(ecdf6h)
mean(data6h)
#metodo dei momenti
library("evd")
mean(data6h) -> m6h
var(data6h) -> v6h
eulergamma <- 0.5772156649
b1.gumbel.6h <- sqrt(6*v6h)/pi
a1.gumbel.6h <- (m6h - b1.gumbel.6h*eulergamma)
z6 <- sort(data6h)
plot(z6,pgumbel(z6,loc=a1.gumbel.6h, scale=b1.gumbel.6h), col="blue", type="l", xlab="h[mm]", ylab="P[H
plot(ecdf6h, add=T)
z6
#metodo dei minimi quadrati
ec6 <- ecdf(sort(data6h))
ec6
z6 <- z6[-56]
z6 <- z6[-55]
Y6 <- -log(-log(ec6(z6)))
Y6
lsfit(z6,Y6) -> fts6
fts6
b2.gumbel.6h <- fts6$coefficients[[2]]^-1
a2.gumbel.6h <- -fts6$coefficients[[1]]*b2.gumbel.6h
#metodo di massima verosimiglianza
library("MASS")
fitdistr(data6h, densfun=dgumbel, start=list(loc=a1.gumbel.6h, scale=b1.gumbel.6h)) -> mlab.6h
mlab.6h
a1.gumbel.6h
b1.gumbel.6h
a2.gumbel.6h
b2.gumbel.6h
a3.gumbel.6h <- 34.5098162
b3.gumbel.6h <- 6.8655925
#test di pearson
plot(z6, pgumbel(z6, loc=a1.gumbel.6h, scale=b1.gumbel.6h), type="l", col="blue")
plot(z6, pgumbel(z6, loc=a2.gumbel.6h, scale=b2.gumbel.6h), type="l", col="purple")
plot(z6, pgumbel(z6, loc=a3.gumbel.6h, scale=b3.gumbel.6h), type="l", col="red")
q <- c(0.2, 0.4, 0.6, 0.8)
qgumbel(q, loc=a1.gumbel.6h, scale=b1.gumbel.6h) -> qi_1.6h
qi_1.6h
data6h
c(0,ec6(qi_1.6h))*56 -> no1_1.6h
c(ec6(qi_1.6h)*56,56) -> no2_1.6h
no2_1.6h - no1_1.6h -> no_1.6h
0.2*length(data6h) -> deltapi_1.6h
X1.6h <- sum((no_1.6h - deltapi_1.6h)^2/deltapi_1.6h)
X1.6h
qgumbel(q, loc=a2.gumbel.6h, scale=b2.gumbel.6h) -> qi_2.6h
qi_2.6h
c(0,ec6(qi_2.6h)*56) -> no1_2.6h
c(ec6(qi_2.6h)*56,56) -> no2_2.6h
no2_2.6h - no1_2.6h -> no_2.6h
0.2*length(data6h) -> deltapi_2.6h
X2.6h <- sum((no_2.6h - deltapi_2.6h)^2/deltapi_2.6h)
X2.6h
qgumbel(q, loc=a3.gumbel.6h, scale=b3.gumbel.6h) -> qi_3.6h
qi_3.6h
c(0,ec6(qi_3.6h))*56 -> no1_3.6h
c(ec6(qi_3.6h)*56,56) -> no2_3.6h
no2_3.6h - no1_3.6h -> no_3.6h
0.2*length(data6h) -> deltapi_3.6h
X3.6h <- sum((no_3.6h - deltapi_3.6h)^2/deltapi_3.6h)
X3.6h
#Osservazione 12 ore
read.table("massimi di precipitazione.txt",header=TRUE) -> data
data[,5] -> data12h
data[,1] -> anni
#rimuovo i valori NA
anni <- anni[!is.na(data12h)]
data12h <- data12h[!is.na(data12h)]
data12h
hist(data12h, breaks=40)
plot(anni,data12h,xlab="anno",ylab="h [mm]",main="Precipitazioni max annuali di durata 12h",col="black",lwd="1",type="l")
ecdf(data12h) -> ecdf12h
ecdf12h
plot(ecdf12h)
mean(data12h)
#metodo dei momenti
library("evd")
mean(data12h) -> m12h
var(data12h) -> v12h
eulergamma <- 0.5772156649
b1.gumbel.12h <- sqrt(6*v12h)/pi
a1.gumbel.12h <- (m12h - b1.gumbel.12h*eulergamma)
z12 <- sort(data12h)
plot(z12,pgumbel(z12,loc=a1.gumbel.12h, scale=b1.gumbel.12h), col="blue", type="l", xlab="h[mm]", ylab="P[H
plot(ecdf12h, add=T)
z12
#metodo dei minimi quadrati
ec12 <- ecdf(sort(data12h))
ec12
z12 <- z12[-56]
Y12 <- -log(-log(ec12(z12)))
Y12
lsfit(z12,Y12) -> fts12
fts12
b2.gumbel.12h <- fts12$coefficients[[2]]^-1
a2.gumbel.12h <- -fts12$coefficients[[1]]*b2.gumbel.12h
#metodo di massima verosimiglianza
library("MASS")
fitdistr(data12h, densfun=dgumbel, start=list(loc=a1.gumbel.12h, scale=b1.gumbel.12h)) -> mlab.12h
mlab.12h
a1.gumbel.12h
b1.gumbel.12h
a2.gumbel.12h
b2.gumbel.12h
a3.gumbel.12h <- 48.039083
b3.gumbel.12h <- 11.613895
#test di pearson
plot(z12, pgumbel(z12, loc=a1.gumbel.12h, scale=b1.gumbel.12h), type="l", col="blue")
plot(z12, pgumbel(z12, loc=a2.gumbel.12h, scale=b2.gumbel.12h), type="l", col="purple")
plot(z12, pgumbel(z12, loc=a3.gumbel.12h, scale=b3.gumbel.12h), type="l", col="red")
q <- c(0.2, 0.4, 0.6, 0.8)
qgumbel(q, loc=a1.gumbel.12h, scale=b1.gumbel.12h) -> qi_1.12h
qi_1.12h
data12h
c(0,ec12(qi_1.12h))*56 -> no1_1.12h
c(ec12(qi_1.12h)*56,56) -> no2_1.12h
no2_1.12h - no1_1.12h -> no_1.12h
0.2*length(data12h) -> deltapi_1.12h
X1.12h <- sum((no_1.12h - deltapi_1.12h)^2/deltapi_1.12h)
X1.12h
qgumbel(q, loc=a2.gumbel.12h, scale=b2.gumbel.12h) -> qi_2.12h
qi_2.12h
c(0,ec12(qi_2.12h)*56) -> no1_2.12h
c(ec12(qi_2.12h)*56,56) -> no2_2.12h
no2_2.12h - no1_2.12h -> no_2.12h
0.2*length(data12h) -> deltapi_2.12h
X2.12h <- sum((no_2.12h - deltapi_2.12h)^2/deltapi_2.12h)
X2.12h
qgumbel(q, loc=a3.gumbel.12h, scale=b3.gumbel.12h) -> qi_3.12h
qi_3.12h
c(0,ec12(qi_3.12h))*56 -> no1_3.12h
c(ec12(qi_3.12h)*56,56) -> no2_3.12h
no2_3.12h - no1_3.12h -> no_3.12h
0.2*length(data12h) -> deltapi_3.12h
X3.12h <- sum((no_3.12h - deltapi_3.12h)^2/deltapi_3.12h)
X3.12h
#Osservazione 24 ore
read.table("massimi di precipitazione.txt",header=TRUE) -> data
data[,6] -> data24h
data[,1] -> anni
#rimuovo i valori NA
anni <- anni[!is.na(data24h)]
data24h <- data24h[!is.na(data24h)]
data24h
hist(data24h, breaks=40)
plot(anni,data24h,xlab="anno",ylab="h [mm]",main="Precipitazioni max annuali di durata 24h",col="black",lwd="1",type="l")
ecdf(data24h) -> ecdf24h
ecdf24h
plot(ecdf24h)
mean(data24h)
#metodo dei momenti
library("evd")
mean(data24h) -> m24h
var(data24h) -> v24h
eulergamma <- 0.5772156649
b1.gumbel.24h <- sqrt(6*v24h)/pi
a1.gumbel.24h <- (m24h - b1.gumbel.24h*eulergamma)
z24 <- sort(data24h)
plot(z24,pgumbel(z24,loc=a1.gumbel.24h, scale=b1.gumbel.24h), col="blue", type="l", xlab="h[mm]", ylab="P[H
plot(ecdf24h, add=T)
z24
#metodo dei minimi quadrati
ec24 <- ecdf(sort(data24h))
ec24
z24 <- z24[-56]
Y24 <- -log(-log(ec24(z24)))
Y24
lsfit(z24,Y24) -> fts24
fts24
b2.gumbel.24h <- fts24$coefficients[[2]]^-1
a2.gumbel.24h <- -fts24$coefficients[[1]]*b2.gumbel.24h
#metodo di massima verosimiglianza
library("MASS")
fitdistr(data24h, densfun=dgumbel, start=list(loc=a1.gumbel.24h, scale=b1.gumbel.24h)) -> mlab.24h
mlab.24h
a1.gumbel.24h
b1.gumbel.24h
a2.gumbel.24h
b2.gumbel.24h
a3.gumbel.24h <- 65.346178
b3.gumbel.24h <- 18.501232
#test di pearson
plot(z24, pgumbel(z24, loc=a1.gumbel.24h, scale=b1.gumbel.24h), type="l", col="blue")
plot(z24, pgumbel(z24, loc=a2.gumbel.24h, scale=b2.gumbel.24h), type="l", col="purple")
plot(z24, pgumbel(z24, loc=a3.gumbel.24h, scale=b3.gumbel.24h), type="l", col="red")
q <- c(0.2, 0.4, 0.6, 0.8)
qgumbel(q, loc=a1.gumbel.24h, scale=b1.gumbel.24h) -> qi_1.24h
qi_1.24h
data24h
c(0,ec24(qi_1.24h))*56 -> no1_1.24h
c(ec24(qi_1.24h)*56,56) -> no2_1.24h
no2_1.24h - no1_1.24h -> no_1.24h
0.2*length(data24h) -> deltapi_1.24h
X1.24h <- sum((no_1.24h - deltapi_1.24h)^2/deltapi_1.24h)
X1.24h
qgumbel(q, loc=a2.gumbel.24h, scale=b2.gumbel.24h) -> qi_2.24h
qi_2.24h
c(0,ec24(qi_2.24h)*56) -> no1_2.24h
c(ec24(qi_2.24h)*56,56) -> no2_2.24h
no2_2.24h - no1_2.24h -> no_2.24h
0.2*length(data24h) -> deltapi_2.24h
X2.24h <- sum((no_2.24h - deltapi_2.24h)^2/deltapi_2.24h)
X2.24h
qgumbel(q, loc=a3.gumbel.24h, scale=b3.gumbel.24h) -> qi_3.24h
qi_3.24h
c(0,ec24(qi_3.24h))*56 -> no1_3.24h
c(ec24(qi_3.24h)*56,56) -> no2_3.24h
no2_3.24h - no1_3.24h -> no_3.24h
0.2*length(data24h) -> deltapi_3.24h
X3.24h <- sum((no_3.24h - deltapi_3.24h)^2/deltapi_3.24h)
X3.24h
#15minuti
plot(z15, pgumbel(z15, loc=a1.gumbel.15m, scale=b1.gumbel.15m), col ="blue", type="l" , xlab="h[mm]", ylab="P[H
lines(z15, pgumbel(z15, loc=a2.gumbel.15m, scale=b2.gumbel.15m), col ="red", type="l" , xlab="h[mm]", ylab="P[H
lines(z15, pgumbel(z15, loc=a3.gumbel.15m, scale=b3.gumbel.15m), col ="green", type="l" , xlab="h[mm]", ylab="P[H
plot(ecdf15m, add=T)
legend(10.25,0.25,c("momenti","minimi quadrati","massima verosimiglianza "),lty=c(1,1,1), lwd=c(2.5,2.5,2.5), col=c("blue","red","green"))
#1h
plot(z1, pgumbel(z1, loc=a1.gumbel.1h, scale=b1.gumbel.1h), col ="blue", type="l" , xlab="h[mm]", ylab="P[H
lines(z1, pgumbel(z1, loc=a2.gumbel.1h, scale=b2.gumbel.1h), col ="red", type="l" , xlab="h[mm]", ylab="P[H
lines(z1, pgumbel(z1, loc=a3.gumbel.1h, scale=b3.gumbel.1h), col ="green", type="l" , xlab="h[mm]", ylab="P[H
plot(ecdf1h, add=T)
legend(19,0.25,c("momenti","minimi quadrati","massima verosimiglianza "),lty=c(1,1,1), lwd=c(2.5,2.5,2.5), col=c("blue","red","green"))
#3h
plot(z3, pgumbel(z3, loc=a1.gumbel.3h, scale=b1.gumbel.3h), col ="blue", type="l" , xlab="h[mm]", ylab="P[H
lines(z3, pgumbel(z3, loc=a2.gumbel.3h, scale=b2.gumbel.3h), col ="red", type="l" , xlab="h[mm]", ylab="P[H
lines(z3, pgumbel(z3, loc=a3.gumbel.3h, scale=b3.gumbel.3h), col ="green", type="l" , xlab="h[mm]", ylab="P[H
plot(ecdf3h, add=T)
legend(26,0.25,c("momenti","minimi quadrati","massima verosimiglianza "),lty=c(1,1,1), lwd=c(2.5,2.5,2.5), col=c("blue","red","green"))
#6h
plot(z6, pgumbel(z6, loc=a1.gumbel.6h, scale=b1.gumbel.6h), col ="blue", type="l" , xlab="h[mm]", ylab="P[H
lines(z6, pgumbel(z6, loc=a2.gumbel.6h, scale=b2.gumbel.6h), col ="red", type="l" , xlab="h[mm]", ylab="P[H
lines(z6, pgumbel(z6, loc=a3.gumbel.6h, scale=b3.gumbel.6h), col ="green", type="l" , xlab="h[mm]", ylab="P[H
plot(ecdf6h, add=T)
legend(35,0.25,c("momenti","minimi quadrati","massima verosimiglianza "),lty=c(1,1,1), lwd=c(2.5,2.5,2.5), col=c("blue","red","green"))
#12h
plot(z12, pgumbel(z12, loc=a1.gumbel.12h, scale=b1.gumbel.12h), col ="blue", type="l" , xlab="h[mm]", ylab="P[H
lines(z12, pgumbel(z12, loc=a2.gumbel.12h, scale=b2.gumbel.12h), col ="red", type="l" , xlab="h[mm]", ylab="P[H
lines(z12, pgumbel(z12, loc=a3.gumbel.12h, scale=b3.gumbel.12h), col ="green", type="l" , xlab="h[mm]", ylab="P[H
plot(ecdf12h, add=T)
legend(54,0.25,c("momenti","minimi quadrati","massima verosimiglianza "),lty=c(1,1,1), lwd=c(2.5,2.5,2.5), col=c("blue","red","green"))
#24h
plot(z24, pgumbel(z24, loc=a1.gumbel.24h, scale=b1.gumbel.24h), col ="blue", type="l" , xlab="h[mm]", ylab="P[H
lines(z24, pgumbel(z24, loc=a2.gumbel.24h, scale=b2.gumbel.24h), col ="red", type="l" , xlab="h[mm]", ylab="P[H
lines(z24, pgumbel(z24, loc=a3.gumbel.24h, scale=b3.gumbel.24h), col ="green", type="l" , xlab="h[mm]", ylab="P[H
plot(ecdf24h, add=T)
legend(76,0.25,c("momenti","minimi quadrati","massima verosimiglianza "),lty=c(1,1,1), lwd=c(2.5,2.5,2.5), col=c("blue","red","green"))
#Analisi complessiva
#Si scelgono le coppie di valori (a,b) corrispondenti a x_i.nh minore
a.15m <- 9.84838642862797
b.15m <- 2.86134571842947
a.1h <- 16.6102029
b.1h <- 4.5843401
a.3h <- 24.7465611
b.3h <- 6.1060564
a.6h <- 34.5756040136424
b.6h <- 6.58846882748384
a.12h <- 47.5046985231325
b.12h <- 11.3393627812782
a.24h <- 64.8047813576029
b.24h <- 18.3000502827714
#Funzione di probabilità
seq(from=1, to=200, by=0.1) -> x
plot(x, pgumbel(x, loc=a.15m, scale=b.15m), type="l", xlab="h[mm]", ylab="P[h]", col="purple")
lines(x, pgumbel(x, loc=a.1h, scale=b.1h), col= "blue")
lines(x, pgumbel(x, loc=a.3h, scale=b.3h), col="green")
lines(x, pgumbel(x, loc=a.6h, scale=b.6h), col="yellow")
lines(x, pgumbel(x, loc=a.12h, scale=b.12h), col="orange")
lines(x, pgumbel(x, loc=a.24h, scale=b.24h), col="red")
text(10,0.97,"15min", cex=0.8)
text(20,0.81,"1h",cex=0.8)
text(30,0.7,"3h",cex=0.8)
text(40,0.55,"6h",cex=0.8)
text(55,0.43,"12h",cex=0.8)
text(70,0.35,"24h",cex=0.8)
legend(150,0.37,c("15 minuti ","1 ora","3 ore","6 ore", "12 ore", "24 ore"),lty=c(1,1,1,1,1,1), lwd=c(2.5,2.5,2.5,2.5,2.5,2.5), col=c("purple","blue","green","yellow","orange","red"))
#Densità di probabilità delle curve di Gumbel
plot(x,dgumbel(x,a.15m,b.15m),type="l",col="purple")
lines(x,dgumbel(x,a.1h,b.1h), type="l", col="blue")
lines(x,dgumbel(x,a.3h,b.3h),type="l",col="green")
lines(x,dgumbel(x,a.6h,b.6h),type="l",col="yellow")
lines(x,dgumbel(x,a.12h,b.12h),type="l",col="orange")
lines(x,dgumbel(x,a.24h,b.24h),type="l",col="red")
text(20,0.12,"15min", cex=0.8)
text(20,0.083,"1h",cex=0.8)
text(30,0.062,"3h",cex=0.8)
text(42,0.055,"6h",cex=0.8)
text(50,0.035,"12h",cex=0.8)
text(70,0.023,"24h",cex=0.8)
legend(140,0.13,c("15 minuti ","1 ora","3 ore","6 ore", "12 ore", "24 ore"),lty=c(1,1,1,1,1,1), lwd=c(2.5,2.5,2.5,2.5,2.5,2.5), col=c("purple","blue","green","yellow","orange","red"))
#Curve di possibilità pluviometrica con tempo di ritorno 5,10,25,50,100 anni
p <- function(x){1-1/x}
p(5)
length(x)
rep(0.8,1991) -> y
length(y)
lines(x,y, col="black")
#5anni
c(qgumbel(0.8,loc=a.15m,scale=b.15m),
qgumbel(0.8,loc=a.1h,scale=b.1h),
qgumbel(0.8,loc=a.3h,scale=b.3h),
qgumbel(0.8,loc=a.6h,scale=b.6h),
qgumbel(0.8,loc=a.12h,scale=b.12h),
qgumbel(0.8,loc=a.24h,scale=b.24h)) -> h5
d <- c(0.25,1,3,6,12,24)
lsfit(log(d),log(h5)) -> ft5
ft5$coefficients[[2]]->n5
n5
exp(ft5$coefficients[[1]])->atp5
atp5
p <- function(x){1-1/x}
p(10)
length(x)
#10anni
c(qgumbel(0.9,loc=a.15m,scale=b.15m),
qgumbel(0.9,loc=a.1h,scale=b.1h),
qgumbel(0.9,loc=a.3h,scale=b.3h),
qgumbel(0.9,loc=a.6h,scale=b.6h),
qgumbel(0.9,loc=a.12h,scale=b.12h),
qgumbel(0.9,loc=a.24h,scale=b.24h)) -> h10
d <- c(0.25,1,3,6,12,24)
lsfit(log(d),log(h10)) -> ft10
ft10$coefficients[[2]]->n10
n10
exp(ft10$coefficients[[1]])->atp10
atp10
p <- function(x){1-1/x}
p(25)
length(x)
#25anni
c(qgumbel(0.96,loc=a.15m,scale=b.15m),
qgumbel(0.96,loc=a.1h,scale=b.1h),
qgumbel(0.96,loc=a.3h,scale=b.3h),
qgumbel(0.96,loc=a.6h,scale=b.6h),
qgumbel(0.96,loc=a.12h,scale=b.12h),
qgumbel(0.96,loc=a.24h,scale=b.24h)) -> h25
d <- c(0.25,1,3,6,12,24)
lsfit(log(d),log(h25)) -> ft25
ft25$coefficients[[2]]->n25
n25
exp(ft25$coefficients[[1]])->atp25
atp25
p <- function(x){1-1/x}
p(50)
length(x)
#50anni
c(qgumbel(0.98,loc=a.15m,scale=b.15m),
qgumbel(0.98,loc=a.1h,scale=b.1h),
qgumbel(0.98,loc=a.3h,scale=b.3h),
qgumbel(0.98,loc=a.6h,scale=b.6h),
qgumbel(0.98,loc=a.12h,scale=b.12h),
qgumbel(0.98,loc=a.24h,scale=b.24h)) -> h50
d <- c(0.25,1,3,6,12,24)
lsfit(log(d),log(h50)) -> ft50
ft50$coefficients[[2]]->n50
n50
exp(ft50$coefficients[[1]])->atp50
atp50
p <- function(x){1-1/x}
p(100)
length(x)
#100anni
c(qgumbel(0.99,loc=a.15m,scale=b.15m),
qgumbel(0.99,loc=a.1h,scale=b.1h),
qgumbel(0.99,loc=a.3h,scale=b.3h),
qgumbel(0.99,loc=a.6h,scale=b.6h),
qgumbel(0.99,loc=a.12h,scale=b.12h),
qgumbel(0.99,loc=a.24h,scale=b.24h)) -> h100
d <- c(0.25,1,3,6,12,24)
lsfit(log(d),log(h100)) -> ft100
ft100$coefficients[[2]]->n100
n100
exp(ft100$coefficients[[1]])->atp100
atp100
#INSIEME piano logaritmico
plot(d,exp(ft5$coefficients[[1]])*d^ft5$coefficients[[2]],type="l", col="purple", xlab="durata [h]",ylab="h [mm]",main="LineeSegnalitrici di Possibilità Pluviometrica")
lines(d,exp(ft10$coefficients[[1]])*d^ft10$coefficients[[2]],type="l",col="blue")
lines(d,exp(ft25$coefficients[[1]])*d^ft25$coefficients[[2]],type="l",col="green")
lines(d,exp(ft50$coefficients[[1]])*d^ft50$coefficients[[2]],type="l",col="orange")
lines(d,exp(ft100$coefficients[[1]])*d^ft100$coefficients[[2]],type="l",col="red")
legend(17,40,c("5 anni","10 anni","25 anni","50 anni","100 anni "),lty=c(1,1,1,1,1), lwd=c(2.5,2.5,2.5,2.5,2.5), col=c("purple","blue","green","orange","red"))
points(d,h5,pch=0)
points(d,h10,pch=1)
points(d,h25,pch=2)
points(d,h50,pch=3)
points(d,h100,pch=4)
#INSIEME piano normale
plot(d,exp(ft5$coefficients[[1]])*d^ft5$coefficients[[2]], type="l", col="purple", xlab="durata [h]", ylab="h [mm]",main="Linee Segnalatrici di Possibilita' Pluviometrica", log="xy")
lines(d,exp(ft10$coefficients[[1]])*d^ft10$coefficients[[2]],type="l",col="blue", log="xy")
lines(d,exp(ft25$coefficients[[1]])*d^ft25$coefficients[[2]],type="l",col="green", log="xy")
lines(d,exp(ft50$coefficients[[1]])*d^ft50$coefficients[[2]],type="l",col="orange", log="xy")
lines(d,exp(ft100$coefficients[[1]])*d^ft100$coefficients[[2]],type="l",col="red", log="xy")
legend(6,27,c("5 anni","10 anni","25 anni","50 anni","100 anni "),lty=c(1,1,1,1,1), lwd=c(2.5,2.5,2.5,2.5,2.5), col=c("purple","blue","green","orange","red"))
points(d,h5,pch=0)
points(d,h10,pch=1)
points(d,h25,pch=2)
points(d,h50,pch=3)
points(d,h100,pch=4)
