R/LatticeDistribution.R

In distr: Object Oriented Implementation of Distributions

################################################################################
#
#                  LatticeDistribution
#
################################################################################

LatticeDistribution <- function(lattice = NULL, supp = NULL, prob = NULL, 
                       .withArith = FALSE, .withSim = FALSE, 
                       DiscreteDistribution = NULL, check = TRUE,
                       Symmetry = NoSymmetry()){
    if(is(lattice,"Lattce")){
       if(width(lattice)>0){
          .finS <- c(TRUE,is.finite(Length(lattice)))
       }else{
          .finS <- c(is.finite(Length(lattice)), TRUE)
       }
    }else .finS <- c(TRUE,TRUE)
    if (is(DiscreteDistribution, "AffLinDiscreteDistribution"))
        {  D <- DiscreteDistribution
           if (is(lattice, "Lattice")) 
             { ### check consistency with support of DiscreteDistribution} 
              if (check){
                 if( !.is.consistent(lattice, support(D), eq.space = FALSE))         
                     stop(paste("Argument 'lattice' is inconsistent to",
                            " the support of argument 'DiscreteDistribution'." , 
                            sep = ""))
              }           
              return(new("AffLinLatticeDistribution", r = D@r, d = D@d,
                          q = D@q, p = D@p, support = D@support, 
                          a = D@a, b = D@b, X0 = D@X0,
                          lattice = lattice, .withArith = .withArith, 
                          .withSim = .withSim, img = D@img,
                          param = D@param, Symmetry = Symmetry,
                          .finSupport = .finS))
              }else{
               if (check){
                   if( !.is.vector.lattice(support(D)))
                       stop(paste("Support of argument 'DiscreteDistribution' ",
                              "is not a lattice.", sep = ""))
               }           
               return(new("AffLinLatticeDistribution", r = D@r, d = D@d, 
                          q = D@q, p = D@p, support = D@support, 
                          lattice = .make.lattice.es.vector(D@support), 
                          a = D@a, b = D@b, X0 = D@X0,
                          .withArith = .withArith, 
                          .withSim = .withSim, img = D@img,
                          param = D@param, Symmetry = Symmetry,
                          .finSupport = .finS))
              }                 
        }

    if (is(DiscreteDistribution, "DiscreteDistribution"))
        {  D <- DiscreteDistribution
           if (is(lattice, "Lattice")) 
             { ### check consistency with support of DiscreteDistribution} 
              if (check){
                  if( !.is.consistent(lattice, support(D), eq.space = FALSE))         
                     stop(paste("Argument 'lattice' is inconsistent to the",
                            " support of argument 'DiscreteDistribution'." , 
                            sep = ""))
              }           
              return(new("LatticeDistribution", r = D@r, d = D@d, 
                          q = D@q, p = D@p, support = D@support, 
                          lattice = lattice, .withArith = .withArith, 
                          .withSim = .withSim, img = D@img,
                          param = D@param, Symmetry = Symmetry,
                          .finSupport = .finS))
              }else{
               if (check){
                   if( !.is.vector.lattice(support(D)))
                     stop(paste("Support of argument 'DiscreteDistribution' is",
                              "not a lattice.", sep = " "))
               }           
 
               return(new("LatticeDistribution", r = D@r, d = D@d, 
                          q = D@q, p = D@p, support = D@support, 
                          lattice = .make.lattice.es.vector(D@support), 
                          .withArith = .withArith, 
                          .withSim = .withSim, img = D@img,
                          param = D@param, Symmetry = Symmetry,
                          .finSupport = .finS))
              }                 
        }

    if (is(lattice, "Lattice") && ! is.null(supp))
       {D <- DiscreteDistribution(supp = supp, prob = prob, 
                                   .withArith = .withArith, 
                                   .withSim = .withSim, Symmetry = Symmetry )
        
        if (check){
            if( !.is.consistent(lattice, supp, eq.space = FALSE))         
                stop("Argument 'lattice' is inconsistent to argument 'supp'.")
        }
        
        return(new("LatticeDistribution", r = r(D), d = d(D), 
                    q = q.l(D), p = p(D), support = supp,
                    lattice = lattice, .withArith = .withArith, 
                    .withSim = .withSim, Symmetry = Symmetry,
                    .finSupport = .finS))
       }

    if (is(lattice, "Lattice"))
       {if (is.finite(Length(lattice)))
             {if (is.null(prob))
                  prob <- rep(1/Length(lattice), Length(lattice))
              if (Length(lattice) == length(prob))
                 {supp <- seq( pivot(lattice), length = Length(lattice), 
                               by = width(lattice))
                  D <- DiscreteDistribution(supp = supp, prob = prob, 
                                             .withArith = .withArith, 
                                             .withSim = .withSim, 
                                             Symmetry = Symmetry )
                  return(new("LatticeDistribution", r = r(D), d = d(D), 
                          q = q.l(D), p = p(D), support = supp,
                          lattice = lattice, .withArith = .withArith, 
                          .withSim = .withSim, Symmetry = Symmetry,
                          .finSupport = .finS))
                  }else{ 
                   #if (check)
                       stop("Lengths of lattice and probabilities differ.")
                   #else return(D)
                   }    
              }else {if (is.null(prob))
                        stop(paste("Insufficient information given to ",
                                   "determine distribution.", sep = ""))
                     else{
                         supp <- seq( pivot(lattice), length = length(prob), 
                                     by = width(lattice))
                         D <- DiscreteDistribution(supp = supp, prob = prob, 
                                                   .withArith = .withArith, 
                                                   .withSim = .withSim, 
                                                   Symmetry = Symmetry )
                         return(new("LatticeDistribution", r = r(D), d = d(D), 
                                q = q.l(D), p = p(D), support = supp,
                                lattice = lattice, .withArith = .withArith, 
                                .withSim = .withSim, Symmetry = Symmetry,
                                .finSupport = .finS))
                        }                  
             }
       }else if (!is.null(supp))
            {if (is.null(prob)) prob <- supp*0+1/length(supp)
             D <- DiscreteDistribution(supp, prob, .withArith = .withArith, 
                                       .withSim = .withSim, Symmetry = Symmetry )
             if (check){
                 if (!.is.vector.lattice (supp))
                     stop("Argument 'supp' given is not a lattice.")
             }    
             return(new("LatticeDistribution", r = D@r, d = D@d, 
                             q = D@q, p = D@p, support = D@support, 
                             lattice = .make.lattice.es.vector(D@support), 
                             .withArith = D@.withArith, 
                             .withSim = D@.withSim, img = D@img,
                             param = D@param, Symmetry = Symmetry,
                            .finSupport = .finS))
            }else 
             stop("Insufficient information given to determine distribution.")
}


