[Robast-commits] r468 - in branches/robast-0.9/pkg: . RobExtremes RobExtremes/R RobExtremes/inst RobExtremes/man RobExtremes/src
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Sun May 20 16:05:21 CEST 2012
Author: ruckdeschel
Date: 2012-05-20 16:05:21 +0200 (Sun, 20 May 2012)
New Revision: 468
Added:
branches/robast-0.9/pkg/RobExtremes/
branches/robast-0.9/pkg/RobExtremes/DESCRIPTION
branches/robast-0.9/pkg/RobExtremes/NAMESPACE
branches/robast-0.9/pkg/RobExtremes/R/
branches/robast-0.9/pkg/RobExtremes/R/AllClass.R
branches/robast-0.9/pkg/RobExtremes/R/AllGeneric.R
branches/robast-0.9/pkg/RobExtremes/R/AllInitialize.R
branches/robast-0.9/pkg/RobExtremes/R/Expectation.R
branches/robast-0.9/pkg/RobExtremes/R/Functionals.R
branches/robast-0.9/pkg/RobExtremes/R/GEV.R
branches/robast-0.9/pkg/RobExtremes/R/GPareto.R
branches/robast-0.9/pkg/RobExtremes/R/GParetoFamily.R
branches/robast-0.9/pkg/RobExtremes/R/Gumbel.R
branches/robast-0.9/pkg/RobExtremes/R/GumbelLocationFamily.R
branches/robast-0.9/pkg/RobExtremes/R/Kurtosis.R
branches/robast-0.9/pkg/RobExtremes/R/LDEstimator.R
branches/robast-0.9/pkg/RobExtremes/R/Pareto.R
branches/robast-0.9/pkg/RobExtremes/R/Skewness.R
branches/robast-0.9/pkg/RobExtremes/R/SnQn.R
branches/robast-0.9/pkg/RobExtremes/R/internals.R
branches/robast-0.9/pkg/RobExtremes/R/interpolLM.R
branches/robast-0.9/pkg/RobExtremes/R/interpolSn.R
branches/robast-0.9/pkg/RobExtremes/R/kMAD.R
branches/robast-0.9/pkg/RobExtremes/R/sysdata.rda
branches/robast-0.9/pkg/RobExtremes/inst/
branches/robast-0.9/pkg/RobExtremes/inst/CITATION
branches/robast-0.9/pkg/RobExtremes/inst/NEWS
branches/robast-0.9/pkg/RobExtremes/inst/TOBEDONE
branches/robast-0.9/pkg/RobExtremes/inst/scripts/
branches/robast-0.9/pkg/RobExtremes/man/
branches/robast-0.9/pkg/RobExtremes/man/0RobExtremes-package.Rd
branches/robast-0.9/pkg/RobExtremes/man/E.Rd
branches/robast-0.9/pkg/RobExtremes/man/GEV-class.Rd
branches/robast-0.9/pkg/RobExtremes/man/GEV.Rd
branches/robast-0.9/pkg/RobExtremes/man/GEVParameter-class.Rd
branches/robast-0.9/pkg/RobExtremes/man/GPareto-class.Rd
branches/robast-0.9/pkg/RobExtremes/man/GPareto.Rd
branches/robast-0.9/pkg/RobExtremes/man/GParetoFamily.Rd
branches/robast-0.9/pkg/RobExtremes/man/GParetoParameter-class.Rd
branches/robast-0.9/pkg/RobExtremes/man/Gumbel-class.Rd
branches/robast-0.9/pkg/RobExtremes/man/Gumbel.Rd
branches/robast-0.9/pkg/RobExtremes/man/GumbelLocationFamily.Rd
branches/robast-0.9/pkg/RobExtremes/man/GumbelParameter-class.Rd
branches/robast-0.9/pkg/RobExtremes/man/InternalReturnClasses-class.Rd
branches/robast-0.9/pkg/RobExtremes/man/LDEstimator.Rd
branches/robast-0.9/pkg/RobExtremes/man/Pareto-class.Rd
branches/robast-0.9/pkg/RobExtremes/man/Pareto.Rd
branches/robast-0.9/pkg/RobExtremes/man/ParetoParameter-class.Rd
branches/robast-0.9/pkg/RobExtremes/man/RobExtremesConstants.Rd
branches/robast-0.9/pkg/RobExtremes/man/Var.Rd
branches/robast-0.9/pkg/RobExtremes/man/internal-interpolateLM.Rd
branches/robast-0.9/pkg/RobExtremes/man/internal-interpolateSn.Rd
branches/robast-0.9/pkg/RobExtremes/man/internalldeHelpers.Rd
branches/robast-0.9/pkg/RobExtremes/man/interpolateSn.Rd
branches/robast-0.9/pkg/RobExtremes/man/kMAD.Rd
branches/robast-0.9/pkg/RobExtremes/man/validParameter-methods.Rd
branches/robast-0.9/pkg/RobExtremes/src/
branches/robast-0.9/pkg/RobExtremes/src/kMad.c
branches/robast-0.9/pkg/RobExtremes/tests/
Log:
RobExtremes: a first version of package committed
Added: branches/robast-0.9/pkg/RobExtremes/DESCRIPTION
===================================================================
--- branches/robast-0.9/pkg/RobExtremes/DESCRIPTION (rev 0)
+++ branches/robast-0.9/pkg/RobExtremes/DESCRIPTION 2012-05-20 14:05:21 UTC (rev 468)
@@ -0,0 +1,19 @@
+Package: RobExtremes
+Version: 0.9
+Date: 2012-05-17
+Title: Optimally robust estimation for extreme value distributions
+Description: Optimally robust estimation for extreme value distributions
+ using S4 classes and methods (based on packages distr, distrEx,
+ distrMod, RobAStBase, and ROptEst)
+Depends: R (>= 2.14.0), methods, distr(>= 2.3), distrEx(>= 2.3), RandVar(>= 0.8),
+ distrMod(>= 2.4), RobAStBase(>= 0.8), ROptEst(>= 0.8), robustbase(>= 0.8-0),
+ evd, actuar
+Author: Peter Ruckdeschel, Matthias Kohl, Nataliya Horbenko
+Maintainer: Peter Ruckdeschel <peter.ruckdeschel at itwm.fraunhofer.de>
+LazyLoad: yes
+ByteCompile: yes
+License: LGPL-3
+URL: http://robast.r-forge.r-project.org/
+LastChangedDate: {$LastChangedDate: 2011-09-30 11:10:33 +0200 (Fr, 30 Sep 2011) $}
+LastChangedRevision: {$LastChangedRevision: 453 $}
+SVNRevision: 439
Added: branches/robast-0.9/pkg/RobExtremes/NAMESPACE
===================================================================
--- branches/robast-0.9/pkg/RobExtremes/NAMESPACE (rev 0)
+++ branches/robast-0.9/pkg/RobExtremes/NAMESPACE 2012-05-20 14:05:21 UTC (rev 468)
@@ -0,0 +1,36 @@
+useDynLib("RobExtremes")
+
+import("methods")
+import("distr")
+import("distrEx")
+import("RandVar")
+import("distrMod")
+import("robustbase")
+import("RobAStBase")
+import("ROptEst")
+
+exportClasses("GumbelParameter",
+ "ParetoParameter",
+ "GParetoParameter",
+ "GEVParameter")
+exportClasses("Gumbel", "Pareto", "GPareto", "GEV")
+exportClasses("GParetoFamily", "GumbelLocationFamily")
+exportMethods("initialize")
+exportMethods("loc", "loc<-",
+ "E", "var", "IQR", "median", "kMAD", "Sn", "Qn")
+exportMethods("validParameter",
+ "location", "location<-",
+ "scale", "scale<-",
+ "shape", "shape<-",
+ "+", "*",
+ "Min", "Min<-",
+ "Range",
+ "E", "var", "IQR", "skewness", "kurtosis",
+ "sd", "median", "mad")
+exportMethods("modifyModel")
+export("EULERMASCHERONICONSTANT","APERYCONSTANT")
+export("Gumbel", "Pareto", "GPareto", "GEV")
+export("GParetoFamily", "GumbelLocationFamily")
+export("LDEstimator", "medkMAD", "medSn", "medQn", "medkMADhybr")
+export("getShapeGrid", "getSnGrid")
+export("loc", "loc<-")
\ No newline at end of file
Added: branches/robast-0.9/pkg/RobExtremes/R/AllClass.R
===================================================================
--- branches/robast-0.9/pkg/RobExtremes/R/AllClass.R (rev 0)
+++ branches/robast-0.9/pkg/RobExtremes/R/AllClass.R 2012-05-20 14:05:21 UTC (rev 468)
@@ -0,0 +1,241 @@
+.onLoad <- function(lib, pkg){
+# require("methods", character = TRUE, quietly = TRUE)
+}
+
+
+.onAttach <- function(library, pkg)
+{
+#unlockBinding(".RobExtremesOptions", asNamespace("RobExtremes"))
+#msga <- gettext("Note: Packages \"e1071\", \"moments\", \"fBasics\" should be attached ")
+#msgb <- gettext("/before/ package \"distrEx\". See distrExMASK().")
+buildStartupMessage(pkg = "RobExtremes", msga="", msgb="",
+ library = library, packageHelp = TRUE
+# MANUAL="http://www.uni-bayreuth.de/departments/math/org/mathe7/DISTR/distr.pdf",
+#VIGNETTE = gettext("Package \"distrDoc\" provides a vignette to this package as well as to several related packages; try vignette(\"distr\").")
+)
+ invisible()
+}
+
+
+.onUnload <- function(libpath)
+{
+ library.dynam.unload("RobExtremes", libpath)
+}
+
+
+
+# parameter of Gumbel distribution
+setClass("GumbelParameter", representation(loc = "numeric",
+ scale = "numeric"),
+ prototype(name = gettext("parameter of a Gumbel distribution"),
+ loc = 0, scale = 1),
+ contains = "Parameter",
+ validity = function(object){
+ if(length(object at scale) != 1)
+ stop("length of 'scale' is not equal to 1")
+ if(length(object at loc) != 1)
+ stop("length of 'loc' is not equal to 1")
+ if(object at scale <= 0)
+ stop("'scale' has to be positive")
+ else return(TRUE)
+ })
+
+# Gumbel distribution
+setClass("Gumbel",
+ prototype = prototype(r = function(n){ rgumbel(n, loc = 0, scale = 1) },
+ d = function(x, log){ dgumbel(x, loc = 0, scale = 1, log = FALSE) },
+ p = function(q, lower.tail = TRUE, log.p = FALSE){
+ p0 <- pgumbel(q, loc = 0, scale = 1, lower.tail = lower.tail)
+ if(log.p) return(log(p0)) else return(p0)
+ },
+ q = function(p, loc = 0, scale = 1, lower.tail = TRUE, log.p = FALSE){
+ ## P.R.: changed to vectorized form
+ p1 <- if(log.p) exp(p) else p
+
+ in01 <- (p1>1 | p1<0)
+ i01 <- .isEqual01(p1)
+ i0 <- (i01 & p1<1)
+ i1 <- (i01 & p1>0)
+ ii01 <- .isEqual01(p1) | in01
+
+ p0 <- p
+ p0[ii01] <- if(log.p) log(0.5) else 0.5
+
+ q1 <- qgumbel(p0, loc = 0, scale = 1,
+ lower.tail = lower.tail)
+ q1[i0] <- if(lower.tail) -Inf else Inf
+ q1[i1] <- if(!lower.tail) -Inf else Inf
+ q1[in01] <- NaN
+
+ return(q1)
+ },
+ img = new("Reals"),
+ param = new("GumbelParameter"),
+ .logExact = FALSE,
+ .lowerExact = TRUE),
+ contains = "AbscontDistribution")
+
+
+###### Pareto distribution by Nataliya Horbenko, ITWM, 18-03-09
+## Class: ParetoParameter
+setClass("ParetoParameter",
+ representation = representation(shape = "numeric",
+ Min = "numeric"
+ ),
+ prototype = prototype(shape = 1, Min = 1, name =
+ gettext("Parameter of a Pareto distribution")
+ ),
+ contains = "Parameter"
+ )
+
+## Class: Pareto distribution
+setClass("Pareto",
+ prototype = prototype(
+ r = function(n){ rpareto1(n, shape = 1, min = 1) },
+ d = function(x, log = FALSE){
+ dpareto1(x, shape = 1, min = 1, log = log)
+ },
+ p = function(q, lower.tail = TRUE, log.p = FALSE ){
+ ppareto1(q, shape = 1, min = 1,
+ lower.tail = lower.tail, log.p = log.p)
+ },
+ q = function(p, lower.tail = TRUE, log.p = FALSE ){
+ ## P.R.: changed to vectorized form
+ p1 <- if(log.p) exp(p) else p
+
+ in01 <- (p1>1 | p1<0)
+ i01 <- .isEqual01(p1)
+ i0 <- (i01 & p1<1)
+ i1 <- (i01 & p1>0)
+ ii01 <- .isEqual01(p1) | in01
+
+ p0 <- p
+ p0[ii01] <- if(log.p) log(0.5) else 0.5
+
+ q1 <- qpareto1(p0, shape = 1, min = 1,
+ lower.tail = lower.tail, log.p = log.p)
+ q1[i0] <- if(lower.tail) -Inf else Inf
+ q1[i1] <- if(!lower.tail) -Inf else Inf
+ q1[in01] <- NaN
+
+ return(q1)
+ },
+ param = new("ParetoParameter"),
+ img = new("Reals"),
+ .logExact = TRUE,
+ .lowerExact = TRUE),
+ contains = "AbscontDistribution"
+ )
+
+## Class: GParetoParameter
+setClass("GParetoParameter",
+ representation = representation(loc = "numeric", scale = "numeric", shape = "numeric"
+ ),
+ prototype = prototype(loc = 0, scale = 1, shape = 0, name =
+ gettext("Parameter of a generalized Pareto distribution")
+ ),
+ contains = "Parameter"
+ )
+## Class: Generalized Pareto distribution
+setClass("GPareto",
+ prototype = prototype(
+ r = function(n){ rgpd(n,loc = 0, scale = 1, shape = 1) },
+ d = function(x, log = FALSE){
+ dgpd(x, loc = 0, scale = 1, shape = 1, log = log)
+ },
+ p = function(q, lower.tail = TRUE, log.p = FALSE ){
+ p0 <- pgpd(q, loc = 0, scale = 1, shape = 1)
+ if(!lower.tail ) p0 <- 1-p0
+ if(log.p) p0 <- log(p0)
+ return(p0)},
+ q = function(p, lower.tail = TRUE, log.p = FALSE ){
+ ## P.R.: changed to vectorized form
+ p1 <- if(log.p) exp(p) else p
+ if(!lower.tail) p1 <- 1-p1
+
+ in01 <- (p1>1 | p1<0)
+ i01 <- .isEqual01(p1)
+ i0 <- (i01 & p1<1)
+ i1 <- (i01 & p1>0)
+ ii01 <- .isEqual01(p1) | in01
+
+ p0 <- p
+ p0[ii01] <- if(log.p) log(0.5) else 0.5
+
+ q1 <- qgpd(p0,loc=0, scale = 1, shape = 1)
+ q1[i0] <- if(lower.tail) -Inf else Inf
+ q1[i1] <- if(!lower.tail) -Inf else Inf
+ q1[in01] <- NaN
+
+ return(q1)
+ },
+ param = new("GParetoParameter"),
+ img = new("Reals"),
+ .withArith = FALSE,
+ .withSim = FALSE,
+ .logExact = TRUE,
+ .lowerExact = TRUE),
+ contains = "AbscontDistribution"
+ )
+
+
+## Class: GEVParameter
+setClass("GEVParameter",
+ representation = representation(loc = "numeric", scale = "numeric", shape = "numeric"
+ ),
+ prototype = prototype(loc = 0, scale = 1, shape = 0.5, name =
+ gettext("Parameter of a generalized extreme value distribution")
+ ),
+ contains = "Parameter"
+ )
+## Class: Generalized extreme value distribution
+setClass("GEV",
+ prototype = prototype(
+ r = function(n){ rgev(n,loc = 0, scale = 1, shape = 0.5) },
+ d = function(x, log = FALSE){
+ dgev(x, loc = 0, scale = 1, shape = 0.5, log = log)
+ },
+ p = function(q, lower.tail = TRUE, log.p = FALSE ){
+ p0 <- pgev(q, loc = 0, scale = 1, shape = 0.5)
+ if(!lower.tail ) p0 <- 1-p0
+ if(log.p) p0 <- log(p0)
+ return(p0)},
+ q = function(p, lower.tail = TRUE, log.p = FALSE ){
+ ## analogous to GPD
+ p1 <- if(log.p) exp(p) else p
+ if(!lower.tail) p1 <- 1-p1
+
+ in01 <- (p1>1 | p1<0)
+ i01 <- .isEqual01(p1)
+ i0 <- (i01 & p1<1)
+ i1 <- (i01 & p1>0)
+ ii01 <- .isEqual01(p1) | in01
+
+ p0 <- p
+ p0[ii01] <- if(log.p) log(0.5) else 0.5
+
+ q1 <- qgev(p0,loc=0, scale = 1, shape = 0.5)
+ q1[i0] <- if(lower.tail) -Inf else Inf
+ q1[i1] <- if(!lower.tail) -Inf else Inf
+ q1[in01] <- NaN
+
+ return(q1)
+ },
+ param = new("GEVParameter"),
+ img = new("Reals"),
+ .withArith = FALSE,
+ .withSim = FALSE,
+ .logExact = TRUE,
+ .lowerExact = TRUE),
+ contains = "AbscontDistribution"
+ )
+## Gumbel location family
+setClass("GumbelLocationFamily",
+ contains = "L2LocationFamily")
+
+## class
+setClass("GParetoFamily",
+ prototype= prototype(withPos = TRUE),
+ contains="L2ScaleShapeUnion")
+
+
Added: branches/robast-0.9/pkg/RobExtremes/R/AllGeneric.R
===================================================================
--- branches/robast-0.9/pkg/RobExtremes/R/AllGeneric.R (rev 0)
+++ branches/robast-0.9/pkg/RobExtremes/R/AllGeneric.R 2012-05-20 14:05:21 UTC (rev 468)
@@ -0,0 +1,22 @@
+
+if(!isGeneric("kMAD")){
+ setGeneric("kMAD", function(x, k, ...) standardGeneric("kMAD"))
+}
+
+
+if(!isGeneric("loc")){
+ setGeneric("loc", function(object) standardGeneric("loc"))
+}
+
+if(!isGeneric("loc<-")){
+ setGeneric("loc<-", function(object, value) standardGeneric("loc<-"))
+}
+
+
+if(!isGeneric("Qn")){
+ setGeneric("Qn", function(x, ...) standardGeneric("Qn"))
+}
+
+if(!isGeneric("Sn")){
+ setGeneric("Sn", function(x, ...) standardGeneric("Sn"))
+}
Added: branches/robast-0.9/pkg/RobExtremes/R/AllInitialize.R
===================================================================
--- branches/robast-0.9/pkg/RobExtremes/R/AllInitialize.R (rev 0)
+++ branches/robast-0.9/pkg/RobExtremes/R/AllInitialize.R 2012-05-20 14:05:21 UTC (rev 468)
@@ -0,0 +1,236 @@
+## initialize method
+setMethod("initialize", "Gumbel",
+ function(.Object, loc = 0, scale = 1) {
+ .Object at img <- Reals()
+ .Object at param <- new("GumbelParameter", loc = loc, scale = scale,
+ name = gettext("parameter of a Gumbel distribution"))
+ .Object at r <- function(n){}
+ body(.Object at r) <- substitute({ rgumbel(n, loc = loc1, scale = scale1) },
+ list(loc1 = loc, scale1 = scale))
+ .Object at d <- function(x, log = FALSE){}
+ body(.Object at d) <- substitute({ dgumbel(x, loc = loc1, scale = scale1, log = log) },
+ list(loc1 = loc, scale1 = scale))
+ .Object at p <- function(q, lower.tail = TRUE, log.p = FALSE){}
+ body(.Object at p) <- substitute({p1 <- pgumbel(q, loc = loc1, scale = scale1, lower.tail = lower.tail)
+ return(if(log.p) log(p1) else p1)},
+ list(loc1 = loc, scale1 = scale))
+ .Object at q <- function(p, loc = loc1, scale = scale1, lower.tail = TRUE, log.p = FALSE){}
+ body(.Object at q) <- substitute({
+ ## P.R.: changed to vectorized form
+ p1 <- if(log.p) exp(p) else p
+
+ in01 <- (p1>1 | p1<0)
+ i01 <- .isEqual01(p1)
+ i0 <- (i01 & p1<1)
+ i1 <- (i01 & p1>0)
+ ii01 <- .isEqual01(p1) | in01
+
+ p0 <- p
+ p0[ii01] <- if(log.p) log(0.5) else 0.5
+
+ q1 <- qgumbel(p0, loc = loc1, scale = scale1,
+ lower.tail = lower.tail)
+ q1[i0] <- if(lower.tail) -Inf else Inf
+ q1[i1] <- if(!lower.tail) -Inf else Inf
+ q1[in01] <- NaN
+
+ return(q1)
+ }, list(loc1 = loc, scale1 = scale))
+ .Object at .withSim <- FALSE
+ .Object at .withArith <- FALSE
+ .Object at .logExact <- FALSE
+ .Object at .lowerExact <- TRUE
+ .Object
+ })
+
+## Class: Pareto distribution
+setMethod("initialize", "Pareto",
+ function(.Object, shape = 1, Min = 1, .withArith = FALSE) {
+ .Object at img <- new("Reals")
+ .Object at param <- new("ParetoParameter", shape = shape, Min = Min)
+ .Object at r <- function(n){}
+ .Object at d <- function(x, log = FALSE){}
+ .Object at p <- function(q, lower.tail = TRUE, log.p = FALSE){}
+ .Object at q <- function(p, lower.tail = TRUE, log.p = FALSE){}
+ body(.Object at r) <- substitute(
+ { rpareto1(n, shape = shapeSub, min = MinSub) },
+ list(shapeSub = shape, MinSub = Min)
+ )
+ body(.Object at d) <- substitute(
+ { dpareto1(x, shape = shapeSub, min = MinSub,
+ log = log) },
+ list(shapeSub = shape, MinSub = Min)
+ )
+ body(.Object at p) <- substitute(
+ { ppareto1(q, shape = shapeSub, min = MinSub,
+ lower.tail = lower.tail, log.p = log.p) },
+ list(shapeSub = shape, MinSub = Min)
+ )
+ body(.Object at q) <- substitute({
+ ## P.R.: changed to vectorized form
+ p1 <- if(log.p) exp(p) else p
+
+ in01 <- (p1>1 | p1<0)
+ i01 <- .isEqual01(p1)
+ i0 <- (i01 & p1<1)
+ i1 <- (i01 & p1>0)
+ ii01 <- .isEqual01(p1) | in01
+
+ p0 <- p
+ p0[ii01] <- if(log.p) log(0.5) else 0.5
+
+ q1 <- qpareto1(p0, shape = shapeSub, min = MinSub,
+ lower.tail = lower.tail, log.p = log.p)
+ q1[i0] <- if(lower.tail) -Inf else Inf
+ q1[i1] <- if(!lower.tail) -Inf else Inf
+ q1[in01] <- NaN
+
+ return(q1)
+ }, list(shapeSub = shape, MinSub = Min))
+ .Object at .withArith <- .withArith
+ .Object at .logExact <- TRUE
+ .Object at .lowerExact <- TRUE
+ .Object
+ })
+
+## Class: Generalized Pareto distribution
+setMethod("initialize", "GPareto",
+ function(.Object, loc = 0, scale = 1, shape = 1) {
+ .Object at img <- new("Reals")
+ .Object at param <- new("GParetoParameter", loc = loc, scale = scale, shape = shape)
+ .Object at r <- function(n){}
+ .Object at d <- function(x, log = FALSE){}
+ .Object at p <- function(q, lower.tail = TRUE, log.p = FALSE){}
+ .Object at q <- function(p, lower.tail = TRUE, log.p = FALSE){}
+ body(.Object at r) <- substitute(
+ { rgpd(n, loc = locSub, scale = scaleSub, shape = shapeSub) },
+ list(locSub = loc, scaleSub = scale, shapeSub = shape)
+ )
+ body(.Object at d) <- substitute(
+ { dgpd(x, loc = locSub, scale = scaleSub, shape = shapeSub,
+ log = log) },
+ list(locSub = loc, scaleSub = scale, shapeSub = shape)
+ )
+ body(.Object at p) <- substitute(
+ { if(!lower.tail && log.p){
+ q0 <- (q-locSub)/scaleSub
+ return(-log(1+shapeSub*q0)/shapeSub)
+ }else{
+ p0 <- pgpd(q, loc = locSub, scale = scaleSub,
+ shape = shapeSub)
+ if(!lower.tail ) p0 <- 1-p0
+ if(log.p) p0 <- log(p0)
+ return(p0)}
+ }, list(locSub = loc, scaleSub = scale,
+ shapeSub = shape)
+ )
+ body(.Object at q) <- substitute({
+ if(!lower.tail && log.p){
+ p1 <- p
+ p1[p<.Machine$double.eps] <- 0.5
+ q0 <- (exp(-shapeSub*p1)-1)/shapeSub*scaleSub + locSub
+ q0[p<.Machine$double.eps] <- NaN
+ return(q0)
+ }else{
+
+ ## P.R.: changed to vectorized form
+ p1 <- if(log.p) exp(p) else p
+
+ in01 <- (p1>1 | p1<0)
+ i01 <- .isEqual01(p1)
+ i0 <- (i01 & p1<1)
+ i1 <- (i01 & p1>0)
+ ii01 <- .isEqual01(p1) | in01
+
+ p0 <- p
+ p0[ii01] <- if(log.p) log(0.5) else 0.5
+ if(!lower.tail) p0 <- 1-p0
+
+ q1 <- qgpd(p0, loc = locSub, scale = scaleSub,
+ shape = shapeSub)
+ q1[i0] <- if(lower.tail) locSub else Inf
+ q1[i1] <- if(!lower.tail) locSub else Inf
+ q1[in01] <- NaN
+
+ return(q1)
+ }
+ }, list(locSub = loc, scaleSub = scale, shapeSub = shape))
+
+ .Object at .withSim <- FALSE
+ .Object at .withArith <- FALSE
+ .Object at .logExact <- TRUE
+ .Object at .lowerExact <- TRUE
+ .Object
+ })
+
+
+## Class: Generalized extreme value distribution
+setMethod("initialize", "GEV",
+ function(.Object, loc = 0, scale = 1, shape = 1) {
+ .Object at img <- new("Reals")
+ .Object at param <- new("GEVParameter", loc = loc, scale = scale, shape = shape)
+ .Object at r <- function(n){}
+ .Object at d <- function(x, log = FALSE){}
+ .Object at p <- function(q, lower.tail = TRUE, log.p = FALSE){}
+ .Object at q <- function(p, lower.tail = TRUE, log.p = FALSE){}
+ body(.Object at r) <- substitute(
+ { rgev(n, loc = locSub, scale = scaleSub, shape = shapeSub) },
+ list(locSub = loc, scaleSub = scale, shapeSub = shape)
+ )
+ body(.Object at d) <- substitute(
+ { dgev(x, loc = locSub, scale = scaleSub, shape = shapeSub, log = log) },
+ list(locSub = loc, scaleSub = scale, shapeSub = shape)
+ )
+ body(.Object at p) <- substitute(
+ { if(lower.tail && log.p){
+ q0 <- (q-locSub)/scaleSub
+ p0 <- -(1+shapeSub*q0)^(-1/shapeSub)
+ p0[q0<(-1)] <- -Inf
+ return(p0)
+ }else{
+ p0 <- pgev(q, loc = locSub, scale = scaleSub, shape = shapeSub,lower.tail=TRUE)
+ if(!lower.tail ) p0 <- 1-p0
+ if(log.p) p0 <- log(p0)
+ return(p0)}
+ }, list(locSub = loc, scaleSub = scale,
+ shapeSub = shape)
+ )
+ body(.Object at q) <- substitute({
+ if(lower.tail && log.p){
+ q0 <-((-p)^(-shapeSub)-1)/shapeSub*scaleSub+locSub
+ #q0[p>0|p< -Inf] <- NaN
+ #q0[.isEqual01(p)& p<1] <- Inf
+ #q0[!is.finite(p)& p<0] <- locSub-scaleSub/shapeSub
+ p0 <- exp(p)
+ q0[p0>1|p0<0] <- NaN
+ q0[(.isEqual01(p) & p0>0)] <- Inf
+ q0[(.isEqual01(p) & p0<1)] <- locSub-scaleSub/shapeSub
+ return(q0)
+ }else{
+ ##higher tolerance for .isEqual01
+ tol=1e-20
+ distroptions(TruncQuantile=tol)
+ p1 <- if(log.p) exp(p) else p
+ in01 <- (p1>1 | p1<0)
+ i01 <- .isEqual01(p1)
+ i0 <- (i01 & p1<1)
+ i1 <- (i01 & p1>0)
+ ii01 <- .isEqual01(p1) | in01
+ p0 <- p
+ p0[ii01] <- if(log.p) log(0.5) else 0.5
+ #if(!lower.tail) p0 <- 1-p0
+ q1 <- qgev(p0, loc = locSub, scale = scaleSub, shape = shapeSub, lower.tail=lower.tail)
+ q1[i0] <- if(lower.tail) locSub-scaleSub/shapeSub else Inf
+ q1[i1] <- if(!lower.tail) locSub-scaleSub/shapeSub else Inf
+ q1[in01] <- NaN
+ return(q1)
+ }
+ }, list(locSub = loc, scaleSub = scale, shapeSub = shape))
+
+ .Object at .withSim <- FALSE
+ .Object at .withArith <- FALSE
+ .Object at .logExact <- TRUE
+ .Object at .lowerExact <- TRUE
+ .Object
+ })
+
Added: branches/robast-0.9/pkg/RobExtremes/R/Expectation.R
===================================================================
--- branches/robast-0.9/pkg/RobExtremes/R/Expectation.R (rev 0)
+++ branches/robast-0.9/pkg/RobExtremes/R/Expectation.R 2012-05-20 14:05:21 UTC (rev 468)
@@ -0,0 +1,144 @@
+
+setMethod("E", signature(object = "Pareto",
+ fun = "missing",
+ cond = "missing"),
+ function(object, low = NULL, upp = NULL, ...){
+ if(!is.null(low)) if(low <= Min(object)) low <- NULL
+ a <- shape(object); b <- Min(object)
+ if(is.null(low) && is.null(upp)){
+ if(a<=1) return(Inf)
+ else return(b*a/(a-1))
+ }
+ else
+ return(E(as(object,"AbscontDistribution"), low=low, upp=upp, ...))
+ })
+
+### source http://mathworld.wolfram.com/ParetoDistribution.html
+
+
+setMethod("E", signature(object = "Gumbel",
+ fun = "missing",
+ cond = "missing"),
+ function(object, low = NULL, upp = NULL, ...){a <- loc(object); b <- scale(object)
+ if(is.null(low) && is.null(upp))
+ return(a- EULERMASCHERONICONSTANT * b)
+ else
+ return(E(as(object,"AbscontDistribution"), low=low, upp=upp, ...))
+ })
+## http://mathworld.wolfram.com/GumbelDistribution.html
+
+setMethod("E", signature(object = "GPareto",
+ fun = "missing",
+ cond = "missing"),
+ function(object, low = NULL, upp = NULL, ...){
+ if(!is.null(low)) if(low <= Min(object)) low <- NULL
+ k <- shape(object); s <- scale(object); mu <- loc(object)
+ if(is.null(low) && is.null(upp)){
+ if(k>=1) return(Inf)
+ else return(mu+s/(1-k))
+ }
+ else
+ return(E(as(object,"AbscontDistribution"), low=low, upp=upp, ...))
+ })
+
+### source http://en.wikipedia.org/wiki/Pareto_distribution
+
+setMethod("E", signature(object = "GPareto",
+ fun = "function",
+ cond = "missing"),
+ function(object, fun, low = NULL, upp = NULL,
+ rel.tol= getdistrExOption("ErelativeTolerance"),
+ lowerTruncQuantile = getdistrExOption("ElowerTruncQuantile"),
+ upperTruncQuantile = getdistrExOption("EupperTruncQuantile"),
+ IQR.fac = max(1e4,getdistrExOption("IQR.fac")), ...
+ ){
+
+ dots <- list(...)
+ dots.withoutUseApply <- dots
+ useApply <- TRUE
+ if(!is.null(dots$useApply)) useApply <- dots$useApply
+ dots.withoutUseApply$useApply <- NULL
+ integrand <- function(x, dfun, ...){ di <- dim(x)
+ y <- exp(x)
+ if(useApply){
+ funy <- sapply(y,fun, ...)
+ dim(y) <- di
+ dim(funy) <- di
+ }else funy <- fun(y,...)
+ return(funy * y * dfun(y)) }
+
+ if(is.null(low)) low <- -Inf
+ if(is.null(upp)) upp <- Inf
+
+ Ib <- .getIntbounds(object, low, upp, lowerTruncQuantile,
+ upperTruncQuantile, IQR.fac)
+ low <- if(Ib["low"]<=0) -Inf else log(Ib["low"])
+ upp <- log(Ib["upp"])
+
+ return(do.call(distrExIntegrate, c(list(f = integrand,
+ lower = low,
+ upper = upp,
+ rel.tol = rel.tol,
+ distr = object, dfun = d(object)), dots.withoutUseApply)))
+
+ })
+
+
+setMethod("E", signature(object = "GEV",
+ fun = "missing",
+ cond = "missing"),
+ function(object, low = NULL, upp = NULL, ...){
+ if(!is.null(low)) if(low <= Min(object)) low <- NULL
+ xi <- shape(object); sigma <- scale(object); mu <- loc(object)
+ if(is.null(low) && is.null(upp)){
+ if (xi==0) return(mu+sigma*EULERMASCHERONICONSTANT)
+ else if(xi>=1) return(Inf)
+ else return(mu+sigma*(gamma(1-xi)-1)/xi)
+ }
+ else
+ return(E(as(object,"AbscontDistribution"), low=low, upp=upp, ...))
+ })
+
+setMethod("E", signature(object = "GEV",
+ fun = "function",
+ cond = "missing"),
+ function(object, fun, low = NULL, upp = NULL,
+ rel.tol= getdistrExOption("ErelativeTolerance"),
+ lowerTruncQuantile = getdistrExOption("ElowerTruncQuantile"),
+ upperTruncQuantile = getdistrExOption("EupperTruncQuantile"),
+ IQR.fac = max(1e4,getdistrExOption("IQR.fac")), ...
+ ){
+
+ dots <- list(...)
+ dots.withoutUseApply <- dots
+ useApply <- TRUE
+ if(!is.null(dots$useApply)) useApply <- dots$useApply
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/robast -r 468
More information about the Robast-commits
mailing list