[Robast-commits] r538 - in branches/robast-0.9/pkg/RobExtremes: R tests tests/TestSuite
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Thu Jan 17 01:41:38 CET 2013
Author: ruckdeschel
Date: 2013-01-17 01:41:38 +0100 (Thu, 17 Jan 2013)
New Revision: 538
Added:
branches/robast-0.9/pkg/RobExtremes/R/comment.txt
Removed:
branches/robast-0.9/pkg/RobExtremes/R/drawPointsIC.R
branches/robast-0.9/pkg/RobExtremes/R/plotScaledIC.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/LDEstimator.R
branches/robast-0.9/pkg/RobExtremes/R/PickandsEstimator.R
branches/robast-0.9/pkg/RobExtremes/R/asvarPickands.R
branches/robast-0.9/pkg/RobExtremes/R/interpolSn.R
branches/robast-0.9/pkg/RobExtremes/R/plotOutlyingness.R
branches/robast-0.9/pkg/RobExtremes/tests/TestSuite.R
branches/robast-0.9/pkg/RobExtremes/tests/TestSuite/TestExpectation.R
Log:
RobExtremes: some intermediate results; still to be debugged
Modified: branches/robast-0.9/pkg/RobExtremes/R/GEVFamily.R
===================================================================
--- branches/robast-0.9/pkg/RobExtremes/R/GEVFamily.R 2013-01-17 00:40:40 UTC (rev 537)
+++ branches/robast-0.9/pkg/RobExtremes/R/GEVFamily.R 2013-01-17 00:41:38 UTC (rev 538)
@@ -29,9 +29,39 @@
## trafo: optional parameter transformation
## start0Est: startEstimator for MLE and MDE --- if NULL HybridEstimator is used;
-GEVFamily <- function(loc = 0, scale = 1, shape = 0.5,trafo = NULL,start0Est = NULL){
- if(is.null(trafo)) trafo = diag(2)
-
+GEVFamily <- function(loc = 0, scale = 1, shape = 0.5,
+ of.interest = c("scale", "shape"),
+ p = NULL, N = NULL, trafo = NULL,
+ start0Est = NULL, withPos = TRUE){
+ if(is.null(trafo)){
+ of.interest <- unique(of.interest)
+ if(length(of.interest) > 2)
+ stop("A maximum number of two parameters resp. parameter transformations may be selected.")
+ if(!all(of.interest %in% c("scale", "shape", "quantile", "expected loss", "expected shortfall")))
+ stop("Parameters resp. transformations of interest have to be selected from: ",
+ "'scale', 'shape', 'quantile', 'expected loss', 'expected shortfall'.")
+
+ ## reordering of of.interest
+ if(("scale" %in% of.interest) && ("scale" != of.interest[1])){
+ of.interest[2] <- of.interest[1]
+ of.interest[1] <- "scale"
+ }
+ if(!("scale" %in% of.interest) && ("shape" %in% of.interest) && ("shape" != of.interest[1])){
+ of.interest[2] <- of.interest[1]
+ of.interest[1] <- "shape"
+ }
+ if(!any(c("scale", "shape") %in% of.interest) && ("quantile" %in% of.interest)
+ && ("quantile" != of.interest[1])){
+ of.interest[2] <- of.interest[1]
+ of.interest[1] <- "quantile"
+ }
+ if(!any(c("scale", "shape", "quantile") %in% of.interest)
+ && ("expected shortfall" %in% of.interest)
+ && ("expected shortfall" != of.interest[1])){
+ of.interest[2] <- of.interest[1]
+ of.interest[1] <- "expected shortfall"
+ }
+ }
theta <- c(loc, scale, shape)
##symmetry
@@ -40,9 +70,129 @@
## parameters
names(theta) <- c("loc", "scale", "shape")
- param <- ParamFamParameter(name = "theta", main = theta[2:3],
+
+ if(!is.null(p)){
+ btq <- substitute({ q <- loc0 + theta[1]*((-log(1-p0))^(-theta[2])-1)/theta[2]
+ names(q) <- "quantile"
+ }, list(loc0 = loc, p0 = p))
+
+ bDq <- substitute({ scale <- theta[1]; shape <- theta[2]
+ D1 <- ((-log(1-p0))^(-shape)-1)/shape
+ D2 <- -scale/shape*(D1 + log(-log(1-p0))*(-log(1-p0))^(-shape))
+ D <- t(c(D1, D2))
+ rownames(D) <- "quantile"; colnames(D) <- NULL
+ D }, list(p0 = p))
+ btes <- substitute({ if(theta[2]>=1L) es <- NA else {
+ pg <- pgamma(-log(p0),1-theta[2], lower.tail = FALSE)
+ es <- theta[1] * gamma(1-theta[2]) * pg / p0 /
+ theta[2] + loc0 }
+ names(es) <- "expected shortfall"
+ es }, list(loc0 = loc, p0 = p))
+ bDes <- substitute({ if(theta[2]>=1L){ D1 <- D2 <- NA} else {
+ scale <- theta[1]; shape <- theta[2]
+ pg <- pgamma(-log(p0), 1-theta[2], lower.tail = FALSE)
+ dd <- ddigamma(-log(p0),1-theta[2])
+ D1 <- gamma(1-theta[2])*pg/p0/theta[2]
+ D21 <- -theta[1]*gamma(1-theta[2])*pg/p0/theta[2]^2
+ D22 < -theta[1]*digamma(1-theta[2])*pg/p0/theta[2]
+ D23 <- theta[1]*dd/p0/theta[2]
+ D2 <- D21+D22+D23}
+ D <- t(c(D1, D2))
+ rownames(D) <- "expected shortfall"
+ colnames(D) <- NULL
+ D }, list(loc0 = loc, p0 = p))
+ }
+ if(!is.null(N)){
+ btel <- substitute({ if(theta[2]>=1L) el <- NA else{
+ el <- N0*(loc0+theta[1]*gamma(1-theta[2])/theta[2])}
+ names(el) <- "expected loss"
+ el }, list(loc0 = loc,N0 = N))
+ bDel <- substitute({ if(theta[2]>=1L){ D1 <- D2 <- NA}else{
+ scale <- theta[1]; shape <- theta[2]
+ D1 <- N0*gamma(1-shape)/shape
+ D2 <- -N0*theta[1]*digamma(1-theta[2])/theta[2]-
+ D1*scale/(1-shape)}
+ D <- t(c(D1, D2))
+ rownames(D) <- "expected loss"
+ colnames(D) <- NULL
+ D }, list(loc0 = loc, N0 = N))
+ }
+ if(is.null(trafo)){
+ tau <- NULL
+ if("scale" %in% of.interest){
+ tau <- function(theta){ th <- theta[1]; names(th) <- "scale"; th}
+ Dtau <- function(theta){ D <- t(c(1, 0)); rownames(D) <- "scale"; D}
+ }
+ if("shape" %in% of.interest){
+ if(is.null(tau)){
+ tau <- function(theta){th <- theta[2]; names(th) <- "shape"; th}
+ Dtau <- function(theta){D <- t(c(0,1));rownames(D) <- "shape";D}
+ }else{
+ tau <- function(theta){th <- theta
+ names(th) <- c("scale", "shape"); th}
+ Dtau <- function(theta){ D <- diag(2);
+ rownames(D) <- c("scale", "shape");D}
+ }
+ }
+ }
+ if("quantile" %in% of.interest){
+ if(is.null(p)) stop("Probability 'p' has to be specified.")
+ if(is.null(tau)){
+ tau <- function(theta){ }; body(tau) <- btq
+ Dtau <- function(theta){ };body(Dtau) <- bDq
+ }else{
+ tau1 <- tau
+ tau <- function(theta){ }
+ body(tau) <- substitute({ btq0; c(tau0(theta), q) },
+ list(btq0=btq, tau0 = tau1))
+ Dtau1 <- Dtau
+ Dtau <- function(theta){}
+ body(Dtau) <- substitute({ bDq0; rbind(Dtau0(theta), D) },
+ list(Dtau0 = Dtau1, bDq0 = bDq))
+ }
+ }
+ if("expected shortfall" %in% of.interest){
+ if(is.null(p)) stop("Probability 'p' has to be specified.")
+ if(is.null(tau)){
+ tau <- function(theta){ }; body(tau) <- btes
+ Dtau <- function(theta){ }; body(Dtau) <- bDes
+ }else{
+ tau1 <- tau
+ tau <- function(theta){ }
+ body(tau) <- substitute({ btes0; c(tau0(theta), es) },
+ list(tau0 = tau1, btes0=btes))
+ Dtau1 <- Dtau
+ Dtau <- function(theta){}
+ body(Dtau) <- substitute({ bDes0; rbind(Dtau0(theta), D) },
+ list(Dtau0 = Dtau1, bDes0=bDes))
+ }
+ }
+ if("expected loss" %in% of.interest){
+ if(is.null(N)) stop("Expected frequency 'N' has to be specified.")
+ if(is.null(tau)){
+ tau <- function(theta){ }; body(tau) <- btel
+ Dtau <- function(theta){ }; body(Dtau) <- bDel
+ }else{
+ tau1 <- tau
+ tau <- function(theta){ }
+ body(tau) <- substitute({ btel0; c(tau0(theta), el) },
+ list(tau0 = tau1, btel0=btel))
+ Dtau1 <- Dtau
+ Dtau <- function(theta){}
+ body(Dtau) <- substitute({ bDel0; rbind(Dtau0(theta), D) },
+ list(Dtau0 = Dtau1, bDel0=bDel))
+ }
+ }
+ trafo <- function(x){ list(fval = tau(x), mat = Dtau(x)) }
+ }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]),
fixed = theta[1],
- trafo = trafo)
+ trafo = trafo, withPosRestr = withPos,
+ .returnClsName ="ParamWithScaleAndShapeFamParameter")
## distribution
distribution <- GEV(loc = loc, scale = scale, shape = shape)
@@ -53,8 +203,9 @@
## Pickand estimator
if(is.null(start0Est)){
- source("kMedMad_Qn_Estimators.R")
- e0 <- EVTEst(x,est="kMAD")
+ #source("kMedMad_Qn_Estimators.R")
+ e0 <- PickandsEstimator(x,ParamFamily=GParetoFamily(loc = theta[1],
+ scale = theta[2], shape = theta[3]))
}else{
if(is(start0Est,"function")){
e1 <- start0Est(x, ...)
@@ -63,40 +214,52 @@
if(!is.null(names(e0)))
e0 <- e0[c("scale", "shape")]
}
- if(any(x < mu-e0["scale"]/e0["shape"])) stop("some data smaller than 'loc-scale/shape' ")
+ if(any(x < mu-e0["scale"]/e0["shape"]))
+ stop("some data smaller than 'loc-scale/shape' ")
names(e0) <- NULL
return(e0)
}
+ ## what to do in case of leaving the parameter domain
+ makeOKPar <- function(theta) {
+ if(withPos){
+ if(!is.null(names(theta)))
+ theta["shape"] <- abs(theta["shape"])
+ else theta[2] <- abs(theta[2])
+ }
+ return(theta)
+ }
+
modifyPar <- function(theta){
theta <- abs(theta)
GEV(loc = loc, scale = theta[1], shape = theta[2])
}
- ## what to do in case of leaving the parameter domain
- makeOKPar <- function(theta) {
- theta <- abs(theta)
- theta[2] <- pmin(theta[2],10)
- return(theta)
- }
## L2-derivative of the distribution
L2deriv.fct <- function(param) {
- beta <- force(main(param)[1])
- xi <- force(main(param)[2])
- mu <- fixed(param)[1]
+ sc <- force(main(param)[1])
+ k <- force(main(param)[2])
+ tr <- fixed(param)[1]
Lambda1 <- function(x) {
y <- x*0
- ind <- x>(mu-beta/xi)
- y[ind] <- -1/beta - xi*(-1/xi-1)*(x[ind]-mu)/beta^2/(1+xi*(x[ind]-mu)/beta) - (x[ind]-mu)*(1+xi*(x[ind]-mu)/beta)^(-1/xi-1)/beta^2
+ ind <- (x > mu-sc/k) # = [later] (x1>0)
+ x <- (x[ind]-tr)/sc
+ x1 <- 1 + k * x
+ y[ind] <- (x*(1-x1^(-1/k))-1)/x1/sc
+# xi*(-1/xi-1)*(x[ind]-mu)/beta^2/(1+xi*(x[ind]-mu)/beta) - (x[ind]-mu)*(1+xi*(x[ind]-mu)/beta)^(-1/xi-1)/beta^2
return(y)
}
Lambda2 <- function(x) {
y <- x*0
- ind <- x>(mu-beta/xi)
- y[ind]<- log(1+xi*(x[ind]-mu)/beta)/xi^2+(-1/xi-1)*(x[ind]-mu)/beta/(1+xi*(x[ind]-mu)/beta) - (1+xi*(x[ind]-mu)/beta)^(-1/xi)*log(1+xi*(x[ind]-mu)/beta)/xi^2 + (1+xi*(x[ind]-mu)/beta)^(-1/xi-1)*(x[ind]-mu)/beta/xi
+ ind <- (x > tr-sc/k) # = [later] (x1>0)
+ x <- (x[ind]-tr)/sc
+ x1 <- 1 + k * x
+ x2 <- x / x1
+ y[ind]<- (1-x1^(-1/k))/k*(log(x1)/k-x2)-x2
+# log(1+xi*(x[ind]-mu)/beta)/xi^2+(-1/xi-1)*(x[ind]-mu)/beta/(1+xi*(x[ind]-mu)/beta) - (1+xi*(x[ind]-mu)/beta)^(-1/xi)*log(1+xi*(x[ind]-mu)/beta)/xi^2 + (1+xi*(x[ind]-mu)/beta)^(-1/xi-1)*(x[ind]-mu)/beta/xi
return(y)
}
## additional centering of scores to increase numerical precision!
@@ -107,33 +270,89 @@
## Fisher Information matrix as a function of parameters
FisherInfo.fct <- function(param) {
- beta <- force(main(param)[1])
- xi <- force(main(param)[2])
- mu <- force(fixed(param)[1])
- fct <- L2deriv.fct(param)
- P <- GEV(loc = mu, scale = beta, shape = xi)
- E11 <- E(P,function(x)fct[[1]](x)*fct[[1]](x))
- E12 <- E(P,function(x)fct[[1]](x)*fct[[2]](x))
- E22 <- E(P,function(x)fct[[2]](x)*fct[[2]](x))
- return(PosSemDefSymmMatrix(matrix(c(E11,E12,E12,E22),2,2)))
+ sc <- force(main(param)[1])
+ k <- force(main(param)[2])
+# tr <- force(fixed(param)[1])
+# fct <- L2deriv.fct(param)
+# P2 <- GPareto(loc = tr, scale = sc, shape = k)
+ G20 <- gamma(2*k)
+ G10 <- gamma(k)
+ G21 <- digamma(2*k)
+ G11 <- digamma(k)
+ G01 <- digamma(1)
+ G02 <- trigamma(1)
+ I11 <- G20*2*k*(k^2-4*k+1)+G10*2*k*(k^2+k-1)+1
+ I11 <- I11/sc^2/k^2
+ I12 <- G20*k*(2*k^2-1)+ G10*k^2*(2-4*k) - k + 3
+ I12 <- I12 + G21*2*k^2 -G11*k^2*(k^2+2*k+1) -G01*k
+ I12 <- I12/sc/k^3
+ I22 <- G20*2*k*(k^2+1)-G10*2*k*(2+3*k+k^2) -3*k^2 +8*k +3
+ I22 <- I22 - G11*2*k^2*(1+k)+G10*2*k*(3-k)+k^2 *G02
+ I22 <- I22 /k^4
+ mat <- PosSemDefSymmMatrix(matrix(c(I11,I12,I12,I22),2,2))
+ dimnames(mat) <- list(scaleshapename,scaleshapename)
+ return(mat)
}
+
+# FisherInfo.fct <- function(param) {
+# beta <- force(main(param)[1])
+# xi <- force(main(param)[2])
+# mu <- force(fixed(param)[1])
+# fct <- L2deriv.fct(param)
+# P <- GEV(loc = mu, scale = beta, shape = xi)
+# E11 <- E(P,function(x)fct[[1]](x)*fct[[1]](x))
+# E12 <- E(P,function(x)fct[[1]](x)*fct[[2]](x))
+# E22 <- E(P,function(x)fct[[2]](x)*fct[[2]](x))
+# return(PosSemDefSymmMatrix(matrix(c(E11,E12,E12,E22),2,2)))
+# }
+
FisherInfo <- FisherInfo.fct(param)
name <- "Generalized Extreme Value Family with positive shape parameter: Frechet Family"
## initializing the GPareto family with components of L2-family
- res <- L2ParamFamily(name = name, param = param,
- distribution = distribution,
- L2deriv.fct = L2deriv.fct,
- FisherInfo.fct = FisherInfo.fct,
- FisherInfo = FisherInfo,
- startPar = startPar,
- makeOKPar = makeOKPar,
- modifyParam = modifyPar,
- .returnClsName = "GEVFamily")
- f.call <- substitute(GEVFamily(loc = loc0, scale = scale0, shape = shape0,trafo = trafo0),
- list(loc0 = loc, scale0 = scale, shape0 = shape,trafo0 = trafo))
- res at fam.call <- f.call
- return(res)
+ L2Fam <- new("GEVFamily")
+ L2Fam at scaleshapename <- scaleshapename
+ L2Fam at name <- name
+ L2Fam at param <- param
+ L2Fam at distribution <- distribution
+ L2Fam at L2deriv.fct <- L2deriv.fct
+ L2Fam at FisherInfo.fct <- FisherInfo.fct
+ L2Fam at FisherInfo <- FisherInfo
+ L2Fam at startPar <- startPar
+ L2Fam at makeOKPar <- makeOKPar
+ L2Fam at modifyParam <- modifyPar
+ L2Fam at L2derivSymm <- FunSymmList(NonSymmetric(), NonSymmetric())
+ L2Fam at L2derivDistrSymm <- DistrSymmList(NoSymmetry(), NoSymmetry())
+
+ L2deriv <- EuclRandVarList(RealRandVariable(L2deriv.fct(param),
+ Domain = Reals()))
+
+ L2Fam at fam.call <- substitute(GEV(loc = loc0, scale = scale0,
+ shape = shape0, of.interest = of.interest0,
+ p = p0, N = N0, trafo = trafo0,
+ withPos = withPos0),
+ list(loc0 = loc, scale0 = scale, shape0 = shape,
+ of.interest0 = of.interest, p0 = p, N0 = N,
+ trafo0 = trafo, withPos0 = withPos))
+
+ L2Fam at LogDeriv <- function(x){
+ x0 <- (x-loc)/scale
+ x1 <- 1 + x0 * shape
+ (shape+1)/scale/x1 + x1^(-1-1/shape)/scale
+ }
+
+ L2Fam at L2deriv <- L2deriv
+
+ L2Fam at L2derivDistr <- imageDistr(RandVar = L2deriv, distr = distribution)
+
+ return(L2Fam)
}
+#ddigamma(t,s) is d/ds \int_t^\infty exp(-x) (-log(x)) x^(-s) dx
+
+ddigamma <- function(t,s){
+ int <- function(x) exp(-x)*(-log(x))*x^(-s)
+ integrate(int, lower=t, upper=Inf)$value
+ }
+
\ No newline at end of file
Modified: branches/robast-0.9/pkg/RobExtremes/R/GParetoFamily.R
===================================================================
--- branches/robast-0.9/pkg/RobExtremes/R/GParetoFamily.R 2013-01-17 00:40:40 UTC (rev 537)
+++ branches/robast-0.9/pkg/RobExtremes/R/GParetoFamily.R 2013-01-17 00:41:38 UTC (rev 538)
@@ -76,169 +76,113 @@
names(theta) <- c("loc", "scale", "shape")
scaleshapename <- c("scale", "shape")
+ if(!is.null(p)){
+ btq <- substitute({ q <- loc0 + theta[1]*((1-p0)^(-theta[2])-1)/theta[2]
+ names(q) <- "quantile"
+ }, list(loc0 = loc, p0 = p))
+
+ bDq <- substitute({ scale <- theta[1]; shape <- theta[2]
+ D1 <- ((1-p0)^(-shape)-1)/shape
+ D2 <- -scale/shape*(D1 + log(1-p0)*(1-p0)^(-shape))
+ D <- t(c(D1, D2))
+ rownames(D) <- "quantile"; colnames(D) <- NULL
+ D }, list(p0 = p))
+ btes <- substitute({ if(theta[2]>=1L) es <- NA else {
+ q <- loc0 + theta[1]*((1-p0)^(-theta[2])-1)/theta[2]
+ es <- (q + theta[1] - theta[2]*loc0)/(1-theta[2])}
+ names(es) <- "expected shortfall"
+ es }, list(loc0 = loc, p0 = p))
+ bDes <- substitute({ if(theta[2]>=1L){ D1 <- D2 <- NA}else{
+ scale <- theta[1]; shape <- theta[2]
+ q <- loc0 + theta[1]*((1-p0)^(-theta[2])-1)/theta[2]
+ dq1 <- ((1-p0)^(-shape)-1)/shape
+ dq2 <- -scale/shape*(dq1 + log(1-p0)*(1-p0)^(-shape))
+ D1 <- (dq1 + 1)/(1-shape)
+ D2 <- (dq2 - loc0)/(1-shape) + (q + scale -
+ loc0*shape)/(1-shape)^2}
+ D <- t(c(D1, D2))
+ rownames(D) <- "expected shortfall"
+ colnames(D) <- NULL
+ D }, list(loc0 = loc, p0 = p))
+ }
+ if(!is.null(N)){
+ btel <- substitute({ if(theta[2]>=1L) el <- NA else {
+ el <- N0*(loc0 + theta[1]/(1-theta[2]))}
+ names(el) <- "expected loss"
+ el }, list(loc0 = loc,N0 = N))
+ bDel <- substitute({ if(theta[2]>=1L){ D1 <- D2 <- NA}else{
+ scale <- theta[1]; shape <- theta[2]
+ D1 <- N0/(1-shape)
+ D2 <- D1*scale/(1-shape)}
+ D <- t(c(D1, D2))
+ rownames(D) <- "expected loss"
+ colnames(D) <- NULL
+ D }, list(loc0 = loc, N0 = N))
+ }
if(is.null(trafo)){
tau <- NULL
if("scale" %in% of.interest){
- tau <- function(theta){
- th <- theta[1]
- names(th) <- "scale"
- th
- }
- Dtau <- function(theta){
- D <- t(c(1, 0))
- rownames(D) <- "scale"
- D
- }
+ tau <- function(theta){ th <- theta[1]; names(th) <- "scale"; th}
+ Dtau <- function(theta){ D <- t(c(1, 0)); rownames(D) <- "scale"; D}
}
if("shape" %in% of.interest){
if(is.null(tau)){
- tau <- function(theta){
- th <- theta[2]
- names(th) <- "shape"
- th
- }
- Dtau <- function(theta){
- D <- t(c(0, 1))
- rownames(D) <- "shape"
- D
- }
+ tau <- function(theta){th <- theta[2]; names(th) <- "shape"; th}
+ Dtau <- function(theta){D <- t(c(0,1));rownames(D) <- "shape";D}
}else{
- tau <- function(theta){
- th <- theta
- names(th) <- c("scale", "shape")
- th
+ tau <- function(theta){th <- theta
+ names(th) <- c("scale", "shape"); th}
+ Dtau <- function(theta){ D <- diag(2);
+ rownames(D) <- c("scale", "shape");D}
}
- Dtau <- function(theta){
- D <- diag(2)
- rownames(D) <- c("scale", "shape")
- D
- }
}
}
if("quantile" %in% of.interest){
if(is.null(p)) stop("Probability 'p' has to be specified.")
if(is.null(tau)){
- tau <- function(theta){ }
- body(tau) <- substitute({ q <- loc0 + theta[1]*((1-p0)^(-theta[2])-1)/theta[2]
- names(q) <- "quantile"
- q },
- list(loc0 = loc, p0 = p))
- Dtau <- function(theta){ }
- body(Dtau) <- substitute({ scale <- theta[1]
- shape <- theta[2]
- D1 <- ((1-p0)^(-shape)-1)/shape
- D2 <- -scale/shape*(D1 + log(1-p0)*(1-p0)^(-shape))
- D <- t(c(D1, D2))
- rownames(D) <- "quantile"
- colnames(D) <- NULL
- D },
- list(p0 = p))
+ tau <- function(theta){ }; body(tau) <- btq
+ Dtau <- function(theta){ };body(Dtau) <- bDq
}else{
tau1 <- tau
tau <- function(theta){ }
- body(tau) <- substitute({ q <- loc0 + theta[1]*((1-p0)^(-theta[2])-1)/theta[2]
- names(q) <- "quantile"
- c(tau0(theta), q) },
- list(tau0 = tau1, loc0 = loc, p0 = p))
+ body(tau) <- substitute({ btq0; c(tau0(theta), q) },
+ list(btq0=btq, tau0 = tau1))
Dtau1 <- Dtau
Dtau <- function(theta){}
- body(Dtau) <- substitute({ scale <- theta[1]
- shape <- theta[2]
- D1 <- ((1-p0)^(-shape)-1)/shape
- D2 <- -scale/shape*(D1 + log(1-p0)*(1-p0)^(-shape))
- D <- t(c(D1, D2))
- rownames(D) <- "quantile"
- colnames(D) <- NULL
- rbind(Dtau0(theta), D) },
- list(Dtau0 = Dtau1, p0 = p))
+ body(Dtau) <- substitute({ bDq0; rbind(Dtau0(theta), D) },
+ list(Dtau0 = Dtau1, bDq0 = bDq))
}
}
if("expected shortfall" %in% of.interest){
if(is.null(p)) stop("Probability 'p' has to be specified.")
if(is.null(tau)){
- tau <- function(theta){ }
- body(tau) <- substitute({ q <- loc0 + theta[1]*((1-p0)^(-theta[2])-1)/theta[2]
- es <- (q + theta[1] - theta[2]*loc0)/(1-theta[2])
- names(es) <- "expected shortfall"
- es },
- list(loc0 = loc, p0 = p))
- Dtau <- function(theta){ }
- body(Dtau) <- substitute({ scale <- theta[1]
- shape <- theta[2]
- q <- loc0 + theta[1]*((1-p0)^(-theta[2])-1)/theta[2]
- dq1 <- ((1-p0)^(-shape)-1)/shape
- dq2 <- -scale/shape*(dq1 + log(1-p0)*(1-p0)^(-shape))
- D1 <- (dq1 + 1)/(1-shape)
- D2 <- (dq2 - loc0)/(1-shape) + (q + scale - loc0*shape)/(1-shape)^2
- D <- t(c(D1, D2))
- rownames(D) <- "expected shortfall"
- colnames(D) <- NULL
- D },
- list(loc0 = loc, p0 = p))
+ tau <- function(theta){ }; body(tau) <- btes
+ Dtau <- function(theta){ }; body(Dtau) <- bDes
}else{
tau1 <- tau
tau <- function(theta){ }
- body(tau) <- substitute({ q <- loc0 + theta[1]*((1-p0)^(-theta[2])-1)/theta[2]
- es <- (q + theta[1] - theta[2]*loc0)/(1-theta[2])
- names(es) <- "expected shortfall"
- c(tau0(theta), es) },
- list(tau0 = tau1, loc0 = loc, p0 = p))
+ body(tau) <- substitute({ btes0; c(tau0(theta), es) },
+ list(tau0 = tau1, btes0=btes))
Dtau1 <- Dtau
Dtau <- function(theta){}
- body(Dtau) <- substitute({ scale <- theta[1]
- shape <- theta[2]
- q <- loc0 + theta[1]*((1-p0)^(-theta[2])-1)/theta[2]
- dq1 <- ((1-p0)^(-shape)-1)/shape
- dq2 <- -scale/shape*(dq1 + log(1-p0)*(1-p0)^(-shape))
- D1 <- (dq1 + 1)/(1-shape)
- D2 <- (dq2 - loc0)/(1-shape) + (q + scale - loc0*shape)/(1-shape)^2
- D <- t(c(D1, D2))
- rownames(D) <- "expected shortfall"
- colnames(D) <- NULL
- rbind(Dtau0(theta), D) },
- list(Dtau0 = Dtau1, loc0 = loc, p0 = p))
+ body(Dtau) <- substitute({ bDes0; rbind(Dtau0(theta), D) },
+ list(Dtau0 = Dtau1, bDes0=bDes))
}
}
if("expected loss" %in% of.interest){
if(is.null(N)) stop("Expected frequency 'N' has to be specified.")
if(is.null(tau)){
- tau <- function(theta){ }
- body(tau) <- substitute({ el <- N0*(loc0 + theta[1]*gamma(1/theta[2]-1)/(theta[2]^2*gamma(1/theta[2]+1)))
- names(el) <- "expected loss"
- el },
- list(loc0 = loc,N0 = N))
- Dtau <- function(theta){ }
- body(Dtau) <- substitute({ scale <- theta[1]
- shape <- theta[2]
- Gneg <- gamma(1/shape-1)
- Gpos <- gamma(1/shape+1)
- D1 <- N0*Gneg/(shape^2*Gpos)
- D2 <- N0*scale*Gneg*(digamma(1/shape+1) - 2*shape - digamma(1/shape-1))/(shape^4*Gpos)
- D <- t(c(D1, D2))
- rownames(D) <- "expected loss"
- colnames(D) <- NULL
- D },
- list(loc0 = loc, N0 = N))
+ tau <- function(theta){ }; body(tau) <- btel
+ Dtau <- function(theta){ }; body(Dtau) <- bDel
}else{
tau1 <- tau
tau <- function(theta){ }
- body(tau) <- substitute({ el <- N0*(loc0 + theta[1]*gamma(1/theta[2]-1)/(theta[2]^2*gamma(1/theta[2]+1)))
- names(el) <- "expected loss"
- c(tau0(theta), el) },
- list(tau0 = tau1, loc0 = loc,N0 = N))
+ body(tau) <- substitute({ btel0; c(tau0(theta), el) },
+ list(tau0 = tau1, btel0=btel))
Dtau1 <- Dtau
Dtau <- function(theta){}
- body(Dtau) <- substitute({ scale <- theta[1]
- shape <- theta[2]
- Gneg <- gamma(1/shape-1)
- Gpos <- gamma(1/shape+1)
- D1 <- N0*Gneg/(shape^2*Gpos)
- D2 <- N0*scale*Gneg*(digamma(1/shape+1) - 2*shape - digamma(1/shape-1))/(shape^4*Gpos)
- D <- t(c(D1, D2))
- rownames(D) <- "expected loss"
- colnames(D) <- NULL
- rbind(Dtau0(theta), D) },
- list(Dtau0 = Dtau1, loc0 = loc, N0 = N))
+ body(Dtau) <- substitute({ bDel0; rbind(Dtau0(theta), D) },
+ list(Dtau0 = Dtau1, bDel0=bDel))
}
}
trafo <- function(x){ list(fval = tau(x), mat = Dtau(x)) }
@@ -257,8 +201,6 @@
startPar <- function(x,...){
tr <- theta[1]
- if(any(x < tr)) stop("some data smaller than 'loc' parameter")
-
## Pickand estimator
if(is.null(start0Est)){
e0 <- estimate(medkMADhybr(x, k=10, ParamFamily=GParetoFamily(loc = theta[1],
@@ -272,6 +214,10 @@
if(!is.null(names(e0)))
e0 <- e0[c("scale", "shape")]
}
+
+ if(any(x < tr-e0["scale"]/e0["shape"]))
+ stop("some data smaller than 'loc-scale/shape' ")
+
names(e0) <- NULL
return(e0)
}
@@ -310,16 +256,18 @@
Lambda1 <- function(x) {
y <- x*0
- x0 <- (x-tr)/sc
- x1 <- x0[x0>0]
- y[x0>0] <- -1/sc + (1+k)/(1+k*x1)*x1/sc
+ ind <- (x > tr-sc/k) # = [later] (x1>0)
+ x <- (x[ind]-tr)/sc
+ x1 <- 1 + k * x
+ y[ind] <- -1/sc + (1+k)/x1*x/sc
return(y)
}
Lambda2 <- function(x) {
y <- x*0
- x0 <- (x-tr)/sc
- x1 <- x0[x0>0]
- y[x0>0] <- log(1+k*x1)/k^2 - (1/k+1)*x1/(1+k*x1)
+ ind <- (x > tr-sc/k) # = [later] (x1>0)
+ x <- (x[ind]-tr)/sc
+ x1 <- 1 + k * x
+ y[ind] <- log(x1)/k^2 - (1/k+1)*x/x1
return(y)
}
## additional centering of scores to increase numerical precision!
@@ -372,7 +320,7 @@
of.interest0 = of.interest, p0 = p, N0 = N,
trafo0 = trafo, withPos0 = withPos))
- L2Fam at LogDeriv <- function(x) (shape+1)/(shape*(scale+(x-loc)))
+ L2Fam at LogDeriv <- function(x) (shape+1)/(scale+shape*(x-loc))
L2Fam at L2deriv <- L2deriv
L2Fam at L2derivDistr <- imageDistr(RandVar = L2deriv, distr = distribution)
Modified: branches/robast-0.9/pkg/RobExtremes/R/LDEstimator.R
===================================================================
--- branches/robast-0.9/pkg/RobExtremes/R/LDEstimator.R 2013-01-17 00:40:40 UTC (rev 537)
+++ branches/robast-0.9/pkg/RobExtremes/R/LDEstimator.R 2013-01-17 00:41:38 UTC (rev 538)
@@ -22,7 +22,7 @@
loc.fctal.0, disp.fctal.0, ParamFamily.0,
loc.est.ctrl.0 = NULL, loc.fctal.ctrl.0=NULL,
disp.est.ctrl.0 = NULL, disp.fctal.ctrl.0=NULL,
- q.lo.0 =0, q.up.0=Inf, log.q.0 =TRUE, ...
+ q.lo.0 =0, q.up.0=Inf, log.q.0 =TRUE, ..., vdbg=FALSE
){
dots <- list(...)
loc.emp <- do.call(loc.est.0, args = .prepend(x.0,loc.est.ctrl.0, dots))
@@ -34,6 +34,7 @@
sc.th <- do.call(disp.fctal.0, args = .prepend(distr.new,disp.fctal.ctrl.0, dots))
val <- if(log.q.0) log(loc.th)-log(sc.th) - q.emp else
loc.th/sc.th-q.emp
+ if(vdbg) print(val)
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/robast -r 538
More information about the Robast-commits
mailing list