setMethod("lattice", "LatticeDistribution", function(object) object@lattice)


## canceling out of lattice points with mass 0
#setAs("LatticeDistribution", "DiscreteDistribution", 
#       def = function(from){
#    cF <- class(from)[1]
#    value <- if (cF!="LatticeDistribution") 
#                 new(cF) else new("DiscreteDistribution")
#    for (what in slotNames("DiscreteDistribution")) 
#         slot(value, what) <- slot(from, what)
#    supp.old <- from@support
#    o.warn <- getOption("warn"); options(warn = -2)
#    d.old <- from@d(from@support)
#    options(warn = o.warn)
#    supp.new <- supp.old[d.old > 0]
#    value@support <- supp.new
#    value
#       }
#)


setAs("AffLinLatticeDistribution","AffLinDiscreteDistribution", 
       def = function(from){
    value <- new("AffLinDiscreteDistribution")
    for (what in slotNames("AffLinDiscreteDistribution")) 
         slot(value, what) <- slot(from, what)
    supp.old <- from@support
    o.warn <- getOption("warn"); options(warn = -2)
    on.exit(options(warn=o.warn))
    d.old <- from@d(from@support)
    supp.new <- supp.old[d.old > 0]
    options(warn = o.warn)
    value@support <- supp.new
    value
       }
)

