[Distr-commits] r489 - branches/distr-2.2/pkg/distrEx/R

[noreply at r-forge.r-project.org]

Author: stamats
Date: 2009-06-30 06:30:50 +0200 (Tue, 30 Jun 2009)
New Revision: 489

Modified:
   branches/distr-2.2/pkg/distrEx/R/Expectation.R
Log:
included new methods for E and GPareto as well as for E and Gammad which increase the precision for the numerical calculation of certain integrals; in particular, integrals involving "log" as integrand.

Modified: branches/distr-2.2/pkg/distrEx/R/Expectation.R
===================================================================
--- branches/distr-2.2/pkg/distrEx/R/Expectation.R	2009-06-30 04:05:59 UTC (rev 488)
+++ branches/distr-2.2/pkg/distrEx/R/Expectation.R	2009-06-30 04:30:50 UTC (rev 489)
@@ -585,6 +585,47 @@
 
 ### source http://mathworld.wolfram.com/GammaDistribution.html
 
+setMethod("E", signature(object = "Gammad",
+                         fun = "function",
+                         cond = "missing"),
+    function(object, fun, low = NULL, upp = NULL,
+             rel.tol= getdistrExOption("ErelativeTolerance"),
+             lowerTruncQuantile = getdistrExOption("ElowerTruncQuantile"),
+             upperTruncQuantile = getdistrExOption("EupperTruncQuantile"),
+             IQR.fac = 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 = "Geom", 
                          fun = "missing", 
                          cond = "missing"),
@@ -755,7 +796,47 @@
 
 ### 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 = 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)))
+
+    })
+
+
 ############################ Expectation for UnivarLebDecDistribution
 ### merged from Expectation_LebDec.R on Apr 15 2009
 setMethod("E", signature(object = "UnivarLebDecDistribution",