Add an activity center image to the README and make the exponential e…

…ncounter model a bit more streamlined

README.md

@@ -1,4 +1,6 @@
 
+![](activity-center.png)
+
 # scr-stan: spatial capture-recapture examples in Stan
 
 <!-- badges: start -->

ch07/exponential-encounter-fn-data-aug.R

@@ -30,3 +30,28 @@ stan_d <- list(M = n_obs + n_aug,
 m_init <- stan_model("ch07/exponential-encounter-fn-s-data-aug.stan")
 m_fit <- sampling(m_init, data = stan_d)
 traceplot(m_fit, pars = c("alpha1", "alpha0", "N"))
+
+
+# Visualize some activity centers --------------------------------
+library(ggplot2)
+library(patchwork)
+
+post <- rstan::extract(m_fit)
+
+plot_center <- function(ind) {
+  ggplot(trap_coords, aes(V1, V2)) + 
+    theme_void() + 
+    stat_density_2d(aes(fill = after_stat(level)), geom = "polygon", 
+                    data = as.data.frame(post$s[, ind, ]), 
+                    color = "white", size = .05) +
+    scale_fill_gradient(low = "white", high = "#b2001d") +
+    theme(legend.position = "none") + 
+    geom_point(alpha = .4, shape = 3, size = .8) + 
+    geom_point(data = trap_coords[y[ind, ] > 0, ]) + 
+    xlim(stan_d$xlim) + 
+    ylim(stan_d$ylim) + 
+    coord_equal()
+}
+p <- wrap_plots(lapply(c(40, 3, 46), plot_center), nrow = 1)
+p
+ggsave("activity-center.png", width = 8, height = 2, dpi = 450)

ch07/exponential-encounter-fn-s-data-aug.stan

@@ -26,6 +26,8 @@ parameters {
 transformed parameters {
   matrix[M, n_trap] logit_p;
   matrix[M, 2] s = append_col(s1, s2);
+  vector[M] lp_if_present;
+  vector[M] log_lik;
 
   {
     matrix[M, n_trap] dist_to_trap;
@@ -40,6 +42,17 @@ transformed parameters {
     log_p = log_inv_logit(alpha0) - alpha1 * dist_to_trap;
     logit_p = log_p - log1m_exp(log_p);
   }
+
+  for (i in 1:M) {
+    lp_if_present[i] = bernoulli_lpmf(1 | psi)
+                       + binomial_logit_lpmf(y[i, ] | n_occasion, logit_p[i, ]);
+    if (sum(y[i, ]) > 0) {
+      log_lik[i] = lp_if_present[i];
+    } else {
+      log_lik[i] = log_sum_exp(lp_if_present[i], bernoulli_lpmf(0 | psi));
+    }
+  }
+
 }
 
 model {
@@ -48,18 +61,7 @@ model {
   alpha1 ~ normal(0, 3);
 
   // likelihood
-  for (i in 1:M) {
-    if (sum(y[i, ]) > 0) {
-      1 ~ bernoulli(psi);
-      y[i, ] ~ binomial_logit_lpmf(n_occasion, logit_p[i, ]);
-    } else {
-      target += log_sum_exp(
-        bernoulli_lpmf(1 | psi) 
-          + binomial_logit_lpmf(y[i, ] | n_occasion, logit_p[i, ]), 
-        bernoulli_lpmf(0 | psi)
-      );
-    }
-  }
+  target += sum(log_lik);
 }
 
 generated quantities {
@@ -72,16 +74,9 @@ generated quantities {
 
     for (i in 1:M) {
       if (sum(y[i, ]) > 0) {
-        lp_present[i] = 0;
         z[i] = 1;
       } else {
-        lp_present[i] = bernoulli_lpmf(1 | psi) 
-                        + binomial_logit_lpmf(y[i, ] | n_occasion, logit_p[i, ])
-                        - log_sum_exp(
-                          bernoulli_lpmf(1 | psi) 
-                          + binomial_logit_lpmf(y[i, ] | n_occasion, logit_p[i, ]), 
-                          bernoulli_lpmf(0 | psi)
-                        );
+        lp_present[i] = lp_if_present[i] - log_lik[i];
         z[i] = bernoulli_rng(exp(lp_present[i]));
       }
     }