[Robast-commits] r732 - in branches/robast-1.0/pkg/RobExtremes: . R man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Thu Apr 3 01:09:09 CEST 2014
Author: ruckdeschel
Date: 2014-04-03 01:09:09 +0200 (Thu, 03 Apr 2014)
New Revision: 732
Added:
branches/robast-1.0/pkg/RobExtremes/R/startEstGEV.R
Modified:
branches/robast-1.0/pkg/RobExtremes/NAMESPACE
branches/robast-1.0/pkg/RobExtremes/R/GEVFamily.R
branches/robast-1.0/pkg/RobExtremes/R/GEVFamilyMuUnknown.R
branches/robast-1.0/pkg/RobExtremes/R/GParetoFamily.R
branches/robast-1.0/pkg/RobExtremes/R/getCVaR.R
branches/robast-1.0/pkg/RobExtremes/man/GEVFamily.Rd
branches/robast-1.0/pkg/RobExtremes/man/GEVFamilyMuUnknown.Rd
branches/robast-1.0/pkg/RobExtremes/man/GEVParameter-class.Rd
branches/robast-1.0/pkg/RobExtremes/man/getCVaR.Rd
Log:
RobExtremes:
+ fixed some issues with help of getVaR, getCVaR, ...
+ new print method for the results of these functions (and a corresponding S3class)
+ GEVFamily and GEVFamilyMuUnknown now have changed default starting estimators
realized in startEstGEV.R : a CvM-MDE with xi varying on a grid...
+ GParetoFamily now handles left endpoint correctly and catches xi < -1/2
+ GEVFamily[MuUnknown] for xi>0 now handles left endpoint correctly and catches xi < -1/2
+ warning for large xi is switched off in GEVFamily[MuUnknown] if called internally
+ double definition of ddigamma eliminated
+ GEV now has a robust starting estimator for mu unknown
Modified: branches/robast-1.0/pkg/RobExtremes/NAMESPACE
===================================================================
--- branches/robast-1.0/pkg/RobExtremes/NAMESPACE 2014-03-23 16:50:05 UTC (rev 731)
+++ branches/robast-1.0/pkg/RobExtremes/NAMESPACE 2014-04-02 23:09:09 UTC (rev 732)
@@ -45,3 +45,4 @@
export("loc", "loc<-", "kMAD", "Sn", "Qn",
"asvarMedkMAD","asvarPickands", "asvarQBCC")
exportMethods("rescaleFunction")
+S3method(print, riskMeasure)
\ No newline at end of file
Modified: branches/robast-1.0/pkg/RobExtremes/R/GEVFamily.R
===================================================================
--- branches/robast-1.0/pkg/RobExtremes/R/GEVFamily.R 2014-03-23 16:50:05 UTC (rev 731)
+++ branches/robast-1.0/pkg/RobExtremes/R/GEVFamily.R 2014-04-02 23:09:09 UTC (rev 732)
@@ -127,8 +127,10 @@
return(FALSE)
if (any(param[1] <= tol))
return(FALSE)
- if (any(param[2] <= tol))
+ if(object at param@withPosRestr) if (any(param[2] <= tol))
return(FALSE)
+ if (any(param[2] <= -1/2))
+ return(FALSE)
return(TRUE)
})
@@ -146,9 +148,10 @@
start0Est = NULL, withPos = TRUE,
withCentL2 = FALSE,
withL2derivDistr = FALSE,
- ..ignoreTrafo = FALSE){
+ ..ignoreTrafo = FALSE,
+ ..withWarningGEV = TRUE){
theta <- c(loc, scale, shape)
- .warningGEVShapeLarge(shape)
+ if(..withWarningGEV).warningGEVShapeLarge(shape)
of.interest <- .pretreat.of.interest(of.interest,trafo)
@@ -234,13 +237,16 @@
## starting parameters
startPar <- function(x,...){
mu <- theta[1]
-
+
## Pickand estimator
if(is.null(start0Est)){
+ ### replaced 20140402: CvMMDE-with xi on Grid
#source("kMedMad_Qn_Estimators.R")
- PF <- GEVFamily(loc = theta[1], scale = theta[2], shape = theta[3])
- e1 <- PickandsEstimator(x,ParamFamily=PF)
- e0 <- estimate(e1)
+ ### replaced 20140402:
+ # PF <- GEVFamily(loc = theta[1], scale = theta[2], shape = theta[3])
+ # e1 <- PickandsEstimator(x,ParamFamily=PF)
+ # e0 <- estimate(e1)
+ e0 <- .getBetaXiGEV(x=x, mu=mu, xiGrid=.getXiGrid(), withPos=withPos)
}else{
if(is(start0Est,"function")){
e1 <- start0Est(x, ...)
@@ -263,9 +269,11 @@
theta <- abs(theta)
}else{
if(!is.null(names(theta))){
+ if(theta["shape"]< (-1/2)) theta["shape"] <- -1/2+1e-4
theta["scale"] <- abs(theta["scale"])
}else{
theta[1] <- abs(theta[1])
+ if(theta[2]< (-1/2)) theta[2] <- -1/2+1e-4
}
}
return(theta)
@@ -273,12 +281,12 @@
modifyPar <- function(theta){
theta <- makeOKPar(theta)
- .warningGEVShapeLarge(theta["shape"])
+ if(..withWarningGEV).warningGEVShapeLarge(theta["shape"])
if(!is.null(names(theta))){
sc <- theta["scale"]
sh <- theta["shape"]
}else{
- theta <- abs(theta)
+ # changed 20140402: theta <- abs(theta)
sc <- theta[1]
sh <- theta[2]
}
@@ -291,7 +299,7 @@
sc <- force(main(param)[1])
k <- force(main(param)[2])
tr <- fixed(param)[1]
- .warningGEVShapeLarge(k)
+ if(..withWarningGEV).warningGEVShapeLarge(k)
Lambda1 <- function(x) {
y <- x*0
@@ -327,7 +335,7 @@
FisherInfo.fct <- function(param) {
sc <- force(main(param)[1])
k <- force(main(param)[2])
- .warningGEVShapeLarge(k)
+ if(..withWarningGEV).warningGEVShapeLarge(k)
G20 <- gamma(2*k)
G10 <- gamma(k)
G11 <- digamma(k)*gamma(k)
Modified: branches/robast-1.0/pkg/RobExtremes/R/GEVFamilyMuUnknown.R
===================================================================
--- branches/robast-1.0/pkg/RobExtremes/R/GEVFamilyMuUnknown.R 2014-03-23 16:50:05 UTC (rev 731)
+++ branches/robast-1.0/pkg/RobExtremes/R/GEVFamilyMuUnknown.R 2014-04-02 23:09:09 UTC (rev 732)
@@ -13,8 +13,10 @@
return(FALSE)
if (any(param[2] <= tol))
return(FALSE)
- if (any(param[3] <= tol))
+ if(object at param@withPosRestr) if (any(param[3] <= tol))
return(FALSE)
+ if (any(param[3] <= -1/2))
+ return(FALSE)
return(TRUE)
})
@@ -32,9 +34,10 @@
start0Est = NULL, withPos = TRUE,
withCentL2 = FALSE,
withL2derivDistr = FALSE,
- ..ignoreTrafo = FALSE){
+ ..ignoreTrafo = FALSE,
+ ..withWarningGEV = TRUE){
theta <- c(loc, scale, shape)
- .warningGEVShapeLarge(shape)
+ if(..withWarningGEV).warningGEVShapeLarge(shape)
of.interest <- .pretreat.of.interest(of.interest,trafo)
@@ -121,24 +124,25 @@
## starting parameters
startPar <- function(x,...){
- mu <- min(x)
-
+
## Pickand estimator
if(is.null(start0Est)){
+ ### replaced 20140402: CvMMDE-with xi on Grid
#source("kMedMad_Qn_Estimators.R")
- PF <- GEVFamily(loc = theta[1], scale = theta[2], shape = theta[3])
- e1 <- PickandsEstimator(x,ParamFamily=PF)
- e0 <- estimate(e1)
+ # PF <- GEVFamily(loc = theta[1], scale = theta[2], shape = theta[3])
+ # e1 <- PickandsEstimator(x,ParamFamily=PF)
+ # e0 <- estimate(e1)
+ e0 <- .getMuBetaXiGEV(x=x, xiGrid=.getXiGrid(), withPos=withPos)
}else{
if(is(start0Est,"function")){
e1 <- start0Est(x, ...)
e0 <- if(is(e1,"Estimate")) estimate(e1) else e1
}else stop("Argument 'start0Est' must be a function or NULL.")
if(!is.null(names(e0)))
- e0 <- e0[c("scale", "shape")]
+ e0 <- e0[c("loc","scale", "shape")]
}
# print(e0); print(str(x)); print(head(summary(x))); print(mu)
- if(any(x < mu-e0["scale"]/e0["shape"]))
+ if(any(x < e0[1]-e0[2]/e0[3]))
stop("some data smaller than 'loc-scale/shape' ")
names(e0) <- NULL
@@ -148,12 +152,14 @@
## what to do in case of leaving the parameter domain
makeOKPar <- function(theta) {
if(withPos){
- theta <- abs(theta)
+ theta[2:3] <- abs(theta[2:3])
}else{
if(!is.null(names(theta))){
+ if(theta["shape"]< (-1/2)) theta["shape"] <- -1/2+1e-4
theta["scale"] <- abs(theta["scale"])
}else{
- theta[1] <- abs(theta[1])
+ theta[2] <- abs(theta[2])
+ if(theta[3]< (-1/2)) theta[3] <- -1/2+1e-4
}
}
return(theta)
@@ -161,14 +167,14 @@
modifyPar <- function(theta){
theta <- makeOKPar(theta)
- .warningGEVShapeLarge(theta["shape"])
+ if(..withWarningGEV).warningGEVShapeLarge(theta["shape"])
if(!is.null(names(theta))){
loc <- theta["loc"]
sc <- theta["scale"]
sh <- theta["shape"]
}else{
loc <- theta[1]
- theta[2:3] <- abs(theta[2:3])
+ #theta[2:3] <- abs(theta[2:3])
sc <- theta[2]
sh <- theta[3]
}
@@ -181,7 +187,7 @@
sc <- force(main(param)[2])
k <- force(main(param)[3])
tr <- force(main(param)[1])
- .warningGEVShapeLarge(k)
+ if(..withWarningGEV).warningGEVShapeLarge(k)
k1 <- k+1
Lambda0 <- function(x) {
@@ -233,7 +239,7 @@
sc <- force(main(param)[2])
k <- force(main(param)[3])
k1 <- k+1
- .warningGEVShapeLarge(k)
+ if(..withWarningGEV).warningGEVShapeLarge(k)
G20 <- gamma(2*k)
G10 <- gamma(k)
G11 <- digamma(k)*gamma(k)
@@ -319,11 +325,3 @@
L2Fam at .withEvalL2derivDistr <- FALSE
return(L2Fam)
}
-
-#ddigamma(t,s) is d/ds \int_0^t exp(-x) x^(s-1) dx
-
-ddigamma <- function(t,s){
- int <- function(x) exp(-x)*(log(x))*x^(s-1)
- integrate(int, lower=0, upper=t)$value
- }
-
\ No newline at end of file
Modified: branches/robast-1.0/pkg/RobExtremes/R/GParetoFamily.R
===================================================================
--- branches/robast-1.0/pkg/RobExtremes/R/GParetoFamily.R 2014-03-23 16:50:05 UTC (rev 731)
+++ branches/robast-1.0/pkg/RobExtremes/R/GParetoFamily.R 2014-04-02 23:09:09 UTC (rev 732)
@@ -12,11 +12,13 @@
param <- main(param)
if (!all(is.finite(param)))
return(FALSE)
- #if (any(param[1] <= tol))
- # return(FALSE)
+ if (any(param[1] <= tol))
+ return(FALSE)
if(object at param@withPosRestr)
if (any(param[2] <= tol))
return(FALSE)
+ if (any(param[2] <= -1/2))
+ return(FALSE)
return(TRUE)
})
@@ -142,8 +144,10 @@
e0 <- e0[c("scale", "shape")]
}
- if(any(x < tr-e0["scale"]/e0["shape"]))
- stop("some data smaller than 'loc-scale/shape' ")
+ if(any(x < tr-.Machine$double.eps))
+ stop("some data smaller than 'loc' ")
+# if(any(x < tr-e0["scale"]/e0["shape"]))
+# stop("some data smaller than 'loc-scale/shape' ")
names(e0) <- NULL
return(e0)
@@ -156,9 +160,11 @@
theta <- abs(theta)
}else{
if(!is.null(names(theta))){
+ if(theta["shape"]< (-1/2)) theta["shape"] <- -1/2+1e-4
theta["scale"] <- abs(theta["scale"])
}else{
theta[1] <- abs(theta[1])
+ if(theta[2]< (-1/2)) theta[2] <- -1/2+1e-4
}
}
return(theta)
Modified: branches/robast-1.0/pkg/RobExtremes/R/getCVaR.R
===================================================================
--- branches/robast-1.0/pkg/RobExtremes/R/getCVaR.R 2014-03-23 16:50:05 UTC (rev 731)
+++ branches/robast-1.0/pkg/RobExtremes/R/getCVaR.R 2014-04-02 23:09:09 UTC (rev 732)
@@ -15,9 +15,25 @@
res <- param(L2Fam)@trafo(estimate(est))
VaR <- res[[1]]
varVaR <- (res[[2]]) %*% asvar(est) %*% t(res[[2]])
- return(c(VaR=VaR,sqrt(varVaR/length(data))))
+ res <- c(VaR,sqrt(varVaR/length(data)))
+ names(res) <- c("Risk","varofRisk")
+ class(res) <- "riskMeasure"
+ res
}
+print.riskMeasure <- function(x, level=NULL, ...){
+ mc <- as.list(match.call(expand.dots=TRUE)[-1])
+ digits <- if(is.null(mc$digits)) 3 else mc$digits
+ if(is.null(level)){
+ cat(" ",signif(x[1],digits),"\n")
+ cat("(",signif(x[2],digits),")\n")
+ }else{qn <- qnorm((level+1)/2)
+ CI <- c(-1,1)*qn*x[2]+x[1]
+ cat(" ",signif(x[1],digits)," [", signif(CI[1],digits), ",",
+ signif(CI[2],digits),"]\n")
+ }
+}
+
getVaR <- function(data, model, level, rob=TRUE)
.getTau(data, model, level, rob, of.interest="quantile", substitute(L2FamC$p <- level))
Added: branches/robast-1.0/pkg/RobExtremes/R/startEstGEV.R
===================================================================
--- branches/robast-1.0/pkg/RobExtremes/R/startEstGEV.R (rev 0)
+++ branches/robast-1.0/pkg/RobExtremes/R/startEstGEV.R 2014-04-02 23:09:09 UTC (rev 732)
@@ -0,0 +1,38 @@
+.getXiGrid <- function(){seq(-0.4,4,by=0.3)}
+
+
+.getBetaXiGEV <- function(x, mu, xiGrid = .getXiGrid(), withPos=TRUE){
+ x0 <- x-mu
+ s0 <- max(x0)-min(x0)
+ crit0 <- Inf
+ fu <- function(x,...) .getBetaXiGEV(x,mu,xiGrid = xiGrid,withPos=withPos)
+ for(xi in xiGrid){
+ funl <- function(sig){
+ mygev1 <- GEV(loc=0,scale=sig,shape=xi)
+ CvMDist(x0,mygev1)
+ }
+ sigCvMMD1 <- optimize(funl, interval=c(1e-5,s0))$minimum
+ print(c("sigma"=sigCvMMD1,"xi"=xi))
+ mygev <- GEVFamily(loc=0,scale=sigCvMMD1,shape=xi, withPos=withPos,
+ start0Est = fu, ..withWarningGEV=FALSE)
+ print(mygev)
+ print(param(mygev))
+ mde0 <- MDEstimator(x0, mygev, distance=CvMDist, startPar=c("scale"=sigCvMMD1,"shape"=xi))
+ print(c("roh"=estimate(mde0)))
+ if(criterion(mde0)<crit0){
+ mdeb <- mde0
+ crit0 <- criterion(mde0)
+ }
+ }
+ es <- estimate(mdeb)
+ names(es) <- c("scale","shape")
+ return(es)
+}
+
+.getMuBetaXiGEV <- function(x, xiGrid = .getXiGrid(), withPos=TRUE){
+ mu <- quantile(x,exp(-1))
+ es <- .getBetaXiGEV(x=x, mu=mu, xiGrid=xiGrid, withPos=withPos)
+ es0 <- c(mu,es)
+ names(es0) <- c("loc","scale","shape")
+ return(es0)
+}
Modified: branches/robast-1.0/pkg/RobExtremes/man/GEVFamily.Rd
===================================================================
--- branches/robast-1.0/pkg/RobExtremes/man/GEVFamily.Rd 2014-03-23 16:50:05 UTC (rev 731)
+++ branches/robast-1.0/pkg/RobExtremes/man/GEVFamily.Rd 2014-04-02 23:09:09 UTC (rev 732)
@@ -9,7 +9,8 @@
\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, ..ignoreTrafo = FALSE)
+ withCentL2 = FALSE, withL2derivDistr = FALSE, ..ignoreTrafo = FALSE,
+ ..withWarningGEV = TRUE)
}
\arguments{
\item{loc}{ real: known/fixed threshold/location parameter }
@@ -29,6 +30,7 @@
\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.}
+ \item{..withWarningGEV}{logical: shall warnings be issued if shape is large?}
}
\details{
The slots of the corresponding L2 differentiable
Modified: branches/robast-1.0/pkg/RobExtremes/man/GEVFamilyMuUnknown.Rd
===================================================================
--- branches/robast-1.0/pkg/RobExtremes/man/GEVFamilyMuUnknown.Rd 2014-03-23 16:50:05 UTC (rev 731)
+++ branches/robast-1.0/pkg/RobExtremes/man/GEVFamilyMuUnknown.Rd 2014-04-02 23:09:09 UTC (rev 732)
@@ -9,7 +9,8 @@
\usage{
GEVFamilyMuUnknown(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)
+ withCentL2 = FALSE, withL2derivDistr = FALSE, ..ignoreTrafo = FALSE,
+ ..withWarningGEV = TRUE)
}
\arguments{
\item{loc}{ real: known/fixed threshold/location parameter }
@@ -29,6 +30,7 @@
\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.}
+ \item{..withWarningGEV}{logical: shall warnings be issued if shape is large?}
}
\details{
The slots of the corresponding L2 differentiable
Modified: branches/robast-1.0/pkg/RobExtremes/man/GEVParameter-class.Rd
===================================================================
--- branches/robast-1.0/pkg/RobExtremes/man/GEVParameter-class.Rd 2014-03-23 16:50:05 UTC (rev 731)
+++ branches/robast-1.0/pkg/RobExtremes/man/GEVParameter-class.Rd 2014-04-02 23:09:09 UTC (rev 732)
@@ -19,11 +19,11 @@
\section{Slots}{
\describe{
\item{\code{loc}}{ real number: location parameter of
- a generalized Pareto distribution. }
+ a GEV distribution. }
\item{\code{scale}}{ real number: scale parameter of
- a generalized Pareto distribution. }
+ a GEV distribution. }
\item{\code{shape}}{ real number: shape parameter of
- a generalized Pareto distribution. }
+ a GEV distribution. }
\item{\code{name}}{ default name is
\dQuote{parameter of a GEV distribution}. }
}
Modified: branches/robast-1.0/pkg/RobExtremes/man/getCVaR.Rd
===================================================================
--- branches/robast-1.0/pkg/RobExtremes/man/getCVaR.Rd 2014-03-23 16:50:05 UTC (rev 731)
+++ branches/robast-1.0/pkg/RobExtremes/man/getCVaR.Rd 2014-04-02 23:09:09 UTC (rev 732)
@@ -1,7 +1,8 @@
\name{getCVaR}
\alias{getVaR}
-\alias{getVaR}
+\alias{getCVaR}
\alias{getEL}
+\alias{print.riskMeasure}
\title{Risk Measures for Scale-Shape Families}
\description{
@@ -12,19 +13,29 @@
getVaR(data, model, level, rob=TRUE)
getCVaR(data, model, level, rob=TRUE)
getEL(data, model, N0, rob=TRUE)
+\method{print}{riskMeasure}(x, level=NULL, ...)
}
\arguments{
\item{data}{data at which to compute the risk measure. }
\item{model}{an object of class \code{"L2ScaleShapeFamily"}.
The parametric family at which to evaluate the risk measure. }
\item{level}{real: probability needed for VaR and CVaR. }
- \item{N0}{real: expected frequency for expected loss }
+ \item{N0}{real: expected frequency for expected loss. }
\item{rob}{logical; if \code{TRUE} (default) the RMXE-parametric estimator is
- used; otherwise the MLE.}
+ used; otherwise the MLE. }
+ \item{x}{an object of (S3-)class \code{"riskmeasure"}. }
+ \item{...}{further arguments for \code{print}. }
}
-\value{A vector of length 2 with first coordinate the respective risk measure
-and the second a corresponding (asymptotic) standard error for the risk
-measure.}
+\value{The risk measures \code{getVaR}, \code{getCVaR}, \code{getEL} return
+an (S3) object of class \code{"riskMeasure"}, i.e., a numeric vector
+of length 2 with components \code{"Risk"} and \code{"varofRisk"}
+containing the respective risk measure
+and a corresponding (asymptotic) standard error for the risk
+measure. To the return class \code{"riskMeasure"},
+there is a particular \code{print}-method; if the corresponding argument
+\code{level} is \code{NULL} (default) the corresponding standard error
+is printed together with the risk measure; otherwise a corresponding
+CLT-based confidence interval for the risk meausre is produced.}
\references{
P. Ruckdeschel, N. Horbenko (2013): Optimally-Robust Estimators in Generalized
Pareto Models. Statistics 47(4), 762--791.
@@ -35,7 +46,7 @@
\seealso{\code{\link{GParetoFamily}}, \code{\link{GEVFamily}}, \code{\link{WeibullFamily}}, \code{\link{GammaFamily}}}
\examples{
set.seed(123)
- GPD <- GPartoFamily(loc=20480, scale=7e4, shape=0.3)
+ GPD <- GParetoFamily(loc=20480, scale=7e4, shape=0.3)
data <- r(GPD)(500)
getVaR(data,GPD,0.99)
getVaR(data,GPD,0.99, rob=FALSE)
More information about the Robast-commits
mailing list