[Robast-commits] r1029 - in branches/robast-1.1/pkg/RobExtremes: R inst inst/scripts
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Mon Jul 23 21:47:56 CEST 2018
Author: ruckdeschel
Date: 2018-07-23 21:47:56 +0200 (Mon, 23 Jul 2018)
New Revision: 1029
Modified:
branches/robast-1.1/pkg/RobExtremes/R/Expectation.R
branches/robast-1.1/pkg/RobExtremes/inst/NEWS
branches/robast-1.1/pkg/RobExtremes/inst/scripts/RobFitsAtRealData.R
Log:
[RobExtremes] branch 1.1
+ Expectation E for DistributionsIntegratingByQuantiles to speed up things
now uses stop.on.error = FALSE and for accuracy splits up the integration range in [0,0.98] and [0.98, upp] as soon as upp>0.98
+ updated script RobFitsAtRealData.R
Modified: branches/robast-1.1/pkg/RobExtremes/R/Expectation.R
===================================================================
--- branches/robast-1.1/pkg/RobExtremes/R/Expectation.R 2018-07-23 19:44:10 UTC (rev 1028)
+++ branches/robast-1.1/pkg/RobExtremes/R/Expectation.R 2018-07-23 19:47:56 UTC (rev 1029)
@@ -57,7 +57,10 @@
dots.withoutUseApply <- dots
useApply <- TRUE
if(!is.null(dots$useApply)) useApply <- dots$useApply
+
dots.withoutUseApply$useApply <- NULL
+ dots.withoutUseApply$stop.on.error <- NULL
+
integrand <- function(x, dfun, ...){ di <- dim(x)
y <- q.l(object)(x)##quantile transformation
if(useApply){
@@ -76,12 +79,29 @@
upp <- p(object)(Ib["upp"])
if(is.nan(low)) low <- 0
if(is.nan(upp)) upp <- 1
- return(do.call(distrExIntegrate, c(list(f = integrand,
+
+ if(upp < 0.98){
+ int <- do.call(distrExIntegrate, c(list(f = integrand,
lower = low,
upper = upp,
- rel.tol = rel.tol,
- distr = object, dfun = dunif), dots.withoutUseApply)))
+ rel.tol = rel.tol, stop.on.error = FALSE,
+ distr = object, dfun = dunif), dots.withoutUseApply))
+ }else{
+ int1 <- do.call(distrExIntegrate, c(list(f = integrand,
+ lower = low,
+ upper = 0.98,
+ rel.tol = rel.tol, stop.on.error = FALSE,
+ distr = object, dfun = dunif), dots.withoutUseApply))
+ int2 <- do.call(distrExIntegrate, c(list(f = integrand,
+ lower = 0.98,
+ upper = upp,
+ rel.tol = rel.tol, stop.on.error = FALSE,
+ distr = object, dfun = dunif), dots.withoutUseApply))
+ int <- int1+int2
+ }
+ return(int)
+
})
setMethod("E", signature(object = "GPareto",
Modified: branches/robast-1.1/pkg/RobExtremes/inst/NEWS
===================================================================
--- branches/robast-1.1/pkg/RobExtremes/inst/NEWS 2018-07-23 19:44:10 UTC (rev 1028)
+++ branches/robast-1.1/pkg/RobExtremes/inst/NEWS 2018-07-23 19:47:56 UTC (rev 1029)
@@ -42,6 +42,8 @@
+ in ParetoFamily.R L2derivDistr was not attached to return value
under the hood:
++ Expectation E for DistributionsIntegratingByQuantiles to speed up things
+ now uses stop.on.error = FALSE and for accuracy splits up the integration range in [0,0.98] and [0.98, upp] as soon as upp>0.98
+ prepared everything for first release on CRAN
+ due to speed gains moved back diagnostic plot examples out of \donttest
+ added MakeIC calls in getStartIC to enhance accuracy
Modified: branches/robast-1.1/pkg/RobExtremes/inst/scripts/RobFitsAtRealData.R
===================================================================
--- branches/robast-1.1/pkg/RobExtremes/inst/scripts/RobFitsAtRealData.R 2018-07-23 19:44:10 UTC (rev 1028)
+++ branches/robast-1.1/pkg/RobExtremes/inst/scripts/RobFitsAtRealData.R 2018-07-23 19:47:56 UTC (rev 1029)
@@ -21,16 +21,27 @@
##
mlEi <- MLEstimator(portpiriei, GEVFamilyMuUnknown(withPos=FALSE))
system.time(MBRi <- roptest(portpiriei, GEVFamilyMuUnknown(withPos=FALSE),risk=MBRRisk()))
+system.time(RMXi <- roptest(portpiriei, GEVFamilyMuUnknown(withPos=FALSE),risk=RMXRRisk()))
+## in fact the precision of the pIC is not too good, but the resp. estimate only differs
+## little to the situation where we enforce IC conditions
+checkIC(pIC(RMXi))
+system.time(RMXiw <- roptest(portpiriei, GEVFamilyMuUnknown(withPos=FALSE),risk=RMXRRisk(),withMakeIC=TRUE))
+checkIC(pIC(RMXiw))
+estimate(RMXi)
+estimate(RMXiw)
+
mlEi
MBRi
estimate(mlEi)
estimate(MBRi)
+estimate(RMXi)
+estimate(RMXiw)
attr(MBRi, "timings")
gev.diag(mlEi)
gev.diag(MBRi)
gev.prof(mlEi, m = 10, 4.1, 5)
gev.profxi(MBRi, -0.3, 0.3)
-plot(MBRi at pIC)
+plot(pIC(MBRi))
## contaminated:
ppfitc <- ismev::gev.fit(portpiriec)
@@ -38,12 +49,15 @@
##
mlEc <- MLEstimator(portpiriec, GEVFamilyMuUnknown(withPos=FALSE))
system.time(MBRc <- roptest(portpiriec, GEVFamilyMuUnknown(withPos=FALSE),risk=MBRRisk()))
+system.time(RMXc <- roptest(portpiriec, GEVFamilyMuUnknown(withPos=FALSE),risk=RMXRRisk()))
mlEc
MBRc
+estimate(mlEi)
+estimate(MBRi)
+estimate(RMXi)
estimate(mlEc)
estimate(MBRc)
-estimate(mlEi)
-estimate(MBRi)
+estimate(RMXc)
gev.diag(mlEc)
gev.diag(MBRc)
gev.prof(mlEc, m = 10, 4.1, 5)
@@ -53,7 +67,7 @@
returnlevelplot(portpiriec,MBRc)
returnlevelplot(portpiriec,MBRc,ylim=c(3.5,5))
-plot(MBRc at pIC)
+plot(pIC(MBRc))
##--------------------------
## GPD
@@ -65,6 +79,7 @@
mlE2i <- MLEstimator(raini, GParetoFamily(loc=10))
gpd.diag(mlE2i)
system.time(MBR2i <- roptest(raini, GParetoFamily(loc=10),risk=MBRRisk()))
+system.time(RMX2i <- roptest(raini, GParetoFamily(loc=10),risk=RMXRRisk()))
mlE2i
MBR2i
estimate(mlE2i)
@@ -85,23 +100,29 @@
mlE2c <- MLEstimator(rainc, GParetoFamily(loc=10))
gpd.diag(mlE2c)
system.time(MBR2c <- roptest(rainc, GParetoFamily(loc=10),risk=MBRRisk()))
+system.time(RMX2c <- roptest(rainc, GParetoFamily(loc=10),risk=RMXRRisk()))
mlE2c
MBR2c
+estimate(mlE2i)
+estimate(MBR2i)
+estimate(RMX2i)
estimate(mlE2c)
estimate(MBR2c)
-estimate(mlE2i)
-estimate(MBR2i)
+estimate(RMX2c)
gpd.diag(mlE2c)
gpd.diag(MBR2c)
+gpd.diag(RMX2c)
gpd.prof(mlE2c, m = 10, 55, 77)
gpd.profxi(mlE2c, -0.02, 0.02)
-plot(MBR2c at pIC)
+plot(pIC(MBR2c))
qqplot(rainc,MBR2c)
qqplot(rainc,MBR2c,ylim=c(5,100))
qqplot(rainc,MBR2c,xlim=c(5,100),ylim=c(5,100),log="xy")
+qqplot(rainc,MBR2c,xlim=c(5,100),ylim=c(5,100),log="xy",cex.pch=8,exp.fadcol.pch = .55,pch=19)
returnlevelplot(raini,MBR2i,MaxOrPot="POT",threshold=0)
+returnlevelplot(raini,MBR2i,MaxOrPot="POT",threshold=0, withLab=TRUE, cex.lbl=0.8)
returnlevelplot(rainc,MBR2c,MaxOrPot="POT",threshold=0)
returnlevelplot(rainc,MBR2c,ylim=c(10,100),MaxOrPot="POT",threshold=0)
#
@@ -118,3 +139,28 @@
## wrong data set
qqplot(portpiriei-10,dI2i)
qqplot(portpiriec,MBR2c)
+
+
+#######################################################
+# synthetic Pareto data
+#######################################################
+set.seed(20180723)
+x <- actuar::rpareto1(100, min=2, shape=0.4)
+xc <- c(x, 1e15,1e12,1e69, 2.001,2.00000001)
+PM <- ParetoFamily(Min=2)
+mlE3i <- MLEstimator(x,PM)
+mlE3c <- MLEstimator(xc,PM)
+qqplot(x, mlE3i, log="xy")
+qqplot(xc, mlE3c, log="xy")
+system.time(MBR3i <- roptest(x, PM,risk=MBRRisk()))
+system.time(RMX3i <- roptest(x, PM,risk=RMXRRisk()))
+system.time(MBR3c <- roptest(xc, PM,risk=MBRRisk()))
+system.time(RMX3c <- roptest(xc, PM,risk=RMXRRisk()))
+estimate(mlE3i)
+estimate(MBR3i)
+estimate(RMX3i)
+estimate(mlE3c)
+estimate(MBR3c)
+estimate(RMX3c)
+plot(pIC(MBR3i))
+plot(pIC(RMX3i))
More information about the Robast-commits
mailing list