[Robast-commits] r666 - in branches/robast-0.9/pkg: RobExtremes/R RobExtremes/man RobExtremesBuffer
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Mon May 27 22:27:32 CEST 2013
Author: ruckdeschel
Date: 2013-05-27 22:27:32 +0200 (Mon, 27 May 2013)
New Revision: 666
Added:
branches/robast-0.9/pkg/RobExtremesBuffer/checkOfInterest.R
Modified:
branches/robast-0.9/pkg/RobExtremes/R/GEVFamily.R
branches/robast-0.9/pkg/RobExtremes/R/GParetoFamily.R
branches/robast-0.9/pkg/RobExtremes/R/WeibullFamily.R
branches/robast-0.9/pkg/RobExtremes/man/GEVFamily.Rd
branches/robast-0.9/pkg/RobExtremes/man/GParetoFamily.Rd
branches/robast-0.9/pkg/RobExtremes/man/WeibullFamily.Rd
branches/robast-0.9/pkg/RobExtremesBuffer/MishaLMScripts.R
Log:
RobExtremes: fixed some bugs in of.interest
RobExtremesBuffer: submitted some checks in checkOfInterest.R
Modified: branches/robast-0.9/pkg/RobExtremes/R/GEVFamily.R
===================================================================
--- branches/robast-0.9/pkg/RobExtremes/R/GEVFamily.R 2013-05-27 13:04:23 UTC (rev 665)
+++ branches/robast-0.9/pkg/RobExtremes/R/GEVFamily.R 2013-05-27 20:27:32 UTC (rev 666)
@@ -145,7 +145,8 @@
p = NULL, N = NULL, trafo = NULL,
start0Est = NULL, withPos = TRUE,
withCentL2 = FALSE,
- withL2derivDistr = FALSE){
+ withL2derivDistr = FALSE,
+ ..ignoreTrafo = FALSE){
theta <- c(loc, scale, shape)
.warningGEVShapeLarge(shape)
@@ -163,6 +164,7 @@
if(!is.null(p)){
btq <- substitute({ q <- loc0 + theta[1]*((-log(p0))^(-theta[2])-1)/theta[2]
names(q) <- "quantile"
+ q
}, list(loc0 = loc, p0 = p))
bDq <- substitute({ scale <- theta[1]; shape <- theta[2]
@@ -208,10 +210,11 @@
D }, list(loc0 = loc, N0 = N))
}
- if(is.null(trafo))
+ fromOfInt <- FALSE
+ if(is.null(trafo)||..ignoreTrafo){fromOfInt <- TRUE
trafo <- .define.tau.Dtau(of.interest, btq, bDq, btes, bDes,
btel, bDel, p, N)
- else if(is.matrix(trafo) & nrow(trafo) > 2)
+ }else if(is.matrix(trafo) & nrow(trafo) > 2)
stop("number of rows of 'trafo' > 2")
####
param <- ParamFamParameter(name = "theta", main = c(theta[2],theta[3]),
@@ -229,8 +232,9 @@
## Pickand estimator
if(is.null(start0Est)){
#source("kMedMad_Qn_Estimators.R")
- e0 <- estimate(PickandsEstimator(x,ParamFamily=GEVFamily(
- loc = theta[1], scale = theta[2], shape = theta[3])))
+ PF <- GEVFamily(loc = theta[1], scale = theta[2], shape = theta[3])
+ e1 <- PickandsEstimator(x,ParamFamily=PF)
+ e0 <- estimate(e1)
}else{
if(is(start0Est,"function")){
e1 <- start0Est(x, ...)
@@ -365,14 +369,25 @@
imageDistr(RandVar = L2deriv, distr = distribution))
}
- L2Fam at fam.call <- substitute(GEVFamily(loc = loc0, scale = scale0,
+ if(fromOfInt){
+ L2Fam at fam.call <- substitute(GEVFamily(loc = loc0, scale = scale0,
shape = shape0, of.interest = of.interest0,
+ p = p0, N = N0,
+ withPos = withPos0, withCentL2 = FALSE,
+ withL2derivDistr = FALSE, ..ignoreTrafo = TRUE),
+ list(loc0 = loc, scale0 = scale, shape0 = shape,
+ of.interest0 = of.interest, p0 = p, N0 = N,
+ withPos0 = withPos))
+ }else{
+ L2Fam at fam.call <- substitute(GEVFamily(loc = loc0, scale = scale0,
+ shape = shape0, of.interest = NULL,
p = p0, N = N0, trafo = trafo0,
withPos = withPos0, withCentL2 = FALSE,
withL2derivDistr = FALSE),
list(loc0 = loc, scale0 = scale, shape0 = shape,
- of.interest0 = of.interest, p0 = p, N0 = N,
- trafo0 = trafo, withPos0 = withPos))
+ p0 = p, N0 = N,
+ withPos0 = withPos, trafo0 = trafo))
+ }
L2Fam at LogDeriv <- function(x){
x0 <- (x-loc)/scale
Modified: branches/robast-0.9/pkg/RobExtremes/R/GParetoFamily.R
===================================================================
--- branches/robast-0.9/pkg/RobExtremes/R/GParetoFamily.R 2013-05-27 13:04:23 UTC (rev 665)
+++ branches/robast-0.9/pkg/RobExtremes/R/GParetoFamily.R 2013-05-27 20:27:32 UTC (rev 666)
@@ -40,7 +40,8 @@
p = NULL, N = NULL, trafo = NULL,
start0Est = NULL, withPos = TRUE,
withCentL2 = FALSE,
- withL2derivDistr = FALSE){
+ withL2derivDistr = FALSE,
+ ..ignoreTrafo = FALSE){
theta <- c(loc, scale, shape)
of.interest <- .pretreat.of.interest(of.interest,trafo)
@@ -57,6 +58,7 @@
if(!is.null(p)){
btq <- substitute({ q <- loc0 + theta[1]*((1-p0)^(-theta[2])-1)/theta[2]
names(q) <- "quantile"
+ q
}, list(loc0 = loc, p0 = p))
bDq <- substitute({ scale <- theta[1]; shape <- theta[2]
@@ -98,10 +100,11 @@
D }, list(loc0 = loc, N0 = N))
}
- if(is.null(trafo))
+ fromOfInt <- FALSE
+ if(is.null(trafo)||..ignoreTrafo){fromOfInt <- TRUE
trafo <- .define.tau.Dtau(of.interest, btq, bDq, btes, bDes,
btel, bDel, p, N)
- else if(is.matrix(trafo) & nrow(trafo) > 2)
+ }else if(is.matrix(trafo) & nrow(trafo) > 2)
stop("number of rows of 'trafo' > 2")
# code .define.tau.Dtau is in file GEVFamily.R
@@ -119,9 +122,11 @@
## Pickand estimator
if(is.null(start0Est)){
- e0 <- estimate(medkMADhybr(x, k=10, ParamFamily=GParetoFamily(loc = theta[1],
- scale = theta[2], shape = theta[3]),
- q.lo = 1e-3, q.up = 15))
+ 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)
}else{
if(is(start0Est,"function")){
e1 <- start0Est(x, ...)
@@ -242,14 +247,26 @@
imageDistr(RandVar = L2deriv, distr = distribution))
}
- L2Fam at fam.call <- substitute(GParetoFamily(loc = loc0, scale = scale0,
+ if(fromOfInt){
+ L2Fam at fam.call <- substitute(GParetoFamily(loc = loc0, scale = scale0,
shape = shape0, of.interest = of.interest0,
+ p = p0, N = N0,
+ withPos = withPos0, withCentL2 = FALSE,
+ withL2derivDistr = FALSE, ..ignoreTrafo = TRUE),
+ list(loc0 = loc, scale0 = scale, shape0 = shape,
+ of.interest0 = of.interest, p0 = p, N0 = N,
+ withPos0 = withPos))
+ }else{
+ L2Fam at fam.call <- substitute(GParetoFamily(loc = loc0, scale = scale0,
+ shape = shape0, of.interest = NULL,
p = p0, N = N0, trafo = trafo0,
withPos = withPos0, withCentL2 = FALSE,
withL2derivDistr = FALSE),
list(loc0 = loc, scale0 = scale, shape0 = shape,
- of.interest0 = of.interest, p0 = p, N0 = N,
- trafo0 = trafo, withPos0 = withPos))
+ p0 = p, N0 = N,
+ withPos0 = withPos, trafo0 = trafo))
+ }
+
L2Fam at LogDeriv <- function(x) (shape+1)/(scale+shape*(x-loc))
L2Fam at L2deriv <- L2deriv
Modified: branches/robast-0.9/pkg/RobExtremes/R/WeibullFamily.R
===================================================================
--- branches/robast-0.9/pkg/RobExtremes/R/WeibullFamily.R 2013-05-27 13:04:23 UTC (rev 665)
+++ branches/robast-0.9/pkg/RobExtremes/R/WeibullFamily.R 2013-05-27 20:27:32 UTC (rev 666)
@@ -38,7 +38,8 @@
p = NULL, N = NULL, trafo = NULL,
start0Est = NULL, withPos = TRUE,
withCentL2 = FALSE,
- withL2derivDistr = FALSE){
+ withL2derivDistr = FALSE,
+ ..ignoreTrafo = FALSE){
theta <- c(scale, shape)
of.interest <- .pretreat.of.interest(of.interest,trafo)
@@ -55,7 +56,8 @@
if(!is.null(p)){
btq <- substitute({ q <- theta[1]*(-log(1-p0))^(1/theta[2])
names(q) <- "quantile"
- }, list(loc0 = loc, p0 = p))
+ q
+ }, list(p0 = p))
bDq <- substitute({ scale <- theta[1]; shape <- theta[2]
lp <- -log(1-p0)
@@ -68,9 +70,9 @@
s1 <- 1+1/theta[2]
pg <- pgamma(-log(p0),s1, lower.tail = FALSE)
g0 <- gamma(s1)
- es <- theta[1] * g0 * pg /(1-p0) + loc0 }
+ es <- theta[1] * g0 * pg /(1-p0)}
names(es) <- "expected shortfall"
- es }, list(loc0 = loc, p0 = p))
+ es }, list(p0 = p))
bDes <- substitute({ if(theta[2]>=1L){ D1 <- D2 <- NA} else {
s1 <- 1+1/theta[2]
pg <- pgamma(-log(p0), s1, lower.tail = FALSE)
@@ -81,12 +83,12 @@
D <- t(c(D1, D2))
rownames(D) <- "expected shortfall"
colnames(D) <- NULL
- D }, list(loc0 = loc, p0 = p))
+ D }, list(p0 = p))
}
if(!is.null(N)){
btel <- substitute({ el <- N0*(theta[1]*gamma(1+1/theta[2]))
names(el) <- "expected loss"
- el }, list(loc0 = loc,N0 = N))
+ el }, list(N0 = N))
bDel <- substitute({ scale <- theta[1]; shape <- theta[2]
s1 <- 1+1/shape
D1 <- N0*gamma(s1)
@@ -94,13 +96,14 @@
D <- t(c(D1, D2))
rownames(D) <- "expected loss"
colnames(D) <- NULL
- D }, list(loc0 = loc, N0 = N))
+ D }, list(N0 = N))
}
- if(is.null(trafo))
+ fromOfInt <- FALSE
+ if(is.null(trafo)||..ignoreTrafo){fromOfInt <- TRUE
trafo <- .define.tau.Dtau(of.interest, btq, bDq, btes, bDes,
btel, bDel, p, N)
- else if(is.matrix(trafo) & nrow(trafo) > 2)
+ }else if(is.matrix(trafo) & nrow(trafo) > 2)
stop("number of rows of 'trafo' > 2")
# code .define.tau.Dtau is in file GEVFamily.R
@@ -118,7 +121,8 @@
## Pickand estimator
if(is.null(start0Est)){
- e0 <- estimate(QuantileBCCEstimator(x))
+ e1 <- QuantileBCCEstimator(x)
+ e0 <- estimate(e1)
}else{
if(is(start0Est,"function")){
e1 <- start0Est(x, ...)
@@ -232,14 +236,25 @@
imageDistr(RandVar = L2deriv, distr = distribution))
}
- L2Fam at fam.call <- substitute(WeibullFamily(scale = scale0,
+ if(fromOfInt){
+ L2Fam at fam.call <- substitute(WeibullFamily(scale = scale0,
shape = shape0, of.interest = of.interest0,
+ p = p0, N = N0,
+ withPos = withPos0, withCentL2 = FALSE,
+ withL2derivDistr = FALSE, ..ignoreTrafo = TRUE),
+ list(scale0 = scale, shape0 = shape,
+ of.interest0 = of.interest, p0 = p, N0 = N,
+ withPos0 = withPos))
+ }else{
+ L2Fam at fam.call <- substitute(WeibullFamily(scale = scale0,
+ shape = shape0, of.interest = NULL,
p = p0, N = N0, trafo = trafo0,
withPos = withPos0, withCentL2 = FALSE,
withL2derivDistr = FALSE),
list(scale0 = scale, shape0 = shape,
- of.interest0 = of.interest, p0 = p, N0 = N,
- trafo0 = trafo, withPos0 = withPos))
+ p0 = p, N0 = N,
+ withPos0 = withPos, trafo0 = trafo))
+ }
L2Fam at LogDeriv <- function(x){ z <- x/scale
log(shape)-log(scale)+(shape-1)*log(z)-shape*z^(shape-1)
Modified: branches/robast-0.9/pkg/RobExtremes/man/GEVFamily.Rd
===================================================================
--- branches/robast-0.9/pkg/RobExtremes/man/GEVFamily.Rd 2013-05-27 13:04:23 UTC (rev 665)
+++ branches/robast-0.9/pkg/RobExtremes/man/GEVFamily.Rd 2013-05-27 20:27:32 UTC (rev 666)
@@ -9,7 +9,7 @@
\usage{
GEVFamily(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)
+ withCentL2 = FALSE, withL2derivDistr = FALSE, ..ignoreTrafo = FALSE)
}
\arguments{
\item{loc}{ real: known/fixed threshold/location parameter }
@@ -28,6 +28,7 @@
when set to \code{TRUE}.}
\item{withL2derivDistr}{logical: shall the distribution of the L2 derivative
be computed? Defaults to \code{FALSE} (to speeds up computations).}
+ \item{..ignoreTrafo}{logical: only used internally in \code{kStepEstimator}; do not change this.}
}
\details{
The slots of the corresponding L2 differentiable
Modified: branches/robast-0.9/pkg/RobExtremes/man/GParetoFamily.Rd
===================================================================
--- branches/robast-0.9/pkg/RobExtremes/man/GParetoFamily.Rd 2013-05-27 13:04:23 UTC (rev 665)
+++ branches/robast-0.9/pkg/RobExtremes/man/GParetoFamily.Rd 2013-05-27 20:27:32 UTC (rev 666)
@@ -9,7 +9,7 @@
\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)
+ withCentL2 = FALSE, withL2derivDistr = FALSE, ..ignoreTrafo = FALSE)
}
\arguments{
\item{loc}{ real: known/fixed threshold/location parameter }
@@ -28,6 +28,7 @@
when set to \code{TRUE}.}
\item{withL2derivDistr}{logical: shall the distribution of the L2 derivative
be computed? Defaults to \code{FALSE} (to speeds up computations).}
+ \item{..ignoreTrafo}{logical: only used internally in \code{kStepEstimator}; do not change this.}
}
\details{
The slots of the corresponding L2 differentiable
Modified: branches/robast-0.9/pkg/RobExtremes/man/WeibullFamily.Rd
===================================================================
--- branches/robast-0.9/pkg/RobExtremes/man/WeibullFamily.Rd 2013-05-27 13:04:23 UTC (rev 665)
+++ branches/robast-0.9/pkg/RobExtremes/man/WeibullFamily.Rd 2013-05-27 20:27:32 UTC (rev 666)
@@ -9,7 +9,7 @@
\usage{
WeibullFamily(scale = 1, shape = 0.5, of.interest = c("scale", "shape"),
p = NULL, N = NULL, trafo = NULL, start0Est = NULL, withPos = TRUE,
- withCentL2 = FALSE, withL2derivDistr = FALSE)
+ withCentL2 = FALSE, withL2derivDistr = FALSE, ..ignoreTrafo = FALSE)
}
\arguments{
\item{scale}{ positive real: scale parameter }
@@ -27,6 +27,7 @@
when set to \code{TRUE}.}
\item{withL2derivDistr}{logical: shall the distribution of the L2 derivative
be computed? Defaults to \code{FALSE} (to speeds up computations).}
+ \item{..ignoreTrafo}{logical: only used internally in \code{kStepEstimator}; do not change this.}
}
\details{
The slots of the corresponding L2 differentiable parameteric family are filled.
Modified: branches/robast-0.9/pkg/RobExtremesBuffer/MishaLMScripts.R
===================================================================
--- branches/robast-0.9/pkg/RobExtremesBuffer/MishaLMScripts.R 2013-05-27 13:04:23 UTC (rev 665)
+++ branches/robast-0.9/pkg/RobExtremesBuffer/MishaLMScripts.R 2013-05-27 20:27:32 UTC (rev 666)
@@ -71,7 +71,7 @@
#Peter
baseDir0 <- "C:/rtest/RobASt"
#Misha
-baseDir0 <- "D:/SVN repositories/robast"
+#baseDir0 <- "D:/SVN repositories/robast"
interpolDir <- "branches/robast-0.9/pkg/RobExtremes/inst/AddMaterial/interpolation"
interpolFile <- "plotInterpol.R"
##
Added: branches/robast-0.9/pkg/RobExtremesBuffer/checkOfInterest.R
===================================================================
--- branches/robast-0.9/pkg/RobExtremesBuffer/checkOfInterest.R (rev 0)
+++ branches/robast-0.9/pkg/RobExtremesBuffer/checkOfInterest.R 2013-05-27 20:27:32 UTC (rev 666)
@@ -0,0 +1,76 @@
+
+require(RobExtremes)
+G0 <- GParetoFamily(shape=0.7,of.interest="quantile", p=0.99)
+IC1 <- optIC(InfRobModel(G0,ContNeighborhood(rad=0.5)), risk=asMSE())
+set.seed(20130527)
+x <- r(G0)(2000)
+res <- kStepEstimator(x, IC1)
+
+G0 <- GParetoFamily(shape=0.7,of.interest="quantile", p=0.99)
+IC0 <- optIC(G0, risk=asCov())
+IC1 <- optIC(InfRobModel(G0,ContNeighborhood(rad=0.5)), risk=asMSE())
+plot(IC1)
+IC2 <- optIC(InfRobModel(G0,ContNeighborhood(rad=0.5)), risk=asBias())
+plot(IC2)
+IC3 <- radiusMinimaxIC(G0,ContNeighborhood(rad=0.5), risk=asMSE())
+IC3
+plot(IC3)
+G0a <- GParetoFamily(shape=0.7,of.interest="expected shortfall", p=0.99)
+IC0a <- optIC(G0a, risk=asCov())
+IC1a <- optIC(InfRobModel(G0a,ContNeighborhood(rad=0.5)), risk=asMSE())
+plot(IC1a)
+IC2a <- optIC(InfRobModel(G0a,ContNeighborhood(rad=0.5)), risk=asBias())
+plot(IC2a)
+IC3a <- radiusMinimaxIC(G0a,ContNeighborhood(rad=0.5), risk=asMSE())
+IC3a
+plot(IC3a)
+set.seed(20130527)
+x <- r(G0)(200)
+res <- kStepEstimator(x, IC1)
+estimate(res)
+confint(res,level=.95)
+q(distribution(G0))(.99)
+set.seed(20130527)
+x <- r(G0)(20000)
+res <- kStepEstimator(x, IC1)
+estimate(res)
+confint(res,level=.95)
+res00 <- MLEstimator(x, G0)
+res0 <- kStepEstimator(x, IC0,steps=20)
+estimate(res0)
+confint(res0,level=.95)
+
+G0 <- GEVFamily(shape=0.7,of.interest="quantile", p=0.99)
+IC0 <- optIC(G0, risk=asCov())
+IC1 <- optIC(InfRobModel(G0,ContNeighborhood(rad=0.5)), risk=asMSE())
+plot(IC1)
+set.seed(20130527)
+x <- r(G0)(1500)
+res <- kStepEstimator(x, IC1)
+estimate(res)
+confint(res,level=.95)
+q(distribution(G0))(.99)
+res00 <- MLEstimator(x, G0)
+res0 <- kStepEstimator(x, IC0,steps=20)
+estimate(res00)
+estimate(res0)
+confint(res0,level=.95)
+q(distribution(G0))(.99)
+
+require(RobExtremes)
+G0 <- WeibullFamily(shape=0.7,of.interest="quantile", p=0.99)
+IC0 <- optIC(G0, risk=asCov())
+IC1 <- optIC(InfRobModel(G0,ContNeighborhood(rad=0.5)), risk=asMSE())
+plot(IC1)
+set.seed(20130527)
+x <- r(G0)(1500)
+res <- kStepEstimator(x, IC1)
+estimate(res)
+confint(res,level=.95)
+q(distribution(G0))(.99)
+res00 <- MLEstimator(x, G0)
+res0 <- kStepEstimator(x, IC0,steps=20)
+estimate(res00)
+estimate(res0)
+confint(res0,level=.95)
+q(distribution(G0))(.99)
More information about the Robast-commits
mailing list