[Distr-commits] r485 - branches/distr-2.2/pkg/distrMod/inst/scripts
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Mon Jun 29 09:55:32 CEST 2009
Author: stamats
Date: 2009-06-29 09:55:32 +0200 (Mon, 29 Jun 2009)
New Revision: 485
Added:
branches/distr-2.2/pkg/distrMod/inst/scripts/distrModExample.R
branches/distr-2.2/pkg/distrMod/inst/scripts/distrModExample1.R
Log:
Added examples which arose from R help question
Added: branches/distr-2.2/pkg/distrMod/inst/scripts/distrModExample.R
===================================================================
--- branches/distr-2.2/pkg/distrMod/inst/scripts/distrModExample.R (rev 0)
+++ branches/distr-2.2/pkg/distrMod/inst/scripts/distrModExample.R 2009-06-29 07:55:32 UTC (rev 485)
@@ -0,0 +1,113 @@
+## Example: skew-normal distribution
+library(distrMod)
+
+fsn <- function(x, loc = 0, scale = 1, shape = 0){
+ u <- (x-loc)/scale
+ dnorm(u)*pnorm(shape*u)/scale
+}
+rsn <- function(n, loc = 0, scale = 1, shape = 0){
+## borrowed from package sn
+ u1 <- rnorm(n)
+ u2 <- rnorm(n)
+ id <- (u2 > shape * u1)
+ u1[id] <- (-u1[id])
+ loc + scale * u1
+}
+
+##Class:parameters
+setClass("SkewNormParameter",
+ representation = representation(loc = "numeric",
+ scale = "numeric",
+ shape = "numeric"),
+ prototype = prototype(loc = 0, scale = 1, shape = 0,
+ name = gettext("Parameter of the Skew Normal distribution")),
+
+ contains = "Parameter",
+ validity = function(object){
+ if(length(object at loc) != 1)
+ stop("length of 'loc' is not equal to 1")
+ if(length(object at scale) != 1)
+ stop("length of 'scale' is not equal to 1")
+ if(length(object at shape) != 1)
+ stop("length of 'shape' is not equal to 1")
+ if(object at scale <= 0)
+ stop("'scale' has to be positive")
+ else return(TRUE)
+ }
+)
+
+##Class: distribution (created using the pdf of the skew normal defined above)
+setClass("SkewNorm",
+ prototype = prototype(d = function(x, log = FALSE){ fsn(x, loc=0, scale=1, shape=0) },
+ r = function(n){ rsn(n, loc = 0, scale = 1, shape = 0) },
+ img = new("Reals"),
+ param = new("SkewNormParameter"),
+ .logExact = FALSE,
+ .lowerExact = FALSE),
+ contains = "AbscontDistribution"
+)
+
+
+## you would need
+#if(!isGeneric("mu")){
+# setGeneric("mu", function(object) standardGeneric("mu"))
+#}
+## before
+#setMethod("mu", "SkewNormParameter", function(object) object at mu)
+
+## loc, scale and shape are already generic
+## accessor methods
+setMethod("loc", "SkewNormParameter", function(object) object at loc)
+setMethod("scale", c("SkewNormParameter", "missing", "missing"), function(x) x at scale)
+setMethod("shape", "SkewNormParameter", function(object) object at shape)
+
+## initialize method
+setMethod("initialize", "SkewNorm",
+ function(.Object, loc = 0, scale = 1, shape = 0) {
+ .Object at img <- Reals()
+ .Object at param <- new("SkewNormParameter", loc = loc, scale = scale, shape = shape,
+ name = gettext("Parameter of a skew norml distribution"))
+ .Object at r <- function(n){}
+ body(.Object at r) <- substitute({ rsn(n, loc = loc1, scale = scale1, shape = shape1) },
+ list(loc1 = loc, scale1 = scale, shape1 = shape))
+ .Object at d <- function(x, log = FALSE){}
+ body(.Object at d) <- substitute({ fsn(x, loc = loc1, scale = scale1, shape = shape1) },
+ list(loc1 = loc, scale1 = scale, shape1 = shape))
+ .withSim <- TRUE
+ r <- .Object at r
+ dpq <- RtoDPQ(r)
+ pfun <- dpq$pfun
+ up1 <- max(rN <- r(10^getdistrOption("RtoDPQ.e")))
+ low1 <- min(r(10^getdistrOption("RtoDPQ.e")))
+ h <- (up1 - low1)/getdistrOption("DefaultNrFFTGridPointsExponent")
+ x <- seq(from = low1, to = up1, by = h)
+ px.l <- pfun(x + 0.5 * h)
+ px.u <- pfun(x + 0.5 * h, lower.tail = FALSE)
+ qfun <- distr:::.makeQNew(x + 0.5 * h, px.l, px.u, FALSE, yL = -Inf, yR = Inf)
+ .Object at .withArith <- FALSE
+ .Object at .logExact <- FALSE
+ .Object at .lowerExact <- FALSE
+ .Object
+ })
+
+SN <- new("SkewNorm", loc = 0, scale = 1, shape = 1)
+
+## Generating Function
+SkewNorm <- function(loc = 0, scale = 1, shape = 0){
+ new("SkewNorm", loc = loc, scale = scale, shape = shape)
+}
+
+## skew normal location family
+theta <- 0
+names(theta) <- "loc"
+SNL <- ParamFamily(name = "Skew normal location family",
+ param = ParamFamParameter(name = "location parameter", main = theta),
+ distribution = SkewNorm(loc = 0, scale = 1, shape = 1), ## scale, shape known!
+ startPar = function(x,...) c(min(x),max(x)),
+ distrSymm <- NoSymmetry(),
+ modifyParam = function(theta){ SkewNorm(loc = theta, scale = 1, shape = 1) })
+SNL
+
+x <- r(SN)(50)
+est <- MLEstimator(x, SNL)
+
Added: branches/distr-2.2/pkg/distrMod/inst/scripts/distrModExample1.R
===================================================================
--- branches/distr-2.2/pkg/distrMod/inst/scripts/distrModExample1.R (rev 0)
+++ branches/distr-2.2/pkg/distrMod/inst/scripts/distrModExample1.R 2009-06-29 07:55:32 UTC (rev 485)
@@ -0,0 +1,105 @@
+## Example: skew-normal distribution
+library(distrMod)
+
+fsn <- function(x, loc = 0, scale = 1, shape = 0){
+ u <- (x-loc)/scale
+ dnorm(u)*pnorm(shape*u)/scale
+}
+rsn <- function(n, loc = 0, scale = 1, shape = 0){
+## borrowed from package sn
+ u1 <- rnorm(n)
+ u2 <- rnorm(n)
+ id <- (u2 > shape * u1)
+ u1[id] <- (-u1[id])
+ loc + scale * u1
+}
+
+##Class:parameters
+setClass("SkewNormParameter",
+ representation = representation(loc = "numeric",
+ scale = "numeric",
+ shape = "numeric"),
+ prototype = prototype(loc = 0, scale = 1, shape = 0,
+ name = gettext("Parameter of the Skew Normal distribution")),
+
+ contains = "Parameter",
+ validity = function(object){
+ if(length(object at loc) != 1)
+ stop("length of 'loc' is not equal to 1")
+ if(length(object at scale) != 1)
+ stop("length of 'scale' is not equal to 1")
+ if(length(object at shape) != 1)
+ stop("length of 'shape' is not equal to 1")
+ if(object at scale <= 0)
+ stop("'scale' has to be positive")
+ else return(TRUE)
+ }
+)
+
+##Class: distribution (created using the pdf of the skew normal defined above)
+setClass("SkewNorm",
+ prototype = prototype(d = function(x, log = FALSE){ fsn(x, loc=0, scale=1, shape=0) },
+ r = function(n){ rsn(n, loc = 0, scale = 1, shape = 0) },
+ img = new("Reals"),
+ param = new("SkewNormParameter"),
+ .logExact = FALSE,
+ .lowerExact = FALSE),
+ contains = "AbscontDistribution"
+)
+
+
+## you would need
+#if(!isGeneric("mu")){
+# setGeneric("mu", function(object) standardGeneric("mu"))
+#}
+## before
+#setMethod("mu", "SkewNormParameter", function(object) object at mu)
+
+## loc, scale and shape are already generic
+## accessor methods
+setMethod("loc", "SkewNormParameter", function(object) object at loc)
+setMethod("scale", c("SkewNormParameter", "missing", "missing"), function(x) x at scale)
+setMethod("shape", "SkewNormParameter", function(object) object at shape)
+
+## initialize method
+setMethod("initialize", "SkewNorm",
+ function(.Object, loc = 0, scale = 1, shape = 0) {
+ .Object at img <- Reals()
+ .Object at param <- new("SkewNormParameter", loc = loc, scale = scale, shape = shape,
+ name = gettext("Parameter of a skew norml distribution"))
+ .Object at r <- function(n){ stop("not yet implemented") }
+ .Object at d <- function(x, log = FALSE){}
+ body(.Object at d) <- substitute({ fsn(x, loc = loc1, scale = scale1, shape = shape1) },
+ list(loc1 = loc, scale1 = scale, shape1 = shape))
+ .withSim <- FALSE
+ .Object at .withArith <- FALSE
+ .Object at .logExact <- FALSE
+ .Object at .lowerExact <- FALSE
+ .Object
+ })
+
+SN <- new("SkewNorm", loc = 0, scale = 1, shape = 1)
+
+## Generating Function
+SkewNorm <- function(loc = 0, scale = 1, shape = 0){
+ new("SkewNorm", loc = loc, scale = scale, shape = shape)
+}
+
+## skew normal location family
+theta <- 0
+names(theta) <- "loc"
+SNL <- ParamFamily(name = "Skew normal location family",
+ param = ParamFamParameter(name = "location parameter", main = theta),
+ distribution = SkewNorm(loc = 0, scale = 1, shape = 1), ## scale, shape known!
+ startPar = function(x,...) c(min(x),max(x)),
+ distrSymm <- NoSymmetry(),
+ modifyParam = function(theta){ SkewNorm(loc = theta, scale = 1, shape = 1) })
+SNL
+
+## some data
+x <- rnorm(50)
+est <- MLEstimator(x, SNL)
+
+## error!
+r(SN)(50)
+
More information about the Distr-commits
mailing list