[Distr-commits] r662 - in branches/distr-2.3/pkg/distrEx: . R man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Fri Jul 9 16:00:49 CEST 2010
Author: horbenko
Date: 2010-07-09 16:00:49 +0200 (Fri, 09 Jul 2010)
New Revision: 662
Added:
branches/distr-2.3/pkg/distrEx/R/GEV.R
branches/distr-2.3/pkg/distrEx/man/GEV.Rd
Modified:
branches/distr-2.3/pkg/distrEx/NAMESPACE
branches/distr-2.3/pkg/distrEx/R/AllClass.R
branches/distr-2.3/pkg/distrEx/R/AllInitialize.R
branches/distr-2.3/pkg/distrEx/R/Expectation.R
branches/distr-2.3/pkg/distrEx/R/Functionals.R
branches/distr-2.3/pkg/distrEx/R/Kurtosis.R
branches/distr-2.3/pkg/distrEx/R/Skewness.R
Log:
Generalized Extreme Value Distribution is ported to distrEx
Modified: branches/distr-2.3/pkg/distrEx/NAMESPACE
===================================================================
--- branches/distr-2.3/pkg/distrEx/NAMESPACE 2010-07-08 17:15:46 UTC (rev 661)
+++ branches/distr-2.3/pkg/distrEx/NAMESPACE 2010-07-09 14:00:49 UTC (rev 662)
@@ -45,7 +45,7 @@
export("LMParameter")
export("DiscreteMVDistribution",
"LMCondDistribution",
- "Gumbel", "Pareto", "GPareto")
+ "Gumbel", "Pareto", "GPareto", "GEV")
export("ConvexContamination")
export("GLIntegrate",
"distrExIntegrate")
Modified: branches/distr-2.3/pkg/distrEx/R/AllClass.R
===================================================================
--- branches/distr-2.3/pkg/distrEx/R/AllClass.R 2010-07-08 17:15:46 UTC (rev 661)
+++ branches/distr-2.3/pkg/distrEx/R/AllClass.R 2010-07-09 14:00:49 UTC (rev 662)
@@ -256,3 +256,56 @@
.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"
+ )
+
Modified: branches/distr-2.3/pkg/distrEx/R/AllInitialize.R
===================================================================
--- branches/distr-2.3/pkg/distrEx/R/AllInitialize.R 2010-07-08 17:15:46 UTC (rev 661)
+++ branches/distr-2.3/pkg/distrEx/R/AllInitialize.R 2010-07-09 14:00:49 UTC (rev 662)
@@ -162,3 +162,74 @@
.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
+ return(-log(1+shapeSub*q0)/shapeSub)
+ }else{
+ p0 <- pgev(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{
+
+ ## analogous to GPD
+ 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)
+ 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
+ })
+
Modified: branches/distr-2.3/pkg/distrEx/R/Expectation.R
===================================================================
--- branches/distr-2.3/pkg/distrEx/R/Expectation.R 2010-07-08 17:15:46 UTC (rev 661)
+++ branches/distr-2.3/pkg/distrEx/R/Expectation.R 2010-07-09 14:00:49 UTC (rev 662)
@@ -843,6 +843,65 @@
})
+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>=1){ return(Inf)}
+ if(xi==0) {return(mu + sigma*(-digamma(1)))}
+ if((xi!=0)&&(xi<1)){ 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
+ dots.withoutUseApply$useApply <- NULL
+ integrand <- function(x, dfun, ...){ di <- dim(x)
+ y <- q(object)##quantile transformation
+ if(useApply){
+ funy <- sapply(y,fun, ...)
+ dim(y) <- di
+ dim(funy) <- di
+ }else funy <- fun(y,...)
+ return(funy) }
+
+# if(is.null(low)) low <- 0
+# if(is.null(upp)) upp <- 1
+#
+# Ib <- .getIntbounds(object, low, upp, #lowerTruncQuantile,
+# upperTruncQuantile, IQR.fac)
+# low <- if(Ib["low"]<=0) -Inf else log(Ib["low"])
+# upp <- log(Ib["upp"])
+
+ low <- 0
+ upp <- 1
+ return(do.call(distrExIntegrate, c(list(f = integrand,
+ lower = low,
+ upper = upp,
+ rel.tol = rel.tol,
+ distr = object, dfun = d(object)), dots.withoutUseApply)))
+
+ })
+
+
+
############################ Expectation for UnivarLebDecDistribution
### merged from Expectation_LebDec.R on Apr 15 2009
setMethod("E", signature(object = "UnivarLebDecDistribution",
Modified: branches/distr-2.3/pkg/distrEx/R/Functionals.R
===================================================================
--- branches/distr-2.3/pkg/distrEx/R/Functionals.R 2010-07-08 17:15:46 UTC (rev 661)
+++ branches/distr-2.3/pkg/distrEx/R/Functionals.R 2010-07-09 14:00:49 UTC (rev 662)
@@ -527,6 +527,23 @@
}})
### source http://en.wikipedia.org/wiki/Pareto_distribution
+
+setMethod("var", signature(x = "GEV"),
+ function(x, ...){
+ dots <- match.call(call = sys.call(sys.parent(1)),
+ expand.dots = FALSE)$"..."
+ fun <- NULL; cond <- NULL; low <- NULL; upp <- NULL
+ if(hasArg(low)) low <- dots$low
+ if(hasArg(upp)) upp <- dots$upp
+ if(hasArg(fun)||hasArg(cond)||!is.null(low)||!is.null(upp))
+ return(var(as(x,"AbscontDistribution"),...))
+ else{ xi <- shape(x); sigma <- scale(x)
+ if(xi>=1/2) return(NA)
+ if(xi==0) return(pi^2/6)
+ if((xi!=0)&&(xi<1/2))return(sigma^2*(gamma(1-2*xi)-gamma(1-xi)^2)/xi^2)
+ }})
+### http://en.wikipedia.org/wiki/Generalized_extreme_value_distribution
+
#################################################################
# some exact medians
#################################################################
@@ -573,6 +590,10 @@
function(x) {k <- shape(x); mu <- loc(x); s <- scale(x)
return(mu + s*(2^k-1)/k)
})
+setMethod("median", signature(x = "GEV"),
+ function(x) {xi <- shape(x); mu <- loc(x); sigma <- scale(x)
+ return(mu + sigma*(log(2)^(-xi)-1)/xi)
+ })
#################################################################
# some exact IQRs
@@ -618,6 +639,10 @@
function(x) {k <- shape(x); s<- scale(x)
return(s/k*4^k*(1-3^(-k)))
})
+setMethod("IQR", signature(x = "GEV"),
+ function(x) {xi <- shape(x); sigma<- scale(x)
+ return(sigma*((log(4/3))^(-xi)-(log(4))^(-xi))/xi)
+ })
#################################################################
# some exact mads
#################################################################
Added: branches/distr-2.3/pkg/distrEx/R/GEV.R
===================================================================
--- branches/distr-2.3/pkg/distrEx/R/GEV.R (rev 0)
+++ branches/distr-2.3/pkg/distrEx/R/GEV.R 2010-07-09 14:00:49 UTC (rev 662)
@@ -0,0 +1,88 @@
+###########################################
+## Class: Generalized Extreme Value Distribution
+##
+## @param: location, scale, shape
+###########################################
+
+## access methods
+setMethod("location", "GEVParameter", function(object) object at loc)
+setMethod("loc", "GEVParameter", function(object) object at loc)
+setMethod("scale", "GEVParameter",
+ function(x, center = TRUE, scale = TRUE) x at scale)
+ ### odd arg-list due to existing function in base package
+setMethod("shape", "GEVParameter", function(object) object at shape)
+
+## replace Methods
+setReplaceMethod("loc", "GEVParameter",
+ function(object, value){ object at loc <- value; object })
+setReplaceMethod("location", "GEVParameter",
+ function(object, value){ object at loc <- value; object })
+setReplaceMethod("scale", "GEVParameter",
+ function(object, value){
+ if(length(value) != 1 || value <= 0)
+ stop("'value' has to be a single positive real number!")
+ object at scale <- value; object})
+setReplaceMethod("shape", "GEVParameter",
+ function(object, value){ object at shape <- value; object})
+
+## wrapped access methods
+setMethod("location", "GEV", function(object) loc(object at param))
+setMethod("loc", "GEV", function(object) loc(object at param))
+setMethod("scale", "GEV",
+ function(x, center = TRUE, scale = TRUE)
+ scale(x at param))
+setMethod("shape", "GEV", function(object) shape(object at param))
+
+
+## wrapped replace methods
+setMethod("loc<-", "GEV", function(object, value)
+ new("GEV", loc = value, scale = scale(object), shape = shape(object)))
+setMethod("location<-", "GEV", function(object, value)
+ new("GEV", loc = value, scale = scale(object), shape = shape(object)))
+setMethod("scale<-", "GEV", function(object, value)
+ new("GEV", loc = loc(object), scale = value, shape = shape(object)))
+setMethod("shape<-", "GEV", function(object, value)
+ new("GEV", loc = loc(object), scale = scale(object), shape = value))
+
+setValidity("GEVParameter", function(object){
+ if(length(loc(object)) != 1)
+ stop("location has to be a numeric of length 1")
+ if(length(scale(object)) != 1)
+ stop("scale has to be a numeric of length 1")
+ if(scale(object) <= 0)
+ stop("scale has to be positive")
+ if(length(shape(object)) != 1)
+ stop("shape has to be a numeric of length 1")
+# if(shape(object) < 0)
+# stop("shape has to be non-negative")
+ else return(TRUE)
+})
+
+## generating function
+GEV <- function(loc = 0, scale = 1, shape = 0, location = loc){
+ if(!missing(loc)&&!missing(location))
+ if(!isTRUE(all.equal(loc,location)))
+ stop("Only one of arguments 'loc' and 'location' may be used.")
+ if(!missing(location)) loc <- location
+ if(abs(shape) < .Machine$double.eps) return(Gumbel(loc=loc,scale=scale))
+ new("GEV", loc = loc, scale = scale, shape = shape) }
+
+## extra methods for GEV distribution
+setMethod("+", c("GEV","numeric"),
+ function(e1, e2){
+ if (length(e2)>1) stop("length of operator must be 1")
+ new("GEV", loc = loc(e1) + e2, scale = scale(e1), shape=shape(e1))
+ })
+
+setMethod("*", c("GEV","numeric"),
+ function(e1, e2){
+ if (length(e2)>1) stop("length of operator must be 1")
+ if (isTRUE(all.equal(e2,0)))
+ return(new("Dirac", location = 0, .withArith = TRUE))
+ GEV <- new("GEV", loc = loc(e1) * abs(e2),
+ scale = scale(e1)*abs(e2), shape=shape(e1))
+ if(e2<0) GEV <- (-1)*as(GEV,"AbscontDistribution")
+ return(GEV)
+ })
+
+
Modified: branches/distr-2.3/pkg/distrEx/R/Kurtosis.R
===================================================================
--- branches/distr-2.3/pkg/distrEx/R/Kurtosis.R 2010-07-08 17:15:46 UTC (rev 661)
+++ branches/distr-2.3/pkg/distrEx/R/Kurtosis.R 2010-07-09 14:00:49 UTC (rev 662)
@@ -460,3 +460,21 @@
}
})
### source Maple ...
+
+setMethod("kurtosis", signature(x = "GEV"),
+ function(x, ...){
+ dots <- match.call(call = sys.call(sys.parent(1)),
+ expand.dots = FALSE)$"..."
+ fun <- NULL; cond <- NULL; low <- NULL; upp <- NULL
+ if(hasArg(low)) low <- dots$low
+ if(hasArg(upp)) upp <- dots$upp
+ if(hasArg(fun)||hasArg(cond)||!is.null(low)||!is.null(upp))
+ return(kurtosis(as(x,"AbscontDistribution"),...))
+ else{
+ xi <- shape(x)
+ if(xi>=1/4) return(NA)
+ else
+ return((gamma(1-4*xi)- 4*gamma(1-xi)*gamma(1-3*xi)+6*gamma(1-2*xi)*gamma(1-xi)^2 - 3*gamma(1-xi)^4)/(gamma(1-2*xi)-gamma(1-xi)^2)^(2))
+ }
+})
+## source http://en.wikipedia.org/wiki/Generalized_extreme_value_distribution
\ No newline at end of file
Modified: branches/distr-2.3/pkg/distrEx/R/Skewness.R
===================================================================
--- branches/distr-2.3/pkg/distrEx/R/Skewness.R 2010-07-08 17:15:46 UTC (rev 661)
+++ branches/distr-2.3/pkg/distrEx/R/Skewness.R 2010-07-09 14:00:49 UTC (rev 662)
@@ -422,3 +422,22 @@
})
### source Maple...
+setMethod("skewness", signature(x = "GPareto"),
+ function(x, ...){
+ dots <- match.call(call = sys.call(sys.parent(1)),
+ expand.dots = FALSE)$"..."
+ fun <- NULL; cond <- NULL; low <- NULL; upp <- NULL
+ if(hasArg(low)) low <- dots$low
+ if(hasArg(upp)) upp <- dots$upp
+ if(hasArg(fun)||hasArg(cond)||!is.null(low)||!is.null(upp))
+ return(skewness(as(x,"AbscontDistribution"),...))
+ else{
+ xi <- shape(x)
+ if(xi>=1/3) return(NA)
+ else
+ return((gamma(1-3*xi)-3*gamma(1-xi)*gamma(1-2*xi) + 2*gamma(1-xi)^3)/(gamma(1-2*xi)-gamma(1-xi)^2)^(3/2))
+ }
+})
+
+### source http://en.wikipedia.org/wiki/Generalized_extreme_value_distribution
+
Added: branches/distr-2.3/pkg/distrEx/man/GEV.Rd
===================================================================
--- branches/distr-2.3/pkg/distrEx/man/GEV.Rd (rev 0)
+++ branches/distr-2.3/pkg/distrEx/man/GEV.Rd 2010-07-09 14:00:49 UTC (rev 662)
@@ -0,0 +1,39 @@
+\name{GEV}
+\alias{GEV}
+
+\title{Generating function for GEV-class}
+\description{
+ Generates an object of class \code{"GEV"}.
+}
+\usage{GEV(loc = 0, scale = 1, shape = 0, location = loc)}
+\arguments{
+ \item{loc}{ real number: location parameter of
+ the GEV distribution. }
+ \item{scale}{ positive real number: scale parameter
+ of the GEV distribution }
+ \item{shape}{ non-negative real number: shape parameter of
+ the GEV distribution. }
+ \item{location}{ real number: location of GEV distribution }
+
+}
+%\details{Generalized Extreme Value Distribution}
+\value{Object of class \code{"GEV"}}
+%\references{}
+\author{Nataliya Horbenko \email{Nataliya.Horbenko at itwm.fraunhofer.de}}
+\note{The class \code{"GEV"} is based on the code provided
+ by the package \pkg{evd} by Alec Stephenson.}
+\seealso{\code{\link{GEV-class}}, \code{\link[evd:gpd]{dgpd}}}
+\examples{
+(P1 <- GEV(loc = 0, scale = 1, shape = 0))
+plot(P1)
+
+E(GEV())
+E(P1, function(x){x^2})
+
+}
+
+\concept{GEV}
+\keyword{distribution}
+\concept{absolutely continuous distribution}
+\concept{GEV distribution}
+\concept{generating function}
More information about the Distr-commits
mailing list