[Distr-commits] r803 - in branches/distr-2.4/pkg/distrMod: . R man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Thu Apr 26 18:06:13 CEST 2012
Author: ruckdeschel
Date: 2012-04-26 18:06:13 +0200 (Thu, 26 Apr 2012)
New Revision: 803
Added:
branches/distr-2.4/pkg/distrMod/R/GParetoFamily.R
branches/distr-2.4/pkg/distrMod/man/GParetoFamily.Rd
Modified:
branches/distr-2.4/pkg/distrMod/NAMESPACE
branches/distr-2.4/pkg/distrMod/R/AllClass.R
branches/distr-2.4/pkg/distrMod/R/AllReturnClasses.R
branches/distr-2.4/pkg/distrMod/man/InternalReturnClasses-class.Rd
branches/distr-2.4/pkg/distrMod/man/validParameter-methods.Rd
Log:
moved GParetoFamily from Nataliya's repostitory for her PhD Thesis to distrMod and wrote Rd file
Modified: branches/distr-2.4/pkg/distrMod/NAMESPACE
===================================================================
--- branches/distr-2.4/pkg/distrMod/NAMESPACE 2012-04-05 00:05:57 UTC (rev 802)
+++ branches/distr-2.4/pkg/distrMod/NAMESPACE 2012-04-26 16:06:13 UTC (rev 803)
@@ -25,7 +25,7 @@
exportClasses("BinomFamily","PoisFamily", "NormLocationFamily",
"GumbelLocationFamily", "NormScaleFamily", "ExpScaleFamily",
"LnormScaleFamily", "GammaFamily", "BetaFamily", "NormLocationScaleFamily",
- "CauchyLocationScaleFamily")
+ "CauchyLocationScaleFamily", "GParetoFamily")
exportClasses("NormType", "QFNorm", "InfoNorm", "SelfNorm")
exportClasses("Estimate", "MCEstimate")
exportClasses("Confint")
@@ -67,7 +67,8 @@
"BinomFamily", "PoisFamily", "NbinomFamily", "NormLocationFamily",
"GumbelLocationFamily", "NormScaleFamily", "ExpScaleFamily",
"LnormScaleFamily", "GammaFamily", "BetaFamily", "NormLocationScaleFamily",
- "CauchyLocationScaleFamily", "NbinomwithSizeFamily", "NbinomMeanSizeFamily" )
+ "CauchyLocationScaleFamily", "NbinomwithSizeFamily", "NbinomMeanSizeFamily",
+ "GParetoFamily")
export("asCov", "trAsCov", "asHampel", "asBias", "asMSE", "asUnOvShoot",
"fiCov", "trFiCov", "fiHampel", "fiMSE", "fiBias", "fiUnOvShoot")
export("positiveBias", "negativeBias", "symmetricBias",
Modified: branches/distr-2.4/pkg/distrMod/R/AllClass.R
===================================================================
--- branches/distr-2.4/pkg/distrMod/R/AllClass.R 2012-04-05 00:05:57 UTC (rev 802)
+++ branches/distr-2.4/pkg/distrMod/R/AllClass.R 2012-04-26 16:06:13 UTC (rev 803)
@@ -237,6 +237,7 @@
fam.call = call("new", "L2LocationScaleFamily")),
contains = "L2LocationScaleUnion")
+
################################################################################
## Norm Classes
################################################################################
Modified: branches/distr-2.4/pkg/distrMod/R/AllReturnClasses.R
===================================================================
--- branches/distr-2.4/pkg/distrMod/R/AllReturnClasses.R 2012-04-05 00:05:57 UTC (rev 802)
+++ branches/distr-2.4/pkg/distrMod/R/AllReturnClasses.R 2012-04-26 16:06:13 UTC (rev 803)
@@ -62,5 +62,7 @@
## Cauchy location scale family
setClass("CauchyLocationScaleFamily",
contains = "L2LocationScaleFamily")
+## class
+setClass("GParetoFamily", contains="L2ParamFamily")
Added: branches/distr-2.4/pkg/distrMod/R/GParetoFamily.R
===================================================================
--- branches/distr-2.4/pkg/distrMod/R/GParetoFamily.R (rev 0)
+++ branches/distr-2.4/pkg/distrMod/R/GParetoFamily.R 2012-04-26 16:06:13 UTC (rev 803)
@@ -0,0 +1,346 @@
+#################################
+##
+## Class: GParetoFamily
+##
+################################
+
+
+## methods
+setMethod("validParameter",signature(object="GParetoFamily"),
+ function(object, param, tol =.Machine$double.eps){
+ if (is(param, "ParamFamParameter"))
+ param <- main(param)
+ if (!all(is.finite(param)))
+ return(FALSE)
+ if (any(param[1] <= tol))
+ return(FALSE)
+ if (any(param[2] <= tol))
+ return(FALSE)
+ return(TRUE)
+ })
+
+
+## generating function
+## loc: known/fixed threshold/location parameter
+## scale: scale parameter
+## shape: shape parameter
+## of.interest: which parameters, transformations are of interest
+## posibilites are: scale, shape, quantile, expected loss, expected shortfall
+## a maximum number of two of these may be selected
+## p: probability needed for quantile and expected shortfall
+## N: expected frequency for expected loss
+## trafo: optional parameter transformation
+## start0Est: startEstimator for MLE and MDE --- if NULL HybridEstimator is used;
+
+GParetoFamily <- function(loc = 0, scale = 1, shape = 0.5,
+ of.interest = c("scale", "shape"),
+ p = NULL, N = NULL, trafo = NULL,
+ start0Est = NULL){
+ 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
+ distrSymm <- NoSymmetry()
+
+ ## parameters
+ names(theta) <- c("loc", "scale", "shape")
+
+ 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) <- 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))
+ }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))
+ 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))
+ }
+ }
+ 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))
+ }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))
+ 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))
+ }
+ }
+ 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))
+ }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))
+ 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))
+ }
+ }
+ 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 = theta[2:3],
+ fixed = theta[1],
+ trafo = trafo)
+
+ ## distribution
+ distribution <- GPareto(loc = loc, scale = scale, shape = shape)
+
+ ## starting parameters
+ startPar <- function(x,...){
+ tr <- theta[1]
+
+ if(any(x < tr)) stop("some data smaller than 'loc' parameter")
+
+ ## Pickand estimator
+ if(is.null(start0Est)){
+ e0 <- HybridEstimator(x, up=300, ...)
+ }else{
+ if(is(start0Est,"function")){
+ e1 <- start0Est(x, ...)
+ e0 <- if(is(e1,"Estimate")) estimate(e1) else e1
+ }
+ if(!is.null(names(e0)))
+ e0 <- e0[c("scale", "shape")]
+ }
+ names(e0) <- NULL
+ return(e0)
+ }
+
+ modifyPar <- function(theta){
+ theta <- abs(theta)
+ GPareto(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) {
+ sc <- force(main(param)[1])
+ k <- force(main(param)[2])
+ tr <- fixed(param)[1]
+
+ 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
+ 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)
+ return(y)
+ }
+ ## additional centering of scores to increase numerical precision!
+ z1 <- E(distribution, fun=Lambda1)
+ z2 <- E(distribution, fun=Lambda2)
+ return(list(function(x){ Lambda1(x)-z1 },function(x){ Lambda2(x)-z2 }))
+ }
+
+ ## Fisher Information matrix as a function of parameters
+ FisherInfo.fct <- function(param) {
+ 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)
+ E11 <- sc^-2
+ E12 <- (sc*(1+k))^-1
+ E22 <- 2/(1+k)
+ return(PosSemDefSymmMatrix(matrix(c(E11,E12,E12,E22)/(1+2*k),2,2)))
+ }
+
+ FisherInfo <- FisherInfo.fct(param)
+ name <- "Generalized Pareto 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 = "GParetoFamily")
+ f.call <- substitute(GParetoFamily(loc = loc0, scale = scale0, shape = shape0,
+ of.interest = of.interest0, p = p0,
+ N = N0, trafo = trafo0),
+ list(loc0 = loc, scale0 = scale, shape0 = shape,
+ of.interest0 = of.interest, p0 = p, N0 = N,
+ trafo0 = trafo))
+ res at fam.call <- f.call
+ return(res)
+}
+
Added: branches/distr-2.4/pkg/distrMod/man/GParetoFamily.Rd
===================================================================
--- branches/distr-2.4/pkg/distrMod/man/GParetoFamily.Rd (rev 0)
+++ branches/distr-2.4/pkg/distrMod/man/GParetoFamily.Rd 2012-04-26 16:06:13 UTC (rev 803)
@@ -0,0 +1,44 @@
+\name{GParetoFamily}
+\alias{GParetoFamily}
+
+\title{Generating function for Generalized Pareto families}
+\description{
+ Generates an object of class \code{"L2ParamFamily"} which
+ represents a Generalized Pareto family.
+}
+\usage{
+GParetoFamily(loc = 0, scale = 1, shape = 0.5, of.interest = c("scale", "shape"),
+ p = NULL, N = NULL, trafo = NULL, start0Est = NULL)
+}
+\arguments{
+ \item{loc}{ real: known/fixed threshold/location parameter }
+ \item{scale}{ positive real: scale parameter }
+ \item{shape}{ positive real: shape parameter }
+ \item{of.interest}{ character: which parameters, transformations are of interest.\cr
+ posibilites are: "scale", "shape", "quantile", "expected loss",
+ "expected shortfall"; a maximum number of two of these may be selected }
+ \item{p}{real or NULL: probability needed for quantile and expected shortfall }
+ \item{N}{real or NULL: expected frequency for expected loss }
+ \item{trafo}{ matrix or NULL: transformation of the parameter }
+ \item{start0Est}{ startEstimator for MLE and MDE --- if NULL HybridEstimator is used }
+}
+\details{
+ The slots of the corresponding L2 differentiable
+ parameteric family are filled.
+}
+\value{Object of class \code{"L2ParamFamily"}}
+\references{
+ Kohl, M. (2005) \emph{Numerical Contributions to
+ the Asymptotic Theory of Robustness}. Bayreuth: Dissertation.}
+\author{Matthias Kohl \email{Matthias.Kohl at stamats.de}\cr
+ Peter Ruckdeschel \email{peter.ruckdeschel at itwm.fraunhofer.de}\cr
+ Nataliya Horbenko \email{nataliya.horbenko at itwm.fraunhofer.de}}
+%\note{}
+\seealso{\code{\link{L2ParamFamily-class}}, \code{\link[distrEx]{GPareto-class}}}
+\examples{
+(G1 <- GParetoFamily())
+FisherInfo(G1)
+checkL2deriv(G1)
+}
+\concept{Generalized Pareto model}
+\keyword{models}
Modified: branches/distr-2.4/pkg/distrMod/man/InternalReturnClasses-class.Rd
===================================================================
--- branches/distr-2.4/pkg/distrMod/man/InternalReturnClasses-class.Rd 2012-04-05 00:05:57 UTC (rev 802)
+++ branches/distr-2.4/pkg/distrMod/man/InternalReturnClasses-class.Rd 2012-04-26 16:06:13 UTC (rev 803)
@@ -10,6 +10,7 @@
\alias{ExpScaleFamily-class}
\alias{LnormScaleFamily-class}
\alias{GammaFamily-class}
+\alias{GParetoFamily-class}
\alias{BetaFamily-class}
\alias{NormLocationScaleFamily-class}
\alias{CauchyLocationScaleFamily-class}
@@ -20,7 +21,7 @@
\section{Described classes}{
In this file we describe classes
-\code{BinomFamily}, \code{PoisFamily}, \code{GammaFamily}
+\code{BinomFamily}, \code{PoisFamily}, \code{GammaFamily}, \code{GParetoFamily},
\code{BetaFamily} ``extending'' (no new slots!) class \code{L2ParamFamily},
classes \code{NormLocationFamily} and \code{GumbelLocationFamily},
``extending'' (no new slots!) class \code{"L2LocationFamily"}, classes
Modified: branches/distr-2.4/pkg/distrMod/man/validParameter-methods.Rd
===================================================================
--- branches/distr-2.4/pkg/distrMod/man/validParameter-methods.Rd 2012-04-05 00:05:57 UTC (rev 802)
+++ branches/distr-2.4/pkg/distrMod/man/validParameter-methods.Rd 2012-04-26 16:06:13 UTC (rev 803)
@@ -9,6 +9,7 @@
\alias{validParameter,BinomFamily-method}
\alias{validParameter,PoisFamily-method}
\alias{validParameter,GammaFamily-method}
+\alias{validParameter,GParetoFamily-method}
\title{ Methods for function validParameter in Package `distrMod' }
@@ -27,6 +28,7 @@
\S4method{validParameter}{BinomFamily}(object, param, tol=.Machine$double.eps)
\S4method{validParameter}{PoisFamily}(object, param, tol=.Machine$double.eps)
\S4method{validParameter}{GammaFamily}(object, param, tol=.Machine$double.eps)
+\S4method{validParameter}{GParetoFamily}(object, param, tol=.Machine$double.eps)
}
\details{
@@ -52,6 +54,9 @@
\item{\code{GammaFamily}}{checks if both parameters are finite by \code{is.finite},
if their length is 1 or 2 (e.g. if one features as nuisance parameter), and if
both are strictly larger than 0 (upto argument \code{tol})}
+ \item{\code{GParetoFamily}}{checks if both parameters are finite by \code{is.finite},
+ if their length is 1 or 2 (e.g. if one features as nuisance parameter), and if
+ both are strictly larger than 0 (upto argument \code{tol})}
}
}
\arguments{
More information about the Distr-commits
mailing list