v1.70 step

- added extract.prior methods for map and map2stan models. See
?extract.prior for details.
- vectorized rmvnorm2
- precis_plot (plot method for precis objects) now understands “add”
argument. Helps with superimposing posterior on prior in these plots.
See example in ?extract.prior.

DESCRIPTION

@@ -1,8 +1,8 @@
 Package: rethinking
 Type: Package
 Title: Statistical Rethinking book package
-Version: 1.64
-Date: 2017-8-8
+Version: 1.70
+Date: 2017-8-15
 Author: Richard McElreath
 Maintainer: Richard McElreath <[email protected]>
 Imports: coda, MASS, mvtnorm, loo

R/denschart.R

@@ -55,6 +55,11 @@ denschart <- function (x, pars,
                 }
             }
             nd <- length(dim(x[[i]]))
+            if ( nd==0 ) {
+                # no dim, so probably a numeric vector
+                # convert to 1D array
+                x[[i]] <- as.array(x[[i]])
+            }
             if ( nd==3 ) {
                 # matrix
                 if ( drop_matrices==FALSE )

R/distributions.r

@@ -376,9 +376,80 @@ dmvnorm2 <- function( x , Mu , sigma , Rho , log=FALSE ) {
     dmvnorm(x, Mu, SIGMA, log=log )
 }
 
+# vectorized version of rmvnorm, with split scale vector and correlation matrix
 rmvnorm2 <- function( n , Mu=rep(0,length(sigma)) , sigma=rep(1,length(Mu)) , Rho=diag(length(Mu)) , method="chol" ) {
-    DS <- diag(sigma)
-    SIGMA <- DS %*% Rho %*% DS
-    rmvnorm( n=n , mean=Mu , sigma=SIGMA , method=method )
+    ldim <- function(x) { 
+        z <- length(dim(x))
+        if ( z==0 ) z <- 1
+        return(z)
+    }
+    if ( ldim(Mu)>1 || ldim(sigma)>1 || ldim(Rho)>1 ) {
+        m <- 0 # dimension of multi_normal
+        mR <- 0 # dim of Rho
+        # vectorization needed
+        # make Mu, sigma, and Rho same implied length
+        K_Mu <- 1
+        if ( ldim(Mu)==2 ) {
+            K_Mu <- dim(Mu)[1] # matrix
+            m <- dim(Mu)[2]
+        } else {
+            m <- length(Mu)
+        }
+
+        K_sigma <- 1
+        if ( ldim(sigma)==2 ) {
+            K_sigma <- dim(sigma)[1] # matrix
+        }
+        K_Rho <- 1
+        if ( ldim(Rho)==3 ) {
+            K_Rho <- dim(Rho)[1] # array (vector of matrices)
+        }
+        mR <- dim(Rho)[2] # should work whether Rho is matrix or array
+
+        # get largest dim
+        #K <- max( K_Mu , K_sigma , K_Rho )
+        K <- n
+
+        # need to fill to same size
+        if ( K_Mu < K ) {
+            Mu_new <- matrix( NA , nrow=K , ncol=m )
+            row_list <- rep( 1:K_Mu , length.out=K )
+            if ( K_Mu==1 )
+                for ( i in 1:K )
+                    Mu_new[i,] <- Mu
+            else
+                for ( i in 1:K )
+                    Mu_new[i,] <- Mu[row_list[i],]
+            Mu <- Mu_new
+        }
+        if ( K_sigma < K ) {
+            sigma_new <- matrix( NA , nrow=K , ncol=m )
+            row_list <- rep( 1:K_sigma , length.out=K )
+            for ( i in 1:K )
+                sigma_new[i,] <- sigma[row_list[i],]
+            sigma <- sigma_new
+        }
+        if ( K_Rho < K ) {
+            Rho_new <- array( NA , dim=c( K , m , m ) )
+            row_list <- rep( 1:K_Rho , length.out=K )
+            for ( i in 1:K )
+                Rho_new[i,,] <- Rho[row_list[i],,]
+            Rho <- Rho_new
+        }
+        # should all be same length now, so can brute force vectorize over and sample
+
+        result <- array( NA , dim=c( K , m ) )
+        for ( i in 1:n ) {
+            DS <- diag(sigma[i,])
+            SIGMA <- DS %*% Rho[i,,] %*% DS
+            result[i,] <- rmvnorm( n=1 , mean=Mu[i,] , sigma=SIGMA , method=method )
+        }
+        return(result)
+
+    } else {
+        DS <- diag(sigma)
+        SIGMA <- DS %*% Rho %*% DS
+        return( rmvnorm( n=n , mean=Mu , sigma=SIGMA , method=method ) )
+    }
 }
 

R/map.r

@@ -527,6 +527,7 @@ map <- function( flist , data , start , method="BFGS" , hessian=TRUE , debug=FAL
             formula_parsed = flist2,
             fminuslogl = fmll,
             links = links )
+    attr(m,"formula_exploded") <- flist # expanded c() constructs
     attr(m,"df") <- length(m@coef)
     attr(m,"veclist") <- veclist
     if (!missing(data)) attr(m,"nobs") = length(data[[1]])

R/map2stan.r

@@ -1961,6 +1961,8 @@ map2stan <- function( flist , data , start , pars , constraints=list() , types=l
             pars = pars,
             formula = flist.orig,
             formula_parsed = fp )
+
+        attr(result,"constraints") <- constraints
 
         attr(result,"df") = length(result@coef)
         if ( DIC==TRUE ) {

R/precis.r

@@ -24,7 +24,7 @@ precis_show <- function( object ) {
 
 setMethod( "show" , "precis" , function(object) precis_show(object) )
 
-precis_plot <- function( x , y , pars , col.ci="black" , xlab="Value" , ... ) {
+precis_plot <- function( x , y , pars , col.ci="black" , xlab="Value" , add=FALSE , xlim=NULL , ... ) {
     #x <- x@output
     if ( !missing(pars) ) {
         x <- x[pars,]
@@ -34,9 +34,13 @@ precis_plot <- function( x , y , pars , col.ci="black" , xlab="Value" , ... ) {
     left <- x[[3]][n:1]
     right <- x[[4]][n:1]
     set_nice_margins()
-    dotchart( mu , labels=rownames(x)[n:1] , xlab=xlab , xlim=c(min(left),max(right)) , ... )
+    if ( is.null(xlim) ) xlim <- c(min(left),max(right))
+    if ( add==FALSE )
+        dotchart( mu , labels=rownames(x)[n:1] , xlab=xlab , xlim=xlim , ... )
+    else
+        points( mu[n:1] , n:1 , ... )
     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) )
+    if ( add==FALSE ) abline( v=0 , lty=1 , col=col.alpha("black",0.15) )
 }
 setMethod( "plot" , "precis" , function(x,y,...) precis_plot(x,y,...) )