[Distr-commits] r951 - in branches/distr-2.6/pkg/distr: R man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Sun Aug 10 15:02:30 CEST 2014
Author: ruckdeschel
Date: 2014-08-10 15:02:30 +0200 (Sun, 10 Aug 2014)
New Revision: 951
Modified:
branches/distr-2.6/pkg/distr/R/internals-qqplot.R
branches/distr-2.6/pkg/distr/R/qqbounds.R
branches/distr-2.6/pkg/distr/man/internals-qqplot.Rd
Log:
[distr]
- qqplot:
+ (already there for a while) gains argument 'debug' to be able to trace where computation of confidence bounds fails;
this 'debug' argument has now been enhanced in the sense that now more detailed information about which
helper function was the culprit and where (uniroot/optimize) this happens
- qqbounds:
+ the respective helper functions .BinomCI.nosym, .BinomCI, .q2kolmogorov capsulate their calls to uniroot / optimize in
individual try-catches (instead of / on top of the one in qqbounds itself) so that in case of several evaluations (e.g. in
pointwise CIs) at least a subset of valid points is produced.
+ the respective helper functions .BinomCI.nosym, .BinomCI, .q2kolmogorov, .q2pw give more detailed information when
a try catch happens
+ In addition the search interval is now more flexible: it is bound away from 0 below and from the upper bound from above, and this
bound moves to the respective lower/upper bounds in up to 20 trials.
+ Also, the search interval for .q2kolmogorov if argument 'exact' is TRUE (i.e., if pKolmogorov2x is called) has been refined.
This all should make the simultaneous intervals more stable.
Modified: branches/distr-2.6/pkg/distr/R/internals-qqplot.R
===================================================================
--- branches/distr-2.6/pkg/distr/R/internals-qqplot.R 2014-08-10 12:52:51 UTC (rev 950)
+++ branches/distr-2.6/pkg/distr/R/internals-qqplot.R 2014-08-10 13:02:30 UTC (rev 951)
@@ -96,21 +96,43 @@
}
-.q2kolmogorov <- function(alpha,n,exact=(n<100)){ ## Kolmogorovstat
+.q2kolmogorov <- function(alpha,n,exact=(n<100), silent0 = TRUE){ ## Kolmogorovstat
if(is.numeric(alpha)) alpha <- as.vector(alpha)
else stop("Level alpha must be numeric.")
if(any(is.na(alpha))) stop("Level alpha must not contain missings.")
if(exact){
- fct <- function(p0){
- ### from ks.test from package stats:
- .pk2(p0,n) -alpha
- }
- res <- uniroot(fct,lower=0,upper=1)$root*sqrt(n)
+ fct <- function(p0){
+ ### from ks.test from package stats:
+ .pk2(p0,n) -alpha
+ }
+ i <- 0
+ oK <- FALSE
+ del <- 0.01
+ while(!oK && i < 20){
+ i <- i + 1
+ res <- try(uniroot(fct,lower=del,upper=3*(1-del)/sqrt(n))$root*sqrt(n), silent=silent0)
+ del <- del / 10
+ if(!is(res, "try-error")) oK <- TRUE
+ }
}else{
- fct <- function(p0){
- ### from ks.test from package stats:
- 1 - .pks2(p0,1e-09)-alpha }
- res <- uniroot(fct,lower=1e-12,upper=sqrt(n))$root
+ fct <- function(p0){
+ ### from ks.test from package stats:
+ 1 - .pks2(p0,1e-09)-alpha }
+ i <- 0
+ oK <- FALSE
+ while(!oK && i < 20){
+ i <- i + 1
+ res <- try(uniroot(fct,lower=del,upper=3*(1-del))$root, silent=silent0)
+ del <- del / 10
+ if(!is(res, "try-error")){
+ oK <- TRUE
+ }else{
+ if(!silent0){
+ cat("Problem in uniroot in .q2kolmogorov:\nSearch bounds were")
+ print(c(lower=del,upper=sqrt(n)*(1-del)))
+ }
+ }
+ }
}
return(res)
}
@@ -128,26 +150,87 @@
}
-.BinomCI <- function(x,p.b,D,n,alpha){
+.BinomCI <- function(x,p.b,D,n,alpha,silent0 = TRUE){
if(length(x)==0) return(NA)
- res <- sapply(1:length(x), function(i) uniroot(.BinomCI.in,
- lower=0, upper=sqrt(n)*max(x[i],n-x[i])+1,
- p.bi = p.b[i], x.i = x[i], del.i = 0,
- D.i = D, n.i = n, alpha.i = alpha, tol = 1e-9)$root)
+ res <- sapply(1:length(x),
+ function(i){
+ .nm.i <- sqrt(n)*max(x[i],n-x[i]) + 1
+ .fct.i <- function(t){
+ .BinomCI.in(t, p.bi = p.b[i], x.i = x[i], del.i = 0,
+ D.i = D, n.i = n, alpha.i = alpha)
+ }
+ ii <- 0
+ oK <- FALSE
+ del <- 0.01
+ while(ii < 20 && !oK){
+ ii <- ii + 1
+ res0 <- try(uniroot(.fct.i, tol = 1e-9, lower = del,
+ upper=.nm.i*(1-del))$root,
+ silent=silent0)
+ del <- del / 10
+ if(!is(res0, "try-error")){
+ oK <- TRUE
+ }else{
+ if(!silent0){
+ cat("Problem in uniroot in .BinomCI:\nSearch bounds were")
+ print(c(lower=del,upper=.nm.i*(1-del)))
+ cat("Further settings:\n")
+ print(c(i=i,p.bi = p.b[i], x.i = x[i], del=del))
+ }
+ res0 <- NA
+ }
+ }
+ return(res0)
+ }
+ )
return(cbind(left=-res, right=res))
}
-.BinomCI.nosym <- function(x,p.b,D,n,alpha){
+.BinomCI.nosym <- function(x,p.b,D,n,alpha,silent0 = TRUE){
if(length(x)==0) return(NA)
- res0 <- sapply(1:length(x), function(i){
- get.t <- function(del.o, p.bi, x.i)
- uniroot(.BinomCI.in,
- lower=0, upper=sqrt(n)*max(x.i,n-x.i)+1,
- p.bi = p.bi, x.i = x.i, del.i=del.o,
- D.i = D, n.i = n, alpha.i = alpha, tol = 1e-9)$root
- res <- optimize(get.t, lower=-sqrt(n)*max(x[i],n-x[i])-1,
- upper = sqrt(n)*max(x[i],n-x[i])+1, p.bi = p.b[i],
- x = x[i], tol = 1e-9)
+ res0 <- sapply(1:length(x),
+ function(i){
+ .nm.i <- max(x[i],n-x[i])*sqrt(n) + 1
+ .fct.i.d <- function(t,del.0){
+ .BinomCI.in(t, p.bi = p.b[i], x.i = x[i], del.i = del.0,
+ D.i = D, n.i = n, alpha.i = alpha)
+ }
+ get.t <- function(del.o){
+ ii <- 0
+ oK <- FALSE
+ del.l <- 0.01
+ while(ii < 20 && !oK){
+ ii <- ii + 1
+ res0 <- try(uniroot(.fct.i.d, tol = 1e-9, lower = del.l,
+ upper = .nm.i*(1-del.l),
+ del.0 = del.o)$root, silent = silent0)
+ del.l <- del.l / 10
+ if(!is(res0, "try-error")){
+ oK <- TRUE
+ }else{
+ if(!silent0){
+ cat("Problem in uniroot in .BinomCI.nosym:\nSearch bounds were")
+ print(c(lower=del.l,upper=.nm.i*(1-del.l)))
+ cat("Further settings:\n")
+ print(c(i=i,p.bi = p.b[i], x.i = x[i], del.i=del.o,
+ del.l=del.l))
+ }
+ res0 <- NA
+ }
+ }
+ return(res0)
+ }
+ res <- try(optimize(get.t, lower=-.nm.i, upper = .nm.i, tol = 1e-9),
+ silent = silent0)
+ if(is(res,"try-error")){
+ if(!silent0){
+ cat("Problem in optimize in .BinomCI.nosym:\nSearch bounds were")
+ print(c(lower=-.nm.i,upper=.nm.i))
+ cat("Further settings:\n")
+ print(c(i=i,p.bi = p.b[i], x.i = x[i]))
+ }
+ return(c(NA,NA))
+ }
t.o <- res$objective
del <- res$minimum
c(left=-t.o+del, right=t.o+del)
@@ -156,10 +239,10 @@
}
-.q2pw <- function(x,p.b,D,n,alpha,exact=(n<100),nosym=FALSE){
+.q2pw <- function(x,p.b,D,n,alpha,exact=(n<100),nosym=FALSE, silent0=TRUE){
if(exact){
fct <- if(nosym) .BinomCI.nosym else .BinomCI
- ro <- fct(x,p.b,D,n,alpha)
+ ro <- fct(x,p.b,D,n,alpha,silent0)
return(ro)
}
pq <- log(p.b)+log(1-p.b)
Modified: branches/distr-2.6/pkg/distr/R/qqbounds.R
===================================================================
--- branches/distr-2.6/pkg/distr/R/qqbounds.R 2014-08-10 12:52:51 UTC (rev 950)
+++ branches/distr-2.6/pkg/distr/R/qqbounds.R 2014-08-10 13:02:30 UTC (rev 951)
@@ -25,8 +25,13 @@
print(l.x)
print(c(alpha,n,exact.sCI))
}
- c.crit <- if(withConf.sim) try(.q2kolmogorov(alpha,n,exact.sCI), silent=TRUE) else NULL
- c.crit.i <- if(withConf.pw) try(.q2pw(x.in,p.r,D,n,alpha,exact.pCI,nosym.pCI),silent=TRUE) else NULL
+
+ silent0 <- !debug
+
+ c.crit <- if(withConf.sim) try(.q2kolmogorov(alpha,n,exact.sCI, silent0),
+ silent=silent0) else NULL
+ c.crit.i <- if(withConf.pw) try(.q2pw(x.in,p.r,D,n,alpha,exact.pCI,nosym.pCI,
+ silent0),silent=silent0) else NULL
#print(cbind(c.crit,c.crit.i))
if(debug){
print(str(c.crit))
Modified: branches/distr-2.6/pkg/distr/man/internals-qqplot.Rd
===================================================================
--- branches/distr-2.6/pkg/distr/man/internals-qqplot.Rd 2014-08-10 12:52:51 UTC (rev 950)
+++ branches/distr-2.6/pkg/distr/man/internals-qqplot.Rd 2014-08-10 13:02:30 UTC (rev 951)
@@ -27,11 +27,11 @@
.makeLenAndOrder(x,ord)
.BinomCI.in(t,p.bi,x.i, del.i=0,D.i,n.i,alpha.i)
-.BinomCI(x,p.b,D,n,alpha)
-.BinomCI.nosym(x,p.b,D,n,alpha)
+.BinomCI(x,p.b,D,n,alpha, silent0 = TRUE)
+.BinomCI.nosym(x,p.b,D,n,alpha, silent0 = TRUE)
-.q2kolmogorov(alpha,n,exact=(n<100))
-.q2pw(x,p.b,D,n,alpha,exact=(n<100),nosym=FALSE)
+.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,
col.pCI, lty.pCI, lwd.pCI, pch.pCI, cex.pCI,
@@ -94,7 +94,9 @@
\item{legend.alpha}{nominal coverage probability}
\item{mcl}{arguments in call as a list}
\item{qqb0}{precomputed return value of \code{qqbounds}}
-\item{debug}{logical; if \code{TRUE} additional output to debug confidence bounds.}
+\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. }
}
\details{
More information about the Distr-commits
mailing list