1.394 - maintenance

- added data(Dissertations) with its man page
- improved numerical stability of dlaplace
- fixed rlaplace to actually work correctly
- added dgnorm, generalized normal type 1 density
- various plotting functions now set narrow margins automatically via
set_nice_margins()
- map2stan uses new vectorized multi_normal syntax for built Stan code,
making sampling substantially faster in some cases
- fixed bug with precis mixing up matrix parameters in ordering. Now
uses postlistprecis() internally to get things right.

DESCRIPTION

@@ -1,8 +1,8 @@
 Package: rethinking
 Type: Package
 Title: Statistical Rethinking book package
-Version: 1.393
-Date: 2014-05-01
+Version: 1.394
+Date: 2014-08-20
 Author: Richard McElreath
 Maintainer: Richard McElreath <[email protected]>
 Depends:

R/distributions.r

@@ -104,15 +104,16 @@ rgampois <- function( n , mu , scale ) {
 # laplace (double exponential)
 dlaplace <- function(x,location=0,lambda=1,log=FALSE) {
     # f(y) = (1/(2b)) exp( -|y-a|/b )
-    l <- (1/(2*lambda))*exp( -abs(x-location)/lambda )
-    if ( log==TRUE ) l <- log(l)
-    l
+    # l <- (1/(2*lambda))*exp( -abs(x-location)/lambda )
+    ll <- -abs(x-location)/lambda - log(2*lambda)
+    if ( log==FALSE ) ll <- exp(ll)
+    ll
 }
 
 rlaplace <- function(n,location=0,lambda=1) {
-    # generate an exponential sample, then generate a random sign, then add location
-    y <- rexp(n=n,rate=lambda) * sample( c(-1,1) , size=n , replace=TRUE )
-    y <- y + location
+    # generate from difference of two random exponentials with rate 1/lambda
+    y <- rexp(n=n,rate=1/lambda) - rexp(n=n,rate=1/lambda)
+    y <- y + location # offset location
     y
 }
 
@@ -133,7 +134,7 @@ dlkjcorr <- function( x , eta=1 , log=TRUE ) {
 rinvwishart <- function(s,df,Sigma,Prec) {
     #s-by-s Inverse Wishart matrix, df degree of freedom, precision matrix
     #Prec. Distribution of W^{-1} for Wishart W with nu=df+s-1 degree of
-    # freedoom, covar martix Prec^{-1}.
+    # freedom, covar martix Prec^{-1}.
     # NOTE  mean of riwish is proportional to Prec
     if ( missing(Prec) ) Prec <- solve(Sigma)
     if (df<=0) stop ("Inverse Wishart algorithm requires df>0")
@@ -221,3 +222,11 @@ rzibinom <- function(n,p_zero,size,prob) {
 dstudent <- function( x , nu , mu , sigma , log=TRUE ) {
 
 }
+
+# generalized normal
+dgnorm <- function( x , mu , alpha , beta , log=FALSE ) {
+    ll <- log(beta) - log( 2*alpha*gamma(1/beta) ) - (abs(x-mu)/alpha)^beta
+    if ( log==FALSE ) ll <- exp(ll)
+    return(ll)
+}
+

R/map2stan-class.r

@@ -283,10 +283,7 @@ tracerplot <- function( object , col=c("slateblue","orange","red","green") , alp
     n_eff <- summary(object@stanfit)$summary[,'n_eff']
 
     # make window
-    mfrow_old <- par("mfrow")
-    on.exit(par(mfrow = mfrow_old))
-    par(mgp = c(1.5, 0.5, 0), mar = c(2.5, 2.5, 2, 1) + 0.1, 
-            tck = -0.02)
+    set_nice_margins()
     par(mfrow=c(n_rows_per_page,n_cols))
 
     # draw traces

R/map2stan.r

@@ -730,7 +730,11 @@ map2stan <- function( flist , data , start , pars , constraints=list() , types=l
             rhstxt <- paste( klist , collapse=" , " )
 
             # add text to model code
-            m_model_txt <- concat( m_model_txt , indent , "for ( j in 1:" , N_txt , " ) " , lhstxt , "[j] ~ " , vprior$density , "( " , rhstxt , " )" , vprior$T_text , ";\n" )
+            # new column vectorized normal and multi_normal
+            if ( vprior$density %in% c("multi_normal","normal") )
+                m_model_txt <- concat( m_model_txt , indent , lhstxt , " ~ " , vprior$density , "( " , rhstxt , " )" , vprior$T_text , ";\n" )
+            else
+                m_model_txt <- concat( m_model_txt , indent , "for ( j in 1:" , N_txt , " ) " , lhstxt , "[j] ~ " , vprior$density , "( " , rhstxt , " )" , vprior$T_text , ";\n" )
 
             # declare each parameter with correct type from template
             outtype <- "vector"

R/plotting.r

@@ -1,11 +1,18 @@
 # plotting
 
+set_nice_margins <- function() {
+    mfrow_old <- par("mfrow")
+    on.exit(par(mfrow = mfrow_old))
+    par(mgp = c(1.5, 0.5, 0), mar = c(2.5, 2.5, 2, 1) + 0.1, tck = -0.02)
+}
+
 dens <- function( x , adj=0.5 , norm.comp=FALSE , main="" , show.HPDI=FALSE , show.zero=FALSE , rm.na=TRUE , add=FALSE , ...) {
     the.class <- class(x)[1]
     if ( the.class=="data.frame" ) {
         # full posterior
         n <- ncol(x)
         cnames <- colnames(x)
+        set_nice_margins()
         par( mfrow=make.grid(n) )
         for ( i in 1:n ) {
             dens( x[,i] , adj=adj , norm.comp=norm.comp , show.HPDI=show.HPDI , show.zero=TRUE , xlab=cnames[i] , ... )
@@ -14,9 +21,10 @@ dens <- function( x , adj=0.5 , norm.comp=FALSE , main="" , show.HPDI=FALSE , sh
         # vector
         if ( rm.na==TRUE ) x <- x[ !is.na(x) ]
         thed <- density(x,adjust=adj)
-        if ( add==FALSE )
+        if ( add==FALSE ) {
+            set_nice_margins()
             plot( thed , main=main , ... )
-        else
+        } else
             lines( thed$x , thed$y , ... )
         if ( show.HPDI != FALSE ) {
             hpd <- HPDI( x , prob=show.HPDI )
@@ -167,6 +175,7 @@ simplehist <- function( x , ylab="Frequency" , xlab="Count" , ycounts=TRUE , adj
             ylim <- c( 0 , max(freqs) )
         if ( is.null(xlim) )
             xlim <- c(min(bins),max(bins))
+        set_nice_margins()
         plot( 0 ~ 0 , col="white" , xlim=xlim , ylim=ylim , ylab=ylab , xlab=xlab , ... )
         for ( i in 1:length(bins) ) {
             if ( freqs[i] > 0 )
@@ -183,6 +192,7 @@ simplehist <- function( x , ylab="Frequency" , xlab="Count" , ycounts=TRUE , adj
         # continuous
        if ( ylab=="Frequency" ) ylab <- "Density"
        if ( xlab=="Count" ) xlab <- "Value"
+       set_nice_margins()
        plot( density( x , adjust=adjust ) , ylab=ylab , xlab=xlab , ... )
     }
 
@@ -274,6 +284,7 @@ mcmcpairs <- function( posterior , cex=0.3 , pch=16 , col=col.alpha("slateblue",
         cy <- sum(range(y))/2
         text( cx , cy , round(k,2) , cex=2*exp(abs(k))/exp(1) )
     }
+    set_nice_margins()
     pairs( posterior , cex=cex , pch=pch , col=col , upper.panel=panel.2d , lower.panel=panel.cor , diag.panel=panel.dens , ... )
 }
 

R/precis.r

@@ -14,7 +14,7 @@ precis.show <- function( object ) {
 }
 setMethod( "show" , "precis" , function(object) precis.show(object) )
 
-precis.plot <- function( x , y , pars , col.ci="black" , ... ) {
+precis.plot <- function( x , y , pars , col.ci="black" , xlab="Value" , ... ) {
     x <- x@output
     if ( !missing(pars) ) {
         x <- x[pars,]
@@ -23,12 +23,54 @@ precis.plot <- function( x , y , pars , col.ci="black" , ... ) {
     mu <- x[n:1,1]
     left <- x[[3]][n:1]
     right <- x[[4]][n:1]
-    dotchart( mu , labels=rownames(x)[n:1] , xlab="Estimate" , xlim=c(min(left),max(right)) , ... )
+    set_nice_margins()
+    dotchart( mu , labels=rownames(x)[n:1] , xlab=xlab , xlim=c(min(left),max(right)) , ... )
     for ( i in 1:length(mu) ) lines( c(left[i],right[i]) , c(i,i) , lwd=2 , col=col.ci )
     abline( v=0 , lty=1 , col=col.alpha("black",0.15) )
 }
 setMethod( "plot" , "precis" , function(x,y,...) precis.plot(x,y,...) )
 
+# function to process a list of posterior samples from extract.samples into a summary table
+# needed because as.data.frame borks the ordering of matrix parameters like varying effects
+postlistprecis <- function( post , prob=0.95 ) {
+    n_pars <- length(post)
+    result <- data.frame( Mean=0 , StdDev=0 , lower=0 , upper=0 )
+    r <- 1
+    for ( k in 1:n_pars ) {
+        dims <- dim( post[[k]] )
+        if ( length(dims)==1 ) {
+            # single parameter
+            hpd <- as.numeric( HPDI( post[[k]] , prob=prob ) )
+            result[r,] <- c( mean(post[[k]]) , sd(post[[k]]) , hpd[1] , hpd[2] )
+            rownames(result)[r] <- names(post)[k]
+            r <- r + 1
+        }
+        if ( length(dims)==2 ) {
+            # vector of parameters
+            # loop over
+            for ( i in 1:dims[2] ) {
+                hpd <- as.numeric( HPDI( post[[k]][,i] , prob=prob ) )
+                result[r,] <- c( mean(post[[k]][,i]) , sd(post[[k]][,i]) , hpd[1] , hpd[2] )
+                rownames(result)[r] <- concat( names(post)[k] , "[" , i , "]" )
+                r <- r + 1
+            }
+        }
+        if ( length(dims)==3 ) {
+            # matrix of parameters
+            for ( i in 1:dims[2] ) {
+                for ( j in 1:dims[3] ) {
+                    hpd <- as.numeric( HPDI( post[[k]][,i,j] , prob=prob ) )
+                    result[r,] <- c( mean(post[[k]][,i,j]) , sd(post[[k]][,i,j]) , hpd[1] , hpd[2] )
+                    rownames(result)[r] <- concat( names(post)[k] , "[" , i , "," , j , "]" )
+                    r <- r + 1
+                }
+            }
+        }
+    }
+    colnames(result)[3:4] <- c(paste("lower", prob), paste("upper", prob))
+    result
+}
+
 precis <- function( model , depth=1 , pars , type.s=FALSE , ci=TRUE , level=0.95 , corr=FALSE , digits=2 , warn=TRUE ) {
     the.class <- class(model)[1]
     found.class <- FALSE
@@ -58,18 +100,18 @@ precis <- function( model , depth=1 , pars , type.s=FALSE , ci=TRUE , level=0.95
         if ( the.class=="data.frame" ) {
             # HPDI from samples
             ci <- t( apply( model , 2 , HPDI , prob=level ) )
+            result <- cbind( result , ci )
         }
         if ( the.class=="map2stan" ) {
             # HPDI from samples
-            post <- as.data.frame( extract.samples(model) )
-            ci <- t( apply( post , 2 , HPDI , prob=level ) )
+            post <- extract.samples(model)
+            result <- postlistprecis( post , prob=level )
         }
         if ( the.class=="stanfit" ) {
             # HPDI from samples
-            post <- as.data.frame( extract.samples(model) )
-            ci <- t( apply( post , 2 , HPDI , prob=level ) )
+            post <- extract.samples(model)
+            result <- postlistprecis( post , prob=level )
         }
-        result <- cbind( result , ci )
     }
     if ( corr==TRUE ) {
         result <- cbind( result , Rho )