Here are the scripts that were “live coded” during the tutorial, all relating to the fish capture re-capture problem.

Original capture re-capture using approximate Bayesian Computation

n_draw <- 100000

# Defining and drawing from the prior distribution
n_fish <- sample(20:250, n_draw, replace = TRUE)

# Defining the generative model
pick_fish <- function(n_fish) {
  fish <- rep(0:1, c(n_fish - 20, 20))
  sum(sample(fish, 20))
}

# Simulating the data
n_marked <- rep(NA, n_draw)
for(i in 1:n_draw) {
  n_marked[i] <- pick_fish(n_fish[i])
}

# Filtering out those parameter values that didn't result in the
# data that we actually observed
post_fish <- n_fish[n_marked == 5]

hist(post_fish)

length(post_fish)
## [1] 8271
# The posterior distribution showing the probability of different number of fish
# (binning here in bins of 20 just make the graph easier to interpret)
barplot(table(cut(post_fish, seq(0, 250, 20))) / length(post_fish), col = "salmon")

Modified model that accounts for that marked fish get “shy”

n_draw <- 100000
n_fish <- sample(20:250, n_draw, replace = TRUE)

# This is the only part of the code that has changed from the original version above.
pick_fish <- function(n_fish) {
  fish <- rep(0:1, c(n_fish - 20, 20))
  prob_pick <- ifelse(fish == 0, 1.0, 0.5)
  sum(sample(fish, 20, prob = prob_pick))
}

n_marked <- rep(NA, n_draw)
for(i in 1:n_draw) {
  n_marked[i] <- pick_fish(n_fish[i])
}

post_fish <- n_fish[n_marked == 5]

hist(post_fish)

length(post_fish)
## [1] 4402
barplot(table(cut(post_fish, seq(0, 250, 20))) / length(post_fish), col = "salmon")

Modified model that accounts for the “expert” opinion of the fisherman

n_draw <- 100000

# This is the only part of the code that has changed from the original version above.
n_fish <- rnbinom(n_draw, mu = 200 - 20, size = 4) + 20
hist(n_fish)

pick_fish <- function(n_fish) {
  fish <- rep(0:1, c(n_fish - 20, 20))
  prob_pick <- ifelse(fish == 0, 1.0, 0.5)
  sum(sample(fish, 20, prob = prob_pick))
}

n_marked <- rep(NA, n_draw)
for(i in 1:n_draw) {
  n_marked[i] <- pick_fish(n_fish[i])
}

post_fish <- n_fish[n_marked == 5]

hist(post_fish)

length(post_fish)
## [1] 1790
barplot(table(cut(post_fish, seq(0, 250, 20))) / length(post_fish), col = "salmon")

A JAGS implementation of the original model

Here the dhyper distribution is used as it implements the same process as the fish picking model.

library(rjags)

data_list <- list(
  n_picked = 20,
  n_unmarked = 20 - 5)

model_string <- "
model {
  n_fish_real   ~ dunif(0, 250)
  n_fish <- round(n_fish_real)
  n_unmarked ~ dhyper(n_fish - n_picked, n_picked, n_picked, 1)
}
"

jags_model <- jags.model(textConnection(model_string), data = data_list)
## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
##    Graph Size: 8
## 
## Initializing model
samples <- coda.samples(jags_model, c("n_fish"), n.iter = 3000)
plot(samples)

s_mat <- as.matrix(samples)

# The probability that there are more than a 100 fish
mean(s_mat[,"n_fish"] > 100)
## [1] 0.4877