Nothing
#------------------------------------------------------------------------------
### internal help function to check whether numeric vectors are integers
#------------------------------------------------------------------------------
.isInteger <- function(x, tol = .Machine$double.eps) abs(as.integer(x)-x)< tol
.isNatural <- function(x, tol = .Machine$double.eps) .isInteger(x, tol) & (x>0)
.isNatural0 <- function(x, tol = .Machine$double.eps) .isInteger(x, tol) & (x>=0)
setAs("numeric","Integer",function(from) new("Integer",as.integer(from)))
#------------------------------------------------------------------------------
### internal help function to check if a vector can be made a lattice
#------------------------------------------------------------------------------
.is.vector.lattice <- function(x)
{ ### is x equally spaced?
all( sapply(diff(x), function(y)
isTRUE(all.equal(y, diff(x)[1],
tolerance = getdistrOption("DistrResolution"),
check.attributes = FALSE))))
}
### internal help function to check consistency of lattice with support:
.is.consistent <- function(lattice, support, eq.space = TRUE)
{
p <- pivot(lattice); l <- Length(lattice); w <- width(lattice)
ds <- diff(support); ms <- min(support); Ms <- max(support)
### is support equally spaced?
if (! .is.vector.lattice(support) && eq.space)
return(FALSE)
### are width of lattice and support consistent
if (! isTRUE(all.equal(min(ds), abs(w), check.attributes = FALSE)))
return(FALSE)
### pivot is left or right endpoint of support
if ( isTRUE(all.equal(ms, p, check.attributes = FALSE)) || isTRUE(all.equal(Ms, p, check.attributes = FALSE)) )
return(TRUE)
if (isTRUE(all.equal(min((support[1]-p)%%w,w-(support[1]-p)%%w),0,
tolerance = getdistrOption("TruncQuantile"),
check.attributes = FALSE)))
return(TRUE)
return(FALSE)
}
## generate a lattice from equally spaced vector
.make.lattice.es.vector <- function(x){
new("Lattice", pivot = x[1], width = x[2]-x[1],
Length = length(x))
}
#------------------------------------------------------------------------------
### internal help function to determin common lattice width
#------------------------------------------------------------------------------
.EuclidAlgo <- function(n1,n2){
r<- 2
m <- round(max(n1,n2))
n <- round(min(n1,n2))
if (n==1) return(1)
while (r>1){
r <- m %% n
m <- n
if(r==1) return(1)
if(r==0) break
n <- r
}
return(n)
}
.getCommonWidth <- function(x1,x2, tol=.Machine$double.eps){
rI <- function(x) .isInteger(x, tol = x * tol)
if(rI(x1)&&rI(x2)) return(.EuclidAlgo(x1,x2))
if(rI(10000*x1) && rI(10000*x2)) return(.EuclidAlgo(10000*x1,10000*x2)/10000)
n <- max(x1,x2); m <- min(x1,x2)
if(rI(n/m)) return(m)
vc <- 1:10000
vecT <- sapply(vc,function(x) rI(x*(n/m)))
vec <- m/vc
# print(min(vc[vecT]))
if(any(vecT)) return((vec/min(vc[vecT]))[1])
return(NULL)
}
#------------------------------------------------------------------------------
### %in% for numerics with tolerance
#------------------------------------------------------------------------------
.inWithTol <- function(x,y,tol=.Machine$double.eps){
sapply(x, function(u) any(abs(u-y)<tol))
}
#------------------------------------------------------------------------------
### brute force convolution
#------------------------------------------------------------------------------
.convDiscrDiscr <- function(e1,e2){
convolutedsupport <- rep(support(e1), each = length(support(e2))) +
support(e2)
gridvalues1 <- d(e1)(support(e1)); gridvalues2 <- d(e2)(support(e2))
convolutedvalues <- rep(gridvalues1, each = length(support(e2))) *
gridvalues2
rm(gridvalues1,gridvalues2)
tmptable <- data.frame(x = convolutedsupport, dx = convolutedvalues)
rm(convolutedsupport,convolutedvalues)
tmp <- tapply(tmptable$dx, tmptable$x, sum)
rm(tmptable)
supp.u <- as.numeric(names(tmp))
prob.u <- as.numeric(tmp)
o <- order(supp.u)
supp <- supp.u[o]
prob <- prob.u[o]
#supp.u <- unique(supp)
len <- length(supp)
if(len > 1){
if (min(diff(supp))< getdistrOption("DistrResolution")){
if (getdistrOption("DistrCollapse")){
erg <- .DistrCollapse(supp, prob,
getdistrOption("DistrResolution"))
if ( len > length(erg$prob) &&
getdistrOption("DistrCollapse.Unique.Warn") )
warning("collapsing to unique support values")
prob <- erg$prob
supp <- erg$supp
}else
stop("grid too narrow --> change DistrResolution")
}
}
rm(tmp, len)
.withSim <- e1@.withSim || e2@.withSim
rfun <- function(n) {}
body(rfun) <- substitute({ f(n) + g(n) },
list(f = e1@r, g = e2@r))
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, p = pfun,
q = qfun, support = supp,
.withSim = .withSim, .withArith = TRUE)
rm(rfun, dfun, qfun, pfun)
if(is(e1@Symmetry,"SphericalSymmetry")&&
is(e2@Symmetry,"SphericalSymmetry"))
object@Symmetry <- SphericalSymmetry(SymmCenter(e1@Symmetry)+
SymmCenter(e2@Symmetry))
object@.finSupport <- e1@.finSupport & e2@.finSupport
object
}
#------------------------------------------------------------------------------
### .fm, .fM, .fM2 functions
#------------------------------------------------------------------------------
.fM <- function(x,f){
owarn <- getOption("warn")
options("warn"=-1)
on.exit(options("warn"=owarn))
if(!.inArgs("log.p", f))
f0 <- f
else f0 <- function(x) f(log(x), log.p = TRUE)
xo <- x
x1 <- (1+xo)/2
i <- 1
while( i < 30)
{ i <- i+1
while(!is.na(f1 <- f0(x1)) && f1 < Inf)
{xo <- x1
x1 <- (1+x1)/2
}
x1 <- (x1+xo)/2}
log(xo)}
.fM2 <- function(x,f){
owarn <- getOption("warn")
options("warn"=-1)
on.exit(options("warn"=owarn))
if(!.inArgs("log.p", f))
f0 <- function(x) f(x, lower.tail = FALSE)
else f0 <- function(x) f(log(x), lower.tail = FALSE, log.p = TRUE)
xo <- x
x1 <- xo/2
i <- 1
while( i < 30)
{ i <- i+1
while(!is.na(f1 <- f0(x1)) && f1 < Inf)
{xo <- x1
x1 <- x1/2
}
x1 <- (x1+xo)/2}
log(xo)}
.fm <- function(x,f){
owarn <- getOption("warn")
options("warn"=-1)
on.exit(options("warn"=owarn))
if(!.inArgs("log.p", f))
f0 <- f
else f0 <- function(x) f(log(x), log.p = TRUE)
xo <- x
x1 <- xo/2
i <- 1
while( i < 30)
{ i <- i+1
while(!is.na(f1 <- f0(x1)) && f1 > -Inf)
{xo <- x1
x1 <- x1/2
}
x1 <- (x1+xo)/2}
log(xo)}
#------------------------------------------------------------------------------
# .presubs : for titles etc
#------------------------------------------------------------------------------
.presubs <- function(inp, frompat, topat){
### replaces in an expression or a string all frompat patterns to topat patterns
logic <- FALSE
inCx <- sapply(inp,
function(inpx){
inC <- deparse(inpx)
l <- length(frompat)
for(i in 1:l)
{ if (is.language(topat[[i]])){
totxt <- deparse(topat[[i]])
totxt <- gsub("expression\\(", "\", ", gsub("\\)$",", \"",totxt))
if (length(grep(frompat[i],inC))) logic <<- TRUE
inC <- gsub(frompat[i],totxt,inC)
}else inC <- gsub(frompat[i], topat[[i]], inC)
}
return(inC)
})
if(length(grep("expression",inCx))>0)
inCx <- gsub("expression\\(", "", gsub("\\)$","",inCx))
if (length(inCx) > 1) {
inCx <- paste(inCx, c(rep(",", length(inCx)-1), ""),
sep = "", collapse = "\"\\n\",")
if ( any(as.logical(c(lapply(inp,is.language)))) | logic )
inCx <- paste("expression(paste(", gsub("\\\\n"," ", inCx), "))", sep ="")
else
inCx <- paste("paste(",inCx,")", sep ="")
}else inCx <- paste("expression(paste(",inCx,"))",sep="")
outC <- eval(parse(text = eval(inCx)))
return(outC)
}
#------------------------------------------------------------------------------
# help check functions
#------------------------------------------------------------------------------
.inArgs <- function(arg, fct)
{as.character(arg) %in% names(formals(fct))}
.isEqual <- function(p0, p1, tol = min( getdistrOption("TruncQuantile")/2,
.Machine$double.eps^.7
))
abs(p0-p1)< tol
.isIn <- function(p0, pmat, tol = min( getdistrOption("TruncQuantile")/2,
.Machine$double.eps^.7
))
{## PR 2018 04 13
## detected by Tuomo.OJALA@3ds.com: the gaps matrix can
## have zero rows -> check this in the following line
if(nrow(pmat)==0) return(FALSE)
list1 <- lapply(1:nrow(pmat), function(x){
(p0+tol > pmat[x,1]) & (p0-tol < pmat[x,2]) })
apply(matrix(unlist(list1), ncol = nrow(pmat)), 1, any)}
.isEqual01<- function(x) .isEqual(x,0)|.isEqual(x,1)
.setEqual <- function(x, y, tol = 1e-7){
### all elements of x equal to some element of y up tol are set to exactly
### the respective element of y
x1 <- round(2*x/tol,0)
y1 <- round(2*y/tol,0)
z <- x
m <- match(x1,y1)
n.ina.m <- !is.na(m)
z[n.ina.m] <- y[m[n.ina.m]]
z
}
.notwithLArg <- function(D) D@.withSim||!.inArgs("lower.tail",p(D))
#------------------------------------------------------------------------------
# other helpers
#------------------------------------------------------------------------------
.getObjName <- function(i = 1)
{Ca <- sys.call(-4); as.character(as.list(Ca[[1]]))[i+1]}
.discretizeP <- function(D, lower, upper, h){
h0 <- 40*(getUp(D)-getLow(D)) /
2^getdistrOption("DefaultNrFFTGridPointsExponent")
if(h > h0 )
warning(paste("Grid for approxfun too wide, ",
"increase DefaultNrFFTGridPointsExponent", sep =""))
x <- seq(from = lower, to = upper, by = h)
if(TRUE){#.notwithLArg(D)){
return(diff(p(D)(x)))
# return((diff(p(D)(x))+diff(rev(p(D)(x,lower=FALSE))))/2)
}else{
M <- q.l(D)(0.5); L <- length(x)
x.l <- x [ x <= M ]; x.u <- x [ x >= M ]
L.l <- length(x.l); L.u <- length(x.u)
if(L.u+L.l-L>0)
return(c(diff(p(D)(x.l, lower.tail = TRUE)),
-diff(p(D)(x.u, lower.tail = FALSE))))
else
return(c(diff(p(D)(x.l, lower.tail = TRUE)),
p(D)(x.u[1])-p(D)(x.l[L.l]),
-diff(p(D)(x.u, lower.tail = FALSE))))
}
}
#------------------------------------------------------------------------------
# .makeD, .makeP, .makeQ
#------------------------------------------------------------------------------
.makeD <- function(object, argList, stand = NULL, fac = NULL)
{ d <- function(x, log = FALSE, ...){}
if(inA <- .inArgs("log", object@d))
argList <- substitute(c(AL, log = log, ...),
list(AL = argList))
else argList <- substitute(c(AL, ...),
list(AL = argList))
myCall <- substitute(d0 <- do.call("@"(objC,d),AL) ,
list(objC = quote(object),
AL = argList))
myBod0 <- myCall
if (!inA) myBod0 <- substitute({myCall
if (log) d0 <- log(d0)})
myBod1 <- myBod0
if (!is.null(stand))
myBod1 <- substitute(
{myBod0
d0 <- if (log) d0 - log(stand) else d0 / stand },
list(stand = stand, myBod0 = myBod0))
myBod <- myBod1
if (!is.null(fac))
myBod <- substitute(
{myBod1
d0 <- if (log) d0 + log(fac) else d0 * fac },
list(fac = fac, myBod1 = myBod1))
body(d) <- substitute({myBod
return(d0)})
return(d)
}
.makeP <- function(object, argList, sign = TRUE, correct = NULL,
fac = NULL, fac2 = NULL)
{
p <- function(q, lower.tail = TRUE, log.p = FALSE, ...){}
siY <- substitute(lower.tail)
siN <- substitute(!lower.tail)
lowS <- if( sign) siY else siN
lowS1 <- if(!sign) siY else siN
if(inA1 <- .inArgs("lower.tail", object@p))
argList <- substitute(c(aL, lower.tail = lowT),
list( lowT = lowS,
aL = argList )
)
else argList <- substitute(aL, list(aL = argList ))
if(inA2 <- (.inArgs("log.p", object@p) && inA1))
argList <- substitute(c(argList, log.p = log.p, ...))
else argList <- substitute(c(argList, ...))
myCall <- substitute(p0 <- facC * do.call("@"(objC,p),AL) + facD,
list(objC = quote(object),
AL = argList,
facC = if(is.null(fac)) 1 else fac,
facD = if(is.null(fac2)) 0 else fac2)
)
if (!inA1)
myBod0 <- substitute({myC
if (lowT) p0 <- 1 - p0},
list(myC=myCall, lowT = lowS1))
else
myBod0 <- substitute(myCall)
if (!sign && !is.null(correct))
myBod1 <- substitute({myC
corC},
list(myC = myBod0, corC = correct))
else myBod1 <- substitute(myBod0)
if (!inA2) myBod <- substitute({myBod1
if (log.p) p0 <- log(p0)
return(p0)} )
else myBod <- substitute({myBod1
return(p0)})
body(p) <- substitute(myBod)
return(p)
}
.makeQ <- function(object, lastCall, sign = TRUE, Cont = TRUE)
## lastCall of form e1 %op% e2
{
siY <- substitute(lower.tail)
siN <- substitute(!lower.tail)
lowS <- if( sign) siY else siN
lowS1 <- if(!sign) siY else siN
q <- function(p, lower.tail = TRUE, log.p = FALSE, ...){}
argList <- alist(p)
if(inA1 <- .inArgs("lower.tail", object@q))
argList1 <- substitute(alist(lower.tail = lowT),
list( lowT = lowS))
else argList1 <- NULL
if(inA2 <- (.inArgs("log.p", object@q) && inA1))
argList1 <- substitute(c(argList1,
alist(log.p = log.p), ...))
else argList1 <- substitute(c(argList1, alist(...)))
argList <- substitute(c(AL,AL1),list(AL=argList,AL1=argList1))
myCall <- substitute({q0 <- do.call("@"(objC,q),AL)
q0 <- lasC
return(q0)},
list(objC = quote(object), AL=argList,
lasC = lastCall))
if (!Cont && !sign){
if (inA1){
dqfS <- substitute({
dq <- getdistrOption("DistrResolution")
if (lower.tail) dq <- -dq
})
}else{
dqfS <- substitute({
dq <- getdistrOption("DistrResolution")
})
}
if (inA2){ indS <- substitute(if(log.p){ ind <- (p>-Inf)&(p<0)
p0[!ind] <- +1 }
else { ind <- (p>0)&(p<1)
p0[!ind] <- -1 } )
}else { indS <- substitute({ ind <- (p>0)&(p<1)
p0[!ind] <- -1 })
}
myBod1 <- substitute({
q01 <- do.call("@"(object,q),AL)
p0 <- do.call("@"(object,p), c(alist(q=q01),AL1))
indC
ind2 <- .isEqual(p,p0)
dqfC # dq * if (with lower.tail-arg & lower.tail == FALSE)
# -1 else 1
p[ind2] <- p[ind2] + dq
my0C
}, list(my0C = substitute(myCall),
AL = substitute(argList), AL1 = substitute(argList1),
indC = substitute(indS), dqfC = substitute(dqfS) )
)
## discrete correction
}else{
myBod1 <- substitute(myCall)
}
if (!inA1) myBod2 <- substitute({if (lowT) p <- 1 - p
myC},
list(lowT = lowS1, myC = myBod1 ))
else myBod2 <- substitute(myBod1)
if (!inA2){
myBod <- substitute({if (log.p) p <- exp(p)
myC },
list(myC = myBod2))
}else{ myBod <- substitute(myBod2) }
body(q) <- substitute(myBod)
return(q)
}
#------------------------------------------------------------------------------
# .plusm, .multm
#------------------------------------------------------------------------------
.plusm <- function(e1, e2, Dclass = "DiscreteDistribution"){
if (length(e2)>1) stop("length of operator must be 1")
if (.isEqual(e2, 0)) return(e1)
if(Dclass %in% c("AffLinDiscreteDistribution",
"AffLinAbscontDistribution"))
if(.isEqual(e1@a,1)&&.isEqual(e1@b+e2,0)) return(e1@X0)
if ((Dclass == "DiscreteDistribution")||
(Dclass == "AffLinDiscreteDistribution"))
supportnew <- e1@support + e2
if ((Dclass == "AbscontDistribution")||
(Dclass == "AffLinAbscontDistribution"))
gapsnew <- if(is.null(e1@gaps)) NULL else e1@gaps + e2
rnew <- function(n, ...){}
body(rnew) <- substitute({ f(n, ...) + g },
list(f = e1@r, g = e2))
dnew <- .makeD(e1, substitute(alist(x = x - e2), list(e2 = e2)))
pnew <- .makeP(e1, substitute(alist(q = q - e2), list(e2 = e2)))
qnew <- .makeQ(e1, substitute(q0 + e2, list(e2 = e2)))
isfe2 <- c((e2>(-Inf)),(e2<Inf))
if (Dclass == "AffLinDiscreteDistribution"){
object <- new("AffLinDiscreteDistribution",
r = rnew, d = dnew, p = pnew,
q = qnew, support = supportnew,
a = e1@a, b = e1@b + e2, X0 = e1@X0,
.withSim = FALSE, .withArith = TRUE,
.logExact = .logExact(e1), .lowerExact = .lowerExact(e1),
.finSupport = e1@.finSupport& isfe2)
rm(supportnew)
}else if (Dclass == "DiscreteDistribution"){
object <- new("AffLinDiscreteDistribution",
r = rnew, d = dnew, p = pnew,
q = qnew, support = supportnew,
a = 1, b = e2, X0 = e1,
.withSim = FALSE, .withArith = TRUE,
.logExact = .logExact(e1), .lowerExact = .lowerExact(e1),
.finSupport = e1@.finSupport&isfe2)
rm(supportnew)
}else if (Dclass == "AffLinAbscontDistribution"){
object <- new("AffLinAbscontDistribution",
r = rnew, d = dnew, p = pnew, q = qnew,
gaps = gapsnew, a = e1@a, b = e1@b + e2,
X0 = e1@X0, .withSim = FALSE, .withArith = TRUE,
.logExact = .logExact(e1), .lowerExact = .lowerExact(e1))
}else if (Dclass == "AbscontDistribution"){
object <- new("AffLinAbscontDistribution", r = rnew,
d = dnew, p = pnew, q = qnew, gaps = gapsnew,
a = 1, b = e2, X0 = e1, .withSim = FALSE,
.withArith = TRUE,
.logExact = .logExact(e1), .lowerExact = .lowerExact(e1))
}
rm(pnew, qnew, dnew, rnew)
object
}
.multm <- function(e1, e2, Dclass = "DiscreteDistribution"){
withg <- getdistrOption("withgaps")
on.exit(distroptions(withgaps=withg))
distroptions(withgaps=FALSE)
if (length(e2)>1) stop("length of operator must be 1")
if (.isEqual(e2, 1)) return(e1)
if (.isEqual(e2, 0)) return(new("Dirac", location = 0))
if(Dclass %in% c("AffLinDiscreteDistribution",
"AffLinAbscontDistribution"))
if(.isEqual(e1@a*e2,1)&&.isEqual(e1@b,0)) return(e1@X0)
rnew <- function(n, ...){}
body(rnew) <- substitute({ f(n, ...) * g },
list(f = e1@r, g = e2))
if (Dclass == "AffLinDiscreteDistribution"){
if(is.finite(e2)){
.finS <- e1@.finSupport
}else{
ep <- .Machine$double.eps
.finS <- c(p(e1)(0)<ep,p(e1)(0)>1-ep)
}
if(e2<0) .finS <- rev(.finS)
supportnew <- e1@support * e2
if (e2 < 0) supportnew <- rev(supportnew)
coR <- substitute({
o.warn <- getOption("warn"); options(warn = -1)
on.exit(options(warn=o.warn))
#
x0 <- .setEqual(q / e2C, support(object))
d0 <- object@d(x = x0)*(x0 %in% support(object))
#
options(warn = o.warn)
if (!lower.tail) d0 <- -d0
p0 <- p0 + d0},
list(e2C = e2)
)
dnew <- .makeD(substitute(e1, list(e1 = e1)),
substitute(alist(x = x / e2), list(e2 = e2)))
pnew <- .makeP(substitute(e1, list(e1 = e1)),
substitute(alist(q = q / e2), list(e2 = e2)),
sign = e2>0, correct = coR)
qnew <- .makeQ(substitute(e1, list(e1 = e1)),
substitute(q0 * e2, list(e2 = e2)),
sign = e2>0, Cont = FALSE)
object <- new(Dclass, r = rnew, d = dnew, p = pnew,
q = qnew, support = supportnew,
a = e1@a * e2, b = e2 * e1@b, X0 = e1@X0,
.withSim = FALSE, .withArith = TRUE,
.logExact = .logExact(e1), .lowerExact = .lowerExact(e1),
.finSupport = .finS)
rm(supportnew)
}else if (Dclass == "DiscreteDistribution"){
if(is.finite(e2)){
.finS <- e1@.finSupport
}else{
ep <- .Machine$double.eps
.finS <- c(p(e1)(0)<ep,p(e1)(0)>1-ep)
}
supportnew <- e1@support * e2
if (e2 < 0) supportnew <- rev(supportnew)
coR <- substitute({
o.warn <- getOption("warn"); options(warn = -1)
on.exit(options(warn=o.warn))
d0 <- object@d(x = q / e2C)
options(warn = o.warn)
if (!lower.tail) d0 <- -d0
p0 <- p0 + d0},
list(e2C = e2)
)
dnew <- .makeD(substitute(e1, list(e1 = e1)),
substitute(alist(x = x / e2), list(e2 = e2)))
pnew <- .makeP(substitute(e1, list(e1 = e1)),
substitute(alist(q = q / e2), list(e2 = e2)),
sign = e2>0, correct = coR)
qnew <- .makeQ(substitute(e1, list(e1 = e1)),
substitute(q0 * e2, list(e2 = e2)),
sign = e2>0, Cont = FALSE)
object <- new("AffLinDiscreteDistribution", r = rnew, d = dnew,
p = pnew, q = qnew, support = supportnew,
a = e2, b = 0, X0 = e1,
.withSim = FALSE, .withArith = TRUE,
.logExact = .logExact(e1), .lowerExact = .lowerExact(e1),
.finSupport = .finS)
rm(supportnew)
}else if (Dclass == "AffLinAbscontDistribution"){
trY <- try(
if(is.null(e1@gaps))
gapsnew <- NULL
else {gapsnew <- e1@gaps
if(nrow(gapsnew)){
if(is.numeric(gapsnew)&&length(gapsnew)>0){
gapsnew <- matrix(gapsnew * e2, ncol=2)
if (e2 < 0) gapsnew <-
gapsnew[rev(seq(nrow(gapsnew))),c(2,1),drop = FALSE]
}}
}, silent=TRUE)
if(is(trY,"try-error")) gapsnew <- NULL
dnew <- .makeD(substitute(e1, list(e1 = e1)),
substitute(alist(x = x / e2), list(e2 = e2)),
stand = abs(e2))
pnew <- .makeP(substitute(e1, list(e1 = e1)),
substitute(alist(q = q / e2), list(e2 = e2)),
sign = e2>0)
qnew <- .makeQ(substitute(e1, list(e1 = e1)),
substitute(q0 * e2, list(e2 = e2)),
sign = e2>0)
object <- new(Dclass, r = rnew, d = dnew, p = pnew, q = qnew,
gaps = gapsnew, a = e1@a * e2, b = e2 * e1@b,
X0 = e1@X0, .withSim = FALSE, .withArith = TRUE,
.logExact = .logExact(e1), .lowerExact = .lowerExact(e1))
}else if (Dclass == "AbscontDistribution"){
trY <- try(
if(is.null(e1@gaps))
gapsnew <- NULL
else {gapsnew <- e1@gaps
if(is.numeric(gapsnew)&&length(gapsnew)>0){
gapsnew <- matrix(gapsnew * e2, ncol=2)
if (e2 < 0) gapsnew <-
gapsnew[rev(seq(nrow(gapsnew))),
c(2,1),drop = FALSE]
}
}, silent=TRUE)
if(is(trY,"try-error")) gapsnew <- NULL
dnew <- .makeD(substitute(e1, list(e1 = e1)),
substitute(alist(x = x / e2), list(e2 = e2)),
stand = abs(e2))
pnew <- .makeP(substitute(e1, list(e1 = e1)),
substitute(alist(q = q / e2), list(e2 = e2)),
sign = e2>0)
qnew <- .makeQ(substitute(e1, list(e1 = e1)),
substitute(q0 * e2, list(e2 = e2)),
sign = e2>0)
object <- new("AffLinAbscontDistribution", r = rnew, d = dnew,
p = pnew, q = qnew, gaps = gapsnew, a = e2,
b = 0, X0 = e1, .withSim = FALSE,
.withArith = TRUE,
.logExact = .logExact(e1), .lowerExact = .lowerExact(e1))
}
rm(pnew, qnew, dnew, rnew)
object
}
#------------------------------------------------------------------------------
# .makeDd, .makePd, .makeQd
#------------------------------------------------------------------------------
.makeDd <- function(x,y, yleft, yright){
intervall <- getdistrOption("DistrResolution") / 2
supp.grid <- c(matrix(rbind(x - intervall, x + intervall), nrow = 1))
prob.grid <- c(matrix(rbind(0, y), nrow = 1), 0)
df0 <- stepfun(x = supp.grid, y = prob.grid)
rm(intervall, supp.grid, prob.grid)
return(df0)
}
.makePd <- function(x,y, yleft, yright){
stepfun(x = x, y = c(yleft, y))
}
.makeQd <- function(x,y, yleft, yright){
force(y)
force(x)
f <- function(u) {
q0 <- sapply(u,
function(z) y[min(sum(x < z-.Machine$double.eps) + 1,
length(y)) ] )
q0[.isEqual(u,0)] <- yleft
q0[.isEqual(u,1)] <- yright
return(q0)}
### individualize local variables
#environment(f) <- new.env()
#assign("y", y, environment(f))
#assign("x", x, environment(f))
#assign("yright", yright, environment(f))
#assign("yleft", yleft, environment(f))
return(f)
}
.makeQc <- function(x,y, yleft, yright){
#f0 <- function(u) {
# q0 <- sapply(u,
# function(z) y[min(sum(x < z-.Machine$double.eps) + 1,
# length(y)) ] )
# return(q0)}
#eps <- getdistrOption("TruncQuantile")
#x00 <- seq(eps, 1-eps, length = getdistrOption("DefaultNrGridPoints"))
#y00 <- f0(x00)
#idx <- cumsum(rle(y00)[[1]]) ### use only unique y's and corresponding
#x0 <- x00[idx] ### maximal x's
#y0 <- y00[idx]
#f1 <- approxfun(x = x0, y = y0, yleft = y0[1], yright = y0[length(y0)])
yleft <- yleft[1]
yright <- yright[1]
isna <- is.na(x)|is.na(y)
x <- x[!isna]
y <- y[!isna]
l0 <- length(unique(x[!.isEqual01(x)]))
if(l0 > 1){
yl <- if(!is.na(yleft) && is.finite(yleft)) yleft else y[1]
yr <- if(!is.na(yright)&& is.finite(yright)) yright else y[length(y)]
f1 <- approxfun(x = x, y = y, yleft = yl, yright = yr)
}else{
i0 <- (1:length(x))[x==unique(x[!.isEqual01(x)])]
y0 <- if(l0 ==1) y[min(i0)] else yleft
f1 <- function(x) return(y0)
}
f <- function(x)
{y1 <- f1(x)
y1[.isEqual(x,0)] <- yleft
y1[.isEqual(x,1)] <- yright
return(y1)
}
return(f)
}
#------------------------------------------------------------------------------
# .makeDNew, .makePNew, .makeQNew
#------------------------------------------------------------------------------
.makeDNew <- function(x, dx, h = NULL, Cont = TRUE, standM = "sum"){
dx <- (dx >= .Machine$double.eps)*dx
if( length(dx) < length(x) ) dx <- c(0,dx)
if (is.null(h)) h <- 1
dx1 <- dx / h
mfun <- if (Cont) approxfun else .makeDd
## density
df1 <- mfun(x = x, y = dx1, yleft = 0, yright = 0)
if (standM == "sum")
stand <- sum(dx)
else {
stand <- try(integrate(df1, -Inf, Inf)$value, TRUE)
if (is(stand,"try-error")){
if(getdistrOption("warn.makeDNew"))
warning("'integrate()' threw an error ---result may be inaccurate.")
stand <- sum(df1(x))*h*(x[2]-x[1])
}
}
dfun <- function(x, log = FALSE)
{if (log)
d0 <- log(df1(x))-log(stand)
else d0 <- df1(x) / stand
return (d0)}
rm(x,dx1,h)
return(dfun)
}
.primefun <- function(f,x, nm = NULL){
h <- diff(x)
l <- length(x)
### 20240127
if(l < 2) stop("Too few grid points for .primefun.")
xm <- (x[-l]+x[-1])/2
fxm <- f(xm)
fx <- f(x)
fxs <- 2 * cumsum(fx) - fx - fx[1]
fxsm <- 4 * cumsum(fxm)
fxx <- c(0, (fxs[-1]+fxsm)* h / 6 )
if (is.null(nm)) nm <- fxx[l]
fx1 <- approxfun(x, fxx, yright = nm, yleft = 0)
ffx <- function(u){
ffy <- fx1(u)
ffy[u > max(x)] <- nm
ffy[u < min(x)] <- 0
return(ffy)
}
return(ffx)
}
.csimpsum <- function(fx){
l <- length(fx)
l2 <- l%/%2
if (l%%2 == 0) {
fx <- c(fx[1:l2],(fx[l2]+fx[l2+1])/2,fx[(l2+1):l])
l <- l+1}
f.even <- fx[seq(l) %% 2 == 0]
f.odd <- fx[seq(l) %% 2 == 1]
fs <- 2 * cumsum(f.odd) - f.odd - f.odd[1]
fsm <- 4 * cumsum(f.even)
ff <- c(0,(fs[2:(l2+1)]+fsm)/3 )
ff
}
.makePNew <- function(x, dx, h = NULL, notwithLLarg = FALSE,
Cont = TRUE, myPf = NULL, pxl = NULL, pxu = NULL){
if (is.null (h)) h <- 0
x.u <- x.l <- x
if (Cont){
mfun <- if (is.null (myPf)) approxfun else myPf
l <- length(x)
### 20240127
#if (l < 2) stop("Too few grid points for .makePNew.")
if ((l%%2==0)&& is.null(myPf)){
l2 <- l/2
if (is.null(pxl))
x.l <- c(x[1:l2],(x[l2]+x[l2+1])/2,x[(l2+1):l])
if (is.null(pxu))
x.u <- c(x[1:l2],(x[l2]+x[l2+1])/2,x[(l2+1):l])
l <- l+1
}
cfun <- .csimpsum
if (is.null(pxl)&& is.null(myPf))
x.l <- x.l[seq(l)%%2==1]
if (is.null(pxu)&& is.null(myPf))
x.u <- x.u[seq(l)%%2==1]
}else {
mfun <- .makePd
cfun <- cumsum
}
p.l <- if(!is.null(pxl)) pxl else cfun(dx)
nm <- max(p.l)
p1.l <- mfun(x = x.l, y = p.l, yleft = 0, yright = nm)
nm <- p1.l(max(x))
if(notwithLLarg){
ifElsePS <- substitute(if (lower.tail) p1.l(q) else 1 - p1.l(q))
}else{
p.u <- if(!is.null(pxu)) pxu else rev(cfun(rev(dx)))
## continuity correction by h/2
if (!Cont) p.u <- c(p.u[-1],0)
p1.u <- mfun(x = x.u, y = p.u, yright = 0, yleft = nm)
rm(p.u)
ifElsePS <- substitute(if (lower.tail) p1.l(q) else p1.u(q))
}
pfun <- function(q, lower.tail = TRUE, log.p = FALSE){}
body(pfun) <- substitute(
{ p0 <- ifElsePC
p0 <- if (log.p) log(p0)-log(nm) else p0/nm
return(p0)
}, list(ifElsePC = ifElsePS))
rm(dx, p.l, notwithLLarg)
return(pfun)
}
.makeQNew <- function(x, px.l, px.u, notwithLLarg = FALSE, yL , yR,
Cont = TRUE){
o.warn <- getOption("warn"); options(warn = -1)
on.exit(options(warn=o.warn))
ix <- .isEqual01(px.l)
mfun <- if (Cont) .makeQc else .makeQd
if(!is.finite(yR)||Cont)
{xx <- px.l[!ix]; yy <- x[!ix]}
else
{xx <- px.l; yy <- x}
q.l0 <- mfun(x = xx, y = yy, yleft = yL, yright = yR)
rm(xx,yy)
if(notwithLLarg){
ifElseQS <- quote(if (lower.tail) q.l0(p01) else q.l0(1-p01))
}else{
# px.u <- rev(px.u);
x <- rev(x)
if (Cont) px.u <- rev(px.u)
ix <- .isEqual01(px.u)
xx <- px.u[!ix]
yy <- if (Cont) x[!ix] else x[rev(!ix)]
q.u <- mfun(x = xx, y = yy, yleft = yR, yright = yL)
rm(xx,yy)
ifElseQS <- quote(if (lower.tail) q.l0(p01) else q.u(p01))
}
options(warn = o.warn)
qfun <- function(p, lower.tail = TRUE, log.p = FALSE){}
body(qfun) <- substitute({
if (log.p) p <- exp(p)
if (any((p < -.Machine$double.eps)|(p > 1+.Machine$double.eps)))
warning(gettextf("q method of %s produced NaN's ", objN))
i01 <- (-.Machine$double.eps<=p)&(p<=1+.Machine$double.eps)
p01 <- p[i01] ## only values in [0,1] are used
q0 <- p*0
q0[!i01] <- NaN
q0[ i01] <- ifElseQC
return(as.numeric(q0))
}, list(ifElseQC = ifElseQS, objN = quote(.getObjName(1))))
return(qfun)
}
#setMethod("+", c("LatticeDistribution","AbscontDistribution"),plusDAC)
#------------------------------------------------------------------------------
# .expm.[dc], .logm.[dc] (exact calculation of exp /log)
# [ similarly possible : other monotone, smooth trafos, e.g. tan / atan ]
#------------------------------------------------------------------------------
.expm.d <- function(e1){
supportnew <- exp(e1@support)
rnew <- function(n, ...){}
body(rnew) <- substitute({ exp(f(n, ...)) },
list(f = e1@r))
dnew <- .makeD(substitute(e1, list(e1 = e1)),
substitute(alist(x = log(x*(x>0)+(x<=0)))),
fac = substitute((x>0)))
pnew <- .makeP(substitute(e1, list(e1 = e1)),
substitute(alist(q = log(q*(q>0)+(q<=0)))),
fac = substitute((q>0)),
fac2 = substitute((q<=0)& !lower.tail))
qnew <- .makeQ(substitute(e1, list(e1 = e1)),
substitute(exp(q0)),
Cont = FALSE)
object <- new("DiscreteDistribution",
r = rnew, d = dnew, p = pnew,
q = qnew, support = supportnew,
.withSim = FALSE, .withArith = TRUE,
.logExact = .logExact(e1), .lowerExact = .lowerExact(e1),
.finSupport = c(TRUE,e1@.finSupport[2]))
rm(supportnew)
rm(pnew, qnew, dnew, rnew)
object
}
.expm.c <- function(e1){
gapsnew <- if(is.null(e1@gaps)) NULL else exp(e1@gaps)
rnew <- function(n, ...){}
body(rnew) <- substitute({ exp(f(n, ...)) },
list(f = e1@r))
dnew0 <- .makeD(substitute(e1, list(e1 = e1)),
substitute(alist(x = log(x*(x>0)+(x<=0)))),
stand = substitute(x+(x==0)),
fac = substitute((x>0)))
# extrapolation for x=0
x0a <- 10^(-(1:10)/4)
f0a <- dnew0(x = x0a, log = FALSE)
f0 <- max(spline(x0a,f0a, xout = 0)$y,0)
lf0 <- log(f0)
dnew <- function(x, log = FALSE, ...){
d0 <- dnew0(x, log = log, ...)
i0 <- .isEqual(x,0)
if(!any(i0)) return(d0)
if(log) d0[i0] <- lf0 else d0[i0] <- f0
return(d0)
}
pnew <- .makeP(substitute(e1, list(e1 = e1)),
substitute(alist(q = log(q*(q>0)+(q<=0)))),
fac = substitute((q>0)),
fac2 = substitute((q<=0)& !lower.tail))
qnew <- .makeQ(substitute(e1, list(e1 = e1)),
substitute(exp(q0)))
object <- AbscontDistribution(
r = rnew, d = dnew, p = pnew,
q = qnew, gaps = gapsnew,
.withSim = FALSE, .withArith = TRUE,
.logExact = .logExact(e1), .lowerExact = .lowerExact(e1))
if(exists("gapsnew")) rm(gapsnew)
rm(pnew, qnew, dnew, rnew)
object
}
.logm.d <- function(e1){
supportnew <- log(e1@support)
rnew <- function(n, ...){}
body(rnew) <- substitute({ log(f(n, ...)) },
list(f = e1@r))
dnew <- .makeD(substitute(e1, list(e1 = e1)),
substitute(alist(x = exp(x))))
pnew <- .makeP(substitute(e1, list(e1 = e1)),
substitute(alist(q = exp(q))))
qnew <- .makeQ(substitute(e1, list(e1 = e1)),
substitute(log(q0)),
Cont = FALSE)
object <- new("DiscreteDistribution",
r = rnew, d = dnew, p = pnew,
q = qnew, support = supportnew,
.withSim = FALSE, .withArith = TRUE,
.logExact = .logExact(e1), .lowerExact = .lowerExact(e1),
.finSupport = c(TRUE,e1@.finSupport[2]))
rm(supportnew)
rm(pnew, qnew, dnew, rnew)
object
}
.logm.c <- function(e1){
gapsnew <- if(is.null(e1@gaps)) NULL else log(e1@gaps)
rnew <- function(n, ...){}
body(rnew) <- substitute({ log(f(n, ...)) },
list(f = e1@r))
dnew <- .makeD(substitute(e1, list(e1 = e1)),
substitute(alist(x = exp(x))),
fac = substitute(exp(x)))
pnew <- .makeP(substitute(e1, list(e1 = e1)),
substitute(alist(q = exp(q))))
qnew <- .makeQ(substitute(e1, list(e1 = e1)),
substitute(log(q0)))
object <- AbscontDistribution(
r = rnew, d = dnew, p = pnew,
q = qnew, gaps = gapsnew,
.withSim = FALSE, .withArith = TRUE,
.logExact = .logExact(e1), .lowerExact = .lowerExact(e1))
if(exists("gapsnew")) rm(gapsnew)
rm(pnew, qnew, dnew, rnew)
object
}
#------------------------------------------------------------------------------
# .P2D, .D2P, .P2Q, .Q2P casting functions
#------------------------------------------------------------------------------
###Functions for AbscontDistribution
#determines slot d from p
.P2D <- function(p, xx, ql, qu, ngrid = getdistrOption("DefaultNrGridPoints"))
{if(missing(xx))
xx <- seq(ql, qu, length = ngrid)
px <- p(xx)
dx <- D1ss(xx,px)
return(.makeDNew(xx, dx, h = NULL, Cont = TRUE, standM = "integrate"))
}
#determines slot q from p
.P2Q <- function(p, xx, ql,qu, ngrid = getdistrOption("DefaultNrGridPoints"),
qL = -Inf, qU = Inf){
if(missing(xx))
{xx <- seq(ql, qu, length = ngrid)
h <- (qu-ql)/ngrid}
else h <- c(diff(xx),(rev(xx)[1]-rev(xx)[2]))
px.l <- p(xx + 0.5*h)
px.u <- p(xx + 0.5*h, lower.tail = FALSE)
return(.makeQNew(xx + 0.5*h, px.l, px.u, FALSE, qL, qU))
}
#determines slot p from d
.D2P <- function(d, xx, ql, qu, ngrid = getdistrOption("DefaultNrGridPoints"))
{if(missing(xx))
{ if(ngrid%%2==0) ngrid <- ngrid+1
xx <- seq(ql, qu, length = ngrid)
h <- (qu-ql)/ngrid}
else h <- c(diff(xx),(rev(xx)[1]-rev(xx)[2]))
dx <- d(xx)
return(.makePNew(xx, dx, h = h))
}
#determines slot p from q
.Q2P <- function(q, ngrid = getdistrOption("DefaultNrGridPoints")){
q.l <- q
ep <- getdistrOption("TruncQuantile")^2
xx0 <- seq(ep, 1-ep, length = ngrid)
qx0 <- q.l(xx0)
qx1 <- unique(qx0)
x1 <- tapply(xx0, qx0, max)
p0 <- approxfun(x = qx1, y = x1, yleft = 0, yright = 1, rule = 2)
p01 <- function(x) {
p11 <- p0(x)
p11[x>max(qx1)] <- 1
p11[x<min(qx1)] <- 0
return(p11)}
return(function(q, lower.tail = TRUE, log.p = FALSE){
p1 <- p01(q)
p1[q==Inf] <- 1
p1[q==-Inf] <- 0
if(!lower.tail) p1 <- 1-p1
if(log.p) p1 <- log(p1)
return(p1)
})
}
#------------------------------------------------------------------------------
# modify slot q for AbscontDistribution if there are gaps
#------------------------------------------------------------------------------
.modifyqgaps <- function(pfun, qfun, gaps, leftright = "left"){
## no modification needed if gaps have no length
if(length(gaps)==0) return(qfun)
if(is.null(dim(gaps))) return(qfun)
finit <- apply(gaps, 1, function(x) all(is.finite(x)&is.numeric(x)))
if(sum(finit)==0) return(qfun)
gaps <- gaps[finit,,drop=FALSE]
# print(gaps)
## p-level of constancy region
if(.inArgs("log.p", pfun)){
lp.gaps <- matrix(pfun(gaps,log.p=TRUE),nrow=nrow(gaps),ncol=2)
if(.inArgs("lower.tail", pfun))
lp.gaps.l <- matrix(pfun(gaps,log.p=TRUE, lower.tail = FALSE),nrow=nrow(gaps),ncol=2)
else
lp.gaps.l <- 1-matrix(pfun(gaps,log.p=TRUE),nrow=nrow(gaps),ncol=2)
}else{
lp.gaps <- log(matrix(pfun(gaps),nrow=nrow(gaps),ncol=2))
if(.inArgs("lower.tail", pfun))
lp.gaps.l <- log(matrix(pfun(gaps,lower.tail = FALSE),nrow=nrow(gaps),ncol=2))
else
lp.gaps.l <- log1p(-matrix(pfun(gaps),nrow=nrow(gaps),ncol=2))
}
# print(lp.gaps)
# print(lp.gaps.l)
## are we heading for left or right continuous quantile fct?
lrpmatch <- pmatch(leftright, table = c("left","right"), nomatch = 1)
## in order to avoid chaining of qgaps modifications:
## place a variable ..q0fun into modified quantile function (after modification)
## which stores the unmodified quantile function
## in the first modification round ..q0fun will not yet exist
## in this situation use qfun instead
qfunN <- qfun
if(!.inArgs("log.p", qfun) || !.inArgs("lower.tail", qfun)){
qfunN <- function(p, lower.tail = TRUE, log.p = FALSE){
if(.inArgs("lower.tail",qfun)) if(log.p) p <- exp(p)
if(.inArgs("log.p", qfun)){
p1 <- if(log.p) exp(p) else p
qval <- if(lower.tail) qfun(p, log.p) else qfun(1-p1, log.p=FALSE)
}else{
if(!.inArgs("lower.tail",qfun)){
qval <- if(lower.tail) qfun(p) else qfun(1-p)
}else{
qval <- qfun(p,lower.tail)
}
}
return(qval)
}
}
qfunE <- environment(qfunN)
qnew <- function(p, lower.tail = TRUE, log.p = FALSE) {}
qnewE <- environment(qnew) <- new.env()
body(qnew) <- substitute({
## .q0fun is the (gaps-)unmodified quantile function
.q0fun <- if(exists("..q0fun", envir=qfunE.)){
get("..q0fun", envir=qfunE.) } else qfun.
q0 <- .q0fun(p, lower.tail = lower.tail, log.p = log.p)
## the gaps-modification: find out which args p coincide
## (numerically, on log scale) with gaps-plevels;
## depending on "leftright" and lower.tail
## set these return values to left or right endpoit of the gap
if(length(lp.gaps.)>0){
i0 <- seq(length=length(p))
lg <- round(3/2-(2*lower.tail-1)*(2*(lrpmatch.==1)-1)/2)
## ==1 if(lower.tail&&leftright==1) or (!lower.tail&&leftright!=1)
## and == 2 otherwise
lpgaps0 <- if(lg==1L) lp.gaps. else lp.gaps.l.
for(i in 1:nrow(lpgaps0)){
i0 <- (log(p)>=lpgaps0[i,1])&(log(p)<=lpgaps0[i,2])
if(length(i0)) q0[i0] <- gaps.[i,lg]
}
}
return(q0)
},list(qfunE. = qfunE, qfun.=qfunN, lp.gaps.=lp.gaps,
lp.gaps.l. = lp.gaps.l[,2:1,drop=FALSE],
lrpmatch. = lrpmatch, gaps. = gaps, ..isEqual = .isEqual)
)
if(exists("..q0fun", envir=qfunE)){
.q0fun <- get("..q0fun", envir=qfunE)
assign("..q0fun", .q0fun, envir = qnewE)
}else{
assign("..q0fun", qfun, envir = qnewE)
}
return(qnew)
}
if(FALSE){ ## old code
.modifyqgaps <- function(pfun, qfun, gaps, leftright = "left"){
if(length(gaps)==0) return(qfun)
p.gaps <- pfun(gaps[,1])
p.gaps.l <- pfun(gaps[,1], lower.tail = FALSE)
dP <- deparse(body(qfun))
syC <- paste(sys.calls())
if(length(grep("q\\.r\\(",syC[length(syC)-2])) == 0)
if(length(grep("q0 <- qfun", dP[2]))>0)
return(qfun)
if(pmatch(leftright, table = c("left","right"), nomatch = 1) == 1){
qnew <- function(p, lower.tail = TRUE, log.p = FALSE) {
q0 <- qfun(p, lower.tail = lower.tail, log.p = log.p)
i0 <- seq(length=length(p))
if(lower.tail){
if(log.p) p.gaps <- log(p.gaps)
for(i in 1:nrow(gaps)){
i1 <- .isEqual(p,p.gaps[i])
i2 <- i0[p<p.gaps[i]]
i3 <- i0[p>p.gaps[i]]
q0[i1] <- gaps[i,1]
if(length(i2)>0)
q0[i2] <- pmin(q0[i2],gaps[i,1])
if(length(i3)>0)
q0[i3] <- pmax(q0[i3],gaps[i,2])
}
}else{
if(log.p) p.gaps.l <- log(p.gaps.l)
for(i in 1:nrow(gaps)){
i1 <- .isEqual(p,p.gaps.l[i])
i2 <- i0[p<p.gaps.l[i]]
i3 <- i0[p>p.gaps.l[i]]
q0[i1] <- gaps[i,2]
if(length(i2)>0)
q0[i2] <- pmax(q0[i2],gaps[i,2])
if(length(i3)>0)
q0[i3] <- pmin(q0[i3],gaps[i,1])
}
}
return(q0)
}
}else{
qnew <- function(p, lower.tail = TRUE, log.p = FALSE) {
q0 <- qfun(p, lower.tail = lower.tail, log.p = log.p)
i0 <- seq(length=length(p))
if(lower.tail){
if(log.p) p.gaps <- log(p.gaps)
for(i in 1:nrow(gaps)){
i1 <- .isEqual(p,p.gaps[i])
i2 <- i0[p<p.gaps[i]]
i3 <- i0[p>p.gaps[i]]
q0[i1] <- gaps[i,2]
if(length(i2)>0)
q0[i2] <- pmin(q0[i2],gaps[i,1])
if(length(i3)>0)
q0[i3] <- pmax(q0[i3],gaps[i,2])
}
}else{
if(log.p) p.gaps.l <- log(p.gaps.l)
for(i in 1:nrow(gaps)){
i1 <- .isEqual(p,p.gaps.l[i])
i2 <- i0[p<p.gaps.l[i]]
i3 <- i0[p>p.gaps.l[i]]
q0[i1] <- gaps[i,1]
if(length(i2)>0)
q0[i2] <- pmax(q0[i2],gaps[i,2])
if(length(i3)>0)
q0[i3] <- pmin(q0[i3],gaps[i,1])
}
}
return(q0)
}
}
return(qnew)
}
}
#------------------------------------------------------------------------------
# issue warnings in show / print as to Arith or print
#------------------------------------------------------------------------------
.IssueWarn <- function(Arith,Sim){
msgA1 <- msgA2 <- msgS1 <- msgS2 <- NULL
if(Arith && getdistrOption("WarningArith")){
msgA1 <- gettext(
"arithmetics on distributions are understood as operations on r.v.'s\n")
msgA2 <- gettext(
"see 'distrARITH()'; for switching off this warning see '?distroptions'")
}
if(Sim && getdistrOption("WarningSim")){
msgS1 <- gettext(
"slots d,p,q have been filled using simulations; ")
msgS2 <- gettext(
"for switching off this warning see '?distroptions'")
}
return(list(msgA=c(msgA1,msgA2), msgS = c(msgS1,msgS2)))
}
#------------------------------------------------------------------------------
# fill a list acc. recycling rules
#------------------------------------------------------------------------------
.List <- function(list0) if(is.list(list0)) list0 else list(list0)
.fillList <- function(list0, len = length(list0)){
if(is.null(list0)) return(vector("list",len))
list0 <- .List(list0)
if(len == length(list0))
return(list0)
i <- 0
ll0 <- length(list0)
li0 <- vector("list",len)
if(ll0){
while(i < len){
j <- 1 + ( i %% ll0)
i <- i + 1
li0[[i]] <- list0[[j]]
}
}
return(li0)
}
#------------------------------------------------------------------------------
# .DistributionAggregate by Jacob van Etten, jacobvanetten@yahoo.com
# on a mail on Feb 27th 2009
#------------------------------------------------------------------------------
.DistrCollapse <- function(support, prob,
eps = getdistrOption("DistrResolution")){
supp <- support
prob <- as.vector(prob)
suppIncr <- diff(c(supp[1]-2*eps,supp)) < eps
groups <- cumsum(!suppIncr)
prob <- as.vector(tapply(prob, groups, sum))
supp <- as.vector(tapply(supp, groups, quantile, probs = 0.5, type = 1))
### in order to get a "support member" take the leftmost median
return(list(supp = supp, prob = prob))
# newDistribution <- DiscreteDistribution(supp=supp,prob=prob)
# return(newDistribution)
}
.panel.mingle <- function(dots, element){
pF <- dots[[element]]
if(is.list(pF)) return(pF)
pFr <- if(typeof(pF)=="symbol") eval(pF) else{
pFc <- as.call(pF)
if(as.list(pFc)[[1]] == "list"){
lis <- vector("list",length(as.list(pFc))-1)
for(i in 1:length(lis)){
lis[[i]] <- pFc[[i+1]]
}
lis
}else pF
}
return(pFr)
}
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.