BayesianMCMC_reg.R
In nsRFA: Non-Supervised Regional Frequency Analysis

### R code from vignette source 'BayesianMCMC_reg.Rnw'

###################################################
### code chunk number 1: BayesianMCMC_reg.Rnw:51-54
###################################################
library(nsRFA)
data(Ardechedata)
ls()


###################################################
### code chunk number 2: BayesianMCMC_reg.Rnw:59-60
###################################################
str(SaintMartin_cont)


###################################################
### code chunk number 3: BayesianMCMC_reg.Rnw:64-65
###################################################
str(SaintMartin_hist)


###################################################
### code chunk number 4: plot_SaintMartin_cont_hist (eval = FALSE)
###################################################
## plot(c(1640,2010), c(0, 8000), type="n", xlab="", ylab="discharge (m3/s)")
##  points(SaintMartin_cont[,"year"], SaintMartin_cont[,"peak"], type="o", pch=21, bg="white")
##  points(SaintMartin_hist$peaks[,"year"], SaintMartin_hist$peaks[,"peak"], pch=16)
##  points(SaintMartin_hist$peaks[,"year"], SaintMartin_hist$peaks[,"peak"], type="h")
##  segments(x0=SaintMartin_hist$thresholds[,"from.yr"], 
##           x1=SaintMartin_hist$thresholds[,"to.yr"],
##           y0=SaintMartin_hist$thresholds[,"threshold"], lty=2)


###################################################
### code chunk number 5: BayesianMCMC_reg.Rnw:79-81
###################################################
par(mar=c(3,3,2,1)+0.03, mgp=c(1.5,0.3,0), tcl=.2, xaxs="i", yaxs="i", las=0)
plot(c(1640,2010), c(0, 8000), type="n", xlab="", ylab="discharge (m3/s)")
 points(SaintMartin_cont[,"year"], SaintMartin_cont[,"peak"], type="o", pch=21, bg="white")
 points(SaintMartin_hist$peaks[,"year"], SaintMartin_hist$peaks[,"peak"], pch=16)
 points(SaintMartin_hist$peaks[,"year"], SaintMartin_hist$peaks[,"peak"], type="h")
 segments(x0=SaintMartin_hist$thresholds[,"from.yr"], 
          x1=SaintMartin_hist$thresholds[,"to.yr"],
          y0=SaintMartin_hist$thresholds[,"threshold"], lty=2)


###################################################
### code chunk number 6: BayesianMCMC_reg.Rnw:92-96
###################################################
xhistSM <- c(-1,SaintMartin_hist$peaks[,"peak"])
seuilSM <- c(7250, 6000, rep(5050, 10), rep(2400, 21))
nbansSM <- c(127, 55, 65, rep(0,9), 71, rep(0, 20))
data.frame(xhistSM, seuilSM, nbansSM)


###################################################
### code chunk number 7: BayesianMCMC_reg.Rnw:105-106
###################################################
set.seed(198)


###################################################
### code chunk number 8: BayesianMCMC_reg.Rnw:108-114
###################################################
fit <- BayesianMCMC(xcont=SaintMartin_cont[,"peak"], 
                    xhist=xhistSM, 
                    seuil=seuilSM, 
                    nbans=nbansSM, 
                    nbpas=1000, nbchaines=3, confint=c(0.05, 0.95), dist="GEV",
                    varparameters0=c(NA, NA, 0.5))


###################################################
### code chunk number 9: BayesianMCMC_reg.Rnw:120-121
###################################################
fit <- BayesianMCMCcont(fit)


###################################################
### code chunk number 10: BayesianMCMC_reg.Rnw:125-126 (eval = FALSE)
###################################################
## plot(fit, which=1:4, ask=TRUE)


###################################################
### code chunk number 11: plot_fit_1
###################################################
par(mar=c(3,3,2,1)+0.03, mgp=c(1.5,0.3,0), tcl=.2, xaxs="r", yaxs="r", las=0)
plot(fit, 1)