setMethod("+", c("LatticeDistribution", "LatticeDistribution"),
function(e1,e2){
            if(length(support(e1))==1) return(e2+support(e1))
            if(length(support(e2))==1) return(e1+support(e2))

### Lattice calculations:

            sup1 <- support(e1)
            sup2 <- support(e2)
            # left and right endpoint of convolution support
            su12.l <- sup1[1]+sup2[1]
            su12.r <- (rev(sup1))[1]+(rev(sup2))[1]

            l1 <- length(sup1)
            l2 <- length(sup2)

            lat1 <- lattice(e1)
            lat2 <- lattice(e2)
            L1 <- Length(lat1)
            L2 <- Length(lat2)
            w1 <- width(lat1)
            w2 <- width(lat2)


            ### take care if lattice is infinite
            L.inf <- !(is.finite(L1)&&is.finite(L2))
            if(L.inf){
               if(is.finite(L2)){
                  if(w1>0)
                     L.lr <- +1
                  else
                     L.lr <- -1   
               }else{    
                  if(is.finite(L1)){
                     if(w2>0)
                        L.lr <- +1
                     else
                        L.lr <- -1   
                  }else{
                     if(w1*w2>0) L.lr <- if(w1>0) +1 else -1
                     if(w1*w2<0) L.lr <- if(abs(su12.l)<abs(su12.r)) +1 else -1
                  }
               }
            }   


            e0 <- NULL
            tol0 <- .distroptions$DistrResolution/1000
            
            ## treat case separately when Discr + Discr is "faster"
            if(l1*l2 < 100){ 
                     d0 <- .convDiscrDiscr(e1,e2)
                     sup0 <- support(d0)
                     md <- min(diff(sup0))
                     sup00 <- seq(from=min(sup0),to=max(sup0),by=md)
                     sup0s <- intersect(sup00,sup0)
                     sup01 <- .inWithTol(sup00, sup0s, tol=tol0)
                     sup10 <- .inWithTol(sup0, sup0s, tol=tol0)
                     if(!all(sup10)) return(d0)
                     pr0 <- sup00*0
                     pr0[sup01] <- (prob(d0))[sup10]
                     pr0 <- pr0/sum(pr0)                              
                     lat <- Lattice(pivot = sup00[1], width = md, 
                                    Length = length(sup00))
                     e0 <- LatticeDistribution(supp = sup00, prob = pr0, 
                                               lattice = lat, check = FALSE)           
                     if(L.inf){
                         wa <- .getCommonWidth(abs(w1),abs(w2), tol=tol0)
                         e0@lattice <- if(L.lr>0){
                            Lattice(pivot = su12.l, width = wa, Length = Inf)
                                }else{ 
                            Lattice(pivot = su12.r, width = -wa, Length = Inf)}
                     }       
               }
 
            ## step 1 common width
              wa <- .getCommonWidth(abs(w1),abs(w2),
                      tol=tol0)


            ## treat case separately when no common support, i.e. when 
            ## w1/w2 is not "rational" enough
            
            if(is.null(wa))  return(.convDiscrDiscr(e1,e2))
            
            
            w0 <- ifelse(w1<0,-1,1) * wa
            pi1 <- pivot(lat1)
            pi2 <- pivot(lat2)                        
            ### Step 2
            supp0 <- seq(by = wa, from = min(sup1-pi1, sup2-pi2),
                                  to = max(sup1-pi1, sup2-pi2))
            s1 <- .inWithTol(supp0,sup1-pi1,tol0)
            s2 <- .inWithTol(supp0,sup2-pi2,tol0)
            d1 <- d2 <- 0*supp0
            d1[s1] <- prob(as(e1,"DiscreteDistribution"))
            d2[s2] <- prob(as(e2,"DiscreteDistribution"))

            L <- length(supp0)
            Ln <- 2^(ceiling(log(L)/log(2))+1)

            ### Step 3
            d1 <- c(d1, numeric(Ln-L))
            d2 <- c(d2, numeric(Ln-L))

            ##STEP 4
            ## computation of DFT
            ftde1 <- fft(d1)
            ftde2 <- fft(d2)

            ## convolution theorem for DFTs
            newd <- (Re(fft(ftde1*ftde2, inverse = TRUE)) / Ln)[1:(2*L+1)]
            newd <- (newd >= .Machine$double.eps^1.5)*newd


            ## reduction to relevant support
            supp1 <- seq(by = wa,
                         from = 2 * min(sup1-pi1, sup2-pi2),
                         to   = 2 * max(sup1-pi1, sup2-pi2))+pi1+pi2

            L1 <- length(supp1)
            newd <- newd[1:L1]

            if (L1 > getdistrOption("DefaultNrGridPoints")){
                rsum.u <- min( sum( rev(cumsum(rev(newd))) >=
                                    getdistrOption("TruncQuantile")/2)+1,
                               length(supp1)
                           )
                rsum.l <- 1 + sum( cumsum(newd) <
                                   getdistrOption("TruncQuantile")/2)
            }else{
                rsum.u <- min( sum( rev(cumsum(rev(newd))) >=
                                    .Machine$double.eps),
                               length(supp1)
                           )
                rsum.l <- 1 + sum( cumsum(newd) < .Machine$double.eps)

            }
            wi1 <- rsum.l:rsum.u
            newd <- newd[wi1]
            newd <- newd/sum(newd)
            supp1 <- supp1[wi1]

            wi2 <- newd > getdistrOption("TruncQuantile")
            supp2 <- supp1[wi2]
            newd2 <- newd[wi2]
            newd2 <- newd2/sum(newd2)

            Symmetry <- NoSymmetry()
            if(is(e1@Symmetry,"SphericalSymmetry")&& 
               is(e2@Symmetry,"SphericalSymmetry"))
               Symmetry <- SphericalSymmetry(SymmCenter(e1@Symmetry)+
                                              SymmCenter(e2@Symmetry))   

            if( length(supp1) >= 2 * length(supp2)){
               res <- DiscreteDistribution(supp = supp2, prob = newd2,
                                           .withArith = TRUE, Symmetry = Symmetry)
               res@.finSupport <- e1@.finSupport & e2@.finSupport
               return(res)
            }else{
               lat <- Lattice(pivot=supp1[1],width=wa, Length=length(supp1))

               e0 <- LatticeDistribution(supp = supp1, prob = newd,
                                         lattice = lat,
                                         .withArith = TRUE, Symmetry = Symmetry,
                                         check = FALSE)
               if(L.inf){
                  e0@lattice <- if(L.lr>0){ 
                            Lattice(pivot = su12.l, width = wa, Length = Inf) 
                         }else{ 
                            Lattice(pivot = su12.r, width = -wa, Length = Inf)}
               }
               e0@.finSupport <- e1@.finSupport & e2@.finSupport
               return(e0)
            }
          })

