[Robast-commits] r1005 - in branches/robast-1.1/pkg/RobExtremes/inst/AddMaterial: . LDE4Pareto
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Fri Jul 20 14:56:24 CEST 2018
Author: ruckdeschel
Date: 2018-07-20 14:56:24 +0200 (Fri, 20 Jul 2018)
New Revision: 1005
Added:
branches/robast-1.1/pkg/RobExtremes/inst/AddMaterial/LDE4Pareto/
branches/robast-1.1/pkg/RobExtremes/inst/AddMaterial/LDE4Pareto/DEstimator.R
Log:
[RobExtremes] branch 1.1 additional material: Code for LDE-type estimators in Pareto case
Added: branches/robast-1.1/pkg/RobExtremes/inst/AddMaterial/LDE4Pareto/DEstimator.R
===================================================================
--- branches/robast-1.1/pkg/RobExtremes/inst/AddMaterial/LDE4Pareto/DEstimator.R (rev 0)
+++ branches/robast-1.1/pkg/RobExtremes/inst/AddMaterial/LDE4Pareto/DEstimator.R 2018-07-20 12:56:24 UTC (rev 1005)
@@ -0,0 +1,357 @@
+## optional file if one wants analoga to LDEstimators in ParetoCase.
+
+.DMatch <- function(x.0,disp.est.0, disp.fctal.0, ParamFamily.0,
+ disp.est.ctrl.0 = NULL, disp.fctal.ctrl.0=NULL,
+ q.lo.0 =0, q.up.0=Inf, log.q.0 =TRUE, ..., vdbg=FALSE
+ ){
+ dots <- list(...)
+
+ disp.emp <- do.call(disp.est.0, args = .prepend(x.0,disp.est.ctrl.0, dots))
+ q.f <- function(xi){
+ distr.new <- ParamFamily.0 at modifyParam(xi)
+ sc.th <- do.call(disp.fctal.0, args = .prepend(distr.new,disp.fctal.ctrl.0, dots))
+ val <- if(log.q.0) log(sc.th) - log(disp.emp) else
+ sc.th-disp.emp
+ if(vdbg) print(val)
+ return(val)
+ }
+ xi.01 <- try(uniroot(q.f,lower=q.lo.0,upper=q.up.0), silent=TRUE)
+ if(is(xi.01, "try-error")) stop("Error in calculating LD-estimator: 'uniroot' did not converge.")
+ xi.0 <- xi.01$root
+ names(xi.0) <- c("shape")
+ distr.new.0 <- ParamFamily.0 at modifyParam(xi)
+ val <- c(xi.0, disp.emp)
+ names(val) <- c("shape", "disp")
+ return(val)
+}
+
+DEstimator <- function(x, disp.est, disp.fctal, ParamFamily,
+ disp.est.ctrl = NULL, disp.fctal.ctrl=NULL,
+ q.lo =1e-3, q.up=15, log.q =TRUE,
+ name, Infos, asvar = NULL, nuis.idx = NULL,
+ trafo = NULL, fixed = NULL, asvar.fct = NULL, na.rm = TRUE,
+ ..., .withEvalAsVar = FALSE, vdbg = FALSE){
+ param0 <- main(param(ParamFamily))
+ name.est <- "DEstimator"
+ es.call <- match.call()
+ if(missing(name))
+ name <- "Some estimator"
+ Dname <- paste("Dispersion:",
+ paste(deparse(substitute(disp.fctal))))
+
+ DMval <- NULL
+ estimator <- function(x,...){
+ DMval <<- .LDMatch(x.0= x,
+ disp.est.0 = disp.est,
+ disp.fctal.0 = disp.fctal,
+ ParamFamily.0 = ParamFamily,
+ disp.est.ctrl.0 = disp.est.ctrl,
+ disp.fctal.ctrl.0 = disp.fctal.ctrl,
+ q.lo.0 = q.lo,
+ q.up.0 = q.up,
+ log.q.0 = log.q, vdbg = vdbg)
+ return(LDMval[1])
+ }
+
+ asvar.fct0 <- asvar.fct
+ asvar.0 <- asvar
+ nuis.idx.0 <- nuis.idx
+ trafo.0 <- trafo
+ if(is.null(fixed)) fixed <- fixed(ParamFamily)
+ fixed.0 <- fixed
+ na.rm.0 <- na.rm
+
+ estimate <- Estimator(x, estimator, name, Infos,
+ asvar = asvar.0, nuis.idx = nuis.idx.0,
+ trafo = trafo.0, fixed = fixed.0,
+ asvar.fct = asvar.fct0,
+ na.rm = na.rm.0, ...,
+ .withEvalAsVar = .withEvalAsVar,
+ ParamFamily = ParamFamily)
+
+ estimate at estimate.call <- es.call
+
+ if(missing(Infos))
+ Infos <- matrix(c("DEstimator", Dname),
+ ncol=2, dimnames=list(character(0), c("method", "message")))
+ else{
+ Infos <- matrix(c(rep("DEstimator", length(Infos)+1), c(Dname,Infos)),
+ ncol = 2)
+ colnames(Infos) <- c("method", "message")
+ }
+ estimate at Infos <- Infos
+
+ estim <- new("DEstimate")
+
+ sln <- names(getSlots(class(estimate)))
+ for( i in 1:length(sln))
+ slot(estim, sln[i]) <- slot(estimate, sln[i])
+ rm(estimate)
+
+ estim at dispersion <- DMval["disp"]
+
+ return(.checkEstClassForParamFamily(ParamFamily,estim))
+}
+
+kMADShapeEstimator <- function(x, ParamFamily, k=1, q.lo =1e-3, q.up=15, nuis.idx = NULL,
+ trafo = NULL, fixed = NULL, asvar.fct = NULL, na.rm = TRUE,
+ ..., .withEvalAsVar = FALSE, vdbg = FALSE){
+ es.call <- match.call()
+ if(missing(k)) k <- 1
+
+ if (is.null(asvar.fct)){asvar.fct <- asvarkMAD
+ asvar <- asvarkMAD(ParamFamily, k=k)}
+
+
+ es <- LDEstimator(x, disp.est = kMAD, disp.fctal = kMAD,
+ ParamFamily = ParamFamily,
+ disp.est.ctrl = list(k=k, na.rm = na.rm),
+ disp.fctal.ctrl=list(k=k),
+ q.lo =q.lo, q.up=q.up, log.q=TRUE,
+ name = "kMADShEst", Infos="kMADShEst",
+ asvar = asvar, nuis.idx = nuis.idx, trafo = trafo, fixed = fixed,
+ asvar.fct = asvar.fct, na.rm = na.rm, ...,
+ .withEvalAsVar = .withEvalAsVar, vdbg = vdbg)
+ es at estimate.call <- es.call
+
+ return(.checkEstClassForParamFamily(ParamFamily,es))
+ }
+
+QnShapeEstimator <- function(x, ParamFamily, q.lo =1e-3, q.up=15, nuis.idx = NULL,
+ trafo = NULL, fixed = NULL, asvar.fct = NULL, na.rm = TRUE,
+ ..., .withEvalAsVar = FALSE){
+ es.call <- match.call()
+ es <- LDEstimator(x, disp.est = Qn, disp.fctal = Qn,
+ ParamFamily = ParamFamily,
+ disp.est.ctrl = list(constant=1,na.rm = na.rm),
+ disp.fctal.ctrl = NULL,
+ q.lo =q.lo, q.up=q.up, log.q=TRUE,
+ name = "QnShEst", Infos="QnShEst",
+ asvar = NULL, nuis.idx = nuis.idx, trafo = trafo, fixed = fixed,
+ asvar.fct = asvar.fct, na.rm = na.rm, ...,
+ .withEvalAsVar = .withEvalAsVar)
+ es at estimate.call <- es.call
+ return(.checkEstClassForParamFamily(ParamFamily,es))
+ }
+
+SnShapeEstimator <- function(x, ParamFamily, q.lo =1e-3, q.up=10, nuis.idx = NULL,
+ trafo = NULL, fixed = NULL, asvar.fct = NULL, na.rm = TRUE,
+ accuracy = 100, ..., .withEvalAsVar = FALSE){
+ es.call <- match.call()
+ es <- LDEstimator(x, disp.est = Sn, disp.fctal = Sn,
+ ParamFamily = ParamFamily,
+ disp.est.ctrl = list(constant=1,na.rm = na.rm),
+ disp.fctal.ctrl = list(accuracy=accuracy),
+ q.lo =q.lo, q.up=q.up, log.q=TRUE,
+ name = "SnShEst", Infos="SnShEst",
+ asvar = NULL, nuis.idx = nuis.idx, trafo = trafo, fixed = fixed,
+ asvar.fct = asvar.fct, na.rm = na.rm, ...,
+ .withEvalAsVar = .withEvalAsVar)
+ es at estimate.call <- es.call
+ return(.checkEstClassForParamFamily(ParamFamily,es))
+ }
+
+kMADhybrShapeEstimator <- function(x, ParamFamily, k=1, q.lo =1e-3, q.up=15,
+ KK=20, nuis.idx = NULL,
+ trafo = NULL, fixed = NULL, asvar.fct = NULL, na.rm = TRUE,
+ ..., .withEvalAsVar = FALSE){
+ i <- 1
+ es <- try(kMADShapeEstimator(x, ParamFamily = ParamFamily, k = k,
+ q.lo = q.lo, q.up = q.up,
+ nuis.idx = nuis.idx, trafo = trafo,
+ fixed = fixed, asvar.fct = asvar.fct, na.rm = na.rm,
+ ..., .withEvalAsVar = FALSE),
+ silent=TRUE)
+ if(! any(is.na(estimate(es))) && !is(es,"try-error"))
+ {return(.checkEstClassForParamFamily(ParamFamily,es))}
+
+ k1 <- 3.23
+ while(i<KK){
+ i <- i + 1
+ es <- try(kMADShapeEstimator(x, k = k1, ParamFamily = ParamFamily,
+ q.lo = q.lo, q.up = q.up,
+ nuis.idx = nuis.idx, trafo = trafo,
+ fixed = fixed, asvar.fct = asvar.fct, na.rm = na.rm,
+ ..., .withEvalAsVar = FALSE), silent=TRUE)
+ k1 <- k1 * 3
+ if(! any(is.na(es)) && !is(es,"try-error"))
+ {if(!missing(asvar.fct)) if(!is.null(asvar.fct)) if(.withEvalAsVar){
+ if(is.call(es at asvar)) es at asvar <- eval(es at asvar)
+ if(is.call(es at untransformed.asvar))
+ es at untransformed.asvar <- eval(es at untransformed.asvar)
+ }
+ return(.checkEstClassForParamFamily(ParamFamily,es))}
+ }
+ return(c("scale"=NA,"shape"=NA))
+}
+
+setMethod("dispersion", "DEstimate", function(object) object at dispersion)
+setMethod("location", "LEstimate", function(object) object at location)
+
+medShapeEstimator <- function(x, ParamFamily, q.lo =1e-3, q.up=10, nuis.idx = NULL,
+ trafo = NULL, fixed = NULL, asvar.fct = NULL, na.rm = TRUE,
+ accuracy = 100, ..., .withEvalAsVar = FALSE){
+ es.call <- match.call()
+ if (is.null(asvar.fct)){asvar.fct <- asvarMedShape
+ asvar <- asvarMedShape(ParamFamily)}
+ Lname <- paste("Location: Median")
+ LMval <- NULL
+ estimator <- function(x,...){
+ LMval <<- .DMatch(x.0= x,
+ disp.est.0 = median,
+ disp.fctal.0 = median,
+ ParamFamily.0 = ParamFamily,
+ disp.est.ctrl.0 = list(na.rm = na.rm),,
+ disp.fctal.ctrl.0 = NULL,
+ q.lo.0 = q.lo,
+ q.up.0 = q.up,
+ log.q.0 = log.q, vdbg = vdbg)
+ return(LMval[1])
+ }
+
+ asvar.fct0 <- asvar.fct
+ asvar.0 <- asvar
+ nuis.idx.0 <- nuis.idx
+ trafo.0 <- trafo
+ if(is.null(fixed)) fixed <- fixed(ParamFamily)
+ fixed.0 <- fixed
+ na.rm.0 <- na.rm
+
+ estimate <- Estimator(x, estimator, name, Infos,
+ asvar = asvar.0, nuis.idx = nuis.idx.0,
+ trafo = trafo.0, fixed = fixed.0,
+ asvar.fct = asvar.fct0,
+ na.rm = na.rm.0, ...,
+ .withEvalAsVar = .withEvalAsVar,
+ ParamFamily = ParamFamily)
+
+ estimate at estimate.call <- es.call
+
+ if(missing(Infos))
+ Infos <- matrix(c("DEstimator", Lname),
+ ncol=2, dimnames=list(character(0), c("method", "message")))
+ else{
+ Infos <- matrix(c(rep("DEstimator", length(Infos)+1), c(Lname,Infos)),
+ ncol = 2)
+ colnames(Infos) <- c("method", "message")
+ }
+ estimate at Infos <- Infos
+
+ estim <- new("LEstimate")
+
+ sln <- names(getSlots(class(estimate)))
+ for( i in 1:length(sln))
+ slot(estim, sln[i]) <- slot(estimate, sln[i])
+ rm(estimate)
+
+ estim at location <- LMval[2]
+ estim at estimate.call <- es.call
+
+ return(.checkEstClassForParamFamily(ParamFamily,estim))
+ }
+
+setClass("DEstimate",
+ representation(dispersion = "numeric"
+ ),
+ prototype(name = "Disp estimate",
+ estimate = numeric(0),
+ samplesize = numeric(0),
+ completecases = logical(0),
+ asvar = NULL,
+ estimate.call = call("{}"),
+ location = 0,
+ dispersion = 1,
+ Infos = matrix(c(character(0),character(0)), ncol=2,
+ dimnames=list(character(0), c("method", "message"))),
+ nuis.idx = NULL,
+ trafo = list(fct = function(x){
+ list(fval = x, mat = matrix(1))},
+ mat = matrix(1))
+ ),
+ contains = "Estimate")
+
+setClass("LEstimate",
+ representation(location = "numeric",
+ ),
+ prototype(name = "Loc estimate",
+ estimate = numeric(0),
+ samplesize = numeric(0),
+ completecases = logical(0),
+ asvar = NULL,
+ estimate.call = call("{}"),
+ location = 0,
+ dispersion = 1,
+ Infos = matrix(c(character(0),character(0)), ncol=2,
+ dimnames=list(character(0), c("method", "message"))),
+ nuis.idx = NULL,
+ trafo = list(fct = function(x){
+ list(fval = x, mat = matrix(1))},
+ mat = matrix(1))
+ ),
+ contains = "Estimate")
+
+setClass("ParetoEstimate", contains="Estimate")
+setClass("ParetoDEstimate", contains=c("DEstimate", "ParetoEstimate"))
+setClass("ParetoLEstimate", contains=c("LEstimate", "ParetoEstimate"))
+setClass("ParetokStepEstimate", contains=c("kStepEstimate", "ParetoEstimate"))
+setClass("ParetoMCEstimate", contains=c("MCEstimate", "ParetoEstimate"))
+
+
+asvarkMADShape <- function(model, k=1){
+
+ xi <- main(model at param)
+ beta <- fixed(model at param)
+
+ M <- kMAD(model at distribution, k=k)
+ m <- q.l(model)(.5)
+
+ x1.0 <- m - M
+ x2.0 <- m + k * M
+
+ ## joint Variance of median and kMAD, see Serfling Mazumder
+ dmm <- d(model)(x1.0)
+ dmp <- d(model)(x2.0)
+ dm <- d(model)(m)
+ alpha <- p(model)(m-M)+p(model)(m+k*M)
+ betA <- dmm-dmp
+ ceta <- dmm+k*dmp
+ eta <- betA^2 + 4*(1-alpha)*betA*dm
+
+ V <- 1/4/ceta^2*(1+eta/dm^2)
+
+ psi_kMad <- function(x){
+ cp <- k*dmp+dmm
+ cm = dmp-dmm
+ return((0.5-((x<=m+k*M)&(x>=m-M)))/cp + cm/cp*((x<=m)-0.5)/dm)
+ }
+
+ L2d <- model at L2deriv[[1]]
+ L_xi.f = function(x) evalRandVar(L2d,x)[2,]
+
+ E12 <- E(distribution(model),fun=function(x) psi_kMad(x) * L_xi.f(x))
+
+ ASV_Med <- PosSemDefSymmMatrix(V/E12^2)
+ dimnames(ASV_Med) <- list(shapename(model),shapename(model))
+ return(ASV_Med)
+}
+
+asvarMedShape <- function(model){
+
+ m <- q.l(model)(.5)
+
+ dm <- d(model)(m)
+
+ V <- 1/4/dm^2
+
+ psi_med <- function(x) (0.5-(x<=m))/dm
+
+ L2d <- model at L2deriv[[1]]
+ L_xi.f = function(x) evalRandVar(L2d,x)[2,]
+
+ E12 <- E(distribution(model),fun=function(x) psi_med(x) * L_xi.f(x))
+ D <- 1/E12
+
+ ASV_Med <- PosSemDefSymmMatrix(V/E12^2)
+ dimnames(ASV_Med) <- list(shapename(model),shapename(model))
+ return(ASV_Med)
+}
More information about the Robast-commits
mailing list