###################################################
### code chunk number 12: BayesianMCMC_reg.Rnw:145-147
###################################################
par(mar=c(3,3,2,1)+0.03, mgp=c(1.5,0.3,0), tcl=.2, xaxs="r", yaxs="r", las=0)
plot(fit, 2)


###################################################
### code chunk number 13: plot_fit_3
###################################################
par(mar=c(3,3,2,1)+0.03, mgp=c(1.5,0.3,0), tcl=.2, xaxs="r", yaxs="r", las=0)
plot(fit, 3)


###################################################
### code chunk number 14: plot_fit_4
###################################################
par(mar=c(3,3,2,1)+0.03, mgp=c(1.5,0.3,0), tcl=.2, xaxs="r", yaxs="r", las=0)
plot(fit, 4)


###################################################
### code chunk number 15: BayesianMCMC_reg.Rnw:184-189 (eval = FALSE)
###################################################
## SaintMartin_cont
## Vogue_cont
## SaintLaurent_cont
## Beauvene_cont
## Chambonas_cont


###################################################
### code chunk number 16: plot_Ardeche_cont (eval = FALSE)
###################################################
## plot(SaintMartin_cont[,"year"], SaintMartin_cont[,"peak"], type="o", pch=21, bg="white", 
##      xlab="", ylab="discharge (m3/s)", main=" Saint Martin on the Ardeche River (2240 km2)")
## plot(Vogue_cont[,"year"], Vogue_cont[,"peak"], type="o", pch=21, bg="white", 
##      xlab="", ylab="discharge (m3/s)", main="Vogue on the Ardeche River (636 km2)")
## plot(SaintLaurent_cont[,"year"], SaintLaurent_cont[,"peak"], type="o", pch=21, bg="white", 
##      xlab="", ylab="discharge (m3/s)", main="Saint Laurent on the Borne River (63 km2)")
## plot(Beauvene_cont[,"year"], Beauvene_cont[,"peak"], type="o", pch=21, bg="white", 
##      xlab="", ylab="discharge (m3/s)", main="Beauvene on the Eyrieux River (392 km2)")
## plot(Chambonas_cont[,"year"], Chambonas_cont[,"peak"], type="o", pch=21, bg="white", 
##      xlab="", ylab="discharge (m3/s)", main="Chambonas on the Chassezac River (507 km2)")


###################################################
### code chunk number 17: BayesianMCMC_reg.Rnw:204-207
###################################################
layout(matrix(1:6, ncol=2))
par(mar=c(3,3,2,1)+0.03, mgp=c(1.5,0.3,0), tcl=.2, xaxs="r", yaxs="r", las=0)
plot(SaintMartin_cont[,"year"], SaintMartin_cont[,"peak"], type="o", pch=21, bg="white", 
     xlab="", ylab="discharge (m3/s)", main=" Saint Martin on the Ardeche River (2240 km2)")
plot(Vogue_cont[,"year"], Vogue_cont[,"peak"], type="o", pch=21, bg="white", 
     xlab="", ylab="discharge (m3/s)", main="Vogue on the Ardeche River (636 km2)")
plot(SaintLaurent_cont[,"year"], SaintLaurent_cont[,"peak"], type="o", pch=21, bg="white", 
     xlab="", ylab="discharge (m3/s)", main="Saint Laurent on the Borne River (63 km2)")
plot(Beauvene_cont[,"year"], Beauvene_cont[,"peak"], type="o", pch=21, bg="white", 
     xlab="", ylab="discharge (m3/s)", main="Beauvene on the Eyrieux River (392 km2)")
plot(Chambonas_cont[,"year"], Chambonas_cont[,"peak"], type="o", pch=21, bg="white", 
     xlab="", ylab="discharge (m3/s)", main="Chambonas on the Chassezac River (507 km2)")