## extra methods
## binary operators

setMethod("+", c("LatticeDistribution", "numeric"),
          function(e1, e2) 
             {L <- lattice(e1)
              pivot(L) <- pivot(L) + e2
              Distr <- as(e1, "DiscreteDistribution") + e2 
              if(is(Distr, "AffLinDistribution"))
                    Distr@X0 <- e1

              Symmetry <- NoSymmetry()
              if(is(e1@Symmetry,"SphericalSymmetry"))
                 Symmetry <- SphericalSymmetry(SymmCenter(e1@Symmetry)+e2)   
              
              res <- LatticeDistribution(lattice = L,
                     DiscreteDistribution = Distr, Symmetry = Symmetry, 
                     check = FALSE)
              res@.finSupport <- e1@.finSupport & c(e2>(-Inf),(e2<Inf))
              return(res)
              })       

setMethod("*", c("LatticeDistribution", "numeric"),
          function(e1, e2) 
             {if (.isEqual(e2,0))
                  return(Dirac( location = 0 ))
              else{
                  L <- lattice(e1)
                  pivot(L) <- pivot(L) * e2
                  width(L) <- width(L) * e2
                  Distr <- as(e1, "DiscreteDistribution") * e2 
                  if(is(Distr, "AffLinDistribution"))
                        Distr@X0 <- e1
            
                  Symmetry <- NoSymmetry()
                  if(is(e1@Symmetry,"SphericalSymmetry"))
                     Symmetry <- SphericalSymmetry(SymmCenter(e1@Symmetry) * e2)   
              
                  res <- LatticeDistribution(lattice = L,
                          DiscreteDistribution = Distr, Symmetry = Symmetry, 
                          check = FALSE)
                  if(is.finite(e2)){
                      res@.finSupport <- e1@.finSupport
                  }else{
                      ep <- .Machine$double.eps
                      res@.finSupport <- c((p(e1)(0)<=  ep),(p(e1)(0)>=1-ep))
                  }
                  if(e2<0) res@.finSupport <- rev(res@.finSupport)
                  return(res)
                }
             }
          )              

setMethod("+", c("AffLinLatticeDistribution", "numeric"),
          function(e1, e2) 
             {L <- lattice(e1)
              pivot(L) <- pivot(L) + e2
              Symmetry <- NoSymmetry()
              if(is(e1@Symmetry,"SphericalSymmetry"))
                 Symmetry <- SphericalSymmetry(SymmCenter(e1@Symmetry) + e2)   
              res <- LatticeDistribution(lattice = L,
                     DiscreteDistribution = 
                        as(e1, "AffLinDiscreteDistribution") + e2,
                        Symmetry = Symmetry, 
                     check = FALSE)                     
              res@.finSupport <- e1@.finSupport & c(e2>(-Inf),(e2<Inf))
              return(res)
              })

setMethod("*", c("AffLinLatticeDistribution", "numeric"),
          function(e1, e2) 
             {if (isTRUE(all.equal(e2,0)))
                  return(Dirac( location = 0 ))
              else     
                { L <- lattice(e1)
                  pivot(L) <- pivot(L) * e2
                  width(L) <- width(L) * e2
                  Symmetry <- NoSymmetry()
                  if(is(e1@Symmetry,"SphericalSymmetry"))
                     Symmetry <- SphericalSymmetry(SymmCenter(e1@Symmetry) * e2)   
                  res <- LatticeDistribution(lattice = L,
                          DiscreteDistribution = 
                             as(e1, "AffLinDiscreteDistribution") * 
                             e2, Symmetry = Symmetry, 
                             check = FALSE)
                  if(is.finite(e2)){
                      res@.finSupport <- e1@.finSupport
                  }else{
                      ep <- .Machine$double.eps
                      res@.finSupport <- c((p(e1)(0)<=  ep),(p(e1)(0)>=1-ep))
                  }
                  if(e2<0) res@.finSupport <- rev(res@.finSupport)
                  return(res)
                }
             }
          )