[Distr-commits] r588 - branches/distr-2.2/pkg
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Fri Sep 11 18:53:14 CEST 2009
Author: ruckdeschel
Date: 2009-09-11 18:53:13 +0200 (Fri, 11 Sep 2009)
New Revision: 588
Modified:
branches/distr-2.2/pkg/qqplot.R
Log:
qqplot.R : asymmetrische, exakte (k?\195?\188rzeste) Konfidenz-Intervalle eingeschlossen
Modified: branches/distr-2.2/pkg/qqplot.R
===================================================================
--- branches/distr-2.2/pkg/qqplot.R 2009-09-10 19:58:15 UTC (rev 587)
+++ branches/distr-2.2/pkg/qqplot.R 2009-09-11 16:53:13 UTC (rev 588)
@@ -72,29 +72,57 @@
return(res)
}
-.BinomCI <- function(x,p.b,D,n,alpha){
- if(length(x)==0) return(NA)
- fct <- function(t,p.bi,x.i){
- p.bi.u <- p(D)(x.i+t/sqrt(n))
- p.bi.l <- p.l(D)(x.i-t/sqrt(n))
- d.r <- if(n*p.bi>floor(n*p.bi)) 0 else
- dbinom(x = n*p.bi, size = n, prob = pmax(p.bi.u,1))
- p.r <- pbinom(q = n*p.bi, size = n, prob = pmin(p.bi.u,1),lower.tail=FALSE)+d.r
- p.l <- pbinom(q = n*p.bi, size = n, prob = pmax(p.bi.l,0),lower.tail=FALSE)
- r <- p.r -p.l - alpha
+.BinomCI.in <- function(t,p.bi,x.i, del.i=0,D.i,n.i,alpha.i){
+ p.bi.u <- p(D.i)(x.i+(t+del.i)/sqrt(n.i))
+ p.bi.l <- p.l(D.i)(x.i-(t-del.i)/sqrt(n.i))
+ d.r <- if(n.i*p.bi>floor(n.i*p.bi)) 0 else
+ dbinom(x = n.i*p.bi, size = n.i, prob = pmax(p.bi.u,1))
+ p.r <- pbinom(q = n.i*p.bi, size = n.i, prob = pmin(p.bi.u,1),lower.tail=FALSE)+d.r
+ p.l <- pbinom(q = n.i*p.bi, size = n.i, prob = pmax(p.bi.l,0),lower.tail=FALSE)
+ r <- p.r -p.l - alpha.i
# print(c(r=r,p=p.bi,x=x.i,p.u=p.bi.u,p.l=p.bi.l,r.r=p.r,r.l=p.l,t=t,np=n*p.bi))
r
}
- sapply(1:length(x), function(i) uniroot(fct, lower=0, upper=sqrt(n)*max(x[i],n-x[i])+1,
- p.bi = p.b[i], x.i = x[i], tol = 1e-9)$root)
+
+
+.BinomCI <- function(x,p.b,D,n,alpha){
+ 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)
+ return(cbind(left=-res, right=res))
}
-.q2pw <- function(x,p.b,D,n,alpha,exact=(n<100)){
- if(exact)
- return(.BinomCI(x,p.b,D,n,alpha))
+.BinomCI.nosym <- function(x,p.b,D,n,alpha){
+ 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)
+ t.o <- res$objective
+ del <- res$minimum
+ c(left=-t.o+del, right=t.o+del)
+ })
+ return(t(res0))
+}
+
+
+.q2pw <- function(x,p.b,D,n,alpha,exact=(n<100),nosym=FALSE){
+ if(exact){
+ fct <- if(nosym) .BinomCI.nosym else .BinomCI
+ ro <- fct(x,p.b,D,n,alpha)
+ return(ro)
+ }
pq <- log(p.b)+log(1-p.b)
if(is(D,"AbscontDistribution")){
dp <- d(D)(x,log=TRUE)
+ dsupp.p <- dsupp.m<-1
}else{
supp.ind <- sapply(x, function(y)
which.min(abs(y-support(D))))
@@ -108,24 +136,22 @@
# print(c(pq[1:3], x[1:3],dsupp.p[1:3],dsupp.m[1:3],m,s))
dp <- log(pnorm((x+dsupp.p/2-m)/s) - pnorm((x-dsupp.m/2-m)/s))
}
- return(exp(pq/2-dp)*(dsupp.p+dsupp.m)/2*qnorm((1+alpha)/2))
+ ro <- exp(pq/2-dp)*(dsupp.p+dsupp.m)/2*qnorm((1+alpha)/2)
+ return(cbind(left=-ro,right=ro))
}
.confqq <- function(x,D,alpha,col.pCI,lty.pCI,lwd.pCI,col.sCI,lty.sCI,lwd.sCI,
- n,exact.sCI=(n<100),exact.pCI=(n<100)){
+ n,exact.sCI=(n<100),exact.pCI=(n<100), nosym.pCI = FALSE){
p.r <- p(D)(x)
p.l <- p.l(D)(x)
- c.crit <- try(#q.r(D)(p(D)(
- .q2kolmogorov(alpha,n,exact.sCI), silent=TRUE)
- #))
- c.crit.i <- try( #q.r(D)(p(D)(
- .q2pw(x,p.r,D,n,alpha,exact.pCI),silent=TRUE)
- #))
+ c.crit <- try(.q2kolmogorov(alpha,n,exact.sCI), silent=TRUE)
+ c.crit.i <- #try(
+ .q2pw(x,p.r,D,n,alpha,exact.pCI,nosym.pCI)#,silent=TRUE)
# print(c.crit.i)
if(!is(c.crit.i, "try-error")){
- lines(x, x+c.crit.i/sqrt(n),
+ lines(x, x+c.crit.i[,"right"]/sqrt(n),
col=col.pCI,lty=lty.pCI,lwd=lwd.pCI)
- lines(x, x-c.crit.i/sqrt(n),
+ lines(x, x+c.crit.i[,"left"]/sqrt(n),
col=col.pCI,lty=lty.pCI,lwd=lwd.pCI)
}
if(!is(c.crit, "try-error")){
@@ -138,9 +164,10 @@
}
if(! is(c.crit,"try-error") || ! is(c.crit.i,"try-error") ){
expression1 <- substitute(
- "symm. pointw."~ex.p~alpha==alpha0~"%- conf. interval",
+ nosym0~"pointw."~ex.p~alpha==alpha0~"%- conf. interval",
list(ex.p = if(exact.pCI) "exact" else "asympt.",
- alpha0 = alpha*100))
+ alpha0 = alpha*100,
+ nosym0 = if(nosym.pCI) "shortest asymm." else "symm"))
expression2 <- substitute(
"simult."~ex.s~alpha==alpha0~"%- conf. interval",
list(ex.s = if(exact.sCI) "exact" else "asympt.",
@@ -179,6 +206,7 @@
mcl$col.lbl <- mcl$cex.lbl <- mcl$adj.lbl <- NULL
mcl$exp.cex2.pch <- mcl$exp.cex2.lbl <- NULL
mcl$exp.fadcol.pch <- mcl$exp.fadcol.lbl <- NULL
+ mcl$nosym.pCI <- NULL
mcl}
## helper into distrMod
@@ -221,10 +249,10 @@
plot.it = TRUE, xlab = deparse(substitute(x)),
ylab = deparse(substitute(y)), ...,
col.IdL = "red", lty.IdL = 2, lwd.IdL = 2,
- alpha.CI = .95, exact.sCI = (n<100), exact.pCI = (n<100),
+ alpha.CI = .95, exact.sCI = (n<100), exact.pCI = (n<100), nosym.pCI = FALSE,
col.pCI = "orange", lty.pCI = 3, lwd.pCI = 2,
col.sCI = "tomato2", lty.sCI = 4, lwd.sCI = 2,
- cex.pch=par("cex"), col.pch = par("col"),
+ cex.pch = par("cex"), col.pch = par("col"),
jit.fac = 0, check.NotInSupport = TRUE,
col.NotInSupport = "red"){
@@ -290,7 +318,8 @@
}
.confqq(xy, y, alpha.CI, col.pCI, lty.pCI, lwd.pCI,
col.sCI, lty.sCI, lwd.sCI,
- length(xc), exact.sCI = exact.sCI, exact.pCI = exact.pCI)
+ length(xc), exact.sCI = exact.sCI, exact.pCI = exact.pCI,
+ nosym.pCI = nosym.pCI)
}
}
return(ret)
@@ -309,7 +338,7 @@
xlab = deparse(substitute(x)), ## x-label
ylab = deparse(substitute(y)), ## y-label
..., ## further parameters
- withLab = TRUE, ## shall observation labels be plotted in
+ withLab = 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
@@ -320,6 +349,7 @@
alpha.CI = .95, ## confidence level
exact.sCI = (n<100), ## shall simultaneous CIs be determined with exact kolmogorov distribution?
exact.pCI = (n<100), ## shall pointwise CIs be determined with exact Binomial 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
@@ -415,7 +445,8 @@
}
.confqq(xy, y, alpha.CI, col.pCI, lty.pCI, lwd.pCI,
col.sCI, lty.sCI, lwd.sCI,
- length(x), exact.sCI = exact.sCI, exact.pCI = exact.pCI)
+ length(x), exact.sCI = exact.sCI, exact.pCI = exact.pCI,
+ nosym.pCI = nosym.pCI)
}
}
return(ret)
More information about the Distr-commits
mailing list