###################################################
### code chunk number 18: BayesianMCMC_reg.Rnw:211-213
###################################################
Ardeche_areas  # km2
Ardeche_ungauged_extremes  # m3/s


###################################################
### code chunk number 19: BayesianMCMC_reg.Rnw:216-217
###################################################
Ardeche_ungauged_extremes_new <- merge(Ardeche_ungauged_extremes, Ardeche_areas)


###################################################
### code chunk number 20: BayesianMCMC_reg.Rnw:232-246
###################################################
xcont <- c(SaintMartin_cont[,2],
           Vogue_cont[,2], 
           SaintLaurent_cont[,2], 
           Beauvene_cont[,2], Chambonas_cont[,2])
scont <- c(rep(2240, length(SaintMartin_cont[,2])), 
           rep(636, length(Vogue_cont[,2])), 
           rep(63, length(SaintLaurent_cont[,2])), 
           rep(392, length(Beauvene_cont[,2])), 
           rep(507, length(Chambonas_cont[,2])))                           
xhist <- Ardeche_ungauged_extremes_new[,"peak"]
xhist[Ardeche_ungauged_extremes_new[,"station"]=="Chambonas"] <- -1
shist <- Ardeche_ungauged_extremes_new[,"area"]
nbans <- Ardeche_ungauged_extremes_new[,"associated.period"]
seuil <- Ardeche_ungauged_extremes_new[,"peak"]


###################################################
### code chunk number 21: BayesianMCMC_reg.Rnw:250-251
###################################################
set.seed(198)


###################################################
### code chunk number 22: BayesianMCMC_reg.Rnw:253-263
###################################################
fitreg <- BayesianMCMCreg(xcont=xcont, 
                          scont=scont,
                          xhist=xhist,
                          shist=shist,
                          nbans=nbans,
                          seuil=seuil,
                          nbpas=1000, nbchaines=3, confint=c(0.05,0.95), dist="GEV",
                          varparameters0=c(NA, NA, NA, 0.5))
fitreg <- BayesianMCMCregcont(fitreg)
fitreg <- BayesianMCMCregcont(fitreg)


###################################################
### code chunk number 23: BayesianMCMC_reg.Rnw:281-282 (eval = FALSE)
###################################################
## plot(fitreg, which=1:4, ask=TRUE)


###################################################
### code chunk number 24: plot_fitreg_1
###################################################
par(mar=c(3,3,2,1)+0.03, mgp=c(1.5,0.3,0), tcl=.2, xaxs="r", yaxs="r", las=0)
plot(fitreg, 1)


###################################################
### code chunk number 25: BayesianMCMC_reg.Rnw:295-297
###################################################
par(mar=c(3,3,2,1)+0.03, mgp=c(1.5,0.3,0), tcl=.2, xaxs="r", yaxs="r", las=0)
plot(fitreg, 2)


###################################################
### code chunk number 26: plot_fitreg_3
###################################################
par(mar=c(3,3,2,1)+0.03, mgp=c(1.5,0.3,0), tcl=.2, xaxs="r", yaxs="r", las=0)
plot(fitreg, 3)


###################################################
### code chunk number 27: plot_fitreg_4
###################################################
par(mar=c(3,3,2,1)+0.03, mgp=c(1.5,0.3,0), tcl=.2, xaxs="r", yaxs="r", las=0)
plot(fitreg, 4)


###################################################
### code chunk number 28: plotBayesianMCMCreg_surfArdeche (eval = FALSE)
###################################################
## plotBayesianMCMCreg_surf(fitreg, surf=unique(Ardeche_areas$area))


###################################################
### code chunk number 29: plotBayesianMCMCreg_surfArdechegrosso
###################################################
layout(matrix(1:20, ncol=4, byrow=TRUE))
par(mar=c(3,3,2,1)+0.03, mgp=c(1.5,0.3,0), tcl=.2, xaxs="r", yaxs="r", las=0)
plotBayesianMCMCreg_surf(fitreg, surf=unique(Ardeche_areas$area))