[Robast-commits] r748 - in branches/robast-1.0/pkg/RobExtremes: R man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Sat Apr 12 15:05:32 CEST 2014
Author: ruckdeschel
Date: 2014-04-12 15:05:31 +0200 (Sat, 12 Apr 2014)
New Revision: 748
Added:
branches/robast-1.0/pkg/RobExtremes/R/startEstGPD.R
Modified:
branches/robast-1.0/pkg/RobExtremes/R/GParetoFamily.R
branches/robast-1.0/pkg/RobExtremes/R/startEstGEV.R
branches/robast-1.0/pkg/RobExtremes/man/GParetoFamily.Rd
Log:
RobExtremes: fixed issue with check.validity as notified by B. Spangl, and extended starting estimator in GParetoFamily
Modified: branches/robast-1.0/pkg/RobExtremes/R/GParetoFamily.R
===================================================================
--- branches/robast-1.0/pkg/RobExtremes/R/GParetoFamily.R 2014-04-09 21:20:47 UTC (rev 747)
+++ branches/robast-1.0/pkg/RobExtremes/R/GParetoFamily.R 2014-04-12 13:05:31 UTC (rev 748)
@@ -41,6 +41,7 @@
of.interest = c("scale", "shape"),
p = NULL, N = NULL, trafo = NULL,
start0Est = NULL, withPos = TRUE,
+ secLevel = 0.7,
withCentL2 = FALSE,
withL2derivDistr = FALSE,
..ignoreTrafo = FALSE){
@@ -127,14 +128,19 @@
## starting parameters
startPar <- function(x,...){
tr <- theta[1]
-
+ n <- length(x)
+ epsn <- min(floor(secLevel*sqrt(n))+1,n)
+
## Pickand estimator
if(is.null(start0Est)){
PF <- GParetoFamily(loc = theta[1],
scale = theta[2], shape = theta[3])
- e1 <- medkMADhybr(c(x), k=10, ParamFamily = PF,
- q.lo = 1e-3, q.up = 15)
- e0 <- estimate(e1)
+ e1 <- try(
+ medkMADhybr(c(x), k=10, ParamFamily = PF,
+ q.lo = 1e-3, q.up = 15), silent =TRUE)
+ if(is(e1,"try-error")){ e0 <- .getBetaXiGPD(x=x, mu=tr,
+ xiGrid=.getXiGrid(), withPos=withPos)
+ }else e0 <- estimate(e1)
}else{
if(is(start0Est,"function")){
e1 <- start0Est(x, ...)
@@ -144,9 +150,9 @@
e0 <- e0[c("scale", "shape")]
}
- if(any(x < tr-.Machine$double.eps))
+ if(quantile(e0[2]*(x-tr), epsn/n)<.Machine$double.eps)
stop("some data smaller than 'loc' ")
- if(e0[2]<0) if(any(x > tr-e0[1]/e0[2]))
+ if(e0[2]<0) if(quantile(x,1-epsn/n) > tr-e0[1]/e0[2])
stop("shape is negative and some data larger than 'loc-scale/shape' ")
# if(any(x < tr-e0["scale"]/e0["shape"]))
# stop("some data smaller than 'loc-scale/shape' ")
Modified: branches/robast-1.0/pkg/RobExtremes/R/startEstGEV.R
===================================================================
--- branches/robast-1.0/pkg/RobExtremes/R/startEstGEV.R 2014-04-09 21:20:47 UTC (rev 747)
+++ branches/robast-1.0/pkg/RobExtremes/R/startEstGEV.R 2014-04-12 13:05:31 UTC (rev 748)
@@ -1,7 +1,8 @@
.getXiGrid <- function(){c(0.1,seq(-0.48,5,by=0.5))}
-.getBetaXiGEV <- function(x, mu, xiGrid = .getXiGrid(), withPos=TRUE, secLevel = 0.7){
+.getBetaXiGEV <- function(x, mu, xiGrid = .getXiGrid(), withPos=TRUE, secLevel = 0.7,
+ .issueIntermediateParams = FALSE){
n <- length(x)
epsn <- min(floor(secLevel*sqrt(n))+1,n)
@@ -9,8 +10,47 @@
x0 <- x-mu
s0 <- max(x0)-min(x0)
crit0 <- Inf
+
fu <- function(x,...) .getBetaXiGEV(x,mu,xiGrid = xiGrid,withPos=withPos)
+ e0 <- NULL
+
+ ### first try (to ensure global consistency): PickandsEstimator
+ try({mygev <- GEVFamily(loc=0,scale=1,shape=0.1, withPos=withPos,
+ ..withWarningGEV=FALSE)
+ e1 <- PickandsEstimator(x,ParamFamily=mygev)
+ if(.issueIntermediateParams){
+ cat("Pickands:\n");print(e1) }
+ e0 <- estimate(e1)}, silent=TRUE)
+
+ if(!is.null(e0)) if(!is(e0,"try-error")){
+ mygev <- GEVFamily(loc=0,scale=e0[1],shape=e0[2], withPos=withPos,
+ start0Est = fu, ..withWarningGEV=FALSE)
+ mde0 <- try(MDEstimator(x0, mygev, distance=CvMDist, startPar=c("scale"=es0[1],"shape"=es0[2])),silent=TRUE)
+ es0 <- c(NA,NA)
+ if(!is(mde0,"try-error")){
+ es <- estimate(mde0)
+ crit1 <- criterion(mde0)
+ if(.issueIntermediateParams){
+ cat("1st candidate:\n", round(es,6), " crit:", round(crit1,6), , " ")
+ }
+ if(quantile(1+es[2]*x0/es[1], epsn/n)>0){
+ mdeb <- mde0
+ crit0 <- crit1
+ es0 <- es
+ if(.issueIntermediateParams){
+ cat("side condition '1+sc/sh (x-mu) >0' fulfilled;\n")
+ }
+ }else{
+ if(.issueIntermediateParams){
+ cat("side condition '1+sc/sh (x-mu) >0' violated;\n")
+ }
+ }
+ }
+ }
+
+ i <- 0
for(xi in xiGrid){
+ i <- i + 1
funl <- function(sig){
mygev1 <- GEV(loc=0,scale=sig,shape=xi)
CvMDist(x0,mygev1)
@@ -25,16 +65,28 @@
if(!is(mde0,"try-error")){
es <- estimate(mde0)
crit1 <- criterion(mde0)
+ if(.issueIntermediateParams){
+ cat("candidate no",i+1, ":\n", round(es,6), " crit:", round(crit1,6), " ")
+ }
if(quantile(1+es[2]*x0/es[1], epsn/n)>0){
+ if(.issueIntermediateParams){
+ cat("side condition '1+sc/sh (x-mu) >0' fulfilled;\n")
+ }
if(crit1<crit0){
mdeb <- mde0
crit0 <- crit1
es0 <- es
}
}else{
+ if(.issueIntermediateParams){
+ cat("side condition '1+sc/sh (x-mu) >0' violated;\n")
+ }
es[1] <- intlo+1e-5
mygev2 <- GEV(loc=0,scale=es[1],shape=es[2])
crit1 <- CvMDist(x0,mygev2)
+ if(.issueIntermediateParams){
+ cat("candidate no",i+1, "(b):\n", round(es,6), " crit:", round(crit1,6), " ")
+ }
if(crit1<crit0){
crit0 <- crit1
es0 <- es
@@ -46,9 +98,12 @@
return(es)
}
-.getMuBetaXiGEV <- function(x, xiGrid = .getXiGrid(), withPos=TRUE, secLevel = 0.7){
+.getMuBetaXiGEV <- function(x, xiGrid = .getXiGrid(), withPos=TRUE, secLevel = 0.7,
+ .issueIntermediateParams = FALSE){
mu <- quantile(x,exp(-1))
- es <- .getBetaXiGEV(x=x, mu=mu, xiGrid=xiGrid, withPos=withPos, secLevel = secLevel)
+ es <- .getBetaXiGEV(x=x, mu=mu, xiGrid=xiGrid, withPos=withPos,
+ secLevel = secLevel,
+ .issueIntermediateParams = .issueIntermediateParams)
es0 <- c(mu,es)
names(es0) <- c("loc","scale","shape")
return(es0)
Added: branches/robast-1.0/pkg/RobExtremes/R/startEstGPD.R
===================================================================
--- branches/robast-1.0/pkg/RobExtremes/R/startEstGPD.R (rev 0)
+++ branches/robast-1.0/pkg/RobExtremes/R/startEstGPD.R 2014-04-12 13:05:31 UTC (rev 748)
@@ -0,0 +1,106 @@
+.getBetaXiGPD <- function(x, mu, xiGrid = .getXiGrid(), withPos=TRUE, secLevel = 0.7,
+ .issueIntermediateParams = FALSE){
+
+ n <- length(x)
+ epsn <- min(floor(secLevel*sqrt(n))+1,n)
+
+ x0 <- x-mu
+ s0 <- max(x0)-min(x0)
+ crit0 <- Inf
+
+ fu <- function(x,...) .getBetaXiGPD(x,mu,xiGrid = xiGrid,withPos=withPos)
+ e0 <- NULL
+
+ ### first try (to ensure global consistency): PickandsEstimator
+ try({mygpd <- GParetoFamily(loc=0,scale=1,shape=0.1, withPos=withPos)
+ e1 <- medkMADhybr(x,ParamFamily=mygpd, k=10)
+ if(.issueIntermediateParams){
+ cat("MedkMAD:\n");print(e1)}
+ e0 <- estimate(e1)}, silent=TRUE)
+
+ if(!is.null(e0)) if(!is(e0,"try-error")){
+ mygpd <- GParetoFamily(loc=0,scale=e0[1],shape=e0[2], withPos=withPos,
+ start0Est = fu)
+ mde0 <- try(MDEstimator(x0, mygpd, distance=CvMDist, startPar=c("scale"=es0[1],"shape"=es0[2])),silent=TRUE)
+ es0 <- c(NA,NA)
+ if(!is(mde0,"try-error")){
+ es <- estimate(mde0)
+ crit1 <- criterion(mde0)
+ if(.issueIntermediateParams){
+ cat("1st candidate:\n", round(es,6), " crit:", round(crit1,6), , " ")
+ }
+ if(quantile(1+es[2]*x0/es[1], epsn/n)>0){
+ mdeb <- mde0
+ crit0 <- crit1
+ es0 <- es
+ if(.issueIntermediateParams){
+ cat("side condition '1+sc/sh (x-mu) >0' fulfilled;\n")
+ }
+ }else{
+ if(.issueIntermediateParams){
+ cat("side condition '1+sc/sh (x-mu) >0' violated;\n")
+ }
+ }
+ }
+ }
+
+ i <- 0
+ for(xi in xiGrid){
+ i <- i + 1
+ funl <- function(sig){
+ mygpd1 <- GPareto(loc=0,scale=sig,shape=xi)
+ CvMDist(x0,mygpd1)
+ }
+ intlo <- quantile(-xi*(x-mu),1-epsn/n)
+ intv <- c(max(1e-5,intlo), s0)
+ sigCvMMD1 <- optimize(funl, interval=intv)$minimum
+ mygpd <- GParetoFamily(loc=0,scale=sigCvMMD1,shape=xi, withPos=withPos,
+ start0Est = fu)
+ mde0 <- try(MDEstimator(x0, mygpd, distance=CvMDist, startPar=c("scale"=sigCvMMD1,"shape"=xi)),silent=TRUE)
+ es0 <- c(NA,NA)
+ if(!is(mde0,"try-error")){
+ es <- estimate(mde0)
+ crit1 <- criterion(mde0)
+ if(.issueIntermediateParams){
+ cat("candidate no",i+1, ":\n", round(es,6), " crit:", round(crit1,6), " ")
+ }
+ if(quantile(1+es[2]*x0/es[1], epsn/n)>0){
+ if(.issueIntermediateParams){
+ cat("side condition '1+sc/sh (x-mu) >0' fulfilled;\n")
+ }
+ if(crit1<crit0){
+ mdeb <- mde0
+ crit0 <- crit1
+ es0 <- es
+ }
+ }else{
+ if(.issueIntermediateParams){
+ cat("side condition '1+sc/sh (x-mu) >0' violated;\n")
+ }
+ es[1] <- intlo+1e-5
+ mygpd2 <- GPareto(loc=0,scale=es[1],shape=es[2])
+ crit1 <- CvMDist(x0,mygpd2)
+ if(.issueIntermediateParams){
+ cat("candidate no",i+1, "(b):\n", round(es,6), " crit:", round(crit1,6), " ")
+ }
+ if(crit1<crit0){
+ crit0 <- crit1
+ es0 <- es
+ }
+ }
+ }
+ }
+ names(es) <- c("scale","shape")
+ return(es)
+}
+
+#.getMuBetaXiGPD <- function(x, xiGrid = .getXiGrid(), withPos=TRUE, secLevel = 0.7,
+# .issueIntermediateParams = FALSE){
+# mu <- quantile(x,exp(-1))
+# es <- .getBetaXiGPD(x=x, mu=mu, xiGrid=xiGrid, withPos=withPos,
+# secLevel = secLevel,
+# .issueIntermediateParams = .issueIntermediateParams)
+# es0 <- c(mu,es)
+# names(es0) <- c("loc","scale","shape")
+# return(es0)
+#}
Modified: branches/robast-1.0/pkg/RobExtremes/man/GParetoFamily.Rd
===================================================================
--- branches/robast-1.0/pkg/RobExtremes/man/GParetoFamily.Rd 2014-04-09 21:20:47 UTC (rev 747)
+++ branches/robast-1.0/pkg/RobExtremes/man/GParetoFamily.Rd 2014-04-12 13:05:31 UTC (rev 748)
@@ -9,7 +9,8 @@
\usage{
GParetoFamily(loc = 0, scale = 1, shape = 0.5, of.interest = c("scale", "shape"),
p = NULL, N = NULL, trafo = NULL, start0Est = NULL, withPos = TRUE,
- withCentL2 = FALSE, withL2derivDistr = FALSE, ..ignoreTrafo = FALSE)
+ secLevel = 0.7, withCentL2 = FALSE, withL2derivDistr = FALSE,
+ ..ignoreTrafo = FALSE)
}
\arguments{
\item{loc}{ real: known/fixed threshold/location parameter }
@@ -23,6 +24,13 @@
\item{trafo}{ matrix or NULL: transformation of the parameter }
\item{start0Est}{ startEstimator --- if \code{NULL} \code{\link{medkMADhybr}} is used }
\item{withPos}{ logical of length 1: Is shape restricted to positive values? }
+ \item{secLevel}{ a numeric of length 1:
+ In the ideal GEV model, for each observastion \eqn{X_i}{Xi}, the expression
+ \eqn{1+\frac{{\rm shape}(X_i-{\rm loc})}{{\rm scale}}}{1+shape(Xi-loc)/scale}
+ must be positive, which in principle could be attacked by a single outlier.
+ Hence for sample size \eqn{n} we allow for \eqn{\varepsilon n}{eps n}
+ violations, interpreting the violations as outliers. Here
+ \eqn{\varepsilon = {\tt secLevel}/\sqrt{n}}{eps = secLevel/sqrt(n)}. }
\item{withCentL2}{logical: shall L2 derivative be centered by substracting
the E()? Defaults to \code{FALSE}, but higher accuracy can be achieved
when set to \code{TRUE}.}
More information about the Robast-commits
mailing list