[Distr-commits] r411 - in branches/distr-2.1/pkg/distrEx: . R chm man src
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Tue Mar 17 17:40:29 CET 2009
Author: ruckdeschel
Date: 2009-03-17 17:40:29 +0100 (Tue, 17 Mar 2009)
New Revision: 411
Added:
branches/distr-2.1/pkg/distrEx/R/AsymTotalVarDist.R
branches/distr-2.1/pkg/distrEx/chm/AsymTotalVarDist.html
branches/distr-2.1/pkg/distrEx/man/AsymTotalVarDist.Rd
Modified:
branches/distr-2.1/pkg/distrEx/NAMESPACE
branches/distr-2.1/pkg/distrEx/R/AllClass.R
branches/distr-2.1/pkg/distrEx/R/AllGeneric.R
branches/distr-2.1/pkg/distrEx/R/AllInitialize.R
branches/distr-2.1/pkg/distrEx/R/HellingerDist.R
branches/distr-2.1/pkg/distrEx/R/Internalfunctions.R
branches/distr-2.1/pkg/distrEx/R/TotalVarDist.R
branches/distr-2.1/pkg/distrEx/chm/00Index.html
branches/distr-2.1/pkg/distrEx/chm/HellingerDist.html
branches/distr-2.1/pkg/distrEx/chm/TotalVarDist.html
branches/distr-2.1/pkg/distrEx/chm/distrEx.chm
branches/distr-2.1/pkg/distrEx/chm/distrEx.hhp
branches/distr-2.1/pkg/distrEx/chm/distrEx.toc
branches/distr-2.1/pkg/distrEx/man/HellingerDist.Rd
branches/distr-2.1/pkg/distrEx/man/TotalVarDist.Rd
branches/distr-2.1/pkg/distrEx/src/distrEx.dll
Log:
introduced new asymmetric total variation distance
(for joint work with Claudio Agostinelli and Matthias Brandl)
TotalVarDist and HellingerDist gain extra arguments
to better control the integration range and exactness
Modified: branches/distr-2.1/pkg/distrEx/NAMESPACE
===================================================================
--- branches/distr-2.1/pkg/distrEx/NAMESPACE 2009-03-17 15:04:15 UTC (rev 410)
+++ branches/distr-2.1/pkg/distrEx/NAMESPACE 2009-03-17 16:40:29 UTC (rev 411)
@@ -20,6 +20,7 @@
exportMethods("ConvexContamination",
"ContaminationSize",
"TotalVarDist",
+ "AsymTotalVarDist",
"KolmogorovDist",
"HellingerDist",
"CvMDist")
Modified: branches/distr-2.1/pkg/distrEx/R/AllClass.R
===================================================================
--- branches/distr-2.1/pkg/distrEx/R/AllClass.R 2009-03-17 15:04:15 UTC (rev 410)
+++ branches/distr-2.1/pkg/distrEx/R/AllClass.R 2009-03-17 16:40:29 UTC (rev 411)
@@ -102,7 +102,22 @@
prototype = prototype(r = function(n){ rgumbel(n, loc = 0, scale = 1) },
d = function(x, ...){ dgumbel(x, loc = 0, scale = 1, ...) },
p = function(q, ...){ pgumbel(q, loc = 0, scale = 1, ...) },
- q = function(p, ...){ qgumbel(p, loc = 0, scale = 1, ...) },
+ q = function(p, ...){
+ mc <- as.list(match.call(call = sys.call())[-1])
+ lower.tail <- mc$lower.tail
+ if(is.null(lower.tail)) lower.tail <- TRUE
+ in01 <- (p>1 | p<0)
+ i01 <- .isEqual01(p)
+ i0 <- (i01 & p<1)
+ i1 <- (i01 & p>0)
+ ii01 <- .isEqual01(p) | in01
+ p0 <- p
+ p0[ii01] <- 0.5
+ q1 <- qgumbel(p0, loc = 0, scale = 1, ...)
+ q[i0] <- if(lower.tail) -Inf else Inf
+ q[i1] <- if(!lower.tail) -Inf else Inf
+ q[in01] <- NaN
+ },
img = new("Reals"),
param = new("GumbelParameter"),
.logExact = TRUE,
Modified: branches/distr-2.1/pkg/distrEx/R/AllGeneric.R
===================================================================
--- branches/distr-2.1/pkg/distrEx/R/AllGeneric.R 2009-03-17 15:04:15 UTC (rev 410)
+++ branches/distr-2.1/pkg/distrEx/R/AllGeneric.R 2009-03-17 16:40:29 UTC (rev 411)
@@ -137,6 +137,11 @@
function(e1, e2, size) standardGeneric("ConvexContamination"))
}
+if(!isGeneric("AsymTotalVarDist")){
+ setGeneric("AsymTotalVarDist",
+ function(e1, e2, ...) standardGeneric("AsymTotalVarDist"))
+}
+
if(!isGeneric("E")){
setGeneric("E", function(object, fun, cond, ...) standardGeneric("E"))
}
Modified: branches/distr-2.1/pkg/distrEx/R/AllInitialize.R
===================================================================
--- branches/distr-2.1/pkg/distrEx/R/AllInitialize.R 2009-03-17 15:04:15 UTC (rev 410)
+++ branches/distr-2.1/pkg/distrEx/R/AllInitialize.R 2009-03-17 16:40:29 UTC (rev 411)
@@ -13,24 +13,29 @@
.Object at p <- function(q, ...){ pgumbel(q, loc = loc1, scale = scale1, ...) }
body(.Object at p) <- substitute({ pgumbel(q, loc = loc1, scale = scale1, ...) },
list(loc1 = loc, scale1 = scale))
- .Object at q <- function(p, ...){
+ .Object at q <- function(p, loc = loc1, scale = scale1, ...){}
+ body(.Object at q) <- substitute({
## P.R.: changed to vectorized form
+ mc <- as.list(match.call(call = sys.call())[-1])
+ lower.tail <- mc$lower.tail
+ if(is.null(lower.tail)) lower.tail <- TRUE
+
+ in01 <- (p>1 | p<0)
+ i01 <- .isEqual01(p)
+ i0 <- (i01 & p<1)
+ i1 <- (i01 & p>0)
+ ii01 <- .isEqual01(p) | in01
+
p0 <- p
- p0[.isEqual01(p)] <- 0.5
- q0 <- qgumbel(p0, loc = loc1, scale = scale1, ...)
- q0[.isEqual01(p)] <- sign(p[.isEqual01(p)]-0.5)*Inf
- return(q0)
- }
- body(.Object at q) <- substitute({
- ## P.R.: changed to vectorized form
- p0 <- p
- p0[.isEqual01(p)] <- 0.5
- q0 <- qgumbel(p0, loc = loc1,
- scale = scale1, ...)
- q0[.isEqual01(p)] <- sign(p[.isEqual01(p)]-0.5)*Inf
- return(q0)
- },
- list(loc1 = loc, scale1 = scale))
+ p0[ii01] <- 0.5
+
+ q1 <- qgumbel(p0, loc = loc1, scale = scale1, ...)
+ 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
Added: branches/distr-2.1/pkg/distrEx/R/AsymTotalVarDist.R
===================================================================
--- branches/distr-2.1/pkg/distrEx/R/AsymTotalVarDist.R (rev 0)
+++ branches/distr-2.1/pkg/distrEx/R/AsymTotalVarDist.R 2009-03-17 16:40:29 UTC (rev 411)
@@ -0,0 +1,354 @@
+###############################################################################
+## Method: AsymAsymTotalVarDist
+## asymmetric total variation distance of two distributions
+##
+## newly introduced by P.R. 13 03 09 --- for use in collaboration with
+## Matthias Brandl, Claudio Agostinelli
+##
+###############################################################################
+setMethod("AsymTotalVarDist", signature(e1 = "AbscontDistribution",
+ e2 = "AbscontDistribution"),
+ function(e1, e2, rho = 1,
+ rel.tol = .Machine$double.eps^0.3, maxiter=1000, Ngrid = 10000,
+ TruncQuantile = getdistrOption("TruncQuantile"),
+ IQR.fac = 15){
+
+ ## if we have to recall this method with a smaller TruncQuantile arg:
+ mc <- as.list(match.call(call = sys.call(sys.parent(1)))[-1])
+ mc$TruncQuantile <- TruncQuantile * 1.8
+
+ #block warnings:
+ o.warn <- getOption("warn"); options(warn = -1)
+ on.exit(options(warn=o.warn))
+
+ ## find sensible lower and upper bounds for integration of q-cp
+ # (a) quantile based
+ low <- min(getLow(e1, eps = TruncQuantile), getLow(e2, eps = TruncQuantile))
+ up <- max(getUp(e1, eps = TruncQuantile), getUp(e2, eps = TruncQuantile))
+ # (b) scale based
+ s0 <- min(IQR(e1),IQR(e2))*IQR.fac
+ low0 <- min(median(e1),median(e2))-s0
+ up0 <- max(median(e1),median(e2))+s0
+ # (a) & (b)
+ low <- max(low,low0); up <- min(up,up0)
+ #
+ # store densities to avoid dispatch
+ d1 <- d(e1); d2 <- d(e2)
+ #
+ ### integration as a function of c:
+ Eip <- function(f, c00)
+ distrExIntegrate(f, lower = low, upper = up,
+ rel.tol = rel.tol, c00 = c00)
+ # positive part
+ integ.p <- function(x,c00) pmax(d2(x)-c00*d1(x),0)
+ # negative part
+ integ.m <- function(x,c00) pmax(c00*d1(x)-d2(x),0)
+
+ ## function zero c(rho) of which is to be found
+ fct <- function(c0, rho = rho){
+ e.p <- Eip(f=integ.p, c00=c0)
+ e.m <- Eip(f=integ.m, c00=c0)
+ e.p*rho - e.m
+ }
+
+ ## find sensible search range for c-values
+ ## goal: range of density quotient d2(x)/d1(x)
+ ## x-range:
+ x.range <- seq(low, up, length=Ngrid/3)
+ x.range <- c(x.range, q(e1)(seq(TruncQuantile,1-TruncQuantile,length=Ngrid/3)))
+ x.range <- c(x.range, q(e2)(seq(TruncQuantile,1-TruncQuantile,length=Ngrid/3)))
+ ## to avoid division by 0:
+ d1x.range <- d10x.range <- d1(x.range)
+ d1x.range <- d1x.range+(d1x.range<1e-20)
+ ## bound dx.range from 0 and from maximal value:
+ d2x.range <- d2(x.range)
+ dx.range <- (d2x.range/d1x.range*(d10x.range>=1e-10)+
+ d2x.range*1e10*(d10x.range<1e-10))*.9999999+1e-10
+ ## gives range for c:
+ low1 <- min(dx.range); up1 <- max(dx.range)
+
+ ## in any case compare with c=1
+ tef <- fct(1, rho = rho)
+ ### if c=1 is already a zero:
+ if(tef == 0){
+ res <- Eip(f=integ.p,c00=1)
+ names(res) <- "asym. total variation distance"
+ return(res)
+ }
+ # else: only have to search in c in [low1;1] resp [1;up1]
+ if(tef < 0) up1 <- 1 else low1 <- 1
+
+ c.rho <- try(uniroot(fct, lower = low1, upper = up1,
+ rho = rho, tol = rel.tol,
+ maxiter = maxiter)$root,
+ silent = TRUE)
+
+ ## if does not give reasonable solution recall function with
+ # smaller TruncQuantile
+ if(!is.numeric(c.rho))
+ return(do.call(getMethod("AsymTotalVarDist",
+ signature(e1 = "AbscontDistribution",
+ e2 = "AbscontDistribution")),
+ args = mc))
+ ## else:
+ res <- Eip(f=integ.p, c00=c.rho)
+ names(res) <- "asym. total variation distance"
+ return(res)
+ })
+
+setMethod("AsymTotalVarDist", signature(e1 = "DiscreteDistribution",
+ e2 = "DiscreteDistribution"),
+ function(e1, e2, rho = 1, ...){
+ o.warn <- getOption("warn"); options(warn = -1)
+ on.exit(options(warn=o.warn))
+ supp <- union(support(e1), support(e2))
+ # store densities to avoid dispatch
+ d1 <- d(e1); d2 <- d(e2)
+
+ d2.range <- d2(supp)
+ d1.range <- d1(supp)
+
+ integ.p <- function(c00) pmax(d2.range-c00*d1.range,0)
+ integ.m <- function(c00) pmax(c00*d1.range-d2.range,0)
+
+ fct <- function(c0, rho = rho){
+ e.p <- sum(integ.p(c0))
+ e.m <- sum(integ.m(c0))
+ e.p*rho-e.m
+ }
+
+ d10.range <- d1.range
+ d1e.range <- d1.range+(d1.range<1e-10)
+ ## bound dx.range from 0 and from maximal value:
+ d.range <- (d2.range/d1e.range*(d10.range>=1e-10)+
+ d2.range*1e10*(d10.range<1e-10))*.99999999+1e-10
+ ## gives range for c:
+ low1 <- min(d.range); up1 <- max(d.range)
+
+ ## in any case compare with c=1
+ tef <- fct(1, rho = rho)
+ ### if c=1 is already a zero:
+ if(tef == 0){
+ res <- sum(integ.p(1))
+ names(res) <- "asym. total variation distance"
+ return(res)
+ }
+ # else: only have to search in c in [low1;1] resp [1;up1]
+ if(tef < 0) up1 <- 1 else low1 <- 1
+
+ c.rho <- uniroot(fct, lower=low1, upper=up1,rho=rho)$root
+ res <- sum(integ.p(c.rho))
+ names(res) <- "asym. total variation distance"
+ return(res)
+ })
+setMethod("AsymTotalVarDist", signature(e1 = "DiscreteDistribution",
+ e2 = "AbscontDistribution"),
+ function(e1, e2, rho = 1, ...){
+ res <- 1
+ names(res) <- "asym. total variation distance"
+
+ return(res)
+ })
+setMethod("AsymTotalVarDist", signature(e1 = "AbscontDistribution",
+ e2 = "DiscreteDistribution"),
+ function(e1, e2, rho = 1, ...){
+ res <- 1
+ names(res) <- "asym. total variation distance"
+
+ return(res)
+ })
+setMethod("AsymTotalVarDist", signature(e1 = "numeric",
+ e2 = "DiscreteDistribution"),
+ function(e1, e2, rho = 1, ...){
+ t1 <- table(e1)
+ d1 <- t1/length(e1)
+ s1 <- as.numeric(names(t1))
+ e11 <- DiscreteDistribution(supp=s1, prob=d1)
+ return(AsymTotalVarDist(e11,e2, rho = rho, ...))
+ })
+
+setMethod("AsymTotalVarDist", signature(e1 = "DiscreteDistribution",
+ e2 = "numeric"),
+ function(e1, e2, rho = 1, ...){
+ return(AsymTotalVarDist(e2, e1, rho= 1 / rho, ...))
+ })
+
+## to avoid trivial distances (distance = 1)
+## abs.cont. distributions may be discretized
+## resp. empirical distributions may be smoothed
+## (by convolution with a normal distribution)
+setMethod("AsymTotalVarDist", signature(e1 = "numeric",
+ e2 = "AbscontDistribution"),
+ function(e1, e2, rho = 1, asis.smooth.discretize = "discretize", n.discr =
+ getdistrExOption("nDiscretize"), low.discr = getLow(e2),
+ up.discr = getUp(e2), h.smooth = getdistrExOption("hSmooth"),
+ rel.tol = .Machine$double.eps^0.3, maxiter=1000, Ngrid = 10000,
+ TruncQuantile = getdistrOption("TruncQuantile"),
+ IQR.fac = 15){
+ .asis.smooth.discretize.distance(e1, e2, asis.smooth.discretize, n.discr,
+ low.discr, up.discr, h.smooth, AsymTotalVarDist, rho = rho,
+ rel.tol = rel.tol, maxiter = maxiter, Ngrid = Ngrid,
+ TruncQuantile = TruncQuantile, IQR.fac = IQR.fac)
+ })
+setMethod("AsymTotalVarDist", signature(e1 = "AbscontDistribution",
+ e2 = "numeric"),
+ function(e1, e2, rho = 1, asis.smooth.discretize = "discretize", n.discr =
+ getdistrExOption("nDiscretize"), low.discr = getLow(e1),
+ up.discr = getUp(e1), h.smooth = getdistrExOption("hSmooth"),
+ rel.tol = .Machine$double.eps^0.3, maxiter=1000, Ngrid = 10000,
+ TruncQuantile = getdistrOption("TruncQuantile"),
+ IQR.fac = 15){
+ return(AsymTotalVarDist(e2, e1, rho= 1/rho,
+ asis.smooth.discretize = asis.smooth.discretize,
+ low.discr = low.discr, up.discr = up.discr, h.smooth = h.smooth,
+ rel.tol = rel.tol, maxiter = maxiter, Ngrid = Ngrid,
+ TruncQuantile = TruncQuantile, IQR.fac = IQR.fac))
+ })
+
+setMethod("AsymTotalVarDist", signature(e1 = "AcDcLcDistribution",
+ e2 = "AcDcLcDistribution"),
+ function(e1, e2, rho = 1,
+ rel.tol = .Machine$double.eps^0.3, maxiter=1000, Ngrid = 10000,
+ TruncQuantile = getdistrOption("TruncQuantile"),
+ IQR.fac = 15){
+ ## if we have to recall this method with a smaller TruncQuantile arg:
+ mc <- as.list(match.call(call = sys.call(sys.parent(1)))[-1])
+ mc$TruncQuantile <- TruncQuantile * 1.8
+
+ #block warnings:
+ o.warn <- getOption("warn"); options(warn = -1)
+ on.exit(options(warn=o.warn))
+
+ if( is(e1,"AbscontDistribution"))
+ e1 <- as(as(e1,"AbscontDistribution"), "UnivarLebDecDistribution")
+ if( is(e2,"AbscontDistribution"))
+ e2 <- as(as(e2,"AbscontDistribution"), "UnivarLebDecDistribution")
+ if(is(e1,"DiscreteDistribution"))
+ e1 <- as(as(e1,"DiscreteDistribution"), "UnivarLebDecDistribution")
+ if(is(e2,"DiscreteDistribution"))
+ e2 <- as(as(e2,"DiscreteDistribution"), "UnivarLebDecDistribution")
+
+ ac1 <- acPart(e1); ac2 <- acPart(e2)
+ ac1d <- ac1 at d; ac2d <- ac2 at d
+ ac1.d <- function(x) ac1d(x)*acWeight(e1)
+ ac2.d <- function(x) ac2d(x)*acWeight(e2)
+
+ dc1 <- discretePart(e1); dc2 <- discretePart(e2)
+ dc1d <- dc1 at d; dc2d <- dc2 at d
+ dc1.d <- function(x) dc1d(x)*discreteWeight(e1)
+ dc2.d <- function(x) dc2d(x)*discreteWeight(e2)
+
+ ### continuous part
+
+ ## find sensible lower and upper bounds for integration of q-cp
+ # (a) quantile based
+ low <- min(getLow(ac1, eps = TruncQuantile), getLow(ac2, eps = TruncQuantile))
+ up <- max(getUp(ac1, eps = TruncQuantile), getUp(ac2, eps = TruncQuantile))
+ # (b) scale based
+ s0 <- min(IQR(ac1),IQR(ac2))*IQR.fac
+ low0 <- min(median(ac1),median(ac2))-s0
+ up0 <- max(median(ac1),median(ac2))+s0
+ # (a) & (b)
+ low <- max(low,low0); up <- min(up,up0)
+ #
+ ### integration as a function of c:
+ Eip <- function(f, c00)
+ distrExIntegrate(f, lower = low, upper = up,
+ rel.tol = rel.tol, c00 = c00)
+ # positive part
+ integ.p.c <- function(x,c00) pmax(ac2.d(x)-c00*ac1.d(x),0)
+ # negative part
+ integ.m.c <- function(x,c00) pmax(c00*ac1.d(x)-ac2.d(x),0)
+
+ ### discrete part
+
+ supp <- union(support(dc1), support(dc2))
+
+ d2.range <- dc2.d(supp)
+ d1.range <- dc1.d(supp)
+
+ integ.p.d <- function(c00) pmax(d2.range-c00*d1.range,0)
+ integ.m.d <- function(c00) pmax(c00*d1.range-d2.range,0)
+
+
+ ## function zero c(rho) of which is to be found
+ fct <- function(c0, rho = rho){
+ e.p.c <- Eip(f=integ.p.c, c00=c0)
+ e.m.c <- Eip(f=integ.m.c, c00=c0)
+ e.p.d <- sum(integ.p.d(c0))
+ e.m.d <- sum(integ.m.d(c0))
+ e.p <- e.p.c + e.p.d
+ e.m <- e.m.c + e.m.d
+ return(e.p*rho - e.m)
+ }
+
+ ## find sensible search range for c-values
+ ## goal: range of density quotient d2(x)/d1(x)
+
+ ### continuous part
+ ## x-range:
+ x.range <- seq(low, up, length=Ngrid/3)
+ x.range <- c(x.range, q(ac1)(seq(TruncQuantile,1-TruncQuantile,length=Ngrid/3)))
+ x.range <- c(x.range, q(ac2)(seq(TruncQuantile,1-TruncQuantile,length=Ngrid/3)))
+ ## to avoid division by 0:
+ d1x.range <- d10x.range <- ac1.d(x.range)
+ d1x.range <- d1x.range+(d1x.range<1e-20)
+ ## bound dx.range from 0 and from maximal value:
+ d2x.range <- ac2.d(x.range)
+ dx.range <- (d2x.range/d1x.range*(d10x.range>=1e-20)+
+ d2x.range*1e20*(d10x.range<1e-20))*.99+1e-5
+
+ ### discrete part
+ d10.range <- d1.range
+ d1e.range <- d1.range+(d1.range<1e-7)
+ ## bound dx.range from 0 and from maximal value:
+ d.range <- (d2.range/d1e.range*(d10.range>=1e-10)+
+ d2.range*1e10*(d10.range<1e-10))*.999999+1e-10
+
+ ## gives range for c:
+ low1 <- min(d.range,dx.range); up1 <- max(d.range,dx.range)
+ ## in any case compare with c=1
+ tef <- fct(1, rho = rho)
+ ### if c=1 is already a zero:
+ if(tef == 0){
+ res <- Eip(f=integ.p.c,c00=1)+sum(integ.p.d(1))
+ names(res) <- "asym. total variation distance"
+ return(res)
+ }
+ # else: only have to search in c in [low1;1] resp [1;up1]
+ if(tef < 0) up1 <- 1 else low1 <- 1
+
+ c.rho <- try(uniroot(fct, lower = low1, upper = up1,
+ rho = rho, tol = rel.tol,
+ maxiter = maxiter)$root,
+ silent = TRUE)
+
+ ## if does not give reasonable solution recall function with
+ # smaller TruncQuantile
+ if(!is.numeric(c.rho)){
+ return(do.call(getMethod("AsymTotalVarDist",
+ signature(e1 = "AcDcLcDistribution",
+ e2 = "AcDcLcDistribution")),
+ args = mc))
+ }
+ res <- Eip(f=integ.p.c, c00=c.rho)+sum(integ.p.d(c.rho))
+ names(res) <- "asym. total variation distance"
+ return(res)
+ })
+
+setMethod("AsymTotalVarDist", signature(e1 = "LatticeDistribution",
+ e2 = "LatticeDistribution"),
+ getMethod("AsymTotalVarDist", signature(e1 = "DiscreteDistribution",
+ e2 = "DiscreteDistribution")))
+
+setMethod("AsymTotalVarDist", signature(e1 = "LatticeDistribution",
+ e2 = "DiscreteDistribution"),
+ getMethod("AsymTotalVarDist", signature(e1 = "DiscreteDistribution",
+ e2 = "DiscreteDistribution")))
+
+setMethod("AsymTotalVarDist", signature(e1 = "DiscreteDistribution",
+ e2 = "LatticeDistribution"),
+ getMethod("AsymTotalVarDist", signature(e1 = "DiscreteDistribution",
+ e2 = "DiscreteDistribution")))
+
+
Modified: branches/distr-2.1/pkg/distrEx/R/HellingerDist.R
===================================================================
--- branches/distr-2.1/pkg/distrEx/R/HellingerDist.R 2009-03-17 15:04:15 UTC (rev 410)
+++ branches/distr-2.1/pkg/distrEx/R/HellingerDist.R 2009-03-17 16:40:29 UTC (rev 411)
@@ -4,27 +4,32 @@
###############################################################################
setMethod("HellingerDist", signature(e1 = "AbscontDistribution",
e2 = "AbscontDistribution"),
- function(e1, e2){
- TruncQuantile <- getdistrOption("TruncQuantile")
- lower1 <- getLow(e1)
- upper1 <- getUp(e1)
- lower2 <- getLow(e2)
- upper2 <- getUp(e2)
- lower <- min(lower1, lower2)
- upper <- max(upper1, upper2)
+ function(e1, e2, rel.tol=.Machine$double.eps^0.3,
+ TruncQuantile = getdistrOption("TruncQuantile"),
+ IQR.fac = 15, ...){
+ ## find sensible lower and upper bounds for integration
+ # (a) quantile based
+ low <- min(getLow(e1, eps = TruncQuantile), getLow(e2, eps = TruncQuantile))
+ up <- max(getUp(e1, eps = TruncQuantile), getUp(e2, eps = TruncQuantile))
+ # (b) scale based
+ s0 <- min(IQR(e1),IQR(e2))*IQR.fac
+ low0 <- min(median(e1),median(e2))-s0
+ up0 <- max(median(e1),median(e2))+s0
+ # (a) & (b)
+ lower <- max(low,low0); upper <- min(up,up0)
o.warn <- getOption("warn"); options(warn = -1)
on.exit(options(warn=o.warn))
integrand <- function(x, dfun1, dfun2){ 0.5*(sqrt(dfun1(x))-sqrt(dfun2(x)))^2 }
res <- distrExIntegrate(integrand, lower = lower, upper = upper,
- dfun1 = d(e1), dfun2 = d(e2), rel.tol=.Machine$double.eps^0.3)
+ dfun1 = d(e1), dfun2 = d(e2), rel.tol = rel.tol)
names(res) <- "Hellinger distance"
return(sqrt(res)) # ^.5 added P.R. 19-12-06
})
setMethod("HellingerDist", signature(e1 = "DiscreteDistribution",
e2 = "DiscreteDistribution"),
- function(e1, e2){
+ function(e1, e2, ...){
o.warn <- getOption("warn"); options(warn = -1)
on.exit(options(warn=o.warn))
supp <- union(support(e1), support(e2))
@@ -36,7 +41,7 @@
})
setMethod("HellingerDist", signature(e1 = "DiscreteDistribution",
e2 = "AbscontDistribution"),
- function(e1, e2){
+ function(e1, e2, ...){
res <- 1
names(res) <- "Hellinger distance"
@@ -44,7 +49,7 @@
})
setMethod("HellingerDist", signature(e1 = "AbscontDistribution",
e2 = "DiscreteDistribution"),
- function(e1, e2){
+ function(e1, e2, ...){
res <- 1
names(res) <- "Hellinger distance"
@@ -53,7 +58,7 @@
## Hellinger distance
setMethod("HellingerDist", signature(e1 = "numeric",
e2 = "DiscreteDistribution"),
- function(e1, e2){
+ function(e1, e2, ... ){
d1 <- table(e1)/length(e1)
d2 <- d(e2)(sort(unique(e1)))
e21 <- setdiff(support(e2), unique(e1))
@@ -65,7 +70,7 @@
})
setMethod("HellingerDist", signature(e1 = "DiscreteDistribution",
e2 = "numeric"),
- function(e1, e2){
+ function(e1, e2, ...){
return(HellingerDist(e2, e1))
})
@@ -77,23 +82,35 @@
e2 = "AbscontDistribution"),
function(e1, e2, asis.smooth.discretize = "discretize", n.discr =
getdistrExOption("nDiscretize"), low.discr = getLow(e2),
- up.discr = getUp(e2), h.smooth = getdistrExOption("hSmooth")){
+ up.discr = getUp(e2), h.smooth = getdistrExOption("hSmooth"),
+ rel.tol=.Machine$double.eps^0.3,
+ TruncQuantile = getdistrOption("TruncQuantile"),
+ IQR.fac = 15, ...){
.asis.smooth.discretize.distance(e1, e2, asis.smooth.discretize, n.discr,
- low.discr, up.discr, h.smooth, HellingerDist)
+ low.discr, up.discr, h.smooth, HellingerDist,
+ rel.tol = rel.tol, TruncQuantile = TruncQuantile,
+ IQR.fac = IQR.fac, ...)
})
setMethod("HellingerDist", signature(e1 = "AbscontDistribution",
e2 = "numeric"),
function(e1, e2, asis.smooth.discretize = "discretize", n.discr =
getdistrExOption("nDiscretize"), low.discr = getLow(e1),
- up.discr = getUp(e1), h.smooth = getdistrExOption("hSmooth")){
+ up.discr = getUp(e1), h.smooth = getdistrExOption("hSmooth"),
+ rel.tol=.Machine$double.eps^0.3,
+ TruncQuantile = getdistrOption("TruncQuantile"),
+ IQR.fac = 15, ...){
return(HellingerDist(e2, e1, asis.smooth.discretize = asis.smooth.discretize,
- low.discr = low.discr, up.discr = up.discr, h.smooth = h.smooth))
+ low.discr = low.discr, up.discr = up.discr, h.smooth = h.smooth,
+ rel.tol = rel.tol, TruncQuantile = TruncQuantile,
+ IQR.fac = IQR.fac, ...))
})
#### new from version 2.0 on: Distance for Mixing Distributions
setMethod("HellingerDist", signature(e1 = "AcDcLcDistribution",
e2 = "AcDcLcDistribution"),
- function(e1, e2){
+ function(e1, e2, rel.tol=.Machine$double.eps^0.3,
+ TruncQuantile = getdistrOption("TruncQuantile"),
+ IQR.fac = 15, ...){
if( is(e1,"AbscontDistribution"))
e1 <- as(as(e1,"AbscontDistribution"), "UnivarLebDecDistribution")
if( is(e2,"AbscontDistribution"))
@@ -110,7 +127,8 @@
dc1d <- dc1 at d; dc2d <- dc2 at d
dc1 at d <- function(x) dc1d(x)*discreteWeight(e1)
dc2 at d <- function(x) dc2d(x)*discreteWeight(e2)
- da2 <- HellingerDist(ac1,ac2)^2
+ da2 <- HellingerDist(ac1,ac2, rel.tol = rel.tol,
+ TruncQuantile = TruncQuantile, IQR.fac = IQR.fac, ...)^2
dd2 <- HellingerDist(dc1,dc2)^2
res <- (da2+dd2)^.5
names(res) <- "Hellinger distance"
Modified: branches/distr-2.1/pkg/distrEx/R/Internalfunctions.R
===================================================================
--- branches/distr-2.1/pkg/distrEx/R/Internalfunctions.R 2009-03-17 15:04:15 UTC (rev 410)
+++ branches/distr-2.1/pkg/distrEx/R/Internalfunctions.R 2009-03-17 16:40:29 UTC (rev 411)
@@ -67,15 +67,15 @@
## modify distributions to avoid trivial distances
.asis.smooth.discretize.distance <- function(x, Distribution, asis.smooth.discretize,
- n.discr, low.discr, up.discr, h.smooth, distance){
+ n.discr, low.discr, up.discr, h.smooth, distance, ...){
ASD <- pmatch(asis.smooth.discretize, c("asis", "smooth", "discretize"), nomatch = 3)
- if (ASD == 1) return(distance(x, Distribution))
+ if (ASD == 1) return(distance(x, Distribution, ...))
if (ASD == 2){
Dx <- .smoothDistr(.empiricalDistribution(x), h = h.smooth)
- return(distance(Dx, Distribution))
+ return(distance(Dx, Distribution, ...))
}
if (ASD == 3){
DD <- .discretizeDistr(D = Distribution, x = x, n = n.discr, lower = low.discr, upper = up.discr)
- return(distance(x, DD))
+ return(distance(x, DD, ...))
}
}
Modified: branches/distr-2.1/pkg/distrEx/R/TotalVarDist.R
===================================================================
--- branches/distr-2.1/pkg/distrEx/R/TotalVarDist.R 2009-03-17 15:04:15 UTC (rev 410)
+++ branches/distr-2.1/pkg/distrEx/R/TotalVarDist.R 2009-03-17 16:40:29 UTC (rev 411)
@@ -4,35 +4,32 @@
###############################################################################
setMethod("TotalVarDist", signature(e1 = "AbscontDistribution",
e2 = "AbscontDistribution"),
- function(e1, e2){
- TruncQuantile <- getdistrOption("TruncQuantile")
- lower1 <- ifelse(!is.finite(q(e1)(0)), q(e1)(TruncQuantile), q(e1)(0))
- upper1 <- ifelse(!is.finite(q(e1)(1)),
- ifelse("lower.tail" %in% names(formals(e1 at q)),
- q(e1)(TruncQuantile, lower.tail = FALSE),
- q(e1)(1-TruncQuantile)),
- q(e1)(1))
- lower2 <- ifelse(!is.finite(q(e2)(0)), q(e2)(TruncQuantile), q(e2)(0))
- upper2 <- ifelse(!is.finite(q(e2)(1)),
- ifelse("lower.tail" %in% names(formals(e2 at q)),
- q(e2)(TruncQuantile, lower.tail = FALSE),
- q(e2)(1-TruncQuantile)),
- q(e2)(1))
- lower <- min(lower1, lower2)
- upper <- max(upper1, upper2)
+ function(e1, e2, rel.tol=.Machine$double.eps^0.3,
+ TruncQuantile = getdistrOption("TruncQuantile"),
+ IQR.fac = 15, ...){
+ ## find sensible lower and upper bounds for integration
+ # (a) quantile based
+ low <- min(getLow(e1, eps = TruncQuantile), getLow(e2, eps = TruncQuantile))
+ up <- max(getUp(e1, eps = TruncQuantile), getUp(e2, eps = TruncQuantile))
+ # (b) scale based
+ s0 <- min(IQR(e1),IQR(e2))*IQR.fac
+ low0 <- min(median(e1),median(e2))-s0
+ up0 <- max(median(e1),median(e2))+s0
+ # (a) & (b)
+ low <- max(low,low0); up <- min(up,up0)
o.warn <- getOption("warn"); options(warn = -1)
on.exit(options(warn=o.warn))
integrand <- function(x, dfun1, dfun2){ 0.5*abs(dfun1(x)-dfun2(x)) }
- res <- distrExIntegrate(integrand, lower = lower, upper = upper,
- dfun1 = d(e1), dfun2 = d(e2), rel.tol=.Machine$double.eps^0.3)
+ res <- distrExIntegrate(integrand, lower = low, upper = up,
+ dfun1 = d(e1), dfun2 = d(e2), rel.tol = rel.tol)
names(res) <- "total variation distance"
return(res)
})
setMethod("TotalVarDist", signature(e1 = "DiscreteDistribution",
e2 = "DiscreteDistribution"),
- function(e1, e2){
+ function(e1, e2, ...){
o.warn <- getOption("warn"); options(warn = -1)
on.exit(options(warn=o.warn))
supp <- union(support(e1), support(e2))
@@ -44,7 +41,7 @@
})
setMethod("TotalVarDist", signature(e1 = "DiscreteDistribution",
e2 = "AbscontDistribution"),
- function(e1, e2){
+ function(e1, e2, ...){
res <- 1
names(res) <- "total variation distance"
@@ -52,7 +49,7 @@
})
setMethod("TotalVarDist", signature(e1 = "AbscontDistribution",
e2 = "DiscreteDistribution"),
- function(e1, e2){
+ function(e1, e2, ...){
res <- 1
names(res) <- "total variation distance"
@@ -60,19 +57,16 @@
})
setMethod("TotalVarDist", signature(e1 = "numeric",
e2 = "DiscreteDistribution"),
- function(e1, e2){
- d1 <- table(e1)/length(e1)
- d2 <- d(e2)(sort(unique(e1)))
- e21 <- setdiff(support(e2), unique(e1))
- d21 <- d(e2)(e21)
- res <- 1/2*(sum(abs(d2-d1))+sum(d21))
- names(res) <- "Total variation distance"
-
- return(res)
+ function(e1, e2, ...){
+ t1 <- table(e1)
+ d1 <- t1/length(e1)
+ s1 <- as.numeric(names(t1))
+ e11 <- DiscreteDistribution(supp=s1, prob=d1)
+ return(TotalVarDist(e11,e2))
})
setMethod("TotalVarDist", signature(e1 = "DiscreteDistribution",
e2 = "numeric"),
- function(e1, e2){
+ function(e1, e2, ...){
return(TotalVarDist(e2, e1))
})
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/distr -r 411
More information about the Distr-commits
mailing list