-
-
Save pboesu/03f5708ef6434557d67fb62889ce02ec to your computer and use it in GitHub Desktop.
Distance sampling in Stan
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
# 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) |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
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); | |
} |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment