version 1.0-6

ChangeLog

@@ -1,3 +1,7 @@
+2019-04-29  1.0-6  - Poisson distribution added.
+                   - A bug in unimodalization of binomial acceptability values
+                     (function binom.blaker.acc, option type="unimod") fixed.
+
 2015-08-20  1.0-5  - So called Vos-Hudson adjustment implemented.
                    - DESCRIPTION corrected.
                    - Imports from stats added to NAMESPACE.

DESCRIPTION

@@ -1,14 +1,14 @@
 Package: BlakerCI
 Type: Package
-Title: Blaker's Binomial Confidence Limits
-Version: 1.0-5
-Date: 2015-08-20
+Title: Blaker's Binomial and Poisson Confidence Limits
+Version: 1.0-6
+Date: 2019-04-29
 Author: Jan Klaschka
 Maintainer: Jan Klaschka <klaschka@cs.cas.cz>
-Description: Fast and accurate calculation of Blaker's binomial confidence limits (and some related stuff).
+Description: Fast and accurate calculation of Blaker's binomial and Poisson confidence limits (and some related stuff).
 License: GPL-3
 LazyLoad: yes
 NeedsCompilation: no
-Packaged: 2015-08-20 15:45:10 UTC; klaschka
+Packaged: 2019-04-29 18:13:58 UTC; klaschka
 Repository: CRAN
-Date/Publication: 2015-08-20 19:52:46
+Date/Publication: 2019-04-29 19:00:03 UTC

MD5

@@ -1,16 +1,23 @@
-9677a0469d8f43c4c0bda58d24473b77 *ChangeLog
-aeff792d2d47edf07c21f9f1f67a3a86 *DESCRIPTION
-622d97f1a2f90799930a4cd04c5a3258 *NAMESPACE
+ba649d16cb495868e0fcc3a5cbb6f399 *ChangeLog
+f448f4bf40774b9145f26d143e63f1e6 *DESCRIPTION
+197c1a93443fb80dd93ce47723c8d5b2 *NAMESPACE
 8f41fbda3173ebf87f028f2479cbd70b *R/binom.blaker.VHadj.acc.R
 7bd74de7a2337bd2e96a08a9ba3f7c94 *R/binom.blaker.VHadj.limits.R
 b81d0a0e4af032ee3c7a781fcce33e29 *R/binom.blaker.VHadj.lower.limit.R
-38f9b537d00dba55ae6a2717721bdbdd *R/binom.blaker.acc.R
+3be9cf0131e8dad1f3198f099041c5e4 *R/binom.blaker.acc.R
 e31a8b5ad8ee890e82fa31bc5777b9f4 *R/binom.blaker.acc.single.p.R
 cfa0c40d1cbd95045c3400212df7a4d6 *R/binom.blaker.limits.R
 83a0de528aa33eb2acae2965b561638b *R/binom.blaker.lower.limit.R
-f65871c4eefefd0a087c3baed147f43d *man/BlakerCI-internal.Rd
-dc046f448f61ede62eb7e81f31ac9965 *man/BlakerCI-package.Rd
-c7739d5ea97fd342138aa25528d5dc4f *man/binom.blaker.VHadj.acc.Rd
+130001600877fa2f419725c810bc0ab9 *R/poisson.blaker.acc.R
+ea80067747eec875a629d10942cf6b4d *R/poisson.blaker.acc.single.p.R
+6514c667cf6c4a336b8f3e7533b54e1a *R/poisson.blaker.limits.R
+66fe922523c11077a0b9b0c444b8372c *R/poisson.blaker.lower.limit.R
+1e8d7966dc7c9ca26fccb91a07472099 *R/poisson.blaker.upper.limit.R
+dceec160fbea30438c11b6863b33dab4 *man/BlakerCI-internal.Rd
+1cc3225d3d2be2fb2792d3ebfec5ab3c *man/BlakerCI-package.Rd
+51cb2de00ff7ea7f45cfd5c497e5f8a8 *man/binom.blaker.VHadj.acc.Rd
 158ec2cc80c7eccf4fb4209e7726b8a2 *man/binom.blaker.VHadj.limits.Rd
