Nothing
# ---------------------------------------------------------------- #
plotpos <- function(x, a=0, orient="xF", ...) {
# INPUT
# x = colonna
# ... = graphical parameters as xlab, ylab, main...
# OUTPUT
# rappresentazione in carta cartesiana della serie
ordinato <- sort(x)
n <- length(ordinato)
i <- 1:n
plotpos <- (i - a)/(n + 1 - 2*a)
#x11(width = 7, height = 7, pointsize = 12)
x=ordinato
F=plotpos
if (orient=="xF") plot(x,F,...)
else if (orient=="Fx") plot(F,x,...)
else stop("plotpos(x, a, orient): orient unknown")
grid()
}
# ---------------------------------------------------------------- #
pointspos <- function(x, a=0, orient="xF", ...) {
# INPUT
# x = colonna
# ... = graphical parameters as xlab, ylab, main...
# OUTPUT
# rappresentazione in carta cartesiana della serie
ordinato <- sort(x)
n <- length(ordinato)
i <- 1:n
plotpos <- (i - a)/(n + 1 - 2*a)
x=ordinato
F=plotpos
if (orient=="xF") points(x,F,...)
else if (orient=="Fx") points(F,x,...)
else stop("pointspos(x, a, orient): orient unknown")
}
# ---------------------------------------------------------------- #
loglogplot <- function(x, a=0, orient="xF", ...) {
# INPUT
# x = colonna
# ... = graphical parameters as xlab, ylab, main...
# OUTPUT
# rappresentazione log-log
ordinato <- sort(x)
n <- length(ordinato)
i <- 1:n
plotpos <- (i - a)/(n + 1 - 2*a)
x=ordinato
F=plotpos
if (orient=="xF") plot(log(x),log(1-F),...)
else if (orient=="Fx") plot(log(1-F),log(x),...)
else stop("loglogplot(x, a, orient): orient unknown")
grid()
}
# ---------------------------------------------------------------- #
loglogpoints <- function(x, a=0, orient="xF", ...) {
ordinato <- sort(x)
n <- length(ordinato)
i <- 1:n
plotpos <- (i - a)/(n + 1 - 2*a)
x=ordinato
F=plotpos
if (orient=="xF") points(log(x),log(1-F),...)
else if (orient=="Fx") points(log(1-F),log(x),...)
else stop("loglogpoints(x, a, orient): orient unknown")
}
# ---------------------------------------------------------------- #
normplot <- function(x, a=0, orient="xF", line=FALSE,...) {
# INPUT
# x = colonna
# line = linea rappresentante la distribuzione
# ... = graphical parameters as xlab, ylab, main...
# OUTPUT
# rappresentazione in carta probabilistica normale della serie
ordinato <- sort(x)
m <- mean(ordinato)
s <- sd(ordinato)
u <- (ordinato-m)/s
n <- length(ordinato)
i <- 1:n
plotpos <- (i - a)/(n + 1 - 2*a)
uteorica <- qnorm(plotpos)
pp <- c(.001, .01, .05, .1, .25, .50, .75, .9, .95, .99, .999)
qq <- qnorm(pp)
#x11(width = 7, height = 7, pointsize = 12)
x=ordinato
F=uteorica
if (orient=="xF") {
plot(x,F,axes=FALSE,...)
axis(1)
axis(2,at=qq,labels=pp)
grid(nx=NULL, ny=NA)
abline(h=qq, col="gray", lty="dotted")
if(line==TRUE) lines(ordinato,u,lty=2)
}
else if (orient=="Fx") {
plot(F,x,axes=FALSE,...)
axis(2)
axis(1,at=qq,labels=pp)
grid(ny=NULL, nx=NA)
abline(v=qq, col="gray", lty="dotted")
if(line==TRUE) lines(u,ordinato,lty=2)
}
else stop("normplot(x, a, orient, line): orient unknown")
box()
}
# ---------------------------------------------------------------- #
normpoints <- function(x, a=0, orient="xF", ...) {
# INPUT
# x = colonna
# line = linea rappresentante la distribuzione
# ... = graphical parameters as xlab, ylab, main...
# OUTPUT
# rappresentazione in carta probabilistica normale della serie
ordinato <- sort(x)
n <- length(ordinato)
i <- 1:n
plotpos <- (i - a)/(n + 1 - 2*a)
uteorica <- qnorm(plotpos)
x=ordinato
F=uteorica
if (orient=="xF") points(x,F,...)
else if (orient=="Fx") points(F,x,...)
else stop("normpoints(x, a, orient): orient unknown")
}
# ---------------------------------------------------------------- #
lognormplot <- function(x, a=0, orient="xF", line=FALSE,...) {
# INPUT
# x = colonna
# line = linea rappresentante la distribuzione
# ... = graphical parameters as xlab, ylab, main...
# OUTPUT
# rappresentazione in carta probabilistica lognormale della serie
ordinato <- sort(x)
logordinato <- log(ordinato)
m <- mean(logordinato)
s <- sd(logordinato)
u <- (logordinato-m)/s
n <- length(ordinato)
i <- 1:n
plotpos <- (i - a)/(n + 1 - 2*a)
uteorica <- qnorm(plotpos)
pp <- c(.001, .01, .05, .1, .25, .50, .75, .9, .95, .99, .999)
qq <- qnorm(pp)
#x11(width = 7, height = 7, pointsize = 12)
x=ordinato
F=uteorica
if (orient=="xF") {
plot(x,F,log="x",axes=FALSE,...)
axis(1)
axis(2,at=qq,labels=pp)
abline(h=qq, col="gray", lty="dotted")
grid(nx=NULL, ny=NA)
if (line==TRUE) lines(ordinato,u,lty=2)
}
else if (orient=="Fx") {
plot(F,x,log="y",axes=FALSE,...)
axis(2)
axis(1,at=qq,labels=pp)
abline(v=qq, col="gray", lty="dotted")
grid(ny=NULL, nx=NA)
if (line==TRUE) lines(u, ordinato, lty=2)
}
else stop("lognormplot(x, a, orient, line): orient unknown")
box()
}
# ---------------------------------------------------------------- #
studentplot <- function(x, df, a=0, orient="xF", line=FALSE,...) {
# INPUT
# df = degrees of freedom
ordinato <- sort(x)
n <- length(ordinato)
i <- 1:n
plotpos <- (i - a)/(n + 1 - 2*a)
uteorica <- qt(plotpos, df)
pp <- c(.001, .01, .05, .1, .25, .50, .75, .9, .95, .99, .999)
qq <- qt(pp, df)
#x11(width = 7, height = 7, pointsize = 12)
x=ordinato
F=uteorica
if (orient=="xF") {
plot(x,F,axes=FALSE,...)
axis(1)
axis(2,at=qq,labels=pp)
grid(nx=NULL, ny=NA)
abline(h=qq, col="gray", lty="dotted")
#if(line==TRUE) lines(ordinato,u,lty=2)
}
else if (orient=="Fx") {
plot(F,x,axes=FALSE,...)
axis(2)
axis(1,at=qq,labels=pp)
grid(ny=NULL, nx=NA)
abline(v=qq, col="gray", lty="dotted")
#if(line==TRUE) lines(u,ordinato,lty=2)
}
else stop("studentplot(x, df, a, orient, line): orient unknown")
box()
}
# ---------------------------------------------------------------- #
studentpoints <- function(x, df, a=0, orient="xF", ...) {
ordinato <- sort(x)
n <- length(ordinato)
i <- 1:n
plotpos <- (i - a)/(n + 1 - 2*a)
uteorica <- qt(plotpos, df)
x=ordinato
F=uteorica
if (orient=="xF") points(x,F,...)
else if (orient=="Fx") points(F,x,...)
else stop("studentpoints(x, df, a, orient): orient unknown")
}
# ---------------------------------------------------------------- #
logisplot <- function(x, a=0, orient="xF", line=FALSE,...) {
ordinato <- sort(x)
mu <- mean(ordinato)
s <- sqrt(3)*sd(ordinato)/pi
adim <- (ordinato - mu)/s
n <- length(ordinato)
i <- 1:n
plotpos <- (i - a)/(n + 1 - 2*a)
uteorica <- qlogis(plotpos)
pp <- c(.001, .01, .05, .1, .25, .50, .75, .9, .95, .99, .999)
qq <- qlogis(pp)
x=ordinato
F=uteorica
if (orient=="xF") {
plot(x,F,axes=FALSE,...)
axis(1)
axis(2,at=qq,labels=pp)
grid(nx=NULL, ny=NA)
abline(h=qq, col="gray", lty="dotted")
if(line==TRUE) lines(ordinato,adim,lty=2)
}
else if (orient=="Fx") {
plot(F,x,axes=FALSE,...)
axis(2)
axis(1,at=qq,labels=pp)
grid(ny=NULL, nx=NA)
abline(v=qq, col="gray", lty="dotted")
if(line==TRUE) lines(adim,ordinato,lty=2)
}
else stop("logisplot(x, a, orient, line): orient unknown")
box()
}
# ---------------------------------------------------------------- #
logispoints <- function(x, a=0, orient="xF", ...) {
ordinato <- sort(x)
n <- length(ordinato)
i <- 1:n
plotpos <- (i - a)/(n + 1 - 2*a)
uteorica <- qlogis(plotpos)
x=ordinato
F=uteorica
if (orient=="xF") points(x,F,...)
else if (orient=="Fx") points(F,x,...)
else stop("logispoints(x, a, orient): orient unknown")
}
# ---------------------------------------------------------------- #
gammaplot <- function(x, shape, a=0, orient="xF", line=FALSE,...) {
# INPUT
# shape of the gamma distribution
ordinato <- sort(x)
n <- length(ordinato)
i <- 1:n
plotpos <- (i - a)/(n + 1 - 2*a)
uteorica <- qgamma(plotpos, shape)
pp <- c(.001, .01, .05, .1, .25, .50, .75, .9, .95, .99, .999)
qq <- qgamma(pp, shape)
#x11(width = 7, height = 7, pointsize = 12)
x=ordinato
F=uteorica
if (orient=="xF") {
plot(x,F,axes=FALSE,...)
axis(1)
axis(2,at=qq,labels=pp)
grid(nx=NULL, ny=NA)
abline(h=qq, col="gray", lty="dotted")
#if(line==TRUE) lines(ordinato,u,lty=2)
}
else if (orient=="Fx") {
plot(F,x,axes=FALSE,...)
axis(2)
axis(1,at=qq,labels=pp)
grid(ny=NULL, nx=NA)
abline(v=qq, col="gray", lty="dotted")
#if(line==TRUE) lines(u,ordinato,lty=2)
}
else stop("gammaplot(x, shape, a, orient, line): orient unknown")
box()
}
# ---------------------------------------------------------------- #
gammapoints <- function(x, shape, a=0, orient="xF", ...) {
ordinato <- sort(x)
n <- length(ordinato)
i <- 1:n
plotpos <- (i - a)/(n + 1 - 2*a)
uteorica <- qgamma(plotpos, shape)
x=ordinato
F=uteorica
if (orient=="xF") points(x,F,...)
else if (orient=="Fx") points(F,x,...)
else stop("gammapoints(x, shape, a, orient): orient unknown")
}
# ---------------------------------------------------------------- #
gumbelplot <- function(x, a=0, orient="xF", line=FALSE,...) {
# INPUT
# x = colonna
# line = linea rappresentante la distribuzione
# ... = graphical parameters as xlab, ylab, main...
# OUTPUT
# rappresentazione in carta probabilistica di Gumbel della serie
ordinato <- sort(x)
m <- mean(ordinato)
s <- sd(ordinato)
beta <- (s*sqrt(6))/pi
mu <- m - 0.5772*beta
adim <- (ordinato-mu)/beta
n <- length(ordinato)
i <- 1:n
plotpos <- (i - a)/(n + 1 - 2*a)
uteorica <- -log(-log(plotpos))
pp <- c(.001, .01, .05, .1, .25, .50, .75, .9, .95, .99, .999)
qq <- -log(-log(pp))
#x11(width = 7, height = 7, pointsize = 12)
x=ordinato
F=uteorica
if (orient=="xF") {
plot(x,F,axes=FALSE,...)
axis(1)
axis(2,at=qq,labels=pp)
abline(h=qq, col="gray", lty="dotted")
grid(nx=NULL, ny=NA)
if (line==TRUE) lines(ordinato,adim,lty=2)
}
else if (orient=="Fx") {
plot(F,x,axes=FALSE,...)
axis(2)
axis(1,at=qq,labels=pp)
abline(v=qq, col="gray", lty="dotted")
grid(ny=NULL, nx=NA)
if (line==TRUE) lines(adim, ordinato, lty=2)
}
else stop("gumbelplot(x, a, orient, line): orient unknown")
box()
}
# ---------------------------------------------------------------- #
gumbelpoints <- function(x, a=0, orient="xF", ...) {
# INPUT
# x = colonna
# line = linea rappresentante la distribuzione
# ... = graphical parameters as xlab, ylab, main...
# OUTPUT
# rappresentazione in carta probabilistica di Gumbel della serie
ordinato <- sort(x)
n <- length(ordinato)
i <- 1:n
plotpos <- (i - a)/(n + 1 - 2*a)
uteorica <- -log(-log(plotpos))
x=ordinato
F=uteorica
if (orient=="xF") points(x,F,...)
else if (orient=="Fx") points(F,x,...)
else stop("gumbelpoints(x, a, orient): orient unknown")
}
# ---------------------------------------------------------------- #
frechetplot <- function(x, a=0, orient="xF", line=FALSE,...) {
ordinato <- sort(x)
logordinato <- log(ordinato)
m <- mean(logordinato)
s <- sd(logordinato)
beta <- (s*sqrt(6))/pi
mu <- m - 0.5772*beta
adim <- (logordinato-mu)/beta
n <- length(ordinato)
i <- 1:n
plotpos <- (i - a)/(n + 1 - 2*a)
uteorica <- -log(-log(plotpos))
pp <- c(.001, .01, .05, .1, .25, .50, .75, .9, .95, .975, .99, .995, .999)
qq <- -log(-log(pp))
#x11(width = 7, height = 7, pointsize = 12)
x=ordinato
F=uteorica
if (orient=="xF") {
plot(x,F,log="x",axes=FALSE,...)
axis(1)
axis(2,at=qq,labels=pp)
abline(h=qq, col="gray", lty="dotted")
grid(nx=NULL, ny=NA)
if (line==TRUE) lines(ordinato,adim,lty=2)
}
else if (orient=="Fx") {
plot(F,x,log="y",axes=FALSE,...)
axis(2)
axis(1,at=qq,labels=pp)
abline(v=qq, col="gray", lty="dotted")
grid(ny=NULL, nx=NA)
if (line==TRUE) lines(adim, ordinato, lty=2)
}
else stop("frechetplot(x, a, orient, line): orient unknown")
box()
}
# ---------------------------------------------------------------- #
weibullplot <- function(x, a=0, orient="xF", line=FALSE,...) {
ordinato <- sort(x)
n <- length(ordinato)
i <- 1:n
plotpos <- (i - a)/(n + 1 - 2*a)
uteorica <- log(-log(1 - plotpos))
pp <- c(.001, .01, .05, .1, .25, .50, .75, .9, .95, .975, .99, .995, .999)
qq <- log(-log(1 - pp))
x=ordinato
F=uteorica
if (orient=="xF") {
plot(x,F,log="x",axes=FALSE,...)
axis(1)
axis(2,at=qq,labels=pp)
abline(h=qq, col="gray", lty="dotted")
grid(nx=NULL, ny=NA)
#if (line==TRUE) lines(ordinato,adim,lty=2)
}
else if (orient=="Fx") {
plot(F,x,log="y",axes=FALSE,...)
axis(2)
axis(1,at=qq,labels=pp)
abline(v=qq, col="gray", lty="dotted")
grid(ny=NULL, nx=NA)
#if (line==TRUE) lines(adim, ordinato, lty=2)
}
else stop("weibullplot(x, a, orient, line): orient unknown")
box()
}
# ---------------------------------------------------------------- #
weibullpoints <- function(x, a=0, orient="xF", ...) {
ordinato <- sort(x)
n <- length(ordinato)
i <- 1:n
plotpos <- (i - a)/(n + 1 - 2*a)
uteorica <- log(-log(1 - plotpos))
x=ordinato
F=uteorica
if (orient=="xF") points(x,F,...)
else if (orient=="Fx") points(F,x,...)
else stop("weibullpoints(x, a, orient): orient unknown")
}
# ---------------------------------------------------------------- #
unifplot <- function (x, a=0, orient="xF", line=FALSE, ...) {
ordinato <- sort(x)
u = (ordinato - min(ordinato))/(max(ordinato) - min(ordinato))
n <- length(ordinato)
i <- 1:n
plotpos <- (i - a)/(n + 1 - 2*a)
uteorica <- qunif(plotpos)
pp <- c(0.01, 0.05, 0.1, 0.25, 0.5, 0.75, 0.9, 0.95, 0.99)
qq <- qunif(pp)
x = ordinato
F = uteorica
if (orient=="xF") {
plot(x, F, axes = FALSE, ...)
axis(1)
axis(2, at = qq, labels = pp)
abline(h = qq, col = "gray", lty = "dotted")
grid(nx = NULL, ny = NA)
if (line == TRUE) lines(ordinato, u, lty = 2)
}
else if (orient=="Fx") {
plot(F, x, axes = FALSE, ...)
axis(2)
axis(1, at = qq, labels = pp)
abline(v = qq, col = "gray", lty = "dotted")
grid(ny = NULL, nx = NA)
if (line == TRUE) lines(u, ordinato, lty = 2)
}
else stop("unifplot(x, a, orient, line): orient unknown")
box()
}
# ---------------------------------------------------------------- #
unifpoints <- function (x, a=0, orient="xF", ...) {
pointspos(x, a=0, orient="xF", ...)
}
# ---------------------------------------------------------------- #
expplot <- function(x, a=0, orient="xF", line=FALSE,...) {
# INPUT
# x = colonna
# line = linea rappresentante la distribuzione
# ... = graphical parameters as xlab, ylab, main...
# OUTPUT
# rappresentazione in carta probabilistica esponenziale della serie
ordinato <- sort(x)
m <- mean(ordinato)
lambda <- 1/m
adim <- lambda*ordinato
n <- length(ordinato)
i <- 1:n
plotpos <- (i - a)/(n + 1 - 2*a)
uteorica <- -log(1 - plotpos)
pp <- c(.1, .25, .50, .75, .9, .95, .975, .99, .995, .999)
qq <- -log(1 - pp)
#x11(width = 7, height = 7, pointsize = 12)
x=ordinato
F=uteorica
if (orient=="xF") {
plot(x,F,axes=FALSE,...)
axis(1)
axis(2,at=qq,labels=pp)
abline(h=qq, col="gray", lty="dotted")
grid(nx=NULL, ny=NA)
if (line==TRUE) lines(ordinato,adim,lty=2)
}
else if (orient=="Fx") {
plot(F,x,axes=FALSE,...)
axis(2)
axis(1,at=qq,labels=pp)
abline(v=qq, col="gray", lty="dotted")
grid(ny=NULL, nx=NA)
if (line==TRUE) lines(adim, ordinato, lty=2)
}
else stop("expplot(x, a, orient, line): orient unknown")
box()
}
# ---------------------------------------------------------------- #
exppoints <- function(x, a=0, orient="xF", ...) {
ordinato <- sort(x)
n <- length(ordinato)
i <- 1:n
plotpos <- (i - a)/(n + 1 - 2*a)
uteorica <- -log(1 - plotpos)
x=ordinato
F=uteorica
if (orient=="xF") points(x,F,...)
else if (orient=="Fx") points(F,x,...)
else stop("gumbelpoints(x, a, orient): orient unknown")
}
# ---------------------------------------------------------------- #
paretoplot <- function(x, a=0, orient="xF", line=FALSE,...) {
ordinato <- sort(x)
logordinato <- log(ordinato)
m <- mean(logordinato)
lambda <- 1/m
adim <- lambda*logordinato
n <- length(logordinato)
i <- 1:n
plotpos <- (i - a)/(n + 1 - 2*a)
uteorica <- -log(1 - plotpos)
pp <- c(.1, .25, .50, .75, .9, .95, .975, .99, .995, .999)
qq <- -log(1 - pp)
#x11(width = 7, height = 7, pointsize = 12)
x=ordinato
F=uteorica
if (orient=="xF") {
plot(x,F,log="x",axes=FALSE,...)
axis(1)
axis(2,at=qq,labels=pp)
abline(h=qq, col="gray", lty="dotted")
grid(nx=NULL, ny=NA)
if (line==TRUE) lines(ordinato,adim,lty=2)
}
else if (orient=="Fx") {
plot(F,x,log="y",axes=FALSE,...)
axis(2)
axis(1,at=qq,labels=pp)
abline(v=qq, col="gray", lty="dotted")
grid(ny=NULL, nx=NA)
if (line==TRUE) lines(adim, ordinato, lty=2)
}
else stop("expplot(x, a, orient, line): orient unknown")
box()
}
# ---------------------------------------------------------------- #
regionalplotpos <- function(x, cod, a=0, orient="xF", ...) {
# INPUT
# x = colonna
# cod = .
# ... = graphical parameters as xlab, ylab, main...
# OUTPUT
# rappresentazione in carta cartesiana della serie
X=x
cod <- factor(cod)
k <- nlevels(cod)
liv <- levels(cod)
x=c(min(x),max(x))
F=c(0,1)
if (orient=="xF") plot(x,F,type="n",...)
else if (orient=="Fx") plot(F,x,type="n",...)
else stop("regionalplotpos(x, cod, a, orient): orient unknown")
for(j in 1:k) {
ordinato <- sort(X[cod==liv[j]])
n <- length(ordinato)
i <- 1:n
plotpos <- (i - a)/(n + 1 - 2*a)
if (orient=="xF") points(ordinato, plotpos, pch=j, col=j)
else if (orient=="Fx") points(plotpos, ordinato, pch=j, col=j)
}
grid()
}
# --------------------------------------------------------------------------- #
regionalnormplot <- function(x, cod, a=0, orient="xF", ...) {
X=x
cod <- factor(cod)
k <- nlevels(cod)
liv <- levels(cod)
ni <- tapply(X,cod,length)
x=c(min(X),max(X))
F=c(qnorm((1-a)/(max(ni)+1-2*a)), qnorm((max(ni)-a)/(max(ni)+1-2*a)))
pp <- c(.001, .01, .05, .1, .25, .50, .75, .9, .95, .99, .999)
qq <- qnorm(pp)
if (orient=="xF") {
plot(x,F,type="n",axes=FALSE,...)
axis(1)
axis(2,at=qq,labels=pp)
for(j in 1:k) {
ordinato <- sort(X[cod==liv[j]])
n <- length(ordinato)
i <- 1:n
plotpos <- (i - a)/(n + 1 - 2*a)
uteorica <- qnorm(plotpos)
points(ordinato,uteorica,pch=j,col=j)
}
abline(h=qq, col="gray", lty="dotted")
grid(nx=NULL, ny=NA)
}
else if (orient=="Fx") {
plot(F,x,type="n",axes=FALSE,...)
axis(2)
axis(1,at=qq,labels=pp)
for(j in 1:k) {
ordinato <- sort(X[cod==liv[j]])
n <- length(ordinato)
i <- 1:n
plotpos <- (i - a)/(n + 1 - 2*a)
uteorica <- qnorm(plotpos)
points(uteorica,ordinato,pch=j,col=j)
}
abline(v=qq, col="gray", lty="dotted")
grid(ny=NULL, nx=NA)
}
else stop("regionalnormplot(x, cod, a, orient): orient unknown")
box()
}
# --------------------------------------------------------------------------- #
regionallognormplot <- function(x, cod, a=0, orient="xF", ...) {
X=x
cod <- factor(cod)
k <- nlevels(cod)
liv <- levels(cod)
ni <- tapply(X,cod,length)
x=c(min(X),max(X))
F=c(qnorm((1-a)/(max(ni)+1-2*a)), qnorm((max(ni)-a)/(max(ni)+1-2*a)))
pp <- c(.001, .01, .05, .1, .25, .50, .75, .9, .95, .99, .999)
qq <- qnorm(pp)
if (orient=="xF") {
plot(x,F,log="x",type="n",axes=FALSE,...)
axis(1)
axis(2,at=qq,labels=pp)
for(j in 1:k) {
ordinato <- sort(X[cod==liv[j]])
n <- length(ordinato)
i <- 1:n
plotpos <- (i - a)/(n + 1 - 2*a)
uteorica <- qnorm(plotpos)
points(ordinato,uteorica,pch=j,col=j)
}
abline(h=qq, col="gray", lty="dotted")
grid(nx=NULL, ny=NA)
}
else if (orient=="Fx") {
plot(F,x,log="y",type="n",axes=FALSE,...)
axis(2)
axis(1,at=qq,labels=pp)
for(j in 1:k) {
ordinato <- sort(X[cod==liv[j]])
n <- length(ordinato)
i <- 1:n
plotpos <- (i - a)/(n + 1 - 2*a)
uteorica <- qnorm(plotpos)
points(uteorica,ordinato, pch=j,col=j)
}
abline(v=qq, col="gray", lty="dotted")
grid(ny=NULL, nx=NA)
}
else stop("regionallognormplot(x, cod, a, orient): orient unknown")
box()
}
# ---------------------------------------------------------------------- #
regionalgumbelplot <- function(x, cod, a=0, orient="xF", ...) {
X=x
cod <- factor(cod)
k <- nlevels(cod)
liv <- levels(cod)
ni <- tapply(X,cod,length)
x=c(min(X),max(X))
F=c(-log(-log((1-a)/(max(ni)+1-2*a))), -log(-log((max(ni)-a)/(max(ni)+1-2*a))))
pp <- c(.001, .01, .05, .1, .25, .50, .75, .9, .95, .99, .999)
qq <- -log(-log(pp))
if (orient=="xF") {
plot(x,F,type="n",axes=FALSE,...)
axis(1)
axis(2,at=qq,labels=pp)
for(j in 1:k) {
ordinato <- sort(X[cod==liv[j]])
n <- length(ordinato)
i <- 1:n
plotpos <- (i - a)/(n + 1 - 2*a)
uteorica <- -log(-log(plotpos))
points(ordinato,uteorica,pch=j,col=j)
}
abline(h=qq, col="gray", lty="dotted")
grid(nx=NULL, ny=NA)
}
else if (orient=="Fx") {
plot(F,x,type="n",axes=FALSE,...)
axis(2)
axis(1,at=qq,labels=pp)
for(j in 1:k) {
ordinato <- sort(X[cod==liv[j]])
n <- length(ordinato)
i <- 1:n
plotpos <- (i - a)/(n + 1 - 2*a)
uteorica <- -log(-log(plotpos))
points(uteorica,ordinato,pch=j,col=j)
}
abline(v=qq, col="gray", lty="dotted")
grid(ny=NULL, nx=NA)
}
else stop("regionalgumbelplot(x, cod, a, orient): orient unknown")
box()
}
# ---------------------------------------------------------------------- #
regionalfrechetplot <- function(x, cod, a=0, orient="xF", ...) {
X=x
cod <- factor(cod)
k <- nlevels(cod)
liv <- levels(cod)
ni <- tapply(X,cod,length)
x=c(min(X),max(X))
F=c(-log(-log((1-a)/(max(ni)+1-2*a))), -log(-log((max(ni)-a)/(max(ni)+1-2*a))))
pp <- c(.001, .01, .05, .1, .25, .50, .75, .9, .95, .99, .999)
qq <- -log(-log(pp))
if (orient=="xF") {
plot(x,F,log="x",type="n",axes=FALSE,...)
axis(1)
axis(2,at=qq,labels=pp)
for(j in 1:k) {
ordinato <- sort(X[cod==liv[j]])
n <- length(ordinato)
i <- 1:n
plotpos <- (i - a)/(n + 1 - 2*a)
uteorica <- -log(-log(plotpos))
points(ordinato,uteorica,pch=j,col=j)
}
abline(h=qq, col="gray", lty="dotted")
grid(nx=NULL, ny=NA)
}
else if (orient=="Fx") {
plot(F,x,log="y",type="n",axes=FALSE,...)
axis(2)
axis(1,at=qq,labels=pp)
for(j in 1:k) {
ordinato <- sort(X[cod==liv[j]])
n <- length(ordinato)
i <- 1:n
plotpos <- (i - a)/(n + 1 - 2*a)
uteorica <- -log(-log(plotpos))
points(uteorica,ordinato,pch=j,col=j)
}
abline(v=qq, col="gray", lty="dotted")
grid(ny=NULL, nx=NA)
}
else stop("regionalfrechetplot(x, cod, a, orient): orient unknown")
box()
}
# ---------------------------------------------------------------------- #
regionalexpplot <- function(x, cod, a=0, orient="xF", ...) {
X=x
cod <- factor(cod)
k <- nlevels(cod)
liv <- levels(cod)
ni <- tapply(X,cod,length)
x=c(min(X),max(X))
F=c(-log(1 - ((1-a)/(max(ni)+1-2*a))), -log(1 - ((max(ni)-a)/(max(ni)+1-2*a))))
pp <- c(.1, .25, .50, .75, .9, .95, .975, .99, .995, .999)
qq <- -log(1 - pp)
if (orient=="xF") {
plot(x,F,type="n",axes=FALSE,...)
axis(1)
axis(2,at=qq,labels=pp)
for(j in 1:k) {
ordinato <- sort(X[cod==liv[j]])
n <- length(ordinato)
i <- 1:n
plotpos <- (i - a)/(n + 1 - 2*a)
uteorica <- -log(1 - plotpos)
points(ordinato,uteorica,pch=j,col=j)
}
abline(h=qq, col="gray", lty="dotted")
grid(nx=NULL, ny=NA)
}
else if (orient=="Fx") {
plot(F,x,type="n",axes=FALSE,...)
axis(2)
axis(1,at=qq,labels=pp)
for(j in 1:k) {
ordinato <- sort(X[cod==liv[j]])
n <- length(ordinato)
i <- 1:n
plotpos <- (i - a)/(n + 1 - 2*a)
uteorica <- -log(1 - plotpos)
points(uteorica,ordinato,pch=j,col=j)
}
abline(v=qq, col="gray", lty="dotted")
grid(ny=NULL, nx=NA)
}
else stop("regionalexpplot(x, cod, a, orient): orient unknown")
box()
}
# ---------------------------------------------------------------------- #
regionalparetoplot <- function(x, cod, a=0, orient="xF", ...) {
X=x
cod <- factor(cod)
k <- nlevels(cod)
liv <- levels(cod)
ni <- tapply(X,cod,length)
x=c(min(X),max(X))
F=c(-log(1 - ((1-a)/(max(ni)+1-2*a))), -log(1 - ((max(ni)-a)/(max(ni)+1-2*a))))
pp <- c(.1, .25, .50, .75, .9, .95, .975, .99, .995, .999)
qq <- -log(1 - pp)
if (orient=="xF") {
plot(x,F,log="x",type="n",axes=FALSE,...)
axis(1)
axis(2,at=qq,labels=pp)
for(j in 1:k) {
ordinato <- sort(X[cod==liv[j]])
n <- length(ordinato)
i <- 1:n
plotpos <- (i - a)/(n + 1 - 2*a)
uteorica <- -log(1 - plotpos)
points(ordinato,uteorica,pch=j,col=j)
}
abline(h=qq, col="gray", lty="dotted")
grid(nx=NULL, ny=NA)
}
else if (orient=="Fx") {
plot(F,x,log="y",type="n",axes=FALSE,...)
axis(2)
axis(1,at=qq,labels=pp)
for(j in 1:k) {
ordinato <- sort(X[cod==liv[j]])
n <- length(ordinato)
i <- 1:n
plotpos <- (i - a)/(n + 1 - 2*a)
uteorica <- -log(1 - plotpos)
points(uteorica,ordinato,pch=j,col=j)
}
abline(v=qq, col="gray", lty="dotted")
grid(ny=NULL, nx=NA)
}
else stop("regionalparetoplot(x, cod, a, orient): orient unknown")
box()
}
# ---------------------------------------------------------------- #
plotposRP <- function(x, a=0, orient="xF", ...) {
ordinato <- sort(x)
n <- length(ordinato)
i <- 1:n
plotpos <- (i - a)/(n + 1 - 2*a)
x=ordinato
T=1/(1 - plotpos)
if (orient=="xF") plot(x,T, log="y", ...)
else if (orient=="Fx") plot(T,x, log="x", ...)
else stop("plotposRP(x, a, orient): orient unknown")
grid(equilogs=FALSE)
}
# ---------------------------------------------------------------- #
pointsposRP <- function(x, a=0, orient="xF", ...) {
ordinato <- sort(x)
n <- length(ordinato)
i <- 1:n
plotpos <- (i - a)/(n + 1 - 2*a)
x=ordinato
T=1/(1 - plotpos)
if (orient=="xF") points(x,T,...)
else if (orient=="Fx") points(T,x,...)
else stop("pointsposRP(x, a, orient): orient unknown")
}
# ---------------------------------------------------------------- #
plotposRPhist <- function(xcont, xhist=NA, infhist=NA, suphist=NA, nbans=NA, seuil=NA, col12=c(1,1), a=0, orient="xF", ...) {
# Method in Maidment pag 18.41-42:
# Jery R. Stedinger, Richard M. Vogel, and Efi Foufoula-Georgiou. Frequency analysis of extreme events. In David R. Maidment, editor, Hand-Book of Hydrology, chapter 18. McGraw-Hill Companies, international edition, 1993.
# see BayesianMCMC for the notation
# col12 = colors for systematic and historical points
if (all(is.na(c(xcont, xhist, infhist, suphist, seuil)))) {
stop("plotposRPhist(xcont, xhist, infhist, suphist, nbans, seuil, ...): no input data")
}
else if (all(is.na(c(xhist, infhist, suphist, seuil)))) {
# systematic data only
n <- length(xcont)
i <- n + 1 - rank(xcont, ties.method="first")
q <- ((i - a)/(n + 1 - 2*a))
T <- 1/q
x <- xcont
if (orient=="xF") plot(x,T, log="y", col=col12[1], ...)
else if (orient=="Fx") plot(T,x, log="x", col=col12[1], ...)
else stop("plotposRPhist(xcont, xhist, infhist, suphist, nbans, seuil, col12, a, orient=, ...): orient unknown")
grid(equilogs=FALSE)
invisible(cbind(T,x))
}
else if (all(is.na(c(infhist, suphist))) & all(!is.na(c(xhist, seuil, nbans)))) {
# historical data known
if (!all(is.na(xcont))) {
xcont_u <- xcont[xcont >= seuil]
xcont_l <- xcont[xcont < seuil]
} else {
xcont_u <- NULL
xcont_l <- NULL
}
xhist_u <- xhist[xhist >= seuil]
xhist_l <- xhist[xhist < seuil]
x_u <- c(xcont_u, xhist_u)
col_x_u <- rep(col12, c(length(xcont_u), length(xhist_u)))
x_l <- c(xcont_l, xhist_l)
col_x_l <- rep(col12, c(length(xcont_l), length(xhist_l)))
r <- length(x_u) # number of times the threshold is exceeded
se <- length(x_l) # number of times the threshold is not exceeded by station-record data
n <- length(xcont) + nbans # length of the all period (historic + recent)
qe <- r/n # exceedence probability of the perception threshold
i <- r + 1 - rank(x_u, ties.method="first")
qi <- qe*((i - a)/(r + 1 - 2*a)) # eq. 18.6.11 in Maidment
T_u <- 1/qi # return periods above threshold
j <- se + 1 - rank(x_l, ties.method="first")
qj <- qe + (1 - qe)*((j - a)/(se + 1 - 2*a)) # eq. 18.6.12 in Maidment
T_l <- 1/qj # return periods below threshold
x <- c(x_l, x_u)
T <- c(T_l, T_u)
colori <- c(col_x_l, col_x_u)
if (orient=="xF") {
plot(x,T, log="y", col=colori, ...)
abline(v=seuil, lty=3)
}
else if (orient=="Fx") {
plot(T,x, log="x", col=colori, ...)
abline(h=seuil, lty=3)
}
else stop("plotposRPhist(xcont, xhist, infhist, suphist, nbans, seuil, col12, a, orient=, ...): orient unknown")
grid(equilogs=FALSE)
invisible(cbind(T,x))
}
else if (all(is.na(c(xhist))) & all(!is.na(c(infhist, suphist, seuil, nbans)))) {
# just ranges for the historicals
xhist <- (infhist + suphist)/2
if (!all(is.na(xcont))) {
xcont_u <- xcont[xcont >= seuil]
xcont_l <- xcont[xcont < seuil]
} else {
xcont_u <- NULL
xcont_l <- NULL
}
xhist_u <- xhist[xhist >= seuil]
xhist_l <- xhist[xhist < seuil]
x_u <- c(xcont_u, xhist_u)
col_x_u <- rep(c(col12[1],NA), c(length(xcont_u), length(xhist_u)))
x_l <- c(xcont_l, xhist_l)
col_x_l <- rep(c(col12[1],NA), c(length(xcont_l), length(xhist_l)))
r <- length(x_u) # number of times the threshold is exceeded
se <- length(x_l) # number of times the threshold is not exceeded
n <- length(xcont) + nbans # length of the all period (historic + recent)
qe <- r/n # exceedence probability of the perception threshold
i <- r + 1 - rank(x_u, ties.method="first")
qi <- qe*((i - a)/(r + 1 - 2*a)) # eq. 18.6.11 in Maidment
T_u <- 1/qi # return periods above threshold
j <- se + 1 - rank(x_l, ties.method="first")
qj <- qe + (1 - qe)*((j - a)/(se + 1 - 2*a)) # eq. 18.6.12 in Maidment
T_l <- 1/qj # return periods below threshold
x <- c(x_l, x_u)
T <- c(T_l, T_u)
Thist <- sort(T[x %in% xhist])[order(xhist)]
colori <- c(col_x_l, col_x_u)
if (orient=="xF") {
plot(x,T, log="y", col=colori, ...)
abline(v=seuil, lty=3)
points(infhist, Thist, pch=24, bg=col12[2], col=col12[2])
points(suphist, Thist, pch=25, bg=col12[2], col=col12[2])
segments(infhist, Thist, suphist, Thist, lty=3, col=col12[2])
}
else if (orient=="Fx") {
plot(T,x, log="x", col=colori, ...)
abline(h=seuil, lty=3)
points(Thist, infhist, pch=24, bg=col12[2], col=col12[2])
points(Thist, suphist, pch=25, bg=col12[2], col=col12[2])
segments(Thist, infhist, Thist, suphist, lty=3, col=col12[2])
}
else stop("plotposRPhist(xcont, xhist, infhist, suphist, nbans, seuil, col12, a, orient=, ...): orient unknown")
grid(equilogs=FALSE)
invisible(cbind(T,x))
}
}
# ---------------------------------------------------------------- #
pointsposRPhist <- function(xcont, xhist=NA, infhist=NA, suphist=NA, nbans=NA, seuil=NA, col12=c(1,1), a=0, orient="xF", ...) {
# Method in Maidment pag 18.41-42:
# Jery R. Stedinger, Richard M. Vogel, and Efi Foufoula-Georgiou. Frequency analysis of extreme events. In David R. Maidment, editor, Hand-Book of Hydrology, chapter 18. McGraw-Hill Companies, international edition, 1993.
# see BayesianMCMC for the notation
# col12 = colors for systematic and historical points
if (all(is.na(c(xcont, xhist, infhist, suphist, seuil)))) {
stop("plotposRPhist(xcont, xhist, infhist, suphist, nbans, seuil, ...): no input data")
}
else if (all(is.na(c(xhist, infhist, suphist, seuil)))) {
# systematic data only
n <- length(xcont)
i <- n + 1 - rank(xcont, ties.method="first")
q <- ((i - a)/(n + 1 - 2*a))
T <- 1/q
x <- xcont
if (orient=="xF") points(x,T, col=col12[1], ...)
else if (orient=="Fx") points(T,x, col=col12[1], ...)
else stop("plotposRPhist(xcont, xhist, infhist, suphist, nbans, seuil, col12, a, orient=, ...): orient unknown")
grid(equilogs=FALSE)
invisible(cbind(T,x))
}
else if (all(is.na(c(infhist, suphist))) & all(!is.na(c(xhist, seuil, nbans)))) {
# historical data known
if (!all(is.na(xcont))) {
xcont_u <- xcont[xcont >= seuil]
xcont_l <- xcont[xcont < seuil]
} else {
xcont_u <- NULL
xcont_l <- NULL
}
xhist_u <- xhist[xhist >= seuil]
xhist_l <- xhist[xhist < seuil]
x_u <- c(xcont_u, xhist_u)
col_x_u <- rep(col12, c(length(xcont_u), length(xhist_u)))
x_l <- c(xcont_l, xhist_l)
col_x_l <- rep(col12, c(length(xcont_l), length(xhist_l)))
r <- length(x_u) # number of times the threshold is exceeded
se <- length(x_l) # number of times the threshold is not exceeded
n <- length(xcont) + nbans # length of the all period (historic + recent)
qe <- r/n # exceedence probability of the perception threshold
i <- r + 1 - rank(x_u, ties.method="first")
qi <- qe*((i - a)/(r + 1 - 2*a)) # eq. 18.6.11 in Maidment
T_u <- 1/qi # return periods above threshold
j <- se + 1 - rank(x_l, ties.method="first")
qj <- qe + (1 - qe)*((j - a)/(se + 1 - 2*a)) # eq. 18.6.12 in Maidment
T_l <- 1/qj # return periods below threshold
x <- c(x_l, x_u)
T <- c(T_l, T_u)
colori <- c(col_x_l, col_x_u)
if (orient=="xF") {
points(x,T, col=colori, ...)
abline(v=seuil, lty=3)
}
else if (orient=="Fx") {
points(T,x, col=colori, ...)
abline(h=seuil, lty=3)
}
else stop("plotposRPhist(xcont, xhist, infhist, suphist, nbans, seuil, col12, a, orient=, ...): orient unknown")
grid(equilogs=FALSE)
invisible(cbind(T,x))
}
else if (all(is.na(c(xhist))) & all(!is.na(c(infhist, suphist, seuil, nbans)))) {
# just ranges for the historicals
xhist <- (infhist + suphist)/2
if (!all(is.na(xcont))) {
xcont_u <- xcont[xcont >= seuil]
xcont_l <- xcont[xcont < seuil]
} else {
xcont_u <- NULL
xcont_l <- NULL
}
xhist_u <- xhist[xhist >= seuil]
xhist_l <- xhist[xhist < seuil]
x_u <- c(xcont_u, xhist_u)
col_x_u <- rep(c(col12[1],NA), c(length(xcont_u), length(xhist_u)))
x_l <- c(xcont_l, xhist_l)
col_x_l <- rep(c(col12[1],NA), c(length(xcont_l), length(xhist_l)))
r <- length(x_u) # number of times the threshold is exceeded
se <- length(x_l) # number of times the threshold is not exceeded
n <- length(xcont) + nbans # length of the all period (historic + recent)
qe <- r/n # exceedence probability of the perception threshold
i <- r + 1 - rank(x_u, ties.method="first")
qi <- qe*((i - a)/(r + 1 - 2*a)) # eq. 18.6.11 in Maidment
T_u <- 1/qi # return periods above threshold
j <- se + 1 - rank(x_l, ties.method="first")
qj <- qe + (1 - qe)*((j - a)/(se + 1 - 2*a)) # eq. 18.6.12 in Maidment
T_l <- 1/qj # return periods below threshold
x <- c(x_l, x_u)
T <- c(T_l, T_u)
Thist <- T[x %in% xhist]
colori <- c(col_x_l, col_x_u)
if (orient=="xF") {
points(x,T, col=colori, ...)
abline(v=seuil, lty=3)
points(infhist, Thist, pch=24, bg=col12[2], col=col12[2])
points(suphist, Thist, pch=25, bg=col12[2], col=col12[2])
segments(infhist, Thist, suphist, Thist, lty=3, col=col12[2])
}
else if (orient=="Fx") {
points(T,x, col=colori, ...)
abline(h=seuil, lty=3)
points(Thist, infhist, pch=24, bg=col12[2], col=col12[2])
points(Thist, suphist, pch=25, bg=col12[2], col=col12[2])
segments(Thist, infhist, Thist, suphist, lty=3, col=col12[2])
}
else stop("plotposRPhist(xcont, xhist, infhist, suphist, nbans, seuil, col12, a, orient=, ...): orient unknown")
grid(equilogs=FALSE)
invisible(cbind(T,x))
}
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.