Distance sampling in Stan

dist.R

# Distance sampling model with data augmentation. 
# Based on section 8.3.1 of the Applied Hierarchical Modeling book by Royle and Kery

library(rstan)

B <- 50

# note that I'm dividing by 10 and adding 1e-6 to put this on a manageable
# scale and prevent values == 0
d_obs <- c(71.93, 26.05, 58.47, 92.35, 163.83, 84.52, 163.83, 157.33,
           22.27, 72.11, 86.99, 50.8, 0, 73.14, 0, 128.56, 163.83, 71.85,
           30.47, 71.07, 150.96, 68.83, 90, 64.98, 165.69, 38.01, 378.21,
           78.15, 42.13, 0, 400, 175.39, 30.47, 35.07, 86.04, 31.69, 200,
           271.89, 26.05, 76.6, 41.04, 200, 86.04, 0, 93.97, 55.13, 10.46,
           84.52, 0, 77.65, 0, 96.42, 0, 64.28, 187.94, 0, 160.7, 150.45,
           63.6, 193.19, 106.07, 114.91, 143.39, 128.56, 245.75, 123.13,
           123.13, 153.21, 143.39, 34.2, 96.42, 259.81, 8.72) / 10 + 1e-6
n_aug <- 200
n_obs <- length(d_obs)
M <- n_obs + n_aug
y <- c(rep(1, n_obs), rep(0, n_aug))

hist(d_obs)

stan_d <- list(M = M, 
               n_obs = n_obs, 
               d_obs = d_obs, 
               B = B, 
               y = y)

m_init <- stan_model("dist.stan")
m_fit <- sampling(m_init, data = stan_d, cores = 4)
traceplot(m_fit)

dist.stan

data {
  int<lower = 1> M;
  int<lower = 0, upper = 1> y[M];
  int n_obs;
  real B;
  vector<lower = 0, upper = B>[n_obs] d_obs;
}

transformed data {
  int n_aug = M - n_obs;
}

parameters {
  real<lower = 0> sigma;
  real<lower = 0, upper = 1> psi;
  vector<lower = 0, upper = B>[n_aug] d_aug;
}

transformed parameters {
  vector[M] logit_p;
  vector[M] log_lik;
  
  {
    vector[M] logp;
    vector[M] lp_if_present;
    vector[M] distances;
    real denom = 2*sigma^2;
    for (i in 1:n_obs) {
      distances[i] = d_obs[i];
    }
    for (i in 1:n_aug) {
      distances[n_obs + i] = d_aug[i];
    }
    
    for (i in 1:M) {
      logp[i] = - distances[i]^2 / denom;
      logit_p[i] = logp[i] - log1m_exp(logp[i]);
    
      lp_if_present[i] = bernoulli_lpmf(1 | psi)
                       + bernoulli_logit_lpmf(y[i] | logit_p[i]);
      if (y[i]) {
        log_lik[i] = lp_if_present[i];
      } else {
        log_lik[i] = log_sum_exp(lp_if_present[i], bernoulli_lpmf(0 | psi));
      }
    }
  }
}

model {
  // priors
  sigma ~ normal(0, 5);
  
  // likelihood
  target += sum(log_lik);
}