-f1a7c77709703911bc5b6306381576d0 *man/binom.blaker.acc.Rd
-0b2164dd6cd538fc2278ca83bb9ce32f *man/binom.blaker.limits.Rd
+918aacadaf9475d9ec6a52fb19d178c2 *man/binom.blaker.acc.Rd
+854505d4d71511eee26102097450420f *man/binom.blaker.limits.Rd
+e4b88ec1eb0097bd426465e7280eff26 *man/poisson.blaker.acc.Rd
+f0918798af33103a52da2a855fa257e6 *man/poisson.blaker.limits.Rd

NAMESPACE

@@ -1,3 +1,5 @@
-importFrom("stats", "pbinom", "qbinom", "qbeta", "pbeta")
+importFrom("stats", "pbinom", "qbinom", "qbeta", "pbeta","ppois","qgamma","qpois")
 export(binom.blaker.lower.limit,binom.blaker.limits,binom.blaker.acc.single.p,binom.blaker.acc,
-       binom.blaker.VHadj.lower.limit,binom.blaker.VHadj.limits,binom.blaker.VHadj.acc)
+       binom.blaker.VHadj.lower.limit,binom.blaker.VHadj.limits,binom.blaker.VHadj.acc,
+       poisson.blaker.acc,poisson.blaker.acc.single.p,poisson.blaker.lower.limit,
+       poisson.blaker.upper.limit,poisson.blaker.limits)

R/binom.blaker.acc.R

