version 0.72

COPYRIGHTS

@@ -0,0 +1,14 @@
+COPYRIGHT STATUS
+----------------
+
+This bulk of this code is
+
+  Copyright (C) 200x Simon Jackman, ... FIXME
+
+The count data regression functionality in R/*hurdle* and R/*zeroinfl* is
+
+  Copyright (C) 2006 Achim Zeileis
+
+All code is subject to the GNU General Public License, Version 2. See
+the file COPYING for the exact conditions under which you may
+redistribute it.

DESCRIPTION

@@ -1,7 +1,7 @@
 Package: pscl
-Version: 0.61
+Version: 0.72
 Type: Package
-Date: 2006-08-15
+Date: 2006-10-24
 Title: Political Science Computational Laboratory, Stanford University
 Author: Simon Jackman <[email protected]>, with contributions from Alex Tahk, Christina Maimone and Jim Fearon
 Maintainer: Simon Jackman <[email protected]>
@@ -11,4 +11,4 @@ Description: Bayesian analysis of item-response theory (IRT) models, roll call a
 SaveImage: yes
 License: GPL version 2 or newer
 URL: http://pscl.stanford.edu
-Packaged: Fri Aug 18 21:28:26 2006; hornik
+Packaged: Wed Oct 25 12:20:26 2006; hornik

HISTORY

@@ -1,3 +1,27 @@
+0.72     * fixed error in bioChemists data found by Bettina Gr�n
+  	 <[email protected]>, variable kids5 was off by 1 unit, now runs
+  	 from min of zero (no kids).
+
+0.71     * fixed bug in betaHPD discovered by John Bullock
+
+0.70     * completely rewritten version of hurdle() and zeroinfl():
+             - new formula interface of type y ~ x | z where y ~ x
+	       specifies the count model and z the inflation/hurdle
+	       regressors.
+	     - re-structured returned value, is now more similar to
+	       "glm" objects
+	     - extended/enhanced extractor functions
+
+0.62    * plot.ideal.1d: better left plot margin, based on max length of
+ 	 legis.name
+	 * plot.ideal.2d: inconsistent testing of presence of beta in ideal object
+ 	   when overlaying cutting planes
+	   * plot.ideal.Rd: more examples (but in \dontrun)
+	 * tracex: bug for 2d ideal objects
+	 * tracex: 2d, make legend lines heavier for showAll
+	 * tracex: for R >= 2.4, change par() to par(no.readonly=TRUE) 
+	 * fixed typo in plot.ideal.Rd
+
 0.61	   * fixed bug in summary.ideal found by Keith Poole (8/8/06)
            * documentation of ideal section on Identification changed
 	   to reflect presence of postProcess function

NAMESPACE

@@ -1,51 +1,74 @@
-useDynLib(pscl)
-
-export(hurdle, odTest,
-       predprob, predprob.glm, predprob.ideal, predprob.zeroinfl,
-       vuong, zeroinfl, ntable, betaHPD) 
-export(densigamma, pigamma, qigamma, rigamma, igammaHDR)
-export(computeMargins,
-       constrain.items, constrain.legis,
-       convertCodes,
-       dropRollCall,
-       dropUnanimous, 
-       extractRollCallObject,
-       ideal, idealToMCMC,
-       readKH,
-       rollcall, summary.rollcall,
-       plot.predict.ideal,
-       plot.ideal,
-       postProcess,
-       tracex,
-       vectorRepresentation)
-
-importFrom(MASS, glm.nb)
-importFrom(stats, logLik)
-
-S3method(coef, zeroinfl)
-S3method(coef, hurdle)
-S3method(dropUnanimous, matrix)
-S3method(dropUnanimous, rollcall)
-S3method(logLik, zeroinfl)
-S3method(logLik, hurdle)
-S3method(plot, ideal)
-S3method(plot, predict.ideal)
-S3method(predict, zeroinfl)
-S3method(predict, ideal)
-S3method(predprob, glm)
-S3method(predprob, zeroinfl)
-S3method(print, zeroinfl)
-S3method(print, hurdle)
-S3method(print, summary.zeroinfl)
-S3method(print, summary.hurdle)
-S3method(print, ideal)
-S3method(print, rollcall)
-S3method(print, predict.ideal)
-S3method(print, summary.ideal)
-S3method(print, summary.rollcall)
-S3method(summary, zeroinfl)
-S3method(summary, hurdle)
-S3method(summary, ideal)
-S3method(summary, rollcall)
+useDynLib("pscl")
 
+export("hurdle", "hurdle.control",
+       "zeroinfl", "zeroinfl.control",
+       "odTest",
+       "predprob", "vuong", "ntable", "betaHPD")
 
+export("densigamma", "pigamma", "qigamma", "rigamma", "igammaHDR")
+
+export("computeMargins",
+       "constrain.items", "constrain.legis",
+       "convertCodes",
+       "dropRollCall",
+       "dropUnanimous", 
+       "extractRollCallObject",
+       "ideal", "idealToMCMC",
+       "readKH",
+       "rollcall", "summary.rollcall",
+       "plot.predict.ideal",
+       "plot.ideal",
+       "postProcess",
+       "tracex",
+       "vectorRepresentation")
+
+importFrom("MASS", "glm.nb")
+importFrom("stats", "logLik")
+
+## methods for class zeroinfl
+  S3method("print", "zeroinfl")
+  S3method("print", "summary.zeroinfl")
+  S3method("summary", "zeroinfl")
+  S3method("coef", "zeroinfl")
+  S3method("vcov", "zeroinfl")
+  S3method("logLik", "zeroinfl")
+  S3method("predict", "zeroinfl")
+  S3method("residuals", "zeroinfl")
+  S3method("fitted", "zeroinfl")
+  S3method("predprob", "zeroinfl")
+  S3method("terms", "zeroinfl")
+  S3method("model.matrix", "zeroinfl")
+
+## methods for class hurdle
+  S3method("print", "hurdle")
+  S3method("print", "summary.hurdle")
+  S3method("summary", "hurdle")
+  S3method("coef", "hurdle")
+  S3method("vcov", "hurdle")
+  S3method("logLik", "hurdle")
+  S3method("predict", "hurdle")
+  S3method("residuals", "hurdle")
+  S3method("fitted", "hurdle")
+  S3method("predprob", "hurdle")
+  S3method("terms", "hurdle")
+  S3method("model.matrix", "hurdle")
+
+## methods for class ideal
+  S3method("plot", "ideal")
+  S3method("plot", "predict.ideal")
+  S3method("predict", "ideal")
+  S3method("predprob", "ideal")
+  S3method("print", "ideal")
+  S3method("print", "predict.ideal")
+  S3method("print", "summary.ideal")
+  S3method("summary", "ideal")
+
+## methods for class rollcall
+  S3method("dropUnanimous", "rollcall")
+  S3method("print", "rollcall")
+  S3method("print", "summary.rollcall")
+  S3method("summary", "rollcall")
+
+## misc methods
+  S3method("dropUnanimous", "matrix")
+  S3method("predprob", "glm")

R/betaHPD.r

@@ -1,50 +1,104 @@
-betaHPD <- function(alpha,beta,p=.95,plot=FALSE){
-
-    if(is.na(p) | is.nan(p) | p > 1 | p < 0)
-        stop("p not between 0 and 1\n")
-
+betaHPD <- function(alpha,beta,p=.95,plot=FALSE,xlim=NULL,debug=FALSE){
+  
+  if(is.na(p) | is.nan(p) | p > 1 | p < 0)
+    stop("p not between 0 and 1\n")
+  
     if(alpha<=1 | beta <=1)
-        stop("betaHPD only implemented for alpha and beta both > 1\n")
+      stop("betaHPD only implemented for alpha and beta both > 1\n")
 
-    func <- function(x0,alpha,beta){
-        y0 <- dbeta(x0,alpha,beta)
-        p0 <- pbeta(x0,alpha,beta)
-        x1 <- qbeta(p0+p,alpha,beta)
-        y1 <- dbeta(x1,alpha,beta)
-        out <- abs(y0-y1)
-        out
-    }
-
-    foo <- try(optimize(f=func,alpha=alpha,beta=beta,
+  ## initialize internal logical flags
+  compute <- TRUE
+  swap <- FALSE
+
+  if(alpha==beta){
+    if(debug)
+      cat("symmetric case, alpha=",alpha,"beta=",beta,"\n")
+    out <- qbeta((1 + c(-1,1)*p)/2,
+                 alpha,beta)
+    compute <- FALSE
+  }
+  if(alpha>beta){
+    swap <- TRUE
+    alphaStar <- beta
+    betaStar <- alpha
+  }
+  else if(beta>alpha){
+    swap <- FALSE
+    alphaStar <- alpha
+    betaStar <- beta
+  }
+  if(debug)
+    cat("swap=",swap,"\n")
+
+  func <- function(x0,alpha,beta){
+    y0 <- dbeta(x0,alpha,beta)
+    p0 <- pbeta(x0,alpha,beta)
+    x1 <- qbeta(p0+p,alpha,beta)
+    y1 <- dbeta(x1,alpha,beta)
+    out <- abs(y0-y1)
+    out
+  }
+
+  if(compute){
+    foo <- try(optimize(f=func,alpha=alphaStar,beta=betaStar,
+                        tol=.Machine$double.eps^(.6),
                         interval=c(.Machine$double.eps,
-                        qbeta(1-p,
-                              alpha,beta))))
+                          qbeta(1-p,
+                                alphaStar,betaStar))))
     if(inherits(foo,"try-error")){
-        warning("optimization in betaHPD failed\n")
-        out <- rep(NA,2)
+      warning("optimization in betaHPD failed\n")
+      out <- rep(NA,2)
     }
     else{
-        out <- c(foo$minimum,
-                 qbeta(pbeta(foo$minimum,alpha,beta)+p,
-                       alpha,beta)
-                 )
-        if(plot){
-            xseq <- seq(min(qbeta(.0001,alpha,beta),out[1]),
-                        max(qbeta(.9999,alpha,beta),out[2]),
-                        length=1000)
-            plot(xseq,dbeta(xseq,alpha,beta),
-                 xlab=expression(theta),
-                 ylab="",
-                 axes=F,
-                 type="n")
-            axis(1)
-            dseq <- seq(out[1],out[2],length=250)
-            fx <- dbeta(dseq,alpha,beta)
-            polygon(x=c(out[1],dseq,rev(dseq)),
-                    y=c(0,fx,rep(0,250)),
-                    border=F,col=gray(.45))
-            lines(xseq,dbeta(xseq,alpha,beta))
-        }
+      if(debug){
+        cat("results of optimization:\n")
+        print(foo)
+      }
+      out <- c(foo$minimum,
+               qbeta(pbeta(foo$minimum,alphaStar,betaStar)+p,
+                     alphaStar,betaStar)
+               )
     }
-    out
+    if(swap){
+      out <- 1-out
+      out <- sort(out)
+      if(debug){
+        cat("swapped back\n")
+        print(out)
+      }
+    }
+  }
+
+  ## plotting
+  if(plot & all(!is.na(out))){
+    xseq <- NULL
+    if(length(xlim)==2 & all(!is.na(xlim))){
+      if(xlim[2]>xlim[1] & xlim[1] >= 0 & xlim[2] <= 1){
+        xseq <- seq(xlim[1]+(.Machine$double.eps^(.25)),
+                    xlim[2]+(.Machine$double.eps^(.25)),
+                    length=1000)
+      }
+    }
+    if(is.null(xseq))
+      xseq <- seq(min(qbeta(.0001,alpha,beta),out[1]),
+                  max(qbeta(.9999,alpha,beta),out[2]),
+                  length=1000)
+
+    plot(xseq,dbeta(xseq,alpha,beta),
+         xlab=expression(theta),
+         ylab="",
+         axes=F,
+         type="n")
+    axis(1)
+
+    ## get polygon for HDR
+    dseq <- seq(out[1],out[2],length=250)
+    fx <- dbeta(dseq,alpha,beta)
+    polygon(x=c(out[1],dseq,rev(dseq)),
+            y=c(0,fx,rep(0,250)),
+            border=F,col=gray(.45))
+    lines(xseq,dbeta(xseq,alpha,beta))
+  }
+
+  out
 }