Nothing
################################################################
# QQ - Plot functions in package distrMod
################################################################
### to be written into the respective MASKING files....
## helper into distrMod
.labelprep <- function(x,y,lab.pts,col.lbs,cex.lbs,adj.lbs,which.lbs,which.Order,order.traf, which.nonlbs){
n <- length(x)
ind0 <- 1:n
# first selection with which.lbs
ind1 <- ind0
if(!is.null(which.lbs)) ind1 <- ind0[ind0%in%which.lbs]
# second selection with which.Order
n1 <- length(ind1)
x1 <- x[ind1]
rk1.0 <- rank(x1)
if(!is.null(order.traf)) rk1 <- rank(order.traf(x1))
rk1 <- n1+1-rk1.0
#
ind2 <- ind1
if(!is.null(which.Order)) ind2 <- ind1[rk1 %in% which.Order]
#
x2 <- x[ind2]
or2.0 <- order(x2, decreasing = TRUE)
#
ind.s <- ind2[or2.0]
#
ind.ns <- ind0[-ind2]
if(length(ind.ns) && !is.null(which.nonlbs))
ind.ns <- ind.ns[ind.ns%in%which.nonlbs]
#
#------------------------------------------------------------------------
x0 <- x[ind.s]
y0 <- x[ind.s]
col.lbs <- col.lbs[ind.s]
lab.pts <- lab.pts[ind.s]
cex.lbs <- cex.lbs[ind.s]
adj.lbs <- adj.lbs[ind.s]
return(list(x0=x0,y0=y0,lab=lab.pts,col=col.lbs,cex=cex.lbs,adj=adj.lbs,ord=ind.s, ns=ind.ns))
}
setMethod("qqplot", signature(x = "ANY",
y = "UnivariateDistribution"),
function(x, ### observations
y, ### distribution
n = length(x), ### number of points to be plotted
withIdLine = TRUE, ### shall line y=x be plotted in
withConf = TRUE, ### shall confidence lines be plotted
withConf.pw = withConf, ### shall pointwise confidence lines be plotted
withConf.sim = withConf, ### shall simultaneous confidence lines be plotted
plot.it = TRUE, ### shall be plotted at all (inherited from stats::qqplot)
datax = FALSE, ### as in qqnorm
xlab = deparse(substitute(x)), ## x-label
ylab = deparse(substitute(y)), ## y-label
..., ## further parameters
width = 10, ## width (in inches) of the graphics device opened
height = 5.5, ## height (in inches) of the graphics device opened}
withSweave = getdistrOption("withSweave"), ## logical: if \code{TRUE}
## (for working with \command{Sweave}) no extra device is opened and height/width are not set
mfColRow = TRUE, ## shall we use panel partition mfrow=c(1,1)?
n.CI = n, ## number of points to be used for CI
with.lab = FALSE, ## shall observation labels be plotted in
lab.pts = NULL, ## observation labels to be used
which.lbs = NULL, ## which observations shall be labelled
which.Order = NULL, ## which of the ordered (remaining) observations shall be labelled
which.nonlbs = NULL, ## which of the non-labelled observations shall be plotted
attr.pre = FALSE, ## do indices refer to order pre or post ordering
order.traf = NULL, ## an optional trafo; by which the observations are ordered (as order(trafo(obs))
col.IdL = "red", ## color for the identity line
lty.IdL = 2, ## line type for the identity line
lwd.IdL = 2, ## line width for the identity line
alpha.CI = .95, ## confidence level
exact.pCI = (n<100), ## shall pointwise CIs be determined with exact Binomial distribution?
exact.sCI = (n<100), ## shall simultaneous CIs be determined with exact kolmogorov distribution?
nosym.pCI = FALSE, ## shall we use (shortest) asymmetric CIs?
col.pCI = "orange", ## color for the pointwise CI
lty.pCI = 3, ## line type for the pointwise CI
lwd.pCI = 2, ## line width for the pointwise CI
pch.pCI = par("pch"),## symbol for points (for discrete mass points) in pointwise CI
cex.pCI = par("cex"),## magnification factor for points (for discrete mass points) in pointwise CI
col.sCI = "tomato2", ## color for the simultaneous CI
lty.sCI = 4, ## line type for the simultaneous CI
lwd.sCI = 2, ## line width for the simultaneous CI
pch.sCI = par("pch"),## symbol for points (for discrete mass points) in simultaneous CI
cex.sCI = par("cex"),## magnification factor for points (for discrete mass points) in simultaneous CI
added.points.CI = TRUE, ## should the CIs be drawn through additional points?
cex.pch = par("cex"),## magnification factor for the plotted symbols (for backward compatibility only, cex.pts in the sequel)
col.pch = par("col"),## color for the plotted symbols (for backward compatibility only, col.pts in the sequel)
cex.pts = 1, ## magnification factor for labelled shown observations
col.pts = par("col"),## color for labelled shown observations
pch.pts = 19, ## symbol for labelled shown observations
cex.npts = 1, ## magnification factor for non-labelled shown observations
col.npts = grey(.5), ## color for non-labelled shown observations
pch.npts = 20, ## symbol for non-labelled shown observations
cex.lbs = par("cex"),## magnification factor for the plotted observation labels
col.lbs = par("col"),## color for the plotted observation labels
adj.lbs = par("adj"),## adj parameter for the plotted observation labels
alpha.trsp = NA, ## alpha transparency to be added afterwards
jit.fac = 0, ## jittering factor used for discrete distributions
jit.tol = .Machine$double.eps, ## tolerance for jittering: if distance
#is smaller than jit.tol, points are considered replicates
check.NotInSupport = TRUE, ## shall we check if all x lie in support(y)
col.NotInSupport = "red", ## if preceding check TRUE color of x if not in support(y)
with.legend = TRUE, ## shall a legend be plotted
legend.bg = "white", ## background for the legend
legend.pos = "topleft", ## position for the legend
legend.cex = 0.8, ## magnification factor for the legend
legend.pref = "", ## prefix for legend text
legend.postf = "", ## postfix for legend text
legend.alpha = alpha.CI, ## nominal level of CI
debug = FALSE, ## shall additional debug output be printed out?
withSubst = TRUE
){ ## return value as in stats::qqplot
mc <- match.call(call = sys.call(sys.parent(1)))
xcc <- as.character(deparse(mc$x))
ycc <- as.character(deparse(mc$y))
if(missing(xlab)){ xlab <- mc$xlab <- xcc}
if(missing(ylab)){ ylab <- mc$ylab <- ycc}
dots <- match.call(call = sys.call(sys.parent(1)),
expand.dots = FALSE)$"..."
args0 <- list(x = x, y = y, n = n, withIdLine = withIdLine,
withConf = withConf, withConf.pw = withConf.pw,
withConf.sim = withConf.sim, plot.it = plot.it, datax = datax,
xlab = xlab, ylab = ylab, width = width, height = height,
withSweave = withSweave, mfColRow = mfColRow,
n.CI = n.CI, with.lab = with.lab, lab.pts = lab.pts,
which.lbs = which.lbs, which.Order = which.Order,
order.traf = order.traf, col.IdL = col.IdL, lty.IdL = lty.IdL,
lwd.IdL = lwd.IdL, alpha.CI = alpha.CI, exact.pCI = exact.pCI,
exact.sCI = exact.sCI, nosym.pCI = nosym.pCI, col.pCI = col.pCI,
lty.pCI = lty.pCI, lwd.pCI = lwd.pCI, pch.pCI = pch.pCI,
cex.pCI = cex.pCI, col.sCI = col.sCI, lty.sCI = lty.sCI,
lwd.sCI = lwd.sCI, pch.sCI = pch.sCI, cex.sCI = cex.sCI,
added.points.CI = added.points.CI, cex.pch = cex.pch,
col.pch = col.pch, cex.lbs = cex.lbs, col.lbs = col.lbs,
adj.lbs = adj.lbs, alpha.trsp = alpha.trsp, jit.fac = jit.fac,
jit.tol = jit.tol, check.NotInSupport = check.NotInSupport,
col.NotInSupport = col.NotInSupport, with.legend = with.legend,
legend.bg = legend.bg, legend.pos = legend.pos,
legend.cex = legend.cex, legend.pref = legend.pref,
legend.postf = legend.postf, legend.alpha = legend.alpha,
debug = debug, withSubst = withSubst)
plotInfo <- list(call = mc, dots=dots, args=args0)
mcl <- as.list(mc)[-1]
force(x)
if(is.null(mcl$datax)) datax <- FALSE
if(!datax){ mcl$ylab <- xlab; mcl$xlab <- ylab}
.mpresubs <- if(withSubst){
function(inx)
.presubs(inx, c("%C", "%A", "%D" ),
c(as.character(class(x)[1]),
as.character(date()),
xcc))
}else function(inx)inx
rank0x <- rank(x)
xj <- sort(x)
if(any(.isReplicated(x, jit.tol))&&jit.fac>0)
xj[.isReplicated(x, jit.tol)] <- jitter(x[.isReplicated(x, jit.tol)], factor=jit.fac)
rank1x <- rank(xj)[rank0x]
ind.x <- seq(along=x)
xj <- sort(xj)
pp <- ppoints(n)
yc <- q.l(y)(pp)
yc.o <- yc
if("support" %in% names(getSlots(class(y))))
yc <- sort(jitter(yc, factor=jit.fac))
alp.v <- .makeLenAndOrder(alpha.trsp,ind.x)
alp.t <- function(x,a1) if(is.na(x)) x else addAlphTrsp2col(x,a1)
alp.f <- if(length(alpha.trsp)==1L && is.na(alpha.trsp))
function(x,a) x else function(x,a) mapply(x,alp.t,a1=a)
if(missing(cex.lbs)) cex0.lbs <- par("cex")
cex0.lbs <- .makeLenAndOrder(cex.lbs,ind.x)
if(missing(adj.lbs)) adj0.lbs <- par("adj")
adj0.lbs <- .makeLenAndOrder(adj.lbs,ind.x)
if(missing(col.lbs)) col0.lbs <- par("col")
col0.lbs <- alp.f(.makeLenAndOrder(col.lbs,ind.x),alp.v)
if(missing(lab.pts)||is.null(lab.pts)) lab0.pts <- ind.x else
lab0.pts <- .makeLenAndOrder(lab.pts,ind.x)
lbprep <- .labelprep(x = x, y = yc.o[rank1x], lab.pts = lab0.pts,
col.lbs = col0.lbs, cex.lbs = cex0.lbs,
adj.lbs = adj0.lbs, which.lbs = which.lbs,
which.Order = which.Order, order.traf = order.traf,
which.nonlbs = which.nonlbs)
n.ns <- length(lbprep$ns)
n.s <- length(lbprep$ord)
shown <- c(lbprep$ord,lbprep$ns)
xs <- x[shown]
ycs <- (yc.o[rank1x])[shown]
ordx <- order(xs)
xso <- xs[ordx]
ycso <- ycs[ordx]
if(missing(cex.pch)) cex.pch <- par("cex")
if(missing(col.pch)) col.pch <- par("col")
if(missing(cex.pts)) cex.pts <- if(missing(cex.pch)) 1 else cex.pch
if(missing(col.pts)) col.pts <- if(missing(col.pch)) par("col") else col.pch
if(missing(pch.pts)) pch.pts <- 19
if(missing(cex.npts)) cex.npts <- 1
if(missing(col.npts)) col.npts <- par("col")
if(missing(pch.npts)) pch.npts <- 20
if(with.lab) lab.pts <- lbprep$lab.pts
if(attr.pre){
if(with.lab){
col.lbs <- lbprep$col.lbs
cex.lbs <- lbprep$cex.lbs
adj.lbs <- lbprep$adj.lbs
}
cex.pts <- .makeLenAndOrder(cex.pts,ind.x)
col.pts <- alp.f(.makeLenAndOrder(col.pts,ind.x),alp.v)
pch.pts <- .makeLenAndOrder(pch.pts,ind.x)
cex.pts <- cex.pts[shown]
col.pts <- col.pts[shown]
pch.pts <- pch.pts[shown]
}else{
ind.s <- 1:n.s
ind.ns <- 1:n.ns
if(with.lab){
if(missing(cex.lbs)) cex.lbs <- par("cex")
cex.lbs <- (.makeLenAndOrder(cex.lbs,ind.s))
if(missing(adj.lbs)) adj.lbs <- par("adj")
adj.lbs <- (.makeLenAndOrder(adj.lbs,ind.s))
if(missing(col.lbs)) col.lbs <- par("col")
col.lbs <- (alp.f(.makeLenAndOrder(col.lbs,ind.s),alp.v[lbprep$ord]))
}
cex.pts <- .makeLenAndOrder(cex.pts,ind.s)
col.pts <- alp.f(.makeLenAndOrder(col.pts,ind.s),alp.v[lbprep$ord])
pch.pts <- .makeLenAndOrder(pch.pts,ind.s)
cex.npts <- .makeLenAndOrder(cex.npts,ind.ns)
col.npts <- alp.f(.makeLenAndOrder(col.npts,ind.ns),alp.v[lbprep$ns])
pch.npts <- .makeLenAndOrder(pch.npts,ind.ns)
col.pts <- c(col.pts,col.npts)
cex.pts <- c(cex.pts,cex.npts)
pch.pts <- c(pch.pts,pch.npts)
}
cex.pts <- cex.pts[ordx]
col.pts <- col.pts[ordx]
pch.pts <- pch.pts[ordx]
if(check.NotInSupport){
xo <- xso #x[ord.x]
nInSupp <- which(xo < q.l(y)(0))
nInSupp <- unique(sort(c(nInSupp,which( xo > q.l(y)(1)))))
if("support" %in% names(getSlots(class(y))))
nInSupp <- unique(sort(c(nInSupp,which( ! xo %in% support(y)))))
if("gaps" %in% names(getSlots(class(y))))
nInSupp <- unique(sort(c(nInSupp,which( .inGaps(xo,gaps(y))))))
if(length(nInSupp)){
# col.pch[nInSupp] <- col.NotInSupport
col.pts[nInSupp] <- col.NotInSupport
if(with.lab)
# col.lbs[ord.x[nInSupp]] <- col.NotInSupport
col.lbs[nInSupp] <- col.NotInSupport
}
}
if(n < length(x)){
with.lab <- FALSE
nos <- length(shown)
idx <- sample(1:nos,size=n,replace=FALSE)
cex.pts <- cex.pts[idx]
col.pts <- col.pts[idx]
pch.pts <- pch.pts[idx]
xso <- xso[idx]
ycso <- ycso[idx]
}
if(datax){
mcl$x <- xso#xj
mcl$y <- ycso #yc
}else{
mcl$y <- xso# xj
mcl$x <- ycso #yc
}
mcl <- .deleteItemsMCL(mcl)
mcl$pch <- pch.pts
mcl$cex <- cex.pts
mcl$col <- col.pts
mcl$xlab <- .mpresubs(mcl$xlab)
mcl$ylab <- .mpresubs(mcl$ylab)
if (!is.null(eval(mcl$main)))
mcl$main <- .mpresubs(eval(mcl$main))
if (!is.null(eval(mcl$sub)))
mcl$sub <- .mpresubs(eval(mcl$sub))
if (!withSweave){
devNew(width = width, height = height)
}
opar <- par("mfrow", no.readonly = TRUE)
if(mfColRow) on.exit(do.call(par, list(mfrow=opar, no.readonly = TRUE)))
if(mfColRow) opar1 <- par(mfrow = c(1,1), no.readonly = TRUE)
ret <- do.call(stats::qqplot, args=mcl)
qq.usr <- par("usr")
if(with.lab&& plot.it){
xlb0 <- if(datax) lbprep$x0 else lbprep$y0
ylb0 <- if(datax) lbprep$y0 else lbprep$x0
text(x = xlb0, y = ylb0, labels = lbprep$lab,
cex = lbprep$cex, col = lbprep$col, adj = lbprep$adj)
}
qqb <- NULL
if(withIdLine){
if(plot.it) abline(0,1,col=col.IdL,lty=lty.IdL,lwd=lwd.IdL)
if(#is(y,"AbscontDistribution")&&
withConf){
xy <- unique(sort(c(x,yc.o)))
if(added.points.CI){
mxy <- min(xy); Mxy <- max(xy)
mnxy <- (mxy+Mxy)/2
sxy <- (Mxy-mxy)/2*1.1
xyn <- seq(mnxy-sxy,mnxy+sxy,length.out=500)
xy <- unique(sort(c(xy,xyn)))
}
xy <- xy[!.NotInSupport(xy,y)]
lxy <- length(xy)
if(is(y,"DiscreteDistribution")){
n0 <- min(n.CI, length(support(y)))
n1 <- max(n0-lxy,0)
if (n1 >0 ){
notyetInXY <- setdiff(support(y), xy)
xy0 <- sample(notyetInXY, n1)
xy <- sort(unique(c(xy,xy0)))
}
}else{
if(lxy < n.CI){
n1 <- (n.CI-lxy)%/%3
xy0 <- seq(min(xy),max(xy),length=n1)
xy1 <- r(y)(n.CI-lxy-n1)
xy <- sort(unique(c(xy,xy0,xy1)))
}
}
qqplotInfo <- list(xy.0=xy, y.0=y, datax = datax,
withConf.pw=withConf.pw,
withConf.sim=withConf.sim,
alpha.CI=alpha.CI ,
col.pCI = col.pCI , lty.pCI = lty.pCI ,
lwd.pCI = lwd.pCI , pch.pCI = pch.pCI,
cex.pCI = cex.pCI ,
col.sCI = col.sCI , lty.sCI = lty.sCI ,
lwd.sCI = lwd.sCI , pch.sCI = pch.sCI,
cex.sCI = cex.sCI ,
n = n ,
exact.sCI = exact.sCI, exact.pCI = exact.pCI,
nosym.pCI = nosym.pCI, with.legend = with.legend,
legend.bg = legend.bg, legend.pos = legend.pos,
legend.cex = legend.cex, legend.pref = legend.pref,
legend.postf = legend.postf, legend.alpha = legend.alpha,
debug = debug,
args.stats.qqplot = mcl,
with.lab = with.lab,
lbprep = lbprep
)
if(plot.it){
qqb <- .confqq(xy, y, datax=datax, withConf.pw, withConf.sim, alpha.CI,
col.pCI, lty.pCI, lwd.pCI, pch.pCI, cex.pCI,
col.sCI, lty.sCI, lwd.sCI, pch.sCI, cex.sCI,
n, exact.sCI = exact.sCI, exact.pCI = exact.pCI,
nosym.pCI = nosym.pCI, with.legend = with.legend,
legend.bg = legend.bg, legend.pos = legend.pos,
legend.cex = legend.cex, legend.pref = legend.pref,
legend.postf = legend.postf, legend.alpha = legend.alpha, debug = debug)
}else{
qqb <- qqbounds(sort(unique(xy)),y,alpha.CI,n,withConf.pw, withConf.sim,
exact.sCI,exact.pCI,nosym.pCI, debug = debug)
}
}
}
qqplotInfo <- c(plotInfo, ret=ret, usr=qq.usr, qqplotInfo=qqplotInfo, qqb=qqb)
class(qqplotInfo) <- c("qqplotInfo","DiagnInfo")
return(invisible(qqplotInfo))
})
## into distrMod
setMethod("qqplot", signature(x = "ANY",
y = "ProbFamily"), function(x, y,
n = length(x), withIdLine = TRUE, withConf = TRUE,
withConf.pw = withConf, withConf.sim = withConf,
plot.it = TRUE, xlab = deparse(substitute(x)),
ylab = deparse(substitute(y)), ...){
mc <- match.call(call = sys.call(sys.parent(1)))
mc1 <- match.call(call = sys.call(sys.parent(1)), expand.dots=FALSE)
mcx <- as.character(deparse(mc$x))
mcy <- as.character(deparse(mc$y))
dots <- mc1$"..."
args0 <- list(x = x, y = y,
n = if(!missing(n)) n else length(x),
withIdLine = withIdLine, withConf = withConf,
withConf.pw = if(!missing(withConf.pw)) withConf.pw else if(!missing(withConf)) withConf else NULL,
withConf.sim = if(!missing(withConf.sim)) withConf.sim else if(!missing(withConf)) withConf else NULL,
plot.it = plot.it, xlab = xlab, ylab = ylab)
plotInfo <- list(call=mc, dots=dots, args=args0)
if(missing(xlab)) mc$xlab <- mcx
if(missing(ylab)) mc$ylab <- mcy
mcl <- as.list(mc)[-1]
mcl$y <- yD <- y@distribution
if(!is(yD,"UnivariateDistribution"))
stop("Not yet implemented.")
retv <- do.call(getMethod("qqplot", signature(x="ANY", y="UnivariateDistribution")),
args=mcl)
retv$call <- retv$dots <- retv$args <- NULL
plotInfo <- c(plotInfo,retv)
class(plotInfo) <- c("qqplotInfo","DiagnInfo")
return(invisible(plotInfo))
})
setMethod("qqplot", signature(x = "ANY",
y = "Estimate"), function(x, y,
n = length(x), withIdLine = TRUE, withConf = TRUE,
withConf.pw = withConf, withConf.sim = withConf,
plot.it = TRUE, xlab = deparse(substitute(x)),
ylab = deparse(substitute(y)), ...){
mc <- match.call(call = sys.call(sys.parent(1)))
mc1 <- match.call(call = sys.call(sys.parent(1)), expand.dots=FALSE)
mcx <- as.character(deparse(mc$x))
mcy <- as.character(deparse(mc$y))
args0 <- list(x=x,y=y,n=n,withIdLine=withIdLine, withConf=withConf,
withConf.pw = if(!missing(withConf.pw)) withConf.pw else if(!missing(withConf)) withConf else NULL,
withConf.sim = if(!missing(withConf.sim)) withConf.sim else if(!missing(withConf)) withConf else NULL,
plot.it = plot.it, xlab = xlab, ylab = ylab)
dots <- mc1$"..."
plotInfo <- list(call=mc, dots=dots, args=args0)
if(missing(xlab)) mc$xlab <- mcx
mcl <- as.list(mc)[-1]
param <- ParamFamParameter(main=untransformed.estimate(y), nuisance=nuisance(y),
fixed=fixed(y))
es.call <- y@estimate.call
nm.call <- names(es.call)
PFam <- NULL
if("ParamFamily" %in% nm.call)
PFam <- eval(as.list(es.call)[["ParamFamily"]])
if(is.null(PFam))
stop("There is no object of class 'ProbFamily' in the call of 'x'")
PFam0 <- modifyModel(PFam, param)
mcl$y <- PFam0
if(missing(ylab)) mcl$ylab <- paste(name(PFam0),gettext("fitted by"), mcy)
retv <- do.call(getMethod("qqplot", signature(x="ANY", y="ProbFamily")),
args=mcl)
retv$call <- retv$dots <- retv$args <- NULL
plotInfo <- c(plotInfo,retv)
class(plotInfo) <- c("qqplotInfo","DiagnInfo")
return(invisible(plotInfo))
})
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.