1.70 increment

- added hex sticker designs, in case people want to print their own
- stub for Cormack-Jolly-Seber map2stan template
- map2stan declares everything in data list, by default. See
declare_all_data argument.
- precis_plot method understands optional labels argument, for custom
row labels.
- extract prior should function better with map() fits, returning
everything as an array

ERRATA.md

@@ -37,6 +37,8 @@ page 331, line 1: "a women" -> "a woman"
 
 page 386, problem 12H1, first paragraph: 'By the year 200' should read 'By the year 2000'.
 
+page 403: The average effect in the P *C interaction model is typed βP but I think should be βPC.
+
 page 435: "FIGURE 14.4 display ... imputed neocortex values in blue ...
 shown by the blue line segments". The imputed values are actually the
 open black dots (and corresponding black line segments) as the caption

R/map2stan-class.r

@@ -449,3 +449,19 @@ stanergy <- function( x , colscheme="blue" , binwidth=NULL , merge_chains=FALSE
     color_scheme_set(colscheme) # we hates the ggplot
     mcmc_nuts_energy( np , merge_chains = merge_chains , binwidth=binwidth )
 }
+
+# function to do mcmc_parcoord to help find divergences
+divergence_tracker <- function( x , no_lp=TRUE , pars , ... ) {
+    if ( class(x)=="map2stan" ) x <- x@stanfit
+    np <- nuts_params(x)
+    draws <- as.array(x)
+    if ( missing(pars) ) {
+        pars <- dimnames(x)$parameters
+        if ( no_lp==TRUE ) {
+            pars <- pars[ -which(pars=="lp__") ]
+        }
+    }
+    require(bayesplot)
+    bayesplot::mcmc_parcoord(draws,np=np,pars=pars,...)
+}
+# divergence_tracker( m5 )

R/map2stan-templates.r

@@ -1020,6 +1020,76 @@ dev <- dev + (-2)*(bernoulli_lpmf(0|PAR1) + binomial_lpmf(OUTCOME|PAR2,PAR3));",
             return(c(new_k[1],theta_txt));
         },
         vectorized = FALSE
+    ),
+    CormackJollySeberK = list(
+        name = "CormackJollySeberK",
+        R_name = "dCJS",
+        stan_name = "increment_log_prob",
+        stan_code = 
+"if (cjs_last[i] > 0) {
+        for (k in (cjs_first[i]+1):cjs_last[i]) {
+            target += (log(PAR1[k-1]));     // i survived from k-1 to k
+            if (OUTCOME[i,k] == 1)
+                target += (log(PAR2[k]));       // i captured at k
+            else
+                target += (log1m(PAR2[k]));     // i not captured at k
+        }
+        target += (log(PAR1[cjs_last[i]]));   // i not seen after last[i]
+    }",
+        stan_dev = "",
+        stan_loglik = "",
+        num_pars = 2,
+        par_names = c("phi","p"), # phi: survival prob; p: detection prob (conditional on survival)
+        par_bounds = c("lower=0,upper=1","lower=0,upper=1"),
+        par_types = c("real","real"),
+        out_type = "matrix", # each outcome is a vector of 0/1 indicating detection history of individual
+        par_map = function(k,e,...) {
+            # to-do
+            # (1) translate histories to first,list,n_captured variables
+            # (2) compute chi parameters in transformed parameters
+            #     chi[k] = Pr[ no capture > k | alive at k ]
+            # (3) generated quantities code for abundance estimate
+            return(k)
+        },
+        # new 'code' slot uses specification of custom code map2stan argument
+        code = list(
+            list( 
+"int<lower=0,upper=cjs_K+1> cjs_first[cjs_I];     // first[i]: ind i first capture
+int<lower=0,upper=cjs_K+1> cjs_last[cjs_I];      // last[i]:  ind i last capture
+int<lower=0,upper=cjs_I> cjs_n_captured[cjs_K];  // n_capt[k]: num aptured at k",
+                block="transformed data", section="declare" ),
+            list(
+"cjs_first = rep_array(cjs_K+1,cjs_I);
+cjs_last = rep_array(0,cjs_I);
+for (i in 1:cjs_I) {
+    for (k in 1:cjs_K) {
+        if (OUTCOME[i,k] == 1) {
+            if (k < cjs_first[i]) cjs_first[i] <- k;
+            if (k > cjs_last[i]) cjs_last[i] <- k;
+        }
+    }
+}
+n_captured <- rep_array(0,K);
+for (i in 1:I)
+    for (k in 1:K)
+        n_captured[k] <- n_captured[k] + X[i,k];",
+                block="transformed data", section="body" ),
+            list(
+"vector<lower=0,upper=1>[cjs_K] chi;   // chi[k]: Pr[no capture > k | alive at k]",
+                block="transformed parameters", section="declare" ),
+            list(
+"{
+    int k;
+    chi[cjs_K] <- 1.0;              
+    k <- cjs_K - 1;
+    while (k > 0) {
+        chi[k] <- (1 - phi[k]) + phi[k] * (1 - p[k+1]) * chi[k+1]; 
+        k <- k - 1;
+    }
+}",
+                block="transformed parameters", section="body" )
+        ),
+        vectorized = FALSE
     )
 )
 

