From 0f6fcc90767a9f2a81fec954679fda6769dedc58 Mon Sep 17 00:00:00 2001 From: Richard McElreath Date: Wed, 3 Dec 2014 16:50:02 -0800 Subject: [PATCH] 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 | 4 +-- R/distributions.r | 12 ++++++++ R/map.r | 5 +++- R/map2stan-templates.r | 6 ++-- R/sim.r | 2 +- man/dbetabinom.Rd | 64 ++++++++++++++++++++++++++++++++++++++++++ man/dgampois.Rd | 57 +++++++++++++++++++++++++++++++++++++ 7 files changed, 143 insertions(+), 7 deletions(-) create mode 100644 man/dbetabinom.Rd create mode 100644 man/dgampois.Rd diff --git a/DESCRIPTION b/DESCRIPTION index 8524638..8686538 100644 --- a/DESCRIPTION +++ b/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 Depends: diff --git a/R/distributions.r b/R/distributions.r index 7cb1852..9080832 100644 --- a/R/distributions.r +++ b/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 ) { diff --git a/R/map.r b/R/map.r index 17288e1..d2b8270 100644 --- a/R/map.r +++ b/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=" " ) ) diff --git a/R/map2stan-templates.r b/R/map2stan-templates.r index ed5977f..693fb39 100644 --- a/R/map2stan-templates.r +++ b/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("",""), @@ -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 diff --git a/R/sim.r b/R/sim.r index 73ef072..6239586 100644 --- a/R/sim.r +++ b/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 ]] diff --git a/man/dbetabinom.Rd b/man/dbetabinom.Rd new file mode 100644 index 0000000..c9b8934 --- /dev/null +++ b/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{ } + diff --git a/man/dgampois.Rd b/man/dgampois.Rd new file mode 100644 index 0000000..b0276d1 --- /dev/null +++ b/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{ } +