[Robast-commits] r1076 - branches/robast-1.1/pkg/RobExtremes/R branches/robast-1.1/pkg/RobExtremes/inst/scripts branches/robast-1.2/pkg/RobExtremes/R branches/robast-1.2/pkg/RobExtremes/inst/scripts pkg/RobExtremes/R pkg/RobExtremes/inst/scripts
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Sun Jul 29 20:06:43 CEST 2018
Author: ruckdeschel
Date: 2018-07-29 20:06:43 +0200 (Sun, 29 Jul 2018)
New Revision: 1076
Modified:
branches/robast-1.1/pkg/RobExtremes/R/internal-getpsi.R
branches/robast-1.1/pkg/RobExtremes/inst/scripts/RobFitsAtRealData.R
branches/robast-1.2/pkg/RobExtremes/R/internal-getpsi.R
branches/robast-1.2/pkg/RobExtremes/inst/scripts/RobFitsAtRealData.R
pkg/RobExtremes/R/internal-getpsi.R
pkg/RobExtremes/inst/scripts/RobFitsAtRealData.R
Log:
[RobExtremes] trunk & branch 1.1 & branch 1.2
internal-getpsi now for MBRE uses the same LM A inside and outside weight (obtained by averaging);
and a centering zi in the weight such that A %*% zi = centering outside
updated Example script; now produces legends as well; inserted many calls to devNew() ...
Modified: branches/robast-1.1/pkg/RobExtremes/R/internal-getpsi.R
===================================================================
--- branches/robast-1.1/pkg/RobExtremes/R/internal-getpsi.R 2018-07-29 15:58:19 UTC (rev 1075)
+++ branches/robast-1.1/pkg/RobExtremes/R/internal-getpsi.R 2018-07-29 18:06:43 UTC (rev 1076)
@@ -18,15 +18,24 @@
.dbeta <- diag(c(beta,1)); .dbeta1 <- diag(c(1/beta,1))
b <- fct[[1]](xi)
- a <- c(.dbeta%*%c(fct[[2]](xi),fct[[3]](xi)))
- aw <- c(.dbeta1%*%c(fct[[4]](xi),fct[[5]](xi)))
+
+ aa <- c(fct[[2]](xi),fct[[3]](xi))
+ zi <- c(fct[[4]](xi),fct[[5]](xi))
am <- mean(c(fct[[7]](xi),fct[[8]](xi)))
- A <- .dbeta%*%matrix(c(fct[[6]](xi),am,am,fct[[9]](xi)),2,2)%*%.dbeta
+ Aa <- matrix(c(fct[[6]](xi),am,am,fct[[9]](xi)),2,2)
am <- mean(c(fct[[11]](xi),fct[[12]](xi)))
- Aw <- matrix(c(fct[[10]](xi),am,am,fct[[13]](xi)),2,2)%*%.dbeta
+ Ai <- matrix(c(fct[[10]](xi),am,am,fct[[13]](xi)),2,2)
+ if(type==".MBRE"){
+ ai <- Ai %*% zi
+ Am <- (Ai+Aa)/2; Ai <- Aa <- Am
+ am <- (ai+aa)/2; ai <- aa <- am
+ zi <- solve(Ai,ai)
+ }
+ a <- c(.dbeta%*%aa)
+ aw <- c(.dbeta1%*%zi)
+ A <- .dbeta%*%Aa%*%.dbeta
+ Aw <- Ai%*%.dbeta
-
-
normt <- NormType()
biast <- symmetricBias()
nb <- ContNeighborhood(radius=0.5)
@@ -73,16 +82,26 @@
.dbeta <- diag(c(beta,beta,1)); .dbeta1 <- diag(c(1/beta,1/beta,1))
b <- fct[[1]](xi)
- a <- c(.dbeta%*%c(fct[[2]](xi),fct[[3]](xi),fct[[4]](xi)))
- aw <- c(.dbeta1%*%c(fct[[5]](xi),fct[[6]](xi),fct[[7]](xi)))
+ aa <- c(fct[[2]](xi),fct[[3]](xi),fct[[4]](xi))
+ zi <- c(fct[[5]](xi),fct[[6]](xi),fct[[7]](xi))
am1 <- mean(c(fct[[9]](xi),fct[[11]](xi)))
am2 <- mean(c(fct[[10]](xi),fct[[14]](xi)))
am3 <- mean(c(fct[[13]](xi),fct[[15]](xi)))
- A <- .dbeta%*%matrix(c(fct[[8]](xi),am1,am2,am1,fct[[12]](xi),am3,am2,am3,fct[[16]](xi)),3,3)%*%.dbeta
+ Aa <- matrix(c(fct[[8]](xi),am1,am2,am1,fct[[12]](xi),am3,am2,am3,fct[[16]](xi)),3,3)
am1 <- mean(c(fct[[18]](xi),fct[[20]](xi)))
am2 <- mean(c(fct[[19]](xi),fct[[23]](xi)))
am3 <- mean(c(fct[[22]](xi),fct[[24]](xi)))
- Aw <- matrix(c(fct[[17]](xi),am1,am2,am1,fct[[21]](xi),am3,am2,am3,fct[[25]](xi)),3,3)%*%.dbeta
+ Ai <- matrix(c(fct[[8]](xi),am1,am2,am1,fct[[17]](xi),am3,am2,am3,fct[[25]](xi)),3,3)
+ if(type==".MBRE"){
+ ai <- Ai %*% zi
+ Am <- (Ai+Aa)/2; Ai <- Aa <- Am
+ am <- (ai+aa)/2; ai <- aa <- am
+ zi <- solve(Ai,ai)
+ }
+ a <- c(.dbeta%*%aa)
+ aw <- c(.dbeta1%*%zi)
+ A <- .dbeta%*%Aa%*%.dbeta
+ Aw <- Ai%*%.dbeta
normt <- NormType()
biast <- symmetricBias()
Modified: branches/robast-1.1/pkg/RobExtremes/inst/scripts/RobFitsAtRealData.R
===================================================================
--- branches/robast-1.1/pkg/RobExtremes/inst/scripts/RobFitsAtRealData.R 2018-07-29 15:58:19 UTC (rev 1075)
+++ branches/robast-1.1/pkg/RobExtremes/inst/scripts/RobFitsAtRealData.R 2018-07-29 18:06:43 UTC (rev 1076)
@@ -26,6 +26,8 @@
grbsi <- groundbeef$serving
grbsc <- c(grbsi,10000,1000)
+
+options("newDevice"=TRUE) ## opens a new plot for each graphic
##--------------------------
## GEV
##--------------------------
@@ -67,21 +69,56 @@
attr(MBRi, "timings")
## our return values can be plugged into ismev-diagnostics:
+devNew()
gev.diag(mlEi)
+devNew()
gev.diag(MBRi)
+devNew()
gev.prof(mlEi, m = 10, 4.1, 5)
+devNew()
gev.profxi(MBRi, -0.3, 0.3)
## this is how the MLEi IC looks like
+devNew()
plot(pIC(mlEi10ALE),portpiriei)
## with additional arguments
+devNew()
plot(pIC(mlEi10ALE),portpiriei,
- cex.lbs=1.8,with.lab=TRUE, which.Order=1:3, c
- ol.pts="red", adj.lbs=c(1.3,1.3), pch.pts=20,pch.npts=21,
+ cex.lbs=1.8,with.lab=TRUE, which.Order=1:3, col.pts="red",
+ adj.lbs=c(1.3,1.3), pch.pts=20,pch.npts=21,
cex.pts=2,cex.npts=.6,lab.pts=c("H","G","I"))
## this is how the MBRi IC looks like
-plot(pIC(MBRi))
+devNew()
+plot(pIC(MBRi),portpiriei,
+ cex.lbs=1.8,with.lab=TRUE, which.Order=1:3, col.pts="red",
+ adj.lbs=c(1.3,1.3), pch.pts=20,pch.npts=21,
+ cex.pts=2,cex.npts=.6,lab.pts=c("H","G","I"))
+## compare the 4 ICs
+devNew()
+comparePlot(obj1=pIC(mlEi10ALE),obj2=pIC(MBRi),obj3=pIC(RMXi),obj4=pIC(mlEiALE),
+ portpiriei, with.legend=TRUE,# lwd=1, lty=1:4,
+ cex.lbs=1.8,with.lab=TRUE, which.Order=1:3, col.pts="red",
+ adj.lbs=c(1.3,1.3), pch.pts=20,pch.npts=21,
+ cex.pts=2,cex.npts=.6,lab.pts=c("H","G","I"))
+
+devNew()
+comparePlot(obj1=pIC(mlEi10ALE),obj2=pIC(MBRi),obj3=pIC(RMXi),
+ with.legend=TRUE, ylim=matrix(c(-1.5,1.5,-1,3,-5,5),2,3))
+
+devNew()
+infoPlot(pIC(MBRi), portpiriei, with.legend=TRUE, lty=1:4, log="y",
+ cex.lbs=1.8,with.lab=TRUE, which.Order=1:3, col.pts="red",
+ adj.lbs=c(1.3,1.3), pch.pts=20,pch.npts=21,
+ cex.pts=2,cex.npts=.6,lab.pts=c("H","G","I"))
+
+devNew()
+infoPlot(pIC(RMXi), portpiriei, with.legend=TRUE, lty=1:4, log="y",
+ cex.lbs=1.8,with.lab=TRUE, which.Order=1:3, col.pts="red",
+ adj.lbs=c(1.3,1.3), pch.pts=20,pch.npts=21,
+ cex.pts=2,cex.npts=.6,lab.pts=c("H","G","I"))
+
## rescale y axes according to different probit scales
+devNew()
plot(pIC(mlEi10ALE),portpiriei, with.lab=TRUE,
which.lbs=c(1:10), attr.pre=FALSE, col.pts="blue", col.npts="green",
scaleY = TRUE, scaleY.fct=list(function(x)pnorm(x,sd=4),
@@ -90,6 +127,7 @@
function(x)qnorm(x,sd=10),function(x)qnorm(x,sd=30)),scaleN=20)
## contaminated:
ppfitc <- ismev::gev.fit(portpiriec)
+devNew()
gev.diag(ppfitc)
##
mlEc <- MLEstimator(portpiriec, GEVFamilyMuUnknown(withPos=FALSE))
@@ -112,18 +150,27 @@
estimate(RMXc)
## diagnostics from ismev applied to our estimators:
+devNew()
gev.diag(mlEc)
+devNew()
gev.diag(MBRc)
+devNew()
gev.prof(mlEc, m = 10, 4.1, 5)
+devNew()
gev.profxi(mlEc, -0.3, 0.3)
## diagnostics from pkg 'distrMod'/'RobAStBase'
+devNew()
qqplot(portpiriec,MBRc)
+devNew()
qqplot(portpiriec,MBRc,ylim=c(3.5,5))
+devNew()
returnlevelplot(portpiriec,MBRc)
+devNew()
returnlevelplot(portpiriec,MBRc,ylim=c(3.5,5))
## here the MBR-IC looks as follows
+devNew()
plot(pIC(MBRc),portpiriec,log="x",xlim=c(1,1.5*max(portpiriec)))
##--------------------------
@@ -141,10 +188,15 @@
MBR2i
estimate(mlE2i)
estimate(MBR2i)
+devNew()
plot(MBR2i at pIC)
+devNew()
gpd.diag(mlE2i)
+devNew()
gpd.diag(MBR2i)
+devNew()
gpd.prof(mlE2i, m = 10, 55, 77)
+devNew()
gpd.profxi(MBR2i, -0.02, 0.02)
GP <- GParetoFamily(scale=rnfiti$mle[1],shape=rnfiti$mle[2],loc=10)
returnlevelplot(rain, GP, MaxOrPOT="POT", xlim=c(1e-1,1e3))
@@ -152,9 +204,11 @@
## contaminated:
rnfitc <- ismev::gpd.fit(c(rain,1000,10000),10)
+devNew()
gpd.diag(rnfitc)
##
mlE2c <- MLEstimator(rainc, GParetoFamily(loc=10))
+devNew()
gpd.diag(mlE2c)
system.time(MBR2c <- roptest(rainc, GParetoFamily(loc=10),risk=MBRRisk()))
system.time(RMX2c <- roptest(rainc, GParetoFamily(loc=10),risk=RMXRRisk()))
@@ -166,39 +220,61 @@
estimate(mlE2c)
estimate(MBR2c)
estimate(RMX2c)
+devNew()
gpd.diag(mlE2c)
+devNew()
gpd.diag(MBR2c)
+devNew()
gpd.diag(RMX2c)
+devNew()
gpd.prof(mlE2c, m = 10, 55, 77)
+devNew()
gpd.profxi(mlE2c, -0.02, 0.02)
+devNew()
plot(pIC(MBR2c))
+devNew()
qqplot(rainc,MBR2c)
+devNew()
qqplot(rainc,MBR2c,ylim=c(5,100))
+devNew()
qqplot(rainc,MBR2c,xlim=c(5,100),ylim=c(5,100),log="xy")
-qqplot(rainc,di,xlim=c(5,100),ylim=c(5,100),log="xy",
+devNew()
+qqplot(rainc,MBR2c,xlim=c(5,100),ylim=c(5,100),log="xy",
cex.pts=2,col.pts="blue",with.lab=TRUE,cex.lbs=.9,which.Order=1:3)
+devNew()
returnlevelplot(raini,MBR2i,MaxOrPot="POT",threshold=0)
+devNew()
returnlevelplot(raini,MBR2i,MaxOrPot="POT",threshold=0, withLab=TRUE, cex.lbl=0.8)
+devNew()
returnlevelplot(rainc,MBR2c,MaxOrPot="POT",threshold=0)
+devNew()
returnlevelplot(rainc,MBR2c,ylim=c(10,100),MaxOrPot="POT",threshold=0)
#
L2F <- eval(MBR2c at pIC@CallL2Fam)
dI2c <- L2F at distribution
-#loc(dI2c) <- 0
+devNew()
qqplot(rainc,dI2c)
rainc.10 <- rainc-10
+devNew()
qqplot(rainc.10,dI2c-10)
-## to be fixed
+devNew()
returnlevelplot(rainc.10,dI2c-10,MaxOrPot="POT",threshold=0)
+
## wrong data set
dI2i <- distribution(eval(MBR2i at pIC@CallL2Fam))
loc(dI2i) <- 0
+devNew()
qqplot(portpiriei-10,dI2i)
+devNew()
qqplot(portpiriec,MBR2c)
+### all points are red
+
## right data set
+devNew()
qqplot(raini-10,dI2i)
+devNew()
qqplot(rainc,MBR2c)
@@ -211,8 +287,11 @@
PM <- ParetoFamily(Min=2)
mlE3i <- MLEstimator(x,PM)
mlE3c <- MLEstimator(xc,PM)
+devNew()
qqplot(x, mlE3i, log="xy")
+devNew()
qqplot(xc, mlE3c, log="xy")
+devNew()
returnlevelplot(x, mlE3i, MaxOrPOT="POT",ylim=c(1,1e5),log="y")
system.time(MBR3i <- roptest(x, PM,risk=MBRRisk()))
@@ -225,7 +304,9 @@
estimate(mlE3c)
estimate(MBR3c)
estimate(RMX3c)
+devNew()
plot(pIC(MBR3i))
+devNew()
plot(pIC(RMX3i))
#######################################################
@@ -248,9 +329,13 @@
estimate(MBR4c)
estimate(OMS4c)
estimate(RMX4c)
+devNew()
plot(pIC(MBR4i))
+devNew()
plot(pIC(RMX4i))
+devNew()
qqplot(grbsi, RMX4i)
+devNew()
qqplot(grbsc, RMX4c, log="xy")
#######################################################
@@ -261,16 +346,26 @@
system.time(mlE5i <- MLEstimator(grbsi, GF))
system.time(OMS5i <- roptest(grbsi, GF,risk=OMSRRisk()))
system.time(RMX5i <- roptest(grbsi, GF,risk=RMXRRisk()))
+system.time(MBR5i <- roptest(grbsi, GF,risk=MBRRisk()))
system.time(mlE5c <- MLEstimator(grbsc, GF))
system.time(OMS5c <- roptest(grbsc, GF,risk=OMSRRisk()))
system.time(RMX5c <- roptest(grbsc, GF,risk=RMXRRisk()))
+system.time(MBR5c <- roptest(grbsc, GF,risk=MBRRisk()))
estimate(mlE5i)
estimate(RMX5i)
estimate(OMS5i)
+estimate(MBR5i)
estimate(mlE5c)
estimate(OMS5c)
estimate(RMX5c)
+estimate(MBR5c)
+devNew()
plot(pIC(OMS5i))
+devNew()
plot(pIC(RMX5i))
+devNew()
+plot(pIC(MBR5i))
+devNew()
qqplot(grbsi, RMX5i)
+devNew()
qqplot(grbsc, RMX5c, log="xy")
Modified: branches/robast-1.2/pkg/RobExtremes/R/internal-getpsi.R
===================================================================
--- branches/robast-1.2/pkg/RobExtremes/R/internal-getpsi.R 2018-07-29 15:58:19 UTC (rev 1075)
+++ branches/robast-1.2/pkg/RobExtremes/R/internal-getpsi.R 2018-07-29 18:06:43 UTC (rev 1076)
@@ -18,15 +18,24 @@
.dbeta <- diag(c(beta,1)); .dbeta1 <- diag(c(1/beta,1))
b <- fct[[1]](xi)
- a <- c(.dbeta%*%c(fct[[2]](xi),fct[[3]](xi)))
- aw <- c(.dbeta1%*%c(fct[[4]](xi),fct[[5]](xi)))
+
+ aa <- c(fct[[2]](xi),fct[[3]](xi))
+ zi <- c(fct[[4]](xi),fct[[5]](xi))
am <- mean(c(fct[[7]](xi),fct[[8]](xi)))
- A <- .dbeta%*%matrix(c(fct[[6]](xi),am,am,fct[[9]](xi)),2,2)%*%.dbeta
+ Aa <- matrix(c(fct[[6]](xi),am,am,fct[[9]](xi)),2,2)
am <- mean(c(fct[[11]](xi),fct[[12]](xi)))
- Aw <- matrix(c(fct[[10]](xi),am,am,fct[[13]](xi)),2,2)%*%.dbeta
+ Ai <- matrix(c(fct[[10]](xi),am,am,fct[[13]](xi)),2,2)
+ if(type==".MBRE"){
+ ai <- Ai %*% zi
+ Am <- (Ai+Aa)/2; Ai <- Aa <- Am
+ am <- (ai+aa)/2; ai <- aa <- am
+ zi <- solve(Ai,ai)
+ }
+ a <- c(.dbeta%*%aa)
+ aw <- c(.dbeta1%*%zi)
+ A <- .dbeta%*%Aa%*%.dbeta
+ Aw <- Ai%*%.dbeta
-
-
normt <- NormType()
biast <- symmetricBias()
nb <- ContNeighborhood(radius=0.5)
@@ -73,16 +82,26 @@
.dbeta <- diag(c(beta,beta,1)); .dbeta1 <- diag(c(1/beta,1/beta,1))
b <- fct[[1]](xi)
- a <- c(.dbeta%*%c(fct[[2]](xi),fct[[3]](xi),fct[[4]](xi)))
- aw <- c(.dbeta1%*%c(fct[[5]](xi),fct[[6]](xi),fct[[7]](xi)))
+ aa <- c(fct[[2]](xi),fct[[3]](xi),fct[[4]](xi))
+ zi <- c(fct[[5]](xi),fct[[6]](xi),fct[[7]](xi))
am1 <- mean(c(fct[[9]](xi),fct[[11]](xi)))
am2 <- mean(c(fct[[10]](xi),fct[[14]](xi)))
am3 <- mean(c(fct[[13]](xi),fct[[15]](xi)))
- A <- .dbeta%*%matrix(c(fct[[8]](xi),am1,am2,am1,fct[[12]](xi),am3,am2,am3,fct[[16]](xi)),3,3)%*%.dbeta
+ Aa <- matrix(c(fct[[8]](xi),am1,am2,am1,fct[[12]](xi),am3,am2,am3,fct[[16]](xi)),3,3)
am1 <- mean(c(fct[[18]](xi),fct[[20]](xi)))
am2 <- mean(c(fct[[19]](xi),fct[[23]](xi)))
am3 <- mean(c(fct[[22]](xi),fct[[24]](xi)))
- Aw <- matrix(c(fct[[17]](xi),am1,am2,am1,fct[[21]](xi),am3,am2,am3,fct[[25]](xi)),3,3)%*%.dbeta
+ Ai <- matrix(c(fct[[8]](xi),am1,am2,am1,fct[[17]](xi),am3,am2,am3,fct[[25]](xi)),3,3)
+ if(type==".MBRE"){
+ ai <- Ai %*% zi
+ Am <- (Ai+Aa)/2; Ai <- Aa <- Am
+ am <- (ai+aa)/2; ai <- aa <- am
+ zi <- solve(Ai,ai)
+ }
+ a <- c(.dbeta%*%aa)
+ aw <- c(.dbeta1%*%zi)
+ A <- .dbeta%*%Aa%*%.dbeta
+ Aw <- Ai%*%.dbeta
normt <- NormType()
biast <- symmetricBias()
Modified: branches/robast-1.2/pkg/RobExtremes/inst/scripts/RobFitsAtRealData.R
===================================================================
--- branches/robast-1.2/pkg/RobExtremes/inst/scripts/RobFitsAtRealData.R 2018-07-29 15:58:19 UTC (rev 1075)
+++ branches/robast-1.2/pkg/RobExtremes/inst/scripts/RobFitsAtRealData.R 2018-07-29 18:06:43 UTC (rev 1076)
@@ -26,6 +26,8 @@
grbsi <- groundbeef$serving
grbsc <- c(grbsi,10000,1000)
+
+options("newDevice"=TRUE) ## opens a new plot for each graphic
##--------------------------
## GEV
##--------------------------
@@ -67,21 +69,56 @@
attr(MBRi, "timings")
## our return values can be plugged into ismev-diagnostics:
+devNew()
gev.diag(mlEi)
+devNew()
gev.diag(MBRi)
+devNew()
gev.prof(mlEi, m = 10, 4.1, 5)
+devNew()
gev.profxi(MBRi, -0.3, 0.3)
## this is how the MLEi IC looks like
+devNew()
plot(pIC(mlEi10ALE),portpiriei)
## with additional arguments
+devNew()
plot(pIC(mlEi10ALE),portpiriei,
- cex.lbs=1.8,with.lab=TRUE, which.Order=1:3, c
- ol.pts="red", adj.lbs=c(1.3,1.3), pch.pts=20,pch.npts=21,
+ cex.lbs=1.8,with.lab=TRUE, which.Order=1:3, col.pts="red",
+ adj.lbs=c(1.3,1.3), pch.pts=20,pch.npts=21,
cex.pts=2,cex.npts=.6,lab.pts=c("H","G","I"))
## this is how the MBRi IC looks like
-plot(pIC(MBRi))
+devNew()
+plot(pIC(MBRi),portpiriei,
+ cex.lbs=1.8,with.lab=TRUE, which.Order=1:3, col.pts="red",
+ adj.lbs=c(1.3,1.3), pch.pts=20,pch.npts=21,
+ cex.pts=2,cex.npts=.6,lab.pts=c("H","G","I"))
+## compare the 4 ICs
+devNew()
+comparePlot(obj1=pIC(mlEi10ALE),obj2=pIC(MBRi),obj3=pIC(RMXi),obj4=pIC(mlEiALE),
+ portpiriei, with.legend=TRUE,# lwd=1, lty=1:4,
+ cex.lbs=1.8,with.lab=TRUE, which.Order=1:3, col.pts="red",
+ adj.lbs=c(1.3,1.3), pch.pts=20,pch.npts=21,
+ cex.pts=2,cex.npts=.6,lab.pts=c("H","G","I"))
+
+devNew()
+comparePlot(obj1=pIC(mlEi10ALE),obj2=pIC(MBRi),obj3=pIC(RMXi),
+ with.legend=TRUE, ylim=matrix(c(-1.5,1.5,-1,3,-5,5),2,3))
+
+devNew()
+infoPlot(pIC(MBRi), portpiriei, with.legend=TRUE, lty=1:4, log="y",
+ cex.lbs=1.8,with.lab=TRUE, which.Order=1:3, col.pts="red",
+ adj.lbs=c(1.3,1.3), pch.pts=20,pch.npts=21,
+ cex.pts=2,cex.npts=.6,lab.pts=c("H","G","I"))
+
+devNew()
+infoPlot(pIC(RMXi), portpiriei, with.legend=TRUE, lty=1:4, log="y",
+ cex.lbs=1.8,with.lab=TRUE, which.Order=1:3, col.pts="red",
+ adj.lbs=c(1.3,1.3), pch.pts=20,pch.npts=21,
+ cex.pts=2,cex.npts=.6,lab.pts=c("H","G","I"))
+
## rescale y axes according to different probit scales
+devNew()
plot(pIC(mlEi10ALE),portpiriei, with.lab=TRUE,
which.lbs=c(1:10), attr.pre=FALSE, col.pts="blue", col.npts="green",
scaleY = TRUE, scaleY.fct=list(function(x)pnorm(x,sd=4),
@@ -90,6 +127,7 @@
function(x)qnorm(x,sd=10),function(x)qnorm(x,sd=30)),scaleN=20)
## contaminated:
ppfitc <- ismev::gev.fit(portpiriec)
+devNew()
gev.diag(ppfitc)
##
mlEc <- MLEstimator(portpiriec, GEVFamilyMuUnknown(withPos=FALSE))
@@ -112,18 +150,27 @@
estimate(RMXc)
## diagnostics from ismev applied to our estimators:
+devNew()
gev.diag(mlEc)
+devNew()
gev.diag(MBRc)
+devNew()
gev.prof(mlEc, m = 10, 4.1, 5)
+devNew()
gev.profxi(mlEc, -0.3, 0.3)
## diagnostics from pkg 'distrMod'/'RobAStBase'
+devNew()
qqplot(portpiriec,MBRc)
+devNew()
qqplot(portpiriec,MBRc,ylim=c(3.5,5))
+devNew()
returnlevelplot(portpiriec,MBRc)
+devNew()
returnlevelplot(portpiriec,MBRc,ylim=c(3.5,5))
## here the MBR-IC looks as follows
+devNew()
plot(pIC(MBRc),portpiriec,log="x",xlim=c(1,1.5*max(portpiriec)))
##--------------------------
@@ -141,10 +188,15 @@
MBR2i
estimate(mlE2i)
estimate(MBR2i)
+devNew()
plot(MBR2i at pIC)
+devNew()
gpd.diag(mlE2i)
+devNew()
gpd.diag(MBR2i)
+devNew()
gpd.prof(mlE2i, m = 10, 55, 77)
+devNew()
gpd.profxi(MBR2i, -0.02, 0.02)
GP <- GParetoFamily(scale=rnfiti$mle[1],shape=rnfiti$mle[2],loc=10)
returnlevelplot(rain, GP, MaxOrPOT="POT", xlim=c(1e-1,1e3))
@@ -152,9 +204,11 @@
## contaminated:
rnfitc <- ismev::gpd.fit(c(rain,1000,10000),10)
+devNew()
gpd.diag(rnfitc)
##
mlE2c <- MLEstimator(rainc, GParetoFamily(loc=10))
+devNew()
gpd.diag(mlE2c)
system.time(MBR2c <- roptest(rainc, GParetoFamily(loc=10),risk=MBRRisk()))
system.time(RMX2c <- roptest(rainc, GParetoFamily(loc=10),risk=RMXRRisk()))
@@ -166,39 +220,61 @@
estimate(mlE2c)
estimate(MBR2c)
estimate(RMX2c)
+devNew()
gpd.diag(mlE2c)
+devNew()
gpd.diag(MBR2c)
+devNew()
gpd.diag(RMX2c)
+devNew()
gpd.prof(mlE2c, m = 10, 55, 77)
+devNew()
gpd.profxi(mlE2c, -0.02, 0.02)
+devNew()
plot(pIC(MBR2c))
+devNew()
qqplot(rainc,MBR2c)
+devNew()
qqplot(rainc,MBR2c,ylim=c(5,100))
+devNew()
qqplot(rainc,MBR2c,xlim=c(5,100),ylim=c(5,100),log="xy")
-qqplot(rainc,di,xlim=c(5,100),ylim=c(5,100),log="xy",
+devNew()
+qqplot(rainc,MBR2c,xlim=c(5,100),ylim=c(5,100),log="xy",
cex.pts=2,col.pts="blue",with.lab=TRUE,cex.lbs=.9,which.Order=1:3)
+devNew()
returnlevelplot(raini,MBR2i,MaxOrPot="POT",threshold=0)
+devNew()
returnlevelplot(raini,MBR2i,MaxOrPot="POT",threshold=0, withLab=TRUE, cex.lbl=0.8)
+devNew()
returnlevelplot(rainc,MBR2c,MaxOrPot="POT",threshold=0)
+devNew()
returnlevelplot(rainc,MBR2c,ylim=c(10,100),MaxOrPot="POT",threshold=0)
#
L2F <- eval(MBR2c at pIC@CallL2Fam)
dI2c <- L2F at distribution
-#loc(dI2c) <- 0
+devNew()
qqplot(rainc,dI2c)
rainc.10 <- rainc-10
+devNew()
qqplot(rainc.10,dI2c-10)
-## to be fixed
+devNew()
returnlevelplot(rainc.10,dI2c-10,MaxOrPot="POT",threshold=0)
+
## wrong data set
dI2i <- distribution(eval(MBR2i at pIC@CallL2Fam))
loc(dI2i) <- 0
+devNew()
qqplot(portpiriei-10,dI2i)
+devNew()
qqplot(portpiriec,MBR2c)
+### all points are red
+
## right data set
+devNew()
qqplot(raini-10,dI2i)
+devNew()
qqplot(rainc,MBR2c)
@@ -211,8 +287,11 @@
PM <- ParetoFamily(Min=2)
mlE3i <- MLEstimator(x,PM)
mlE3c <- MLEstimator(xc,PM)
+devNew()
qqplot(x, mlE3i, log="xy")
+devNew()
qqplot(xc, mlE3c, log="xy")
+devNew()
returnlevelplot(x, mlE3i, MaxOrPOT="POT",ylim=c(1,1e5),log="y")
system.time(MBR3i <- roptest(x, PM,risk=MBRRisk()))
@@ -225,7 +304,9 @@
estimate(mlE3c)
estimate(MBR3c)
estimate(RMX3c)
+devNew()
plot(pIC(MBR3i))
+devNew()
plot(pIC(RMX3i))
#######################################################
@@ -248,9 +329,13 @@
estimate(MBR4c)
estimate(OMS4c)
estimate(RMX4c)
+devNew()
plot(pIC(MBR4i))
+devNew()
plot(pIC(RMX4i))
+devNew()
qqplot(grbsi, RMX4i)
+devNew()
qqplot(grbsc, RMX4c, log="xy")
#######################################################
@@ -261,16 +346,26 @@
system.time(mlE5i <- MLEstimator(grbsi, GF))
system.time(OMS5i <- roptest(grbsi, GF,risk=OMSRRisk()))
system.time(RMX5i <- roptest(grbsi, GF,risk=RMXRRisk()))
+system.time(MBR5i <- roptest(grbsi, GF,risk=MBRRisk()))
system.time(mlE5c <- MLEstimator(grbsc, GF))
system.time(OMS5c <- roptest(grbsc, GF,risk=OMSRRisk()))
system.time(RMX5c <- roptest(grbsc, GF,risk=RMXRRisk()))
+system.time(MBR5c <- roptest(grbsc, GF,risk=MBRRisk()))
estimate(mlE5i)
estimate(RMX5i)
estimate(OMS5i)
+estimate(MBR5i)
estimate(mlE5c)
estimate(OMS5c)
estimate(RMX5c)
+estimate(MBR5c)
+devNew()
plot(pIC(OMS5i))
+devNew()
plot(pIC(RMX5i))
+devNew()
+plot(pIC(MBR5i))
+devNew()
qqplot(grbsi, RMX5i)
+devNew()
qqplot(grbsc, RMX5c, log="xy")
Modified: pkg/RobExtremes/R/internal-getpsi.R
===================================================================
--- pkg/RobExtremes/R/internal-getpsi.R 2018-07-29 15:58:19 UTC (rev 1075)
+++ pkg/RobExtremes/R/internal-getpsi.R 2018-07-29 18:06:43 UTC (rev 1076)
@@ -18,15 +18,24 @@
.dbeta <- diag(c(beta,1)); .dbeta1 <- diag(c(1/beta,1))
b <- fct[[1]](xi)
- a <- c(.dbeta%*%c(fct[[2]](xi),fct[[3]](xi)))
- aw <- c(.dbeta1%*%c(fct[[4]](xi),fct[[5]](xi)))
+
+ aa <- c(fct[[2]](xi),fct[[3]](xi))
+ zi <- c(fct[[4]](xi),fct[[5]](xi))
am <- mean(c(fct[[7]](xi),fct[[8]](xi)))
- A <- .dbeta%*%matrix(c(fct[[6]](xi),am,am,fct[[9]](xi)),2,2)%*%.dbeta
+ Aa <- matrix(c(fct[[6]](xi),am,am,fct[[9]](xi)),2,2)
am <- mean(c(fct[[11]](xi),fct[[12]](xi)))
- Aw <- matrix(c(fct[[10]](xi),am,am,fct[[13]](xi)),2,2)%*%.dbeta
+ Ai <- matrix(c(fct[[10]](xi),am,am,fct[[13]](xi)),2,2)
+ if(type==".MBRE"){
+ ai <- Ai %*% zi
+ Am <- (Ai+Aa)/2; Ai <- Aa <- Am
+ am <- (ai+aa)/2; ai <- aa <- am
+ zi <- solve(Ai,ai)
+ }
+ a <- c(.dbeta%*%aa)
+ aw <- c(.dbeta1%*%zi)
+ A <- .dbeta%*%Aa%*%.dbeta
+ Aw <- Ai%*%.dbeta
-
-
normt <- NormType()
biast <- symmetricBias()
nb <- ContNeighborhood(radius=0.5)
@@ -73,16 +82,26 @@
.dbeta <- diag(c(beta,beta,1)); .dbeta1 <- diag(c(1/beta,1/beta,1))
b <- fct[[1]](xi)
- a <- c(.dbeta%*%c(fct[[2]](xi),fct[[3]](xi),fct[[4]](xi)))
- aw <- c(.dbeta1%*%c(fct[[5]](xi),fct[[6]](xi),fct[[7]](xi)))
+ aa <- c(fct[[2]](xi),fct[[3]](xi),fct[[4]](xi))
+ zi <- c(fct[[5]](xi),fct[[6]](xi),fct[[7]](xi))
am1 <- mean(c(fct[[9]](xi),fct[[11]](xi)))
am2 <- mean(c(fct[[10]](xi),fct[[14]](xi)))
am3 <- mean(c(fct[[13]](xi),fct[[15]](xi)))
- A <- .dbeta%*%matrix(c(fct[[8]](xi),am1,am2,am1,fct[[12]](xi),am3,am2,am3,fct[[16]](xi)),3,3)%*%.dbeta
+ Aa <- matrix(c(fct[[8]](xi),am1,am2,am1,fct[[12]](xi),am3,am2,am3,fct[[16]](xi)),3,3)
am1 <- mean(c(fct[[18]](xi),fct[[20]](xi)))
am2 <- mean(c(fct[[19]](xi),fct[[23]](xi)))
am3 <- mean(c(fct[[22]](xi),fct[[24]](xi)))
- Aw <- matrix(c(fct[[17]](xi),am1,am2,am1,fct[[21]](xi),am3,am2,am3,fct[[25]](xi)),3,3)%*%.dbeta
+ Ai <- matrix(c(fct[[8]](xi),am1,am2,am1,fct[[17]](xi),am3,am2,am3,fct[[25]](xi)),3,3)
+ if(type==".MBRE"){
+ ai <- Ai %*% zi
+ Am <- (Ai+Aa)/2; Ai <- Aa <- Am
+ am <- (ai+aa)/2; ai <- aa <- am
+ zi <- solve(Ai,ai)
+ }
+ a <- c(.dbeta%*%aa)
+ aw <- c(.dbeta1%*%zi)
+ A <- .dbeta%*%Aa%*%.dbeta
+ Aw <- Ai%*%.dbeta
normt <- NormType()
biast <- symmetricBias()
Modified: pkg/RobExtremes/inst/scripts/RobFitsAtRealData.R
===================================================================
--- pkg/RobExtremes/inst/scripts/RobFitsAtRealData.R 2018-07-29 15:58:19 UTC (rev 1075)
+++ pkg/RobExtremes/inst/scripts/RobFitsAtRealData.R 2018-07-29 18:06:43 UTC (rev 1076)
@@ -26,6 +26,8 @@
grbsi <- groundbeef$serving
grbsc <- c(grbsi,10000,1000)
+
+options("newDevice"=TRUE) ## opens a new plot for each graphic
##--------------------------
## GEV
##--------------------------
@@ -67,21 +69,56 @@
attr(MBRi, "timings")
## our return values can be plugged into ismev-diagnostics:
+devNew()
gev.diag(mlEi)
+devNew()
gev.diag(MBRi)
+devNew()
gev.prof(mlEi, m = 10, 4.1, 5)
+devNew()
gev.profxi(MBRi, -0.3, 0.3)
## this is how the MLEi IC looks like
+devNew()
plot(pIC(mlEi10ALE),portpiriei)
## with additional arguments
+devNew()
plot(pIC(mlEi10ALE),portpiriei,
- cex.lbs=1.8,with.lab=TRUE, which.Order=1:3, c
- ol.pts="red", adj.lbs=c(1.3,1.3), pch.pts=20,pch.npts=21,
+ cex.lbs=1.8,with.lab=TRUE, which.Order=1:3, col.pts="red",
+ adj.lbs=c(1.3,1.3), pch.pts=20,pch.npts=21,
cex.pts=2,cex.npts=.6,lab.pts=c("H","G","I"))
## this is how the MBRi IC looks like
-plot(pIC(MBRi))
+devNew()
+plot(pIC(MBRi),portpiriei,
+ cex.lbs=1.8,with.lab=TRUE, which.Order=1:3, col.pts="red",
+ adj.lbs=c(1.3,1.3), pch.pts=20,pch.npts=21,
+ cex.pts=2,cex.npts=.6,lab.pts=c("H","G","I"))
+## compare the 4 ICs
+devNew()
+comparePlot(obj1=pIC(mlEi10ALE),obj2=pIC(MBRi),obj3=pIC(RMXi),obj4=pIC(mlEiALE),
+ portpiriei, with.legend=TRUE,# lwd=1, lty=1:4,
+ cex.lbs=1.8,with.lab=TRUE, which.Order=1:3, col.pts="red",
+ adj.lbs=c(1.3,1.3), pch.pts=20,pch.npts=21,
+ cex.pts=2,cex.npts=.6,lab.pts=c("H","G","I"))
+
+devNew()
+comparePlot(obj1=pIC(mlEi10ALE),obj2=pIC(MBRi),obj3=pIC(RMXi),
+ with.legend=TRUE, ylim=matrix(c(-1.5,1.5,-1,3,-5,5),2,3))
+
+devNew()
+infoPlot(pIC(MBRi), portpiriei, with.legend=TRUE, lty=1:4, log="y",
+ cex.lbs=1.8,with.lab=TRUE, which.Order=1:3, col.pts="red",
+ adj.lbs=c(1.3,1.3), pch.pts=20,pch.npts=21,
+ cex.pts=2,cex.npts=.6,lab.pts=c("H","G","I"))
+
+devNew()
+infoPlot(pIC(RMXi), portpiriei, with.legend=TRUE, lty=1:4, log="y",
+ cex.lbs=1.8,with.lab=TRUE, which.Order=1:3, col.pts="red",
+ adj.lbs=c(1.3,1.3), pch.pts=20,pch.npts=21,
+ cex.pts=2,cex.npts=.6,lab.pts=c("H","G","I"))
+
## rescale y axes according to different probit scales
+devNew()
plot(pIC(mlEi10ALE),portpiriei, with.lab=TRUE,
which.lbs=c(1:10), attr.pre=FALSE, col.pts="blue", col.npts="green",
scaleY = TRUE, scaleY.fct=list(function(x)pnorm(x,sd=4),
@@ -90,6 +127,7 @@
function(x)qnorm(x,sd=10),function(x)qnorm(x,sd=30)),scaleN=20)
## contaminated:
ppfitc <- ismev::gev.fit(portpiriec)
+devNew()
gev.diag(ppfitc)
##
mlEc <- MLEstimator(portpiriec, GEVFamilyMuUnknown(withPos=FALSE))
@@ -112,18 +150,27 @@
estimate(RMXc)
## diagnostics from ismev applied to our estimators:
+devNew()
gev.diag(mlEc)
+devNew()
gev.diag(MBRc)
+devNew()
gev.prof(mlEc, m = 10, 4.1, 5)
+devNew()
gev.profxi(mlEc, -0.3, 0.3)
## diagnostics from pkg 'distrMod'/'RobAStBase'
+devNew()
qqplot(portpiriec,MBRc)
+devNew()
qqplot(portpiriec,MBRc,ylim=c(3.5,5))
+devNew()
returnlevelplot(portpiriec,MBRc)
+devNew()
returnlevelplot(portpiriec,MBRc,ylim=c(3.5,5))
## here the MBR-IC looks as follows
+devNew()
plot(pIC(MBRc),portpiriec,log="x",xlim=c(1,1.5*max(portpiriec)))
##--------------------------
@@ -141,10 +188,15 @@
MBR2i
estimate(mlE2i)
estimate(MBR2i)
+devNew()
plot(MBR2i at pIC)
+devNew()
gpd.diag(mlE2i)
+devNew()
gpd.diag(MBR2i)
+devNew()
gpd.prof(mlE2i, m = 10, 55, 77)
+devNew()
gpd.profxi(MBR2i, -0.02, 0.02)
GP <- GParetoFamily(scale=rnfiti$mle[1],shape=rnfiti$mle[2],loc=10)
returnlevelplot(rain, GP, MaxOrPOT="POT", xlim=c(1e-1,1e3))
@@ -152,9 +204,11 @@
## contaminated:
rnfitc <- ismev::gpd.fit(c(rain,1000,10000),10)
+devNew()
gpd.diag(rnfitc)
##
mlE2c <- MLEstimator(rainc, GParetoFamily(loc=10))
+devNew()
gpd.diag(mlE2c)
system.time(MBR2c <- roptest(rainc, GParetoFamily(loc=10),risk=MBRRisk()))
system.time(RMX2c <- roptest(rainc, GParetoFamily(loc=10),risk=RMXRRisk()))
@@ -166,39 +220,61 @@
estimate(mlE2c)
estimate(MBR2c)
estimate(RMX2c)
+devNew()
gpd.diag(mlE2c)
+devNew()
gpd.diag(MBR2c)
+devNew()
gpd.diag(RMX2c)
+devNew()
gpd.prof(mlE2c, m = 10, 55, 77)
+devNew()
gpd.profxi(mlE2c, -0.02, 0.02)
+devNew()
plot(pIC(MBR2c))
+devNew()
qqplot(rainc,MBR2c)
+devNew()
qqplot(rainc,MBR2c,ylim=c(5,100))
+devNew()
qqplot(rainc,MBR2c,xlim=c(5,100),ylim=c(5,100),log="xy")
-qqplot(rainc,di,xlim=c(5,100),ylim=c(5,100),log="xy",
+devNew()
+qqplot(rainc,MBR2c,xlim=c(5,100),ylim=c(5,100),log="xy",
cex.pts=2,col.pts="blue",with.lab=TRUE,cex.lbs=.9,which.Order=1:3)
+devNew()
returnlevelplot(raini,MBR2i,MaxOrPot="POT",threshold=0)
+devNew()
returnlevelplot(raini,MBR2i,MaxOrPot="POT",threshold=0, withLab=TRUE, cex.lbl=0.8)
+devNew()
returnlevelplot(rainc,MBR2c,MaxOrPot="POT",threshold=0)
+devNew()
returnlevelplot(rainc,MBR2c,ylim=c(10,100),MaxOrPot="POT",threshold=0)
#
L2F <- eval(MBR2c at pIC@CallL2Fam)
dI2c <- L2F at distribution
-#loc(dI2c) <- 0
+devNew()
qqplot(rainc,dI2c)
rainc.10 <- rainc-10
+devNew()
qqplot(rainc.10,dI2c-10)
-## to be fixed
+devNew()
returnlevelplot(rainc.10,dI2c-10,MaxOrPot="POT",threshold=0)
+
## wrong data set
dI2i <- distribution(eval(MBR2i at pIC@CallL2Fam))
loc(dI2i) <- 0
+devNew()
qqplot(portpiriei-10,dI2i)
+devNew()
qqplot(portpiriec,MBR2c)
+### all points are red
+
## right data set
+devNew()
qqplot(raini-10,dI2i)
+devNew()
qqplot(rainc,MBR2c)
@@ -211,8 +287,11 @@
PM <- ParetoFamily(Min=2)
mlE3i <- MLEstimator(x,PM)
mlE3c <- MLEstimator(xc,PM)
+devNew()
qqplot(x, mlE3i, log="xy")
+devNew()
qqplot(xc, mlE3c, log="xy")
+devNew()
returnlevelplot(x, mlE3i, MaxOrPOT="POT",ylim=c(1,1e5),log="y")
system.time(MBR3i <- roptest(x, PM,risk=MBRRisk()))
@@ -225,7 +304,9 @@
estimate(mlE3c)
estimate(MBR3c)
estimate(RMX3c)
+devNew()
plot(pIC(MBR3i))
+devNew()
plot(pIC(RMX3i))
#######################################################
@@ -248,9 +329,13 @@
estimate(MBR4c)
estimate(OMS4c)
estimate(RMX4c)
+devNew()
plot(pIC(MBR4i))
+devNew()
plot(pIC(RMX4i))
+devNew()
qqplot(grbsi, RMX4i)
+devNew()
qqplot(grbsc, RMX4c, log="xy")
#######################################################
@@ -261,16 +346,26 @@
system.time(mlE5i <- MLEstimator(grbsi, GF))
system.time(OMS5i <- roptest(grbsi, GF,risk=OMSRRisk()))
system.time(RMX5i <- roptest(grbsi, GF,risk=RMXRRisk()))
+system.time(MBR5i <- roptest(grbsi, GF,risk=MBRRisk()))
system.time(mlE5c <- MLEstimator(grbsc, GF))
system.time(OMS5c <- roptest(grbsc, GF,risk=OMSRRisk()))
system.time(RMX5c <- roptest(grbsc, GF,risk=RMXRRisk()))
+system.time(MBR5c <- roptest(grbsc, GF,risk=MBRRisk()))
estimate(mlE5i)
estimate(RMX5i)
estimate(OMS5i)
+estimate(MBR5i)
estimate(mlE5c)
estimate(OMS5c)
estimate(RMX5c)
+estimate(MBR5c)
+devNew()
plot(pIC(OMS5i))
+devNew()
plot(pIC(RMX5i))
+devNew()
+plot(pIC(MBR5i))
+devNew()
qqplot(grbsi, RMX5i)
+devNew()
qqplot(grbsc, RMX5c, log="xy")
More information about the Robast-commits
mailing list