@@ -26,7 +26,7 @@ binom.blaker.acc <- function(x,n,p,type=c("orig","unimod"),acc.tol=1e-10,...) {
 #     "highlight" the leftmost (below x/n) or rightmost (above x/n)
 #     p element.
       ind <- (ind1 & (q1 > c(-Inf,q1[1:(m-1)]))) |
-             (!ind1 & (q1 > c(q1[2:m],-Inf)))
+             (!ind1 & (q1 < c(q1[2:m],Inf)))
 #     Calculate the "unimodalized" version of the acceptability function
 #     just at the "highlighted" points, and leave the rest
 #     to cummax().

R/poisson.blaker.acc.R

@@ -0,0 +1,45 @@
+poisson.blaker.acc <- function(x,p,type=c("orig","unimod"),acc.tol=1e-10,...) {
+  if (x < 0) stop("Parameter x = ",x, " wrong!")
+  if (acc.tol <= 0) stop("Numerical tolerance ",acc.tol," nonpositive!")
+  type <- match.arg(type)
+  if (type != "orig" && type != "unimod") stop("Unknown type ",type,"!")
+  m <- length(p)
+  if (m < 1) stop("Empty vector p!")
+  if (m < 2) {
+    acc <- poisson.blaker.acc.single.p(x,p,type=type,acc.tol=acc.tol,...)
+    return(acc)
+  }
+  else {
+    if (max(p[2:m]-p[1:(m-1)]) <= 0) stop("Vector p not increasing!")
+#   First, regardless of type ("orig"/"unimod"),
+#   calculate "ordinary" acceptabilities.
+    aq <- sapply(p,poisson.blaker.acc.single.p,x=x,acc.tol=acc.tol,output="both",...=...)
+    acc <- aq[1,]
+    if (type == "orig") {
+      return(acc)
+    }
+    if (type == "unimod") {
+      q1 <- aq[2,]
+      p.hat <- x
+      ind1 <- p <= p.hat
+#     In each interval of continuity of the acceptability function,
+#     "highlight" the leftmost (below x) or rightmost (above x)
+#     p element.
+      ind <- (ind1 & (q1 > c(-Inf,q1[1:(m-1)]))) |
+             (!ind1 & (q1 < c(q1[2:m],Inf)))
+#     Calculate the "unimodalized" version of the acceptability function
+#     just at the "highlighted" points, and leave the rest
+#     to cummax().
+#     (The amount of slow iterative calculations is thus minimized.)
+      acc[ind] <- sapply(p[ind],poisson.blaker.acc.single.p,x=x,type=type,acc.tol=acc.tol,...=...)
+      m1 <- length(which(ind1))
+      if (m1 > 1) {
+        acc[1:m1] <- cummax(acc[1:m1])
+      }
+      if (m1 < m - 1) {
+        acc[m:(m1+1)] <- cummax(acc[m:(m1+1)])
+      }
+      return(acc)
+    }
+  }
+}

R/poisson.blaker.acc.single.p.R

@@ -0,0 +1,92 @@
+poisson.blaker.acc.single.p <- function(x,p,type="orig",acc.tol=1e-10,output="acc",maxiter=100) {
+
+# "Ordinary" acceptability at p.
+if (p <= x) {
+  p1.p <- ppois(x-1,p,lower.tail=FALSE)
+  q1.p <- qpois(p1.p,p)-1
+  a1.p <- min(p1.p + ppois(q1.p,p),1)
+} else {
+  p1.p <- ppois(x,p,lower.tail=TRUE)
+  q1.p <- qpois(1-p1.p,p)+1
+  a1.p <- min(p1.p + ppois(q1.p-1,p,lower.tail=FALSE),1)
+}
+
+# "Unimodalization"
+  if (type == "unimod" && q1.p >= 0) { 
+    if (p <= x) { 
+      upper <- p
+      a1.upp <- a1.p
+      lower <- 0
+      a1.low <- 1   
+#
+      iter <- 0
+#   In 1.0-4, ... >= acc.tol replaced with ... > acc.tol
+      while (a1.low > a1.p && a1.low - a1.upp > acc.tol) { 
+        iter <- iter + 1
+        if (iter > maxiter) {  
+          warning("Convergence not attained after ",maxiter, 
+                             " iterations for x = ",x,", p = ",p,
+                             ", and acceptability tolerance limit of ",acc.tol)
+          break
+        } 
+        mid <- (lower+upper)/2
+        p1.mid <- ppois(x-1,mid,lower.tail=FALSE)
+        p2.mid <- ppois(q1.p,mid)
+        a1.mid <- p1.mid + p2.mid
+        if (p1.mid < p2.mid) { 
+          lower <- mid
+          a1.low <- a1.mid
+        } 
+        else { 
+          upper <- mid
+          a1.upp <- a1.mid
+        } 
+      } 
+    } else {   
+## Unimodalisation for p > x
+      lower <- p
+      a1.low <- a1.p
+      upper <- q1.p 
+      a1.upp <- 1   
+
+      iter <- 0
+#   In 1.0-4, ... 
+      while (a1.upp > a1.p && a1.upp - a1.low > acc.tol) { 
+        iter <- iter + 1
+        if (iter > maxiter) {  
+          warning("Convergence not attained after ",maxiter, 
+                             " iterations for x = ",x,", p = ",p,
+                             ", and acceptability tolerance limit of ",acc.tol)
+          break
+        } 
+        mid <- (lower+upper)/2
+        p1.mid <- ppois(x,mid,lower.tail=TRUE)
+        p2.mid <- ppois(q1.p-1,mid,lower.tail=FALSE)
+        a1.mid <- p1.mid + p2.mid
+        if (p1.mid < p2.mid) { 
+          upper <- mid
+          a1.upp <- a1.mid
+        } 
+        else { 
+          lower <- mid
+          a1.low <- a1.mid
+        } 
+      } 
+
+    } 
+    a1.p <- max(a1.p,a1.low)
+  } 
+
+
+  if (output == "acc") {
+    return(a1.p)
+  }
+  else if (output == "both") {
+    return(c(a1.p,q1.p))
+  }
+  else if (output == "q1") {
+    return(q1.p)
+  }
+
+
+}

R/poisson.blaker.limits.R

@@ -0,0 +1,8 @@
+poisson.blaker.limits <- function(x,level=.95,tol=1e-10,...) {
+  if (x < 0) stop("Parameter x = ",x, " wrong!")
+  if (level <= 0 || level >= 1) stop("Confidence level ",level," out of (0, 1)!")
+  if (tol <= 0) stop("Numerical tolerance ",tol," nonpositive!")
+  lower <- poisson.blaker.lower.limit(x,level,tol,...)
+  upper <- poisson.blaker.upper.limit(x,level,tol,...)
+  return(c(lower,upper))
+}

R/poisson.blaker.lower.limit.R

@@ -0,0 +1,34 @@
+poisson.blaker.lower.limit <- function(x,level,tol=1e-10,maxiter=100) {
+  if (x <=0) return(0)
+  if (x>0) {
+    alpha <- 1-level
+#   Clopper-Pearson limit (CPL)
+    lower <- qgamma(alpha/2,x,1)
+    p1 <- ppois(x-1,lower,lower.tail=FALSE)
+    q1.cp <- qpois(p1,lower)-1
+    upper <- x
+    iter <- 0
+    while (upper-lower >= tol) {
+      iter <- iter+1
+      if (iter > maxiter) {
+        warning("Tolerance limit of ",tol, 
+                             " not attained after ",maxiter, 
+                             " iterations for x = ",x)
+        break
+      }
+      mid <- (lower+upper)/2
+      p1 <- ppois(x-1,mid,lower.tail=FALSE)
+#   Blaker's limit is below the midpoint if either
+#   (i)  acceptability at mid > alpha, or
+#   (ii) acceptability function has a discontinuity between
+#        the midpoint and CPL (first test).
+      if (p1 >= ppois(q1.cp+1,mid) || p1 + ppois(q1.cp,mid) > alpha) {
+        upper <- mid
+      }
+      else {
+        lower <- mid
+      }
+    }
+  return(lower)
+  }
+}  

R/poisson.blaker.upper.limit.R

@@ -0,0 +1,32 @@
+poisson.blaker.upper.limit <- function(x,level,tol=1e-10,maxiter=100) {
+    alpha <- 1-level
+#   Clopper-Pearson limit (CPL)
+    upper <- qgamma(1-alpha/2,x+1,1)
+    p1 <- ppois(x,upper,lower.tail=TRUE)
+    q1.cp <- qpois(1-p1,upper)+1
+    lower <- x
+    iter <- 0
+    while (upper-lower >= tol) {
+      iter <- iter+1
+      if (iter > maxiter) {
+        warning("Tolerance limit of ",tol, 
+                             " not attained after ",maxiter, 
+                             " iterations for x = ",x)
+        break
+      }
+      mid <- (lower+upper)/2
+      p1 <- ppois(x,mid,lower.tail=TRUE)
+#   Blaker's limit is below the midpoint if either
+#   (i)  acceptability at mid > alpha (NEW!! orig: >=), or
+#   (ii) acceptability function has a discontinuity between
+#        the midpoint and CPL (first test).
+      if (p1 >= ppois(q1.cp-2,mid,lower.tail=FALSE) || p1 + ppois(q1.cp-1,mid,lower.tail=FALSE) > alpha) {
+        lower <- mid
+      }
+      else {
+        upper <- mid
+      }
+    }
+  return(upper)
+
+}  

man/BlakerCI-internal.Rd

@@ -3,28 +3,43 @@
 \alias{binom.blaker.lower.limit}
 \alias{binom.blaker.acc.single.p}
 \alias{binom.blaker.VHadj.lower.limit}
+\alias{poisson.blaker.lower.limit}
+\alias{poisson.blaker.upper.limit}
+\alias{poisson.blaker.acc.single.p}
+
 %- Also NEED an '\alias' for EACH other topic documented here.
 \title{
 Internal functions, not expected to be called by the user.
 }
 \description{
+For binomial distribution: 
 Calculation of the lower Blaker's confidence limit
 as defined by Blaker (\code{binom.blaker.lower.limit}),
 or with so called Vos-Hudson adjustment (\code{binom.blaker.VHadj.lower.limit});
 a single acceptability value, optionally \dQuote{unimodalized}
 (\code{binom.blaker.acc.single.p}).
+
+For Poisson distribution:
+Calculation of the lower and upper Blaker's confidence limits
+(\code{poisson.blaker.lower.limit}, \code{poisson.blaker.upper.limit});
+a single acceptability value, optionally \dQuote{unimodalized}
+(\code{poisson.blaker.acc.single.p}).
 }
 \usage{
 binom.blaker.lower.limit(x, n, level, tol = 1e-10, maxiter=100)
 binom.blaker.VHadj.lower.limit(x,n,level,tol=1e-10,maxiter=100,
       nmax=n+1000,int.eps=1e-10)
 binom.blaker.acc.single.p(x, n, p, type = "orig", acc.tol = 1e-10,
       output = "acc", maxiter=100)
+poisson.blaker.lower.limit(x, level, tol = 1e-10, maxiter=100)
+poisson.blaker.upper.limit(x, level, tol = 1e-10, maxiter=100)
+poisson.blaker.acc.single.p(x, p, type = "orig", acc.tol = 1e-10,
+      output = "acc", maxiter=100)
 }
 %- maybe also 'usage' for other objects documented here.
 \arguments{
   \item{x}{
-      number of successes.
+      number of successes (binomial case), or events (Poisson case).
 }
   \item{n}{
       number of trials.
@@ -36,12 +51,12 @@ binom.blaker.acc.single.p(x, n, p, type = "orig", acc.tol = 1e-10,
       numerical tolerance (for the confidence limit).
 }
   \item{p}{
-      point (single value between 0 and 1) where to calculate the acceptability.
+      point (binomial or Poisson parameter value) where to calculate the acceptability.
 }
   \item{type}{
       \code{"orig"}, or \code{"unimod"} -- either unmodified, or unimodalized 
       acceptability \cr
-      (see \code{\link{binom.blaker.acc}}).      
+      (see \code{\link{binom.blaker.acc}}, \code{\link{poisson.blaker.acc}}).      
 }
   \item{acc.tol}{
       numerical tolerance (for the acceptability values when \code{type = "unimod"}). 
@@ -55,9 +70,11 @@ binom.blaker.acc.single.p(x, n, p, type = "orig", acc.tol = 1e-10,
 }
   \item{maxiter}{
       Maximum number of interval halving iterations during the search
-      for a confidence limit (\code{binom.blaker.lower.limit}) or
-      a discontinuity point of the acceptability function 
-      (\code{binom.blaker.acc.single.p} with \code{type = "unimod"}). 
+      for a confidence limit (\code{binom.blaker.lower.limit}, 
+      \code{poisson.blaker.lower.limit}, \code{poisson.blaker.upper.limit}), 
+      or a discontinuity point of the acceptability function 
+      (\code{binom.blaker.acc.single.p}, or \code{poisson.blaker.acc.single.p}  
+      with \code{type = "unimod"}). 
       When the required accuracy is not reached in \code{maxiter} steps -- 
       typically when too small \code{tol} or \code{acc.tol}
       exceeds capabilities of machine arithmetic -- last step's result
@@ -105,7 +122,8 @@ Jan Klaschka \email{klaschka@cs.cas.cz}
 
 \seealso{
 \code{\link{binom.blaker.limits}}, \code{\link{binom.blaker.VHadj.limits}},\cr
-\code{\link{binom.blaker.acc}}, \code{\link{binom.blaker.VHadj.acc}}.
+\code{\link{binom.blaker.acc}}, \code{\link{binom.blaker.VHadj.acc}}, 
+\code{\link{poisson.blaker.limits}}, \code{\link{poisson.blaker.acc}}.
 }
 %\examples{
 %}