R/map2stan.r

@@ -11,7 +11,7 @@
 ##################
 # map2stan itself
 
-map2stan <- function( flist , data , start , pars , constraints=list() , types=list() , sample=TRUE , iter=2000 , warmup=floor(iter/2) , chains=1 , debug=FALSE , verbose=FALSE , WAIC=TRUE , cores=1 , rng_seed , rawstanfit=FALSE , control=list(adapt_delta=0.95) , add_unique_tag=TRUE , code , log_lik=FALSE , DIC=FALSE , ... ) {
+map2stan <- function( flist , data , start , pars , constraints=list() , types=list() , sample=TRUE , iter=2000 , warmup=floor(iter/2) , chains=1 , debug=FALSE , verbose=FALSE , WAIC=TRUE , cores=1 , rng_seed , rawstanfit=FALSE , control=list(adapt_delta=0.95) , add_unique_tag=TRUE , code , log_lik=FALSE , DIC=FALSE , declare_all_data=TRUE , ... ) {
 
     if ( missing(rng_seed) ) rng_seed <- sample( 1:1e5 , 1 )
     set.seed(rng_seed)
@@ -1465,7 +1465,10 @@ map2stan <- function( flist , data , start , pars , constraints=list() , types=l
             if ( class(d[[var$var]])=="integer" ) type <- "int"
             if ( class(d[[var$var]])=="matrix" ) {
                 type <- "matrix"
-                if ( is.null(var$N) ) var$N <- dim(d[[var$var]])
+                if ( is.null(var$N) ) 
+                    var$N <- dim(d[[var$var]])
+                else if ( length(var$N)<2 )
+                    var$N <- dim(d[[var$var]])
             }
             # coerce outcome type
             if ( !is.null(var$type) ) type <- var$type
@@ -1806,10 +1809,10 @@ map2stan <- function( flist , data , start , pars , constraints=list() , types=l
 
         if ( is.null(previous_stanfit) ) {
             #fit <- try(stan( model_code=model_code , model_name=modname , data=d , init=initlist , iter=iter , warmup=warmup , chains=chains , cores=cores , pars=pars , seed=rng_seed , ... ))
-            fit <- try(stan( model_code=model_code , data=d , init=initlist , iter=iter , warmup=warmup , chains=chains , cores=cores , pars=pars , seed=rng_seed , ... ))
+            fit <- try(stan( model_code=model_code , data=d , init=initlist , iter=iter , warmup=warmup , chains=chains , cores=cores , pars=pars , seed=rng_seed , control=control , ... ))
         }
         else
-            fit <- try(stan( fit=previous_stanfit , model_name=modname , data=d , init=initlist , iter=iter , warmup=warmup , chains=chains , cores=cores , pars=pars , seed=rng_seed , ... ))
+            fit <- try(stan( fit=previous_stanfit , model_name=modname , data=d , init=initlist , iter=iter , warmup=warmup , chains=chains , cores=cores , pars=pars , seed=rng_seed , control=control , ... ))
 
         if ( class(fit)=="try-error" ) {
                 # something went wrong in at least one chain

R/plotting.r

@@ -303,7 +303,7 @@ shade <- function( object , lim , label=NULL , col=col.alpha("black",0.15) , bor
     }
 }
 
-mcmcpairs <- function( posterior , cex=0.3 , pch=16 , col=col.alpha("slateblue",0.2) , n=1000 , adj=1 , ... ) {
+mcmcpairs_old <- function( posterior , cex=0.3 , pch=16 , col=col.alpha("slateblue",0.2) , n=1000 , adj=1 , ... ) {
     panel.dens <- function(x, ...) {
         usr <- par("usr"); on.exit(par(usr))
         par(usr = c(usr[1:2], 0, 1.5) )
@@ -329,6 +329,35 @@ mcmcpairs <- function( posterior , cex=0.3 , pch=16 , col=col.alpha("slateblue",
     pairs( posterior , cex=cex , pch=pch , col=col , upper.panel=panel.2d , lower.panel=panel.cor , diag.panel=panel.dens , ... )
 }
 
+mcmcpairs <- function( x , alpha=0.7 , cex=0.7 , pch=16 , adj=1 , pars , n=500 , ... ) {
+    # x should be samples
+    panel.dens <- function(x, ...) {
+        usr <- par("usr"); on.exit(par(usr))
+        par(usr = c(usr[1:2], 0, 1.5) )
+        h <- density(x,adj=adj)
+        y <- h$y
+        y <- y/max(y)
+        abline( v=0 , col="gray" , lwd=0.5 )
+        lines( h$x , y )
+    }
+    panel.2d <- function( x , y , ... ) {
+        i <- sample( 1:length(x) , size=n )
+        abline( v=0 , col="gray" , lwd=0.5 )
+        abline( h=0 , col="gray" , lwd=0.5 )
+        dcols <- densCols( x[i] , y[i] )
+        dcols <- sapply( dcols , function(k) col.alpha(k,alpha) )
+        points( x[i] , y[i] , col=dcols , ... )
+    }
+    panel.cor <- function( x , y , ... ) {
+        k <- cor( x , y )
+        cx <- sum(range(x))/2
+        cy <- sum(range(y))/2
+        text( cx , cy , round(k,2) , cex=2*exp(abs(k))/exp(1) )
+    }
+    set_nice_margins()
+    pairs( x , cex=cex , pch=pch , upper.panel=panel.2d , lower.panel=panel.cor , diag.panel=panel.dens , ... )
+}
+
 # wickham's histospark
 # https://github.com/hadley/precis/blob/master/R/histospark.R
 sparks <- c("\u2581", "\u2582", "\u2583", "\u2585", "\u2587")

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" , add=FALSE , xlim=NULL , ... ) {
+precis_plot <- function( x , y , pars , col.ci="black" , xlab="Value" , add=FALSE , xlim=NULL , labels=rownames(x)[n:1] , ... ) {
     #x <- x@output
     if ( !missing(pars) ) {
         x <- x[pars,]
@@ -36,7 +36,7 @@ precis_plot <- function( x , y , pars , col.ci="black" , xlab="Value" , add=FALS
     set_nice_margins()
     if ( is.null(xlim) ) xlim <- c(min(left),max(right))
     if ( add==FALSE )
-        dotchart( mu , labels=rownames(x)[n:1] , xlab=xlab , xlim=xlim , ... )
+        dotchart( mu , labels=labels , 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 )

R/z_extract_prior.r

@@ -106,6 +106,12 @@ function( fit , n=1000 , pars , ... ) {
         result <- new_result
     }
 
+    # make sure each entry is an array, even if 1D
+    for ( i in 1:length(result) ) {
+        if ( is.null(dim(result[[i]])) )
+            result[[i]] <- as.array(result[[i]])
+    }
+
     return(result)
 })
 
@@ -204,7 +210,7 @@ function( fit , n=1000 , ... ) {
                 # vector prior like: a[id] ~ normal( 0 , sigma_a )
                 dims <- fit@stanfit@par_dims[[ pnames[[1]] ]]
                 x <- sapply( 1:dims , function(k) sample_from_prior( dname , args , n ) )
-                result[[ pnames[[1]] ]] <- x
+                result[[ pnames[[1]] ]] <- as.array(x)
             } else {
                 # multivariate like: c(a,b)[id] ~ multi_normal( 0 , Sigma )
                 # sample and then shunt into correct named slots
@@ -233,6 +239,12 @@ function( fit , n=1000 , ... ) {
     }#i
     result <- new_result
 
+    # make sure each entry is an array, even if 1D
+    for ( i in 1:length(result) ) {
+        if ( is.null(dim(result[[i]])) )
+            result[[i]] <- as.array(result[[i]])
+    }
+
     return(result)
 
 })

R/z_link-map2stan.r

@@ -60,6 +60,7 @@ function( fit , data , n=1000 , post , refresh=0.1 , replace=list() , flatten=TR
 
     # number of samples
     n_samples <- dim( post[[1]] )[1]
+    if ( is.null(n_samples) ) n_samples <- length(post[[1]])
     if ( n == 0 ) n <- n_samples # special flag for all samples in fit
     if ( n_samples < n ) n <- n_samples