1.46 - incremental update

- man pages for dbetabinom and dgampois
- added rbetabinom function
- bug fix: map detects linear models like "mu <- 0" now
- map2stan gamma-Poisson template uses Stan's neg_binomial_2
distribution now (so less parameter transformation required)
- postcheck displays beta-binomial like binomial now

DESCRIPTION

@@ -1,8 +1,8 @@
 Package: rethinking
 Type: Package
 Title: Statistical Rethinking book package
-Version: 1.45
-Date: 2014-10-20
+Version: 1.46
+Date: 2014-11-3
 Author: Richard McElreath
 Maintainer: Richard McElreath <[email protected]>
 Depends:

R/distributions.r

@@ -92,6 +92,18 @@ dbetabinom <- function (x, size, prob, theta, shape1, shape2, log = FALSE)
     else exp(v)
 }
 
+rbetabinom <- function( n , size, prob, theta, shape1, shape2 ) {
+    if (missing(prob) && !missing(shape1) && !missing(shape2)) {
+        prob <- shape1/(shape1 + shape2)
+        theta <- shape1 + shape2
+    }
+    # first generate beta probs
+    p <- rbeta2( n , prob , theta )
+    # then generate binomial counts
+    y <- rbinom( n , size , p )
+    return(y)
+}
+
 # gamma-poisson functions
 
 dgamma2 <- function( x , mu , scale , log=FALSE ) {

R/map.r

@@ -127,7 +127,10 @@ map <- function( flist , data , start , method="BFGS" , hessian=TRUE , debug=FAL
         RHS <- f[[3]]
         LHS <- f[[2]]
         fname <- as.character( RHS[[1]] )
-        if ( fname=="+" | fname=="*" | fname=="-" | fname=="/" | fname=="%*%" | fname %in% invlink.names ) {
+        flag_monad_linear_model <- FALSE
+        if ( length(RHS)==1 & class(RHS[[1]])=="numeric" )
+            flag_monad_linear_model <- TRUE
+        if ( fname=="+" | fname=="*" | fname=="-" | fname=="/" | fname=="%*%" | fname %in% invlink.names | flag_monad_linear_model==TRUE ) {
             # linear model formula with no density (but maybe invlink) function
             # return a list with parameter name in [[1]] and text of RHS in [[2]]
             thetext <- list( as.character(LHS) , paste( deparse(RHS) , collapse=" " ) )

R/map2stan-templates.r

@@ -516,7 +516,7 @@ map2stan.templates <- list(
     GammaPoisson = list(
         name = "GammaPoisson",
         R_name = "dgampois",
-        stan_name = "neg_binomial",
+        stan_name = "neg_binomial_2",
         num_pars = 2,
         par_names = c("alpha","beta"),
         par_bounds = c("<lower=0>","<lower=0>"),
@@ -525,8 +525,8 @@ map2stan.templates <- list(
         par_map = function(k,...) {
             mu_name <- k[[1]];
             scale_name <- k[[2]];
-            k[[1]] <- concat(mu_name,"/",scale_name);
-            k[[2]] <- concat("1/",scale_name);
+            k[[1]] <- concat(mu_name);
+            k[[2]] <- concat(scale_name);
             return(k);
         },
         vectorized = TRUE

R/sim.r

@@ -274,7 +274,7 @@ postcheck <- function( fit , prob=0.9 , window=20 , n=1000 , col=rangi2 , ... )
 
     # check for aggregated binomial context
     mumax <- max(c(as.numeric(mu.PI)))
-    if ( ymax > 1 & mumax <= 1 & dname=="dbinom" ) {
+    if ( ymax > 1 & mumax <= 1 & dname %in% c("dbinom","dbetabinom") ) {
         # probably aggregated binomial
         size_var <- as.character(lik[[3]][[2]])
         size_var <- fit@data[[ size_var ]]

man/dbetabinom.Rd

@@ -0,0 +1,64 @@
+\name{dbetabinom}
+\alias{dbetabinom}\alias{rbetabinom}
+\title{Beta-binomial probability density}
+\description{
+  Functions for computing density and producing random samples from a beta-binomial probability distribution.
+}
+\usage{
+dbetabinom( x , size , prob , theta , shape1 , shape2 , log=FALSE )
+rbetabinom( n , size , prob , theta , shape2 , shape2 )
+}
+%- maybe also 'usage' for other objects documented here.
+\arguments{
+  \item{x}{Integer values to compute probabilies of}
+  \item{size}{Number of trials}
+  \item{prob}{Average probability of beta distribution}
+  \item{theta}{Dispersion of beta distribution}
+  \item{shape1}{First shape parameter of beta distribution (alpha)}
+  \item{shape2}{Second shape parameter of beta distribution (beta)}
+  \item{log}{If \code{TRUE}, returns log-probability instead of probability}
+  \item{n}{Number of random observations to sample}
+}
+\details{
+  These functions provide density and random number calculations for beta-binomial observations. The \code{dbetabinom} code is based on Ben Bolker's original code in the \code{emdbook} package. 
+  
+  Either \code{prob} and \code{theta} OR \code{shape1} and \code{shape2} must be provided. The two parameterizations are related by shape1 = prob * theta, shape2 = (1-prob) * theta.
+  
+  The \code{rbetabinom} function generates random beta-binomial observations by using both \code{\link{rbeta}} and \code{\link{rbinom}}. It draws \code{n} values from a beta distribution. Then for each, it generates a random binomial observation.
+}
+\references{}
+\author{Richard McElreath}
+\seealso{}
+\examples{
+data(reedfrogs)
+reedfrogs$pred_yes <- ifelse( reedfrogs$pred=="pred" , 1 , 0 )
+
+# map model fit
+# note exp(log_theta) to constrain theta to positive reals
+m <- map(
+    alist(
+        surv ~ dbetabinom( density , p , exp(log_theta) ),
+        logit(p) <- a + b*pred_yes,
+        a ~ dnorm(0,10),
+        b ~ dnorm(0,1),
+        log_theta ~ dnorm(1,10)
+    ),
+    data=reedfrogs )
+
+# map2stan model fit
+# constraint on theta is passed via contraints list
+m.stan <- map2stan(
+    alist(
+        surv ~ dbetabinom( density , p , theta ),
+        logit(p) <- a + b*pred_yes,
+        a ~ dnorm(0,10),
+        b ~ dnorm(0,1),
+        theta ~ dcauchy(0,1)
+    ),
+    data=reedfrogs,
+    constraints=list(theta="lower=0") )
+}
+% Add one or more standard keywords, see file 'KEYWORDS' in the
+% R documentation directory.
+\keyword{ }
+

man/dgampois.Rd

@@ -0,0 +1,57 @@
+\name{dgampois}
+\alias{dgampois}\alias{rgampois}
+\title{Gamma-Poisson probability density}
+\description{
+  Functions for computing density and producing random samples from a gamma-Poisson (negative-binomial) probability distribution.
+}
+\usage{
+dgampois( x , mu , scale , log=FALSE )
+rgampois( n , mu , scale )
+}
+%- maybe also 'usage' for other objects documented here.
+\arguments{
+  \item{x}{Integer values to compute probabilies of}
+  \item{mu}{Mean of gamma distribution}
+  \item{scale}{Scale parameter of gamma distribution}
+  \item{log}{If \code{TRUE}, returns log-probability instead of probability}
+  \item{n}{Number of random observations to sample}
+}
+\details{
+  These functions provide density and random number calculations for gamma-Poisson observations. These functions use \code{dnbinom} and \code{rnbinom} internally, but convert the parameters from the \code{mu} and \code{scale} form. The \code{size} parameter of the negative-binomial is defined by \code{mu/scale}, and the \code{prob} parameter of the negative-binomial is the same as \code{1/(1+scale)}.
+}
+\references{}
+\author{Richard McElreath}
+\seealso{}
+\examples{
+data(Hurricanes)
+
+# map model fit
+# note exp(log_scale) to constrain scale to positive reals
+m <- map(
+    alist(
+        deaths ~ dgampois( mu , exp(log_scale) ),
+        log(mu) <- a + b*femininity,
+        a ~ dnorm(0,100),
+        b ~ dnorm(0,1),
+        log_scale ~ dnorm(1,10)
+    ),
+    data=Hurricanes )
+
+# map2stan model fit
+# constraint on scale is passed via contraints list
+m.stan <- map2stan(
+    alist(
+        deaths ~ dgampois( mu , scale ),
+        log(mu) <- a + b*femininity,
+        a ~ dnorm(0,100),
+        b ~ dnorm(0,1),
+        scale ~ dcauchy(0,2)
+    ),
+    data=Hurricanes,
+    constraints=list(scale="lower=0"),
+    start=list(scale=2) )
+}
+% Add one or more standard keywords, see file 'KEYWORDS' in the
+% R documentation directory.
+\keyword{ }
+