[Robast-commits] r541 - in branches/robast-0.9/pkg/RobExtremes: . R man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Fri Jan 18 20:07:09 CET 2013
Author: ruckdeschel
Date: 2013-01-18 20:07:09 +0100 (Fri, 18 Jan 2013)
New Revision: 541
Added:
branches/robast-0.9/pkg/RobExtremes/R/QBCC.R
branches/robast-0.9/pkg/RobExtremes/man/GEVFamily.Rd
branches/robast-0.9/pkg/RobExtremes/man/QuantileBCCEstimator.Rd
branches/robast-0.9/pkg/RobExtremes/man/asvarQBCC.Rd
Modified:
branches/robast-0.9/pkg/RobExtremes/NAMESPACE
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/PickandsEstimator.R
branches/robast-0.9/pkg/RobExtremes/R/WeibullFamily.R
branches/robast-0.9/pkg/RobExtremes/R/asvarPickands.R
branches/robast-0.9/pkg/RobExtremes/man/0RobExtremes-package.Rd
branches/robast-0.9/pkg/RobExtremes/man/GParetoFamily.Rd
branches/robast-0.9/pkg/RobExtremes/man/PickandsEstimator.Rd
branches/robast-0.9/pkg/RobExtremes/man/WeibullFamily.Rd
Log:
RobExtremes: implemented the Quantile estimator of Boudt Caliskan Croux for the Weibull family as starting estimator (together with asvarQBCC..); updated Metrika reference. Corrected some errors in
the scores function for the Weibull family; explicit terms for the FI in the Weibull model. PickandsEstimator got a corrected asVar (due to the new computation for beta, increasing the bdp;
mentioned the increase of BDP in the Rd file.
Modified: branches/robast-0.9/pkg/RobExtremes/NAMESPACE
===================================================================
--- branches/robast-0.9/pkg/RobExtremes/NAMESPACE 2013-01-17 21:23:32 UTC (rev 540)
+++ branches/robast-0.9/pkg/RobExtremes/NAMESPACE 2013-01-18 19:07:09 UTC (rev 541)
@@ -28,5 +28,7 @@
export("Gumbel", "Pareto", "GPareto", "GEV")
export("GParetoFamily", "GumbelLocationFamily", "WeibullFamily")
export("LDEstimator", "medkMAD", "medSn", "medQn", "medkMADhybr")
-export("getShapeGrid", "getSnGrid", "PickandsEstimator")
-export("loc", "loc<-", "kMAD", "Sn", "Qn", "asvarMedkMAD","asvarPickands")
\ No newline at end of file
+export("getShapeGrid", "getSnGrid",
+ "PickandsEstimator","QuantileBCCEstimator")
+export("loc", "loc<-", "kMAD", "Sn", "Qn",
+ "asvarMedkMAD","asvarPickands", "asvarQBCC")
\ No newline at end of file
Modified: branches/robast-0.9/pkg/RobExtremes/R/GEVFamily.R
===================================================================
--- branches/robast-0.9/pkg/RobExtremes/R/GEVFamily.R 2013-01-17 21:23:32 UTC (rev 540)
+++ branches/robast-0.9/pkg/RobExtremes/R/GEVFamily.R 2013-01-18 19:07:09 UTC (rev 541)
@@ -81,20 +81,20 @@
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 }
+ pg <- pgamma(-log(p0),1-theta[2], lower.tail = TRUE)
+ es <- theta[1] * (gamma(1-theta[2]) * pg/ (1-p0) - 1 )/
+ 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)
+ pg <- pgamma(-log(p0), 1-theta[2], lower.tail = TRUE)
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}
+ g0 <- gamma(1-theta[2])
+ D1 <- (g0*pg/(1-p0)-1)/theta[2]
+ D21 <- theta[1]*D1/theta[2]
+ D22 <- theta[1]*dd/(1-p0)/theta[2]
+ D2 <- -D21+D22)}
D <- t(c(D1, D2))
rownames(D) <- "expected shortfall"
colnames(D) <- NULL
@@ -221,15 +221,28 @@
## 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])
+ theta <- abs(theta)
+ }else{
+ if(!is.null(names(theta))){
+ theta["scale"] <- abs(theta["scale"])
+ }else{
+ theta[1] <- abs(theta[1])
+ }
}
return(theta)
}
modifyPar <- function(theta){
- theta <- abs(theta)
+ theta <- makeOKPar(theta)
+
+ if(!is.null(names(theta))){
+ sc <- theta["scale"]
+ sh <- theta["shape"]
+ }else{
+ theta <- abs(theta)
+ sc <- theta[1]
+ sh <- theta[2]
+ }
GEV(loc = loc, scale = theta[1], shape = theta[2])
}
@@ -272,8 +285,8 @@
G20 <- gamma(2*k)
G10 <- gamma(k)
G11 <- digamma(k)*gamma(k)
- G01 <- digamma(1)
- G02 <- trigamma(1)+digamma(1)^2
+ G01 <- -0.57721566490153 # digamma(1)
+ G02 <- 1.9781119906559 #trigamma(1)+digamma(1)^2
x0 <- (k+1)^2*2*k
I11 <- G20*x0-2*G10*k*(k+1)+1
I11 <- I11/sc^2/k^2
@@ -284,7 +297,6 @@
I22 <- I22 - G11*2*k^2*(k+1) + G01*2*k*(1+k)+k^2 *G02
I22 <- I22 /k^4
mat <- PosSemDefSymmMatrix(matrix(c(I11,I12,I12,I22),2,2))
- mat <- PosSemDefSymmMatrix(matrix(c(I11,I12,I12,I22),2,2))
dimnames(mat) <- list(scaleshapename,scaleshapename)
return(mat)
}
@@ -337,6 +349,6 @@
ddigamma <- function(t,s){
int <- function(x) exp(-x)*(-log(x))*x^(-s)
- integrate(int, lower=t, upper=Inf)$value
+ integrate(int, lower=0, upper=t)$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 21:23:32 UTC (rev 540)
+++ branches/robast-0.9/pkg/RobExtremes/R/GParetoFamily.R 2013-01-18 19:07:09 UTC (rev 541)
@@ -225,9 +225,13 @@
## 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])
+ theta <- abs(theta)
+ }else{
+ if(!is.null(names(theta))){
+ theta["scale"] <- abs(theta["scale"])
+ }else{
+ theta[1] <- abs(theta[1])
+ }
}
return(theta)
}
@@ -255,7 +259,7 @@
Lambda1 <- function(x) {
y <- x*0
- ind <- (x > tr-sc/k) # = [later] (x1>0)
+ ind <- (x > tr) #
x <- (x[ind]-tr)/sc
x1 <- 1 + k * x
y[ind] <- -1/sc + (1+k)/x1*x/sc
@@ -263,7 +267,7 @@
}
Lambda2 <- function(x) {
y <- x*0
- ind <- (x > tr-sc/k) # = [later] (x1>0)
+ ind <- (x > tr) #
x <- (x[ind]-tr)/sc
x1 <- 1 + k * x
y[ind] <- log(x1)/k^2 - (1/k+1)*x/x1
Modified: branches/robast-0.9/pkg/RobExtremes/R/PickandsEstimator.R
===================================================================
--- branches/robast-0.9/pkg/RobExtremes/R/PickandsEstimator.R 2013-01-17 21:23:32 UTC (rev 540)
+++ branches/robast-0.9/pkg/RobExtremes/R/PickandsEstimator.R 2013-01-18 19:07:09 UTC (rev 541)
@@ -10,21 +10,41 @@
I <- ms[2]
m <- ms[1]
xi <- log((I-m)/m)/log(alpha)
+
+ ##############
+ ###
+ ### the scale estimate we use, i.e. beta <- xi*m/(alpha^xi-1)
+ ### differs from the one given in the original reference
+ ### ---which was: beta <- xi * m^2 /(I-2*m) leading to xi-dependent bdp---
+ ### the one chosen here avoids taking differences I - 2m hence does not
+ ### require I > 2m; this leads to (functional) breakdown point
+ ### min(a1,1-a2,a2-a1)
+ ### note that this value is independent of xi !!
+ ### for GPD the optimal choice of alpha is 2 leading to bdp 1/4
+ ### for GEVD the optimal choice of alpha is 2.248 leading to bdp 0.180
+ ### (and the standard choice alpha = 2 leads to bdp 0.172)
+ ### for comparison: with the original scale estimate, at xi=0.7, this
+ ### gives optimal bdp's 0.060 (GEVD) and 0.070 (GPD),
+ ### resp. bdp's 0.048 (GEVD) and 0.064 (GPD) for alpha = 2
+ ###
+ #############
+
beta <- xi*m/(alpha^xi-1)
+ ###
theta <- c(beta,xi)
names(theta) <- c("scale","shape")
return(theta)
}
PickandsEstimator <- function(x, alpha = 2, ParamFamily=GParetoFamily(),
- name, Infos, asvar = NULL, nuis.idx = NULL,
- trafo = NULL, fixed = NULL, asvar.fct = NULL, na.rm = TRUE,
+ name, Infos, nuis.idx = NULL,
+ trafo = NULL, fixed = NULL, na.rm = TRUE,
...){
isGP <- is(ParamFamily,"GParetoFamily")
if(!(isGP|is(ParamFamily,"GEVFamily")))
stop("Pickands estimator only available for GPD and GEVD.")
es.call <- match.call()
- if(missing(alpha)) alpha <- 2
+ if(missing(alpha)) alpha <- if(isGP) 2 else 2.248
if(length(alpha)>1 || any(!is.finite(alpha)) || any(alpha<=1))
stop("'alpha' has to be a numeric > 1 of length 1.")
@@ -34,7 +54,6 @@
asvar.fct.0 <- function(L2Fam=ParamFamily, param){
asvarPickands(model=L2Fam, alpha = alpha)}
- asvar <- asvarPickands(model=ParamFamily, alpha = alpha)
nuis.idx.0 <- nuis.idx
trafo.0 <- trafo
fixed.0 <- fixed
@@ -42,7 +61,7 @@
.mPick <- function(x) .PickandsEstimator(x,alpha=alpha, GPD.l=isGP)
estimate <- Estimator(x, .mPick, name, Infos,
- asvar.fct = asvar.fct0, asvar = asvar,
+ asvar.fct = asvar.fct0 asvar = asvar,
nuis.idx = nuis.idx.0, trafo = trafo.0,
fixed = fixed.0, na.rm = na.rm.0, ...)
##->
Added: branches/robast-0.9/pkg/RobExtremes/R/QBCC.R
===================================================================
--- branches/robast-0.9/pkg/RobExtremes/R/QBCC.R (rev 0)
+++ branches/robast-0.9/pkg/RobExtremes/R/QBCC.R 2013-01-18 19:07:09 UTC (rev 541)
@@ -0,0 +1,64 @@
+#####################################################################
+# Robust starting estimator for Weibull model
+# acc. Boudt Caliskan Croux (2011)
+# for p1 < p2
+# has bdp min(p1,1-p2,p2-p1), maximal for p1=p2=1/3
+#
+#####################################################################
+.QBCC <- function(x, p1 = 1/3, p2 =1/3){
+ if(p1>=p2) {p<-p1; p1 <- p2; p2 <- p}
+ l1 <- -log(p1); l2 <- -log(p2)
+ ms <- quantile(x,c(p1,p2))
+ names(ms) <- NULL
+ Q2 <- ms[2]
+ Q1 <- ms[1]
+ xi <- (log(Q2)-log(Q1))/(log(l2)-log(l1))
+ beta <- Q2/l2^(1/xi)
+ ###
+ theta <- c(beta,xi)
+ names(theta) <- c("scale","shape")
+ return(theta)
+}
+
+QuantileBCCEstimator <- function(x, p1=1/3, p2=1/3, ParamFamily=WeibullFamily(),
+ name, Infos, nuis.idx = NULL,
+ trafo = NULL, fixed = NULL, na.rm = TRUE,
+ ...){
+ if(!(is(ParamFamily,"WeibullFamily")))
+ stop("Pickands estimator only available for Weibull")
+ es.call <- match.call()
+ if(length(p1)>1 || any(!is.finite(p1)) || p1<=0 || p1>=1)
+ stop("'p1' has to be in [0,1] and of length 1.")
+ if(length(p2)>1 || any(!is.finite(p2)) || p2<=0 || p2>=1 || abs(p1-p2)< 1e-8)
+ stop("'p2' has to be in [0,1] and of length 1 and distinct of 'p1'.")
+
+ if(missing(name))
+ name <- "QuantileBCCEstimator"
+
+
+ asvar.fct.0 <- function(L2Fam=ParamFamily, param){
+ asvarQBCC(model=L2Fam, p1 = p1, p2 = p2)}
+ nuis.idx.0 <- nuis.idx
+ trafo.0 <- trafo
+ fixed.0 <- fixed
+ na.rm.0 <- na.rm
+
+ .mQBCC <- function(x) .QBCC(x,p1=p1,p2=p2)
+ estimate <- Estimator(x, .mQBCC, name, Infos,
+ asvar.fct = asvar.fct0 asvar = asvar,
+ nuis.idx = nuis.idx.0, trafo = trafo.0,
+ fixed = fixed.0, na.rm = na.rm.0, ...)
+ estimate at estimate.call <- es.call
+
+ if(missing(Infos))
+ Infos <- matrix(c(name, ""),
+ ncol=2, dimnames=list(character(0), c("method", "message")))
+ else{
+ Infos <- matrix(c(rep(name, length(Infos)+1), c("",Infos)),
+ ncol = 2)
+ colnames(Infos) <- c("method", "message")
+ }
+ estimate at Infos <- Infos
+
+ return(estimate)
+}
Modified: branches/robast-0.9/pkg/RobExtremes/R/WeibullFamily.R
===================================================================
--- branches/robast-0.9/pkg/RobExtremes/R/WeibullFamily.R 2013-01-17 21:23:32 UTC (rev 540)
+++ branches/robast-0.9/pkg/RobExtremes/R/WeibullFamily.R 2013-01-18 19:07:09 UTC (rev 541)
@@ -74,171 +74,122 @@
## parameters
names(theta) <- c("scale", "shape")
+ if(!is.null(p)){
+ btq <- substitute({ q <- theta[1]*(-log(1-p0))^(1/theta[2])
+ names(q) <- "quantile"
+ }, list(loc0 = loc, p0 = p))
+
+ bDq <- substitute({ scale <- theta[1]; shape <- theta[2]
+ lp <- -log(1-p0)
+ D1 <- lp^(1/shape)
+ D2 <- -scale/shape^2 *lp^(1/shape)*log(lp)
+ D <- t(c(D1, D2))
+ rownames(D) <- "quantile"; colnames(D) <- NULL
+ D }, list(p0 = p))
+ btes <- substitute({ if(theta[2]>=1L) es <- NA else {
+ s1 <- 1+1/theta[2]
+ pg <- pgamma(-log(p0),s1, lower.tail = FALSE)
+ g0 <- gamma(s1)
+ es <- theta[1] * g0 * pg /(1-p0) + loc0 }
+ names(es) <- "expected shortfall"
+ es }, list(loc0 = loc, 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)
+ g0 <- gamma(s1)
+ dd <- digamma(s1)*g0 - ddigamma(-log(p0),s1)
+ D1 <- g0 * pg / (1-p0)
+ D2 <- theta[1] * dd /(1-p0)}
+ D <- t(c(D1, D2))
+ rownames(D) <- "expected shortfall"
+ colnames(D) <- NULL
+ D }, list(loc0 = loc, 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))
+ bDel <- substitute({ scale <- theta[1]; shape <- theta[2]
+ s1 <- 1+1/shape
+ D1 <- N0*gamma(s1)
+ D2 <- -N0*theta[1]*digamma(s1)*gamma(s1)/shape^2
+ 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
- }
- Dtau <- function(theta){
- D <- diag(2)
- rownames(D) <- c("scale", "shape")
- D
- }
+ 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) <- substitute({ q <- exp(shape^(-1)*log(-log(1-p))+log(scale))
- names(q) <- "quantile"
- q },
- list(p0 = p))
- Dtau <- function(theta){ }
- body(Dtau) <- substitute({ scale <- theta[1]
- shape <- theta[2]
- D1 <- exp(log(-log(1-p))/shape^2+log(scale))/scale
- D2 <- (-log(-log(1-p)))*exp(log(-log(1-p))/shape^2+log(scale))/shape^2
- 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 <- exp(theta[2]^(-1)*log(-log(1-p))+log(theta[1]))
- names(q) <- "quantile"
- c(tau0(theta), q) },
- list(tau0 = tau1, 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 <- exp(log(-log(1-p))/shape^2+log(scale))/scale
- D2 <- (-log(-log(1-p)))*exp(log(-log(1-p))/shape^2+log(scale))/shape^2
- 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 <- exp(theta[2]^(-1)*log(-log(1-p))+log(theta[1]))
- es <- E(q,upp=Inf,low=p0)/(1-p0)
- names(es) <- "expected shortfall"
- es },
- list(p0 = p))
- Dtau <- function(theta){ }
- body(Dtau) <- substitute({ scale <- theta[1]
- shape <- theta[2]
- q <- exp(shape^(-1)*log(-log(1-p))+log(scale))
- dq1 <- exp(log(-log(1-p))/shape^2+log(scale))/scale
- dq2 <- (-log(-log(1-p)))*exp(log(-log(1-p))/shape^2+log(scale))/shape^2
- D1 <- E(dq1,upp=Inf,low=p0)/(1-p0)
- D2 <- E(dq2,upp=Inf,low=p0)/(1-p0)
- D <- t(c(D1, D2))
- rownames(D) <- "expected shortfall"
- colnames(D) <- NULL
- D },
- list(p0 = p))
+ tau <- function(theta){ }; body(tau) <- btes
+ Dtau <- function(theta){ }; body(Dtau) <- bDes
}else{
tau1 <- tau
tau <- function(theta){ }
- body(tau) <- substitute({ q <- exp(theta[2]^(-1)*log(-log(1-p))+log(theta[1]))
- es <- E(q,upp=Inf,low=p0)/(1-p0)
- names(es) <- "expected shortfall"
- c(tau0(theta), es) },
- list(tau0 = tau1, 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 <- exp(shape^(-1)*log(-log(1-p))+log(scale))
- dq1 <- exp(log(-log(1-p))/shape^2+log(scale))/scale
- dq2 <- (-log(-log(1-p)))*exp(log(-log(1-p))/shape^2+log(scale))/shape^2
- D1 <- E(dq1,upp=Inf,low=p0)/(1-p0)
- D2 <- E(dq2,upp=Inf,low=p0)/(1-p0)
- D <- t(c(D1, D2))
- rownames(D) <- "expected shortfall"
- colnames(D) <- NULL
- rbind(Dtau0(theta), D) },
- list(Dtau0 = Dtau1, 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*gamma(1+1/shape)*scale
- names(el) <- "expected loss"
- el },
- list(N0 = N))
- Dtau <- function(theta){ }
- body(Dtau) <- substitute({ scale <- theta[1]
- shape <- theta[2]
- D1 <- N0*gamma(1/xi+1)
- D2 <- -N0*digamma(1/xi+1)*beta/xi^2
- D <- t(c(D1, D2))
- rownames(D) <- "expected loss"
- colnames(D) <- NULL
- D },
- list(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*gamma(1+1/shape)*scale
- names(el) <- "expected loss"
- c(tau0(theta), el) },
- list(tau0 = tau1, 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]
- D1 <- N0*gamma(1/xi+1)
- D2 <- -N0*digamma(1/xi+1)*beta/xi^2
- D <- t(c(D1, D2))
- rownames(D) <- "expected loss"
- colnames(D) <- NULL
- rbind(Dtau0(theta), D) },
- list(Dtau0 = Dtau1, 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)) }
}else{
if(is.matrix(trafo) & nrow(trafo) > 2) stop("number of rows of 'trafo' > 2")
}
+
+
param <- ParamFamParameter(name = "theta", main = c(theta[1],theta[2]),
trafo = trafo, withPosRestr = withPos,
.returnClsName ="ParamWithScaleAndShapeFamParameter")
@@ -252,8 +203,7 @@
## Pickand estimator
if(is.null(start0Est)){
- e0 <- estimate(medkMADhybr(x, k=10, ParamFamily=WeibullFamily(scale = theta[1], shape = theta[2]),
- q.lo = 1e-3, q.up = 15))
+ e0 <- estimate(QuantileBCCEstimator(x))
}else{
if(is(start0Est,"function")){
e1 <- start0Est(x, ...)
@@ -270,9 +220,13 @@
## 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])
+ theta <- abs(theta)
+ }else{
+ if(!is.null(names(theta))){
+ theta["scale"] <- abs(theta["scale"])
+ }else{
+ theta[1] <- abs(theta[1])
+ }
}
return(theta)
}
@@ -300,13 +254,15 @@
Lambda1 <- function(x) {
y <- x*0
x1 <- x[x>0]
- y[x>0] <- -1/sc-(sh-1)/sc+sh*x1/sc^2###
+ z <- x1/sc
+ y[x>0] <- sh/sc*(z^sh-1)###
return(y)
}
Lambda2 <- function(x) {
y <- x*0
x1 <- x[x>0]
- y[x>0] <- 1/sh+log(x1/sc)-x1/sc###
+ z <- x1/sc
+ y[x>0] <- 1/sh-log(z)*(z^sh-1)###
return(y)
}
## additional centering of scores to increase numerical precision!
@@ -316,16 +272,22 @@
}
## Fisher Information matrix as a function of parameters
+ ## Fisher Information matrix as a function of parameters
FisherInfo.fct <- function(param) {
sc <- force(main(param)[1])
- sh <- force(main(param)[2])
- Lambda <- L2deriv.fct(param)
- E11 <- E(distribution,fun=function(x)Lambda[[1]](x)^2)
- E12 <- E(distribution,fun=function(x)Lambda[[1]](x)*Lambda[[2]](x))
- E22 <- E(distribution,fun=function(x)Lambda[[2]](x)^2)
- return(PosSemDefSymmMatrix(matrix(c(E11,E12,E12,E22),2,2)))
+ k <- force(main(param)[2])
+ g1 <- 0.42278433509847 # 1+digamma(1)
+ g2 <- 1.8236806608529 # 1+trigamma(1)+digamma(1)^2+2*digamma(1)
+ I11 <- k^2/sc^2
+ I12 <- -g1/sc
+ I22 <- g2/k^2
+ mat <- PosSemDefSymmMatrix(matrix(c(I11,I12,I12,I22),2,2))
+ dimnames(mat) <- list(scaleshapename,scaleshapename)
+ return(mat)
}
+
+
FisherInfo <- FisherInfo.fct(param)
name <- "Weibull Family"
@@ -354,7 +316,9 @@
of.interest0 = of.interest, p0 = p, N0 = N,
trafo0 = trafo, withPos0 = withPos))
- L2Fam at LogDeriv <- function(x) log(shape)-log(scale)+(shape-1)*log(x/scale)-shape*x/scale###
+ L2Fam at LogDeriv <- function(x){ z <- x/scale
+ log(shape)-log(scale)+(shape-1)*log(z)-shape*z^(shape-1)
+ }###
L2Fam at L2deriv <- L2deriv
L2Fam at L2derivDistr <- imageDistr(RandVar = L2deriv, distr = distribution)
Modified: branches/robast-0.9/pkg/RobExtremes/R/asvarPickands.R
===================================================================
--- branches/robast-0.9/pkg/RobExtremes/R/asvarPickands.R 2013-01-17 21:23:32 UTC (rev 540)
+++ branches/robast-0.9/pkg/RobExtremes/R/asvarPickands.R 2013-01-18 19:07:09 UTC (rev 541)
@@ -1,12 +1,12 @@
asvarPickands <- function(model, alpha=2){
- isGP <- is(ParamFamily,"GParetoFamily")
- if(!(isGP|is(ParamFamily,"GEVFamily")))
+ isGP <- is(model,"GParetoFamily")
+ if(!(isGP|is(model,"GEVFamily")))
stop("Pickands estimator only available for GPD and GEVD.")
scshn <- scaleshapename(model)
- par0 <- main(model at param)[scshn]
- beta <- par0[1]; xi <- par0[2]
+# par0 <- main(model at param)[scshn]
+# beta <- par0[1]; xi <- par0[2]
if(isGP){
al1 <- 1-1/alpha
@@ -15,17 +15,36 @@
al1 <- exp(-1/alpha)
al2 <- exp(-1/alpha^2)
}
+
M2 <- q(model)(al1)
M4 <- q(model)(al2)
+ xi <- log((M4-M2)/M2)/log(alpha)
+ qu <- 1/(alpha^xi-1)
+ beta <- xi * M2 * qu
+ # d/dMi xi = h11, h12
+ # d/dMi beta = (d/dMi M2) * xi* qu + (d/dMi xi) * M2 * qu +(d/dMi qu)* M2 * xi
+ # d/dMi M2 = 1 for i=2, = 0 for i=4
+ # d/dMi qu = d/dxi qu * d/dMi xi
+ # d/dxi qu = - qu^2 * alpha^xi * log(alpha)
+ # => d/dMi beta = (i==2)*xi*qu +
+ # (M2 * qu - beta * qu * alpha^xi * log(alpha)) * (h11,h12)
+ # dqu = M2 * qu + beta (d/dxi qu /qu)
+ dqu <- M2 * qu - beta * qu * alpha^xi * log(alpha)
h11 <- -M4/(M2*(M4-M2))/log(alpha)
h12 <- 1/(M4-M2)/log(alpha)
- t1 <- 2*M2*(M4-M2)/(M4-2*M2)^2
- t2 <- -M2^2/(M4-2*M2)^2
- h21 <- h11*M2^2/(M4-2*M2) + t1*log((M4-M2)/M2)/log(alpha)
- h22 <- h12*M2^2/(M4-2*M2) + t2*log((M4-M2)/M2)/log(alpha)
+ h21 <- h11*dqu + xi*qu
+ h22 <- h12*dqu
+ ### corresponding terms for original definition for beta, i.e.
+ ## beta <- xi * M2^2/(M4-2*M2)
+ #t1 <- 2*M2*(M4-M2)/(M4-2*M2)^2
+ #t2 <- -M2^2/(M4-2*M2)^2
+ #h21 <- h11*M2^2/(M4-2*M2) + t1*log((M4-M2)/M2)/log(alpha)
+ #h22 <- h12*M2^2/(M4-2*M2) + t2*
+
+
C <- matrix(c(h21,h22,h11,h12),2,2)
# f1 <- (1-al1)^(1+xi)/beta
@@ -35,11 +54,17 @@
# GES <- max(colSums(Werte^2)^.5)
# GES
+ if(isGP){
s11 <- al1*(1-al1)^(-1-2*xi)
s12 <- al1*(1-al1)^(-1-xi)*(1-al2)^(-xi)
s21 <- s12
s22 <- al2*(1-al2)^(-1-2*xi)
-
+ }else{
+ s11 <- al1^(-1)*(1-al1)*(-log(al1))^(-2-2*xi)
+ s12 <- al2^(-1)*(1-al2)*(log(al1)*log(al2))^(-1-1*xi)
+ s21 <- s12
+ s22 <- al2^(-1)*(1-al2)*(-log(al2))^(-2-2*xi)
+ }
S <- beta^2*matrix(c(s11,s12,s21,s22),2,2)
ASV_Pick <- t(C) %*% S %*% (C)
@@ -48,14 +73,53 @@
return(ASV_Pick)
}
+asvarQBCC <- function(model, p1 = 1/3, p2= 1/3){
+ if(!(is(model,"WeibullFamily")))
+ stop("QuantileBCC estimator only available for Weibull.")
+ scshn <- scaleshapename(model)
+ if(p1>=p2) {p<-p1; p1 <- p2; p2 <- p}
+ qm <- q(model)
+ Q1 <- qm(p1)
+ Q2 <- qm(p2)
+ l1 <- -log(p1); l2 <- -log(p2)
+ lq <- 1/(log(l2)-log(l1))
+ xi <- (log(Q2)-log(Q1))*lq
+ beta <- Q2*l2^(-1/xi)
+ dqu <- beta * log(l2) /xi^2
+ h11 <- -lq/Q1
+ h12 <- lq/Q2
+ h21 <- h11*dqu
+ h22 <- h12*dqu + l2^(-1/xi)
+ C <- matrix(c(h21,h22,h11,h12),2,2)
+ dm <- d(model)
+ s11 <- p1*(1-p1)/dm(Q1)^2
+ s12 <- p1*(1-p2)/dm(Q1)/dm(Q2)
+ s21 <- s12
+ s22 <- p2*(1-p2)/dm(Q2)^2
+ S <- matrix(c(s11,s12,s21,s22),2,2)
+ ASV <- t(C) %*% S %*% (C)
+ ASV <- PosSemDefSymmMatrix(ASV)
+ dimnames(ASV) <- list(scshn,scshn)
+ return(ASV)
+}
+
+
+
+
+
+
+
+
+
+
Modified: branches/robast-0.9/pkg/RobExtremes/man/0RobExtremes-package.Rd
===================================================================
--- branches/robast-0.9/pkg/RobExtremes/man/0RobExtremes-package.Rd 2013-01-17 21:23:32 UTC (rev 540)
+++ branches/robast-0.9/pkg/RobExtremes/man/0RobExtremes-package.Rd 2013-01-18 19:07:09 UTC (rev 541)
@@ -83,6 +83,8 @@
|>|>|>"L2GroupParamFamily" (from distrMod)
|>|>|>|>"L2ScaleShapeUnion" (from distrMod)
|>|>|>|>|>"GParetoFamily" [*]
+|>|>|>|>|>"GEVFamily" [*]
+|>|>|>|>|>"WeibullFamily" [*]
|>|>|>|>"L2LocationScaleUnion" /VIRTUAL/ (from distrMod)
|>|>|>|>|>"L2LocationFamily" (from distrMod)
|>|>|>|>|>|>"GumbelLocationFamily" [*]
@@ -196,9 +198,9 @@
Pareto Models. ArXiv 1005.1476. To appear at \emph{Statistics}.
DOI: 10.1080/02331888.2011.628022. \cr
-P. Ruckdeschel, N. Horbenko (2011): Yet another breakdown point notion:
-EFSBP --illustrated at scale-shape models. ArXiv 1005.1480. To appear at Metrika.
-DOI: 10.1007/s00184-011-0366-4.
+Ruckdeschel, P. and Horbenko, N. (2012): Yet another breakdown point notion:
+EFSBP --illustrated at scale-shape models. \emph{Metrika}, \bold{75}(8),
+1025--1047.
%a more detailed manual for \pkg{distr}, \pkg{distrSim}, \pkg{distrTEst}, and \pkg{RobExtremes} may be downloaded from
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/robast -r 541
More information about the Robast-commits
mailing list