Nothing
###############################################################################
# Methods for Discrete Distributions
###############################################################################
## (c) Matthias Kohl: revised P.R. 030707
DiscreteDistribution <- function(supp, prob, .withArith = FALSE,
.withSim = FALSE, .lowerExact = TRUE, .logExact = FALSE,
.DistrCollapse =
getdistrOption("DistrCollapse"),
.DistrCollapse.Unique.Warn =
getdistrOption("DistrCollapse.Unique.Warn"),
.DistrResolution = getdistrOption("DistrResolution"),
Symmetry = NoSymmetry()){
if(!is.numeric(supp))
stop("'supp' is no numeric vector")
if(any(!is.finite(supp))) # admit +/- Inf?
stop("infinite or missing values in supp")
len <- length(supp)
if(missing(prob)){
prob <- rep(1/len, len)
}else{
if(len != length(prob))
stop("'supp' and 'prob' must have equal lengths")
if(any(!is.finite(prob)))
stop("infinite or missing values in prob")
if(!identical(all.equal(sum(prob), 1,
tolerance = 2*getdistrOption("TruncQuantile")), TRUE))
stop("sum of 'prob' has to be (approximately) 1")
if(!all(prob >= 0))
stop("'prob' contains values < 0")
}
o <- order(supp)
supp <- supp[o]
prob <- prob[o]
rm(o)
if(.DistrCollapse){
if (len>1 && min(diff(supp))< .DistrResolution){
erg <- .DistrCollapse(supp, prob, .DistrResolution)
if (len>length(erg$prob) && .DistrCollapse.Unique.Warn)
warning("collapsing to unique support values")
prob <- erg$prob
supp <- erg$supp
}
}else{
usupp <- unique(supp)
if(length(usupp) < len){
if(.DistrCollapse.Unique.Warn)
warning("collapsing to unique support values")
prob <- as.vector(tapply(prob, supp, sum))
supp <- sort(usupp)
len <- length(supp)
rm(usupp)
}
if(len > 1){
if(min(diff(supp))< .DistrResolution)
stop("grid too narrow --> change DistrResolution")
}
}
rm(len)
if(length(supp) == 1){
rfun <- function(n){
rep(supp, n)
}
}else{
rfun <- function(n){
sample(x = supp, size = n, replace = TRUE, prob = prob)
}
}
dfun <- .makeDNew(supp, prob, Cont = FALSE)
pfun <- .makePNew(supp, prob, .withSim, Cont = FALSE)
qfun <- .makeQNew(supp, cumsum(prob), rev(cumsum(rev(prob))),
.withSim, min(supp), max(supp), Cont = FALSE)
object <- new("DiscreteDistribution", r = rfun, d = dfun, q = qfun, p=pfun,
support = supp, .withArith = .withArith, .withSim = .withSim,
.lowerExact = .lowerExact, .logExact = .logExact, Symmetry = Symmetry)
}
setMethod("support", "DiscreteDistribution", function(object) object@support)
### left continuous cdf
setMethod("p.l", "DiscreteDistribution", function(object){
if (.inArgs("lower.tail", p(object))){
function(q, lower.tail = TRUE, log.p = FALSE){
px <- p(object)(q, lower.tail = lower.tail)
o.warn <- getOption("warn");
on.exit(options(warn=o.warn))
options(warn = -2)
dx <- d(object)(.setEqual(q, support(object)))
options(warn = o.warn)
px0 <- pmax(px + if(lower.tail) -dx else dx,0)
if (log.p) px0 <- log(px0)
return(px0)
}
}else{
function(q, lower.tail = TRUE, log.p = FALSE){
px <- p(object)(q)
o.warn <- getOption("warn")
on.exit(options(warn=o.warn))
options(warn = -2)
dx <- d(object)(.setEqual(q, support(object)))
options(warn = o.warn)
px0 <- pmax(if(lower.tail) px - dx else 1 - px + dx, 0)
if (log.p) px0 <- log(px0)
return(px0)
}
}
})
### right continuous quantile function
setMethod("q.r", "DiscreteDistribution", function(object){
if (.inArgs("log.p", q.l(object))){
if (.inArgs("lower.tail", q.l(object))){
function(p, lower.tail = TRUE, log.p = FALSE){
s <- support(object)
psx <- p(object)(s, lower.tail = lower.tail,
log.p = log.p)
ps0 <- .setEqual(p, psx)
o.warn <- getOption("warn"); options(warn = -2)
on.exit(options(warn=o.warn))
qx0 <- q.l(object)(ps0, lower.tail = lower.tail,
log.p = log.p)
options(warn = o.warn)
m <- match(ps0, psx)
n.ina.m <- !is.na(m)
if(any(n.ina.m))
{ M.n.ina.m <- m[n.ina.m]
qx0[n.ina.m] <- (support(object))[pmin(M.n.ina.m+1,
length(s))]
}
if(any(is.nan(qx0)))
warning("NaN's produced")
return(qx0)
}
}else{
function(p, lower.tail = TRUE, log.p = FALSE){
s <- support(object)
psx <- p(object)(s, log.p = log.p)
if (lower.tail) p <- 1 - p
ps0 <- .setEqual(p, psx)
o.warn <- getOption("warn"); options(warn = -2)
on.exit(options(warn=o.warn))
qx0 <- q.l(object)(ps0, lower.tail = lower.tail,
log.p = log.p)
options(warn = o.warn)
m <- match(ps0, psx)
n.ina.m <- !is.na(m)
if(any(n.ina.m))
{ M.n.ina.m <- m[n.ina.m]
qx0[n.ina.m] <- (support(object))[pmin(M.n.ina.m+1,
length(s))]
}
if(any(is.nan(qx0)))
warning("NaN's produced")
return(qx0)
}
}
}else{
if (.inArgs("lower.tail", q.l(object))){
function(p, lower.tail = TRUE, log.p = FALSE){
if (log.p) p <- exp(p)
s <- support(object)
psx <- p(object)(s, lower.tail = lower.tail)
ps0 <- .setEqual(p, psx)
o.warn <- getOption("warn"); options(warn = -2)
on.exit(options(warn=o.warn))
qx0 <- q.l(object)(ps0, lower.tail = lower.tail,
log.p = log.p)
options(warn = o.warn)
m <- match(ps0, psx)
n.ina.m <- !is.na(m)
if(any(n.ina.m))
{ M.n.ina.m <- m[n.ina.m]
qx0[n.ina.m] <- (support(object))[pmin(M.n.ina.m+1,
length(s))]
}
if(any(is.nan(qx0)))
warning("NaN's produced")
return(qx0)
}
}else{
function(p, lower.tail = TRUE, log.p = FALSE){
if (log.p) p <- exp(p)
s <- support(object)
psx <- p(object)(s)
if (lower.tail) p <- 1 - p
ps0 <- .setEqual(p, psx)
o.warn <- getOption("warn"); options(warn = -2)
on.exit(options(warn=o.warn))
qx0 <- q.l(object)(ps0, lower.tail = lower.tail,
log.p = log.p)
options(warn = o.warn)
m <- match(ps0, psx)
n.ina.m <- !is.na(m)
if(any(n.ina.m))
{ M.n.ina.m <- m[n.ina.m]
qx0[n.ina.m] <- (support(object))[pmin(M.n.ina.m+1,
length(s))]
}
if(any(is.nan(qx0)))
warning("NaN's produced")
}
}
}
})
## Convolution Discrete Distributions
setMethod("+", c("DiscreteDistribution","DiscreteDistribution"),
function(e1,e2){
if(length(support(e1))==1) return(e2+support(e1))
if(length(support(e2))==1) return(e1+support(e2))
e1.L <- as(e1, "LatticeDistribution")
e2.L <- as(e2, "LatticeDistribution")
if(is(e1.L, "LatticeDistribution") & is(e2.L, "LatticeDistribution"))
{w1 <- width(lattice(e1.L))
w2 <- width(lattice(e2.L))
W <- sort(abs(c(w1,w2)))
if (abs(abs(w1)-abs(w2))<getdistrOption("DistrResolution") ||
W[2] %% W[1] < getdistrOption("DistrResolution") ){
res <- e1.L + e2.L
res@.finSupport <- e1.L@.finSupport&e2.L@.finSupport
return(res)
}
}
res <- .convDiscrDiscr(e1,e2)
res@.finSupport <- e1@.finSupport&e2@.finSupport
return(res)
})
setMethod("+", c("Dirac","DiscreteDistribution"),
function(e1,e2){e2+location(e1)})
## binary operators for discrete distributions
setMethod("*", c("DiscreteDistribution","numeric"),
function(e1, e2) { Distr <- .multm(e1,e2, "DiscreteDistribution")
if(is(Distr, "AffLinDistribution"))
Distr@X0 <- e1
if(is(e1@Symmetry,"SphericalSymmetry"))
Distr@Symmetry <-
SphericalSymmetry(SymmCenter(e1@Symmetry)*e2)
if(is.finite(e2)){
Distr@.finSupport <- e1@.finSupport
}else{
ep <- .Machine$double.eps
Distr@.finSupport <- c(p(e1)(0)<ep,p(e1)(0)>1-ep)
}
if(e2<0) Distr@.finSupport <- rev(Distr@.finSupport)
return(Distr)
})
setMethod("+", c("DiscreteDistribution","numeric"),
function(e1, e2) { Distr <- .plusm(e1,e2, "DiscreteDistribution")
if(is(Distr, "AffLinDistribution"))
Distr@X0 <- e1
if(is(e1@Symmetry,"SphericalSymmetry"))
Distr@Symmetry <-
SphericalSymmetry(SymmCenter(e1@Symmetry)+e2)
isfe2 <- c(e2 >(-Inf), e2<Inf)
Distr@.finSupport <- e1@.finSupport & isfe2
return(Distr)
})
setMethod("*", c("AffLinDiscreteDistribution","numeric"),
function(e1, e2) {
Distr <- .multm(e1,e2, "AffLinDiscreteDistribution")
if(is(e1@Symmetry,"SphericalSymmetry"))
Distr@Symmetry <-
SphericalSymmetry(SymmCenter(e1@Symmetry)*e2)
if(is.finite(e2)){
Distr@.finSupport <- e1@.finSupport
}else{
ep <- .Machine$double.eps
Distr@.finSupport <- c(p(e1)(0)<ep,p(e1)(0)>1-ep)
}
if(e2<0) Distr@.finSupport <- rev(Distr@.finSupport)
return(Distr)
})
setMethod("+", c("AffLinDiscreteDistribution","numeric"),
function(e1, e2) {
Distr <- .plusm(e1,e2, "AffLinDiscreteDistribution")
if(is(e1@Symmetry,"SphericalSymmetry"))
Distr@Symmetry <-
SphericalSymmetry(SymmCenter(e1@Symmetry)*e2)
isfe2 <- c(e2 >(-Inf), e2<Inf)
Distr@.finSupport <- e1@.finSupport & isfe2
return(Distr)
})
## Group Math for discrete distributions
setMethod("Math", "DiscreteDistribution",
function(x){
rnew <- function(n, ...){}
body(rnew) <- substitute({ f(g(n, ...)) },
list(f = as.name(.Generic), g = x@r))
object <- new("DiscreteDistribution", r = rnew,
.withSim = TRUE, .withArith = TRUE)
object@.finSupport <- x@.finSupport&NA
object
})
setMethod("Math", "Dirac",
function(x){ loc <- location(x)
lc <- callGeneric(loc)
object <- Dirac(lc)
object@.finSupport <- x@.finSupport&NA
object
})
## exact: abs for discrete distributions
setMethod("abs", "DiscreteDistribution",function(x){
rnew <- function(n, ...){}
body(rnew) <- substitute({ abs(g(n, ...)) },
list(g = x@r))
xx <- x
supportnew <- support(x)
isSym0 <- FALSE
if(is(Symmetry(x),"SphericalSymmetry"))
if(.isEqual(SymmCenter(Symmetry(x)),0))
isSym0 <- TRUE
if(isSym0){
supportnew <- supportnew[supportnew>=0]
.lowerExact = .lowerExact(x)
dxlog <- if("log" %in% names(formals(d(x))))
quote({dx <- d(xx)(x, log = TRUE)})
else quote({dx <- log(d(xx)(x))})
pxlog <- if("log.p" %in% names(formals(p(x))) &&
"lower.tail" %in% names(formals(p(x))))
quote({p(x)(q, lower.tail = FALSE, log.p = TRUE)})
else
quote({log(1-p(x)(q))})
qxlog <- if("lower.tail" %in% names(formals(q.l(x))))
quote({qx <- if(lower.tail)
q.l(x)((1+p1)/2)
else
q.l(x)(p1/2,lower.tail=FALSE)})
else
quote({qx <- q.l(x)(if(lower.tail) (1+p1)/2 else 1-p1/2)})
if("lower.tail" %in% names(formals(q.l(x)))&&
"log.p" %in% names(formals(q.l(x))))
qxlog <- quote({qx <- if(lower.tail) q.l(x)((1+p1)/2)
else
q.l(x)(if(log.p)p-log(2)
else p1/2,lower.tail=FALSE,log.p=log.p)})
dnew <- function(x, log = FALSE){}
body(dnew) <- substitute({
dxlog0
dx[x>0] <- dx+log(2)
if (!log) dx <- exp(dx)
dx[x<0] <- if(log) -Inf else 0
return(dx)
}, list(dxlog0 = dxlog))
pnew <- function(q, lower.tail = TRUE, log.p = FALSE){}
body(pnew) <- substitute({
if (!lower.tail){
px <- (log(2) + pxlog0)*(q>=0)
if(!log.p) px <- exp(px)
}else{
px <- pmax(2 * p(x)(q) - 1,0)
if(log.p) px <- log(px)
}
return(px)
}, list(pxlog0 = pxlog))
qnew <- function(p, lower.tail = TRUE, log.p = FALSE){}
body(qnew) <- substitute({
p1 <- if(log.p) exp(p) else p
qxlog0
qx[p1<0] <- NaN
if (any((p1 < -.Machine$double.eps)|(p1 > 1+.Machine$double.eps)))
warning(gettextf("q method of %s produced NaN's ", objN))
return(qx)
}, list(qxlog0 = qxlog, objN= quote(.getObjName(1))))
}else{
if (.isEqual(p.l(x)(0),0)) return(x)
supportnew <- sort(unique(abs(supportnew)))
dnew <- function(x, log = FALSE){
o.warn <- getOption("warn"); options(warn = -1)
on.exit(options(warn=o.warn))
dx <- (x>=0) * d(xx)(x) + (x>0) * d(xx)(-x)
options(warn = o.warn)
if (log) dx <- log(dx)
return(dx)
}
pxlow <- if("lower.tail" %in% names(formals(p(x))))
substitute({p(x)(q, lower=FALSE)})
else
substitute({1-p(x)(q)})
pnew <- function(q, lower.tail = TRUE, log.p = FALSE){}
body(pnew) <- substitute({
px <- if (lower.tail)
(q>=0) * (p(x)(q) - p.l(x)(-q))
else pxlow0 + p.l(x)(-q)
if (log.p) px <- log(px)
return(px)
}, list(pxlow0=pxlow))
prob <- dnew(supportnew)
qnew <- .makeQNew(supportnew, cumsum(prob),
rev(cumsum(rev(prob))), notwithLLarg = x@.withSim,
min(supportnew), max(supportnew), Cont = FALSE)
}
object <- new("DiscreteDistribution", r = rnew, p = pnew,
q = qnew, d = dnew, support = supportnew,
.withSim = x@.withSim, .withArith = TRUE,
.lowerExact = .lowerExact(x))
object@.finSupport <- c(TRUE, all(x@.finSupport))
object
})
## exact: eps for discrete distributions
setMethod("exp", "DiscreteDistribution",
function(x){ obj <- .expm.d(x)
obj@.finSupport <- c(TRUE, x@.finSupport[2])
obj
}
)
### preliminary to export special functions
if (getRversion()>='2.6.0'){
setMethod("log", "DiscreteDistribution",
function(x, base = exp(1)) {
xs <- as.character(deparse(match.call(
call = sys.call(sys.parent(1)))$x))
ep <- getdistrOption("TruncQuantile")
basl <- log(base)
if(p(x)(0)>ep)
stop(gettextf("log(%s) is not well-defined with positive probability ", xs))
else{
obj <- .logm.d(x)/basl
obj@.finSupport <- c(TRUE, x@.finSupport[2])
return(obj)
}})
setMethod("log", "Dirac",
function(x, base = exp(1)){
xs <- as.character(deparse(match.call(
call = sys.call(sys.parent(1)))$x))
loc <- location(x)
ep <- getdistrOption("TruncQuantile")
basl <- log(base)
if(loc < ep)
stop(gettextf("log(%s) is not well-defined with positive probability ", xs))
Dirac(log(loc)/basl)})
setMethod("log10", "DiscreteDistribution",
function(x) log(x = x, base = 10))
setMethod("sign", "DiscreteDistribution",
function(x){
d0 <- d(x)(0)
DiscreteDistribution(supp=c(-1,0,1),
prob=c(p(x)(-getdistrOption("TruncQuantile")),
d0,
p(x)(getdistrOption("TruncQuantile"), lower=FALSE)))
})
setMethod("digamma", "DiscreteDistribution",
function(x){
px0 <- p(x)(0)
if(px0>0) stop("argument of 'digamma' must be concentrated on positive values")
rnew <- function(n, ...){}
body(rnew) <- substitute({ digamma(g(n, ...)) }, list(g = x@r))
object <- DiscreteDistribution(
supp=digamma(support(x)),
prob=prob(x), .withArith = TRUE)
object@.finSupport <- c(TRUE, x@.finSupport[2])
object
})
setMethod("lgamma", "DiscreteDistribution",
function(x){
rnew = function(n, ...){}
body(rnew) <- substitute({ lgamma(g(n, ...)) }, list(g = x@r))
object <- new("DiscreteDistribution", r = rnew,
.withSim = TRUE, .withArith = TRUE)
object@.finSupport <- c(TRUE, x@.finSupport[2])
object
})
setMethod("gamma", "DiscreteDistribution",
function(x){
rnew = function(n, ...){}
body(rnew) <- substitute({ gamma(g(n, ...)) }, list(g = x@r))
object <- new("DiscreteDistribution", r = rnew,
.withSim = TRUE, .withArith = TRUE)
object@.finSupport <- c(TRUE, x@.finSupport[2])
object
})
setMethod("sqrt", "DiscreteDistribution",
function(x) x^0.5)
}
setMethod("prob", "DiscreteDistribution",
function(object) {sp <- object@support
pr <- object@d(sp)
names(pr) <- paste(sp)
return(pr)
})
## Replace Methods
setReplaceMethod("prob", "DiscreteDistribution",
function(object, value){
return(DiscreteDistribution(supp = object@support,
prob = value,
.withArith = object@.withArith,
.withSim = object@.withSim,
.lowerExact = .lowerExact(object),
.logExact = .logExact(object)))}
)
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.