[Distr-commits] r970 - in branches/distr-2.6/pkg: distr/R distr/man distrMod/R distrMod/man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Mon Oct 6 16:25:00 CEST 2014
Author: ruckdeschel
Date: 2014-10-06 16:25:00 +0200 (Mon, 06 Oct 2014)
New Revision: 970
Modified:
branches/distr-2.6/pkg/distr/R/internals-qqplot.R
branches/distr-2.6/pkg/distr/man/internals-qqplot.Rd
branches/distr-2.6/pkg/distrMod/R/0distrModUtils.R
branches/distr-2.6/pkg/distrMod/R/qqplot.R
branches/distr-2.6/pkg/distrMod/R/returnlevelplot.R
branches/distr-2.6/pkg/distrMod/man/qqplot.Rd
branches/distr-2.6/pkg/distrMod/man/returnlevelplot.Rd
Log:
[distr+distrMod] revised qqplot and returnlevelplot
Modified: branches/distr-2.6/pkg/distr/R/internals-qqplot.R
===================================================================
--- branches/distr-2.6/pkg/distr/R/internals-qqplot.R 2014-10-05 21:06:31 UTC (rev 969)
+++ branches/distr-2.6/pkg/distr/R/internals-qqplot.R 2014-10-06 14:25:00 UTC (rev 970)
@@ -270,7 +270,7 @@
-.confqq <- function(x,D, datax = TRUE, withConf.pw = TRUE,
+.confqq <- function(x,D, datax = FALSE, withConf.pw = TRUE,
withConf.sim = TRUE, alpha,
col.pCI, lty.pCI, lwd.pCI, pch.pCI, cex.pCI,
col.sCI, lty.sCI, lwd.sCI, pch.sCI, cex.sCI,
@@ -278,7 +278,7 @@
with.legend = TRUE, legend.bg = "white",
legend.pos = "topleft", legend.cex = 0.8,
legend.pref = "", legend.postf = "",
- legend.alpha = alpha, qqb0=NULL, debug = FALSE){
+ legend.alpha = alpha, qqb0=NULL, transf0=NULL, debug = FALSE){
x <- sort(unique(x))
if("gaps" %in% names(getSlots(class(D))))
@@ -293,37 +293,38 @@
x.in <- x[SI.in]
x.c <- x.in[SI.c]
x.d <- x.in[!SI.c]
-
+ tx.c <- if(is.null(transf0)) x.c else transf0(x.c)
+ tx.d <- if(is.null(transf0)) x.d else transf0(x.d)
qqb <- if(is.null(qqb0)) qqbounds(x,D,alpha,n,withConf.pw, withConf.sim,
exact.sCI,exact.pCI,nosym.pCI, debug) else qqb0
- qqb$crit <- qqb$crit[SI.in,]
+ qqb$crit <- if(is.null(qqb0)) qqb$crit[SI.in,]
if(qqb$err["pw"]){
if(sum(SI.c)>0){
if(datax){
- lines(x.c, qqb$crit[SI.c,"pw.right"],
+ lines(tx.c, qqb$crit[SI.c,"pw.right"],
col=col.pCI,lty=lty.pCI,lwd=lwd.pCI)
- lines(x.c, qqb$crit[SI.c,"pw.left"],
+ lines(tx.c, qqb$crit[SI.c,"pw.left"],
col=col.pCI,lty=lty.pCI,lwd=lwd.pCI)
}else{
- lines(qqb$crit[SI.c,"pw.right"], x.c,
+ lines(qqb$crit[SI.c,"pw.right"], tx.c,
col=col.pCI,lty=lty.pCI,lwd=lwd.pCI)
- lines(qqb$crit[SI.c,"pw.left"], x.c,
+ lines(qqb$crit[SI.c,"pw.left"], tx.c,
col=col.pCI,lty=lty.pCI,lwd=lwd.pCI)
}
}
if(sum(!SI.c)>0){
if(datax){
- points(x.d, qqb$crit[!SI.c,"pw.right"],
+ points(tx.d, qqb$crit[!SI.c,"pw.right"],
col=col.pCI, pch=pch.pCI, cex = cex.pCI)
- points(x.d, qqb$crit[!SI.c,"pw.left"],
+ points(tx.d, qqb$crit[!SI.c,"pw.left"],
col=col.pCI, pch=pch.pCI, cex = cex.pCI)
}else{
- points(qqb$crit[!SI.c,"pw.right"], x.d,
+ points(qqb$crit[!SI.c,"pw.right"], tx.d,
col=col.pCI, pch=pch.pCI, cex = cex.pCI)
- points(qqb$crit[!SI.c,"pw.left"], x.d,
+ points(qqb$crit[!SI.c,"pw.left"], tx.d,
col=col.pCI, pch=pch.pCI, cex = cex.pCI)
}
}
@@ -331,27 +332,27 @@
if(qqb$err["sim"]){
if(sum(SI.c)>0){
if(datax){
- lines(x.c, qqb$crit[SI.c,"sim.right"],
+ lines(tx.c, qqb$crit[SI.c,"sim.right"],
col=col.sCI,lty=lty.sCI,lwd=lwd.sCI)
- lines(x.c, qqb$crit[SI.c,"sim.left"],
+ lines(tx.c, qqb$crit[SI.c,"sim.left"],
col=col.sCI,lty=lty.sCI,lwd=lwd.sCI)
}else{
- lines(qqb$crit[SI.c,"sim.right"], x.c,
+ lines(qqb$crit[SI.c,"sim.right"], tx.c,
col=col.sCI,lty=lty.sCI,lwd=lwd.sCI)
- lines(qqb$crit[SI.c,"sim.left"], x.c,
+ lines(qqb$crit[SI.c,"sim.left"], tx.c,
col=col.sCI,lty=lty.sCI,lwd=lwd.sCI)
}
}
if(sum(!SI.c)>0){
if(datax){
- points(x.d, qqb$crit[!SI.c,"sim.right"],
+ points(tx.d, qqb$crit[!SI.c,"sim.right"],
col=col.sCI, pch=pch.sCI, cex = cex.sCI)
- points(x.d, qqb$crit[!SI.c,"sim.left"],
+ points(tx.d, qqb$crit[!SI.c,"sim.left"],
col=col.sCI, pch=pch.sCI, cex = cex.sCI)
}else{
- points(qqb$crit[!SI.c,"sim.right"], x.d,
+ points(qqb$crit[!SI.c,"sim.right"], tx.d,
col=col.sCI, pch=pch.sCI, cex = cex.sCI)
- points(qqb$crit[!SI.c,"sim.left"], x.d,
+ points(qqb$crit[!SI.c,"sim.left"], tx.d,
col=col.sCI, pch=pch.sCI, cex = cex.sCI)
}
}
Modified: branches/distr-2.6/pkg/distr/man/internals-qqplot.Rd
===================================================================
--- branches/distr-2.6/pkg/distr/man/internals-qqplot.Rd 2014-10-05 21:06:31 UTC (rev 969)
+++ branches/distr-2.6/pkg/distr/man/internals-qqplot.Rd 2014-10-06 14:25:00 UTC (rev 970)
@@ -33,14 +33,14 @@
.q2kolmogorov(alpha,n,exact=(n<100), silent0 = TRUE)
.q2pw(x,p.b,D,n,alpha,exact=(n<100),nosym=FALSE, silent0 = TRUE)
-.confqq(x,D, datax=TRUE, withConf.pw = TRUE, withConf.sim = TRUE, alpha,
+.confqq(x,D, datax=FALSE, withConf.pw = TRUE, withConf.sim = TRUE, alpha,
col.pCI, lty.pCI, lwd.pCI, pch.pCI, cex.pCI,
col.sCI, lty.sCI, lwd.sCI, pch.sCI, cex.sCI,
n,exact.sCI=(n<100),exact.pCI=(n<100),
nosym.pCI = FALSE, with.legend = TRUE,
legend.bg = "white", legend.pos = "topleft",
legend.cex = 0.8, legend.pref = "", legend.postf = "",
- legend.alpha = alpha, qqb0 = NULL, debug = FALSE)
+ legend.alpha = alpha, qqb0 = NULL, transf0=NULL, debug = FALSE)
.deleteItemsMCL(mcl)
.distrExInstalled
@@ -94,6 +94,7 @@
\item{legend.alpha}{nominal coverage probability}
\item{mcl}{arguments in call as a list}
\item{qqb0}{precomputed return value of \code{qqbounds}}
+\item{transf0}{optional transformation of x-values (by default \code{NULL} and then ignored)}
\item{debug}{logical; if \code{TRUE} additional output to debug confidence bounds. }
\item{silent0}{logical; it is used as argument \code{silent} in \code{try}-catches
within this function. }
Modified: branches/distr-2.6/pkg/distrMod/R/0distrModUtils.R
===================================================================
--- branches/distr-2.6/pkg/distrMod/R/0distrModUtils.R 2014-10-05 21:06:31 UTC (rev 969)
+++ branches/distr-2.6/pkg/distrMod/R/0distrModUtils.R 2014-10-06 14:25:00 UTC (rev 970)
@@ -454,7 +454,7 @@
with.legend = TRUE, legend.bg = "white",
legend.pos = "topleft", legend.cex = 0.8,
legend.pref = "", legend.postf = "",
- legend.alpha = alpha, qqb0=NULL, debug = FALSE){
+ legend.alpha = alpha, qqb0=NULL, transf0=NULL, debug = FALSE){
x <- sort(unique(x))
if("gaps" %in% names(getSlots(class(D))))
@@ -469,6 +469,8 @@
x.in <- x[SI.in]
x.c <- x.in[SI.c]
x.d <- x.in[!SI.c]
+ tx.c <- if(is.null(transf0)) x.c else transf0(x.c)
+ tx.d <- if(is.null(transf0)) x.d else transf0(x.d)
qqb <- if(is.null(qqb0)) qqbounds(x,D,alpha,n,withConf.pw, withConf.sim,
@@ -478,56 +480,56 @@
if(qqb$err["pw"]){
if(sum(SI.c)>0){
- if(datax){
- lines(x.c, qqb$crit[SI.c,"pw.right"],
+ if(!datax){
+ lines(tx.c, qqb$crit[SI.c,"pw.right"],
col=col.pCI,lty=lty.pCI,lwd=lwd.pCI)
- lines(x.c, qqb$crit[SI.c,"pw.left"],
+ lines(tx.c, qqb$crit[SI.c,"pw.left"],
col=col.pCI,lty=lty.pCI,lwd=lwd.pCI)
}else{
- lines(qqb$crit[SI.c,"pw.right"], x.c,
+ lines(qqb$crit[SI.c,"pw.right"], tx.c,
col=col.pCI,lty=lty.pCI,lwd=lwd.pCI)
- lines(qqb$crit[SI.c,"pw.left"], x.c,
+ lines(qqb$crit[SI.c,"pw.left"], tx.c,
col=col.pCI,lty=lty.pCI,lwd=lwd.pCI)
}
}
if(sum(!SI.c)>0){
- if(datax){
- points(x.d, qqb$crit[!SI.c,"pw.right"],
+ if(!datax){
+ points(tx.d, qqb$crit[!SI.c,"pw.right"],
col=col.pCI, pch=pch.pCI, cex = cex.pCI)
- points(x.d, qqb$crit[!SI.c,"pw.left"],
+ points(tx.d, qqb$crit[!SI.c,"pw.left"],
col=col.pCI, pch=pch.pCI, cex = cex.pCI)
}else{
- points(qqb$crit[!SI.c,"pw.right"], x.d,
+ points(qqb$crit[!SI.c,"pw.right"], tx.d,
col=col.pCI, pch=pch.pCI, cex = cex.pCI)
- points(qqb$crit[!SI.c,"pw.left"], x.d,
+ points(qqb$crit[!SI.c,"pw.left"], tx.d,
col=col.pCI, pch=pch.pCI, cex = cex.pCI)
}
}
}
if(qqb$err["sim"]){
if(sum(SI.c)>0){
- if(datax){
- lines(x.c, qqb$crit[SI.c,"sim.right"],
+ if(!datax){
+ lines(tx.c, qqb$crit[SI.c,"sim.right"],
col=col.sCI,lty=lty.sCI,lwd=lwd.sCI)
- lines(x.c, qqb$crit[SI.c,"sim.left"],
+ lines(tx.c, qqb$crit[SI.c,"sim.left"],
col=col.sCI,lty=lty.sCI,lwd=lwd.sCI)
}else{
- lines(qqb$crit[SI.c,"sim.right"], x.c,
+ lines(qqb$crit[SI.c,"sim.right"], tx.c,
col=col.sCI,lty=lty.sCI,lwd=lwd.sCI)
- lines(qqb$crit[SI.c,"sim.left"], x.c,
+ lines(qqb$crit[SI.c,"sim.left"], tx.c,
col=col.sCI,lty=lty.sCI,lwd=lwd.sCI)
}
}
if(sum(!SI.c)>0){
- if(datax){
- points(x.d, qqb$crit[!SI.c,"sim.right"],
+ if(!datax){
+ points(tx.d, qqb$crit[!SI.c,"sim.right"],
col=col.sCI, pch=pch.sCI, cex = cex.sCI)
- points(x.d, qqb$crit[!SI.c,"sim.left"],
+ points(tx.d, qqb$crit[!SI.c,"sim.left"],
col=col.sCI, pch=pch.sCI, cex = cex.sCI)
}else{
- points(qqb$crit[!SI.c,"sim.right"], x.d,
+ points(qqb$crit[!SI.c,"sim.right"], tx.d,
col=col.sCI, pch=pch.sCI, cex = cex.sCI)
- points(qqb$crit[!SI.c,"sim.left"], x.d,
+ points(qqb$crit[!SI.c,"sim.left"], tx.d,
col=col.sCI, pch=pch.sCI, cex = cex.sCI)
}
}
Modified: branches/distr-2.6/pkg/distrMod/R/qqplot.R
===================================================================
--- branches/distr-2.6/pkg/distrMod/R/qqplot.R 2014-10-05 21:06:31 UTC (rev 969)
+++ branches/distr-2.6/pkg/distrMod/R/qqplot.R 2014-10-06 14:25:00 UTC (rev 970)
@@ -42,7 +42,7 @@
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 = TRUE, ### as in qqnorm
+ datax = FALSE, ### as in qqnorm
xlab = deparse(substitute(x)), ## x-label
ylab = deparse(substitute(y)), ## y-label
..., ## further parameters
@@ -211,7 +211,7 @@
}
if(plot.it){
- qqb <- .confqq(xy, y, datax, withConf.pw, withConf.sim, alpha.CI,
+ 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,
Modified: branches/distr-2.6/pkg/distrMod/R/returnlevelplot.R
===================================================================
--- branches/distr-2.6/pkg/distrMod/R/returnlevelplot.R 2014-10-05 21:06:31 UTC (rev 969)
+++ branches/distr-2.6/pkg/distrMod/R/returnlevelplot.R 2014-10-06 14:25:00 UTC (rev 970)
@@ -16,6 +16,7 @@
datax = FALSE, ### as in qqnorm
MaxOrPOT = c("Max","POT"), ### used for block maxima or points over threshold
npy = 365, ### number of observations per year
+ threshold = if(is(y,"GPareto")) NA else 0,
xlab = deparse(substitute(x)), ## x-label
ylab = deparse(substitute(y)), ## y-label
main = "",
@@ -85,6 +86,12 @@
mcl$added.points.CI <- NULL
force(x)
+ thresh0 <- threshold
+ if(is(y,"GPareto")){
+ if(is.na(threshold)) thresh0 <- location(y)
+ y <- y - thresh0
+ x <- x + thresh0
+ }
xj <- sort(x)
if(any(.isReplicated(x))&&jit.fac>0)
@@ -237,9 +244,9 @@
}
}
- qqb <- qqbounds(sort(unique(xy)),y,alpha.CI,n,withConf.pw, withConf.sim,
- exact.sCI,exact.pCI,nosym.pCI, debug = debug)
- qqb$crit <- p2rl(qqb$crit)
+ #qqb <- qqbounds(sort(unique(xy)),y,alpha.CI,n,withConf.pw, withConf.sim,
+ # exact.sCI,exact.pCI,nosym.pCI, debug = debug)
+ #qqb$crit <- p2rl(qqb$crit)
if(plot.it){
qqb <- .confqq(xy, y, datax, withConf.pw, withConf.sim, alpha.CI,
col.pCI, lty.pCI, lwd.pCI, pch.pCI, cex.pCI,
@@ -249,7 +256,7 @@
legend.bg = legend.bg, legend.pos = legend.pos,
legend.cex = legend.cex, legend.pref = legend.pref,
legend.postf = legend.postf, legend.alpha = legend.alpha,
- qqb0=qqb, debug = debug)
+ qqb0=NULL, transf0=p2rl, debug = debug)
}
}}
return(invisible(c(ret,qqb)))
Modified: branches/distr-2.6/pkg/distrMod/man/qqplot.Rd
===================================================================
--- branches/distr-2.6/pkg/distrMod/man/qqplot.Rd 2014-10-05 21:06:31 UTC (rev 969)
+++ branches/distr-2.6/pkg/distrMod/man/qqplot.Rd 2014-10-06 14:25:00 UTC (rev 970)
@@ -6,7 +6,7 @@
\S4method{qqplot}{ANY,UnivariateDistribution}(x,y,
n = length(x), withIdLine = TRUE,
withConf = TRUE, withConf.pw = withConf, withConf.sim = withConf,
- plot.it = TRUE, datax = TRUE, xlab = deparse(substitute(x)),
+ plot.it = TRUE, datax = FALSE, xlab = deparse(substitute(x)),
ylab = deparse(substitute(y)),
..., width = 10, height = 5.5, withSweave = getdistrOption("withSweave"),
mfColRow = TRUE, n.CI = n, withLab = FALSE, lab.pts = NULL, which.lbs = NULL,
@@ -188,8 +188,9 @@
set.seed(123)
x <- rnorm(40,mean=15,sd=30)
qqplot(x, Chisq(df=15))
-qqplot(x, NormLocationScaleFamily())
-mlE <- MLEstimator(x, NormLocationScaleFamily())
+NF <- NormLocationScaleFamily(mean=15, sd=30)
+qqplot(x, NF)
+mlE <- MLEstimator(x, NF)
qqplot(x, mlE)
}
\keyword{hplot}
Modified: branches/distr-2.6/pkg/distrMod/man/returnlevelplot.Rd
===================================================================
--- branches/distr-2.6/pkg/distrMod/man/returnlevelplot.Rd 2014-10-05 21:06:31 UTC (rev 969)
+++ branches/distr-2.6/pkg/distrMod/man/returnlevelplot.Rd 2014-10-06 14:25:00 UTC (rev 970)
@@ -7,6 +7,7 @@
n = length(x), withIdLine = TRUE,
withConf = TRUE, withConf.pw = withConf, withConf.sim = withConf,
plot.it = TRUE, datax = FALSE, MaxOrPOT = c("Max","POT"), npy = 365,
+ threshold = if(is(y,"GPareto")) NA else 0,
xlab = deparse(substitute(x)),
ylab = deparse(substitute(y)),
main = "",
@@ -62,6 +63,9 @@
must be one of "Max" (default) or "POT".
You can specify just the initial letter.}
\item{npy}{number of observations per year/block.}
+\item{threshold}{numerical; in case of \code{MaxOrPot=="POT"}, this captures
+ the (removed) threshold. If it is \code{NA}, it is reconstructed
+ from the distribution \code{y}.}
\item{main}{Main title}
\item{xlab}{x-label}
\item{ylab}{y-label}
@@ -198,7 +202,23 @@
}
\examples{
returnlevelplot(r(Norm(15,sqrt(30)))(40), Chisq(df=15))
-returnlevelplot(r(Norm(15,sqrt(30)))(40), NormLocationFamily())
+### more could be seen after installing RobExtremes and ismev
+#
+# require(RobExtremes)
+# require(ismev)
+#
+# data(portpirie)
+# ppfit <- gev.fit(portpirie[,2])
+# DD2 <- GEVFamily(scale=ppfit$mle[2],shape=ppfit$mle[3],loc=ppfit$mle[1])
+# erg <- returnlevelplot(portpirie[,2], DD2, datax=FALSE)
+# erg <- returnlevelplot(portpirie[,2], DD2, datax=TRUE)
+#
+# data(rain)
+# opfit <- gpd.fit(rain,10)
+# DD1 <- GParetoFamily(scale=opfit$mle[1],shape=opfit$mle[2],loc=10)
+# erg <- returnlevelplot(rain, DD1, datax=FALSE, MaxOrPOT="POT", npy=365, xlim=c(1e-1,1e3))
+# erg <- returnlevelplot(rain, DD1, datax=TRUE, MaxOrPOT="POT", npy=365, xlim=c(1e-1,1e3))
+
}
\keyword{hplot}
\keyword{distribution}
More information about the Distr-commits
mailing list