R/TotalVarDist.R

In distrEx: Extensions of Package 'distr'

###############################################################################
## Method: TotalVarDist
## total variation distance of two distributions
###############################################################################
setMethod("TotalVarDist", signature(e1 = "AbscontDistribution", 
                                    e2 = "AbscontDistribution"),
    function(e1, e2, rel.tol=.Machine$double.eps^0.3, 
             TruncQuantile = getdistrOption("TruncQuantile"), 
             IQR.fac = 15, ..., diagnostic = FALSE){
        ## 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 = low, upper = up, 
                    dfun1 = d(e1), dfun2 = d(e2), rel.tol = rel.tol, ...,
                    diagnostic = diagnostic)
        names(res) <- "total variation distance"
        if(diagnostic){
           diagn <- attr(res,"diagnostic")
           diagn[["call"]] <- match.call()
           class(diagn)<- "DiagnosticClass"
           attr(res,"diagnostic") <- diagn
        }

        return(res)
    })
setMethod("TotalVarDist", signature(e1 = "DiscreteDistribution",
                                    e2 = "DiscreteDistribution"),
    function(e1, e2, ...){
        o.warn <- getOption("warn"); options(warn = -1)
        on.exit(options(warn=o.warn))
        supp <- union(support(e1), support(e2))

        res <- 0.5*sum(abs(d(e1)(supp)-d(e2)(supp)))
        names(res) <- "total variation distance"

        return(res)
    })

## new PR 08-09-16
setMethod("TotalVarDist", signature(e1 = "DiscreteMVDistribution",
                                    e2 = "DiscreteMVDistribution"),
    function(e1, e2, ...){
        o.warn <- getOption("warn"); options(warn = -1)
        on.exit(options(warn=o.warn))
        ## replace univariate line  supp <- union(support(e1), support(e2))   by

        supp <- unique(rbind(support(e1), support(e2)))

        res <- 0.5*sum(abs(d(e1)(supp)-d(e2)(supp)))
        names(res) <- "total variation distance"

        return(res)
    })

setMethod("TotalVarDist", signature(e1 = "DiscreteDistribution",
                                    e2 = "AbscontDistribution"),
    function(e1, e2, ...){
        res <- 1
        names(res) <- "total variation distance"

        return(res)
    })
setMethod("TotalVarDist", signature(e1 = "AbscontDistribution",
                                    e2 = "DiscreteDistribution"),
    function(e1, e2, ...){ 
        res <- 1
        names(res) <- "total variation distance"

        return(res)
    })
setMethod("TotalVarDist", signature(e1 = "numeric",
                                    e2 = "DiscreteDistribution"),
    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, ...){
        return(TotalVarDist(e2, e1))
    })

## 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("TotalVarDist", signature(e1 = "numeric",
                                    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"),
             rel.tol=.Machine$double.eps^0.3, 
             TruncQuantile = getdistrOption("TruncQuantile"), 
             IQR.fac = 15, ..., diagnostic = FALSE){
     res <- .asis.smooth.discretize.distance(e1, e2, asis.smooth.discretize, n.discr,
                 low.discr, up.discr, h.smooth, TotalVarDist, 
                 rel.tol = rel.tol, TruncQuantile = TruncQuantile, 
                 IQR.fac = IQR.fac, ..., diagnostic = diagnostic)
     if(diagnostic){
           diagn <- attr(res,"diagnostic")
           diagn[["call"]] <- match.call()
           class(diagn)<- "DiagnosticClass"
           attr(res,"diagnostic") <- diagn
     }
     return(res)
     })
setMethod("TotalVarDist", 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"),
             rel.tol=.Machine$double.eps^0.3, 
             TruncQuantile = getdistrOption("TruncQuantile"), 
             IQR.fac = 15, ..., diagnostic = FALSE){
        return(TotalVarDist(e2, e1, asis.smooth.discretize = asis.smooth.discretize, 
                  low.discr = low.discr, up.discr = up.discr, h.smooth = h.smooth,
                  rel.tol=rel.tol, 
                  TruncQuantile = TruncQuantile, 
                 IQR.fac = IQR.fac, ..., diagnostic = diagnostic ))
    })
#### new from version 2.0 on: Distance for Mixing Distributions
setMethod("TotalVarDist",  signature(e1 = "AcDcLcDistribution",
                                     e2 = "AcDcLcDistribution"),
           function(e1, e2, rel.tol=.Machine$double.eps^0.3, 
                        TruncQuantile = getdistrOption("TruncQuantile"), 
                        IQR.fac = 15, ..., diagnostic = FALSE){
           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@d; ac2d <- ac2@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@d; dc2d <- dc2@d
              dc1@d <- function(x) dc1d(x)*discreteWeight(e1)
              dc2@d <- function(x) dc2d(x)*discreteWeight(e2)
              res <- TotalVarDist(ac1,ac2, rel.tol = rel.tol, 
                        TruncQuantile = TruncQuantile, IQR.fac = IQR.fac, ...) + 
                     TotalVarDist(dc1,dc2, ..., diagnostic = diagnostic)
              names(res) <- "total variation distance"
              if(diagnostic){
                 diagn <- attr(res,"diagnostic")
                 diagn[["call"]] <- match.call()
                 class(diagn)<- "DiagnosticClass"
                 attr(res,"diagnostic") <- diagn
              }
              res
              })

setMethod("TotalVarDist", signature(e1 = "LatticeDistribution", 
                                         e2 = "LatticeDistribution"),
    getMethod("TotalVarDist", signature(e1 = "DiscreteDistribution", 
                                         e2 = "DiscreteDistribution")))

setMethod("TotalVarDist", signature(e1 = "LatticeDistribution", 
                                         e2 = "DiscreteDistribution"),
    getMethod("TotalVarDist", signature(e1 = "DiscreteDistribution", 
                                         e2 = "DiscreteDistribution")))

setMethod("TotalVarDist", signature(e1 = "DiscreteDistribution", 
                                         e2 = "LatticeDistribution"),
    getMethod("TotalVarDist", signature(e1 = "DiscreteDistribution", 
                                         e2 = "DiscreteDistribution")))