(based on R. Baathâ€™s UseR! 2015 tutorial) # Model 1 â€“ Simple

``````N.draw = 500000 # number of replicates
x = 500 # number of bills marked in the initial capture
y = 300 # number of bills sampled in the second capture
w = 127 # number of marked bills in the sample
upper.limit = 1500 # maximum (theoretical) number of bills
bin.width = 50 # for plotting the posterio
# Defining and drawing from the prior distribution
N.bills = sample (x:1500, N.draw, replace=TRUE) # number of marked bills in each experiment
barplot(table(cut(N.bills, seq(x, upper.limit, bin.width)))/ length(N.bills), col = "gray")``````

``````# Defining the generative model
pick.bills <- function(N.bills) {
bills <- rep(0:1, c(N.bills - x, x)) # 0 for un-marked, 1 for marked in the inital capture
sum(sample(bills, y)) # sampling y bills in the second capture
}
# Simulating the data
N.marked <- rep(NA, N.draw)
for(i in 1:N.draw) {
N.marked[i] <- pick.bills(N.bills[i])
}
# Filtering out parameter values N.marked that were not equal to w
post.bills <- N.bills[N.marked == w]
# Plotting the posterior distribution
length(post.bills)``````
``[1] 4699``
``barplot(table(cut(post.bills, seq(x,upper.limit,bin.width))) / length(post.bills), col = "blue")``

``````# Statistics
min(post.bills)``````
``[1] 996``
``mean(post.bills)``
``[1] 1194``
``median(post.bills)``
``[1] 1191``
``max(post.bills)``
``[1] 1450``

# Model 2 â€“ Marked Bills Are Brittle

``````N.draw = 500000 # number of replicates
x = 500 # number of bills marked in the initial capture
y = 300 # number of bills sampled in the second capture
w = 127 # number of marked bills in the sample
u = 0.9 # probability that marked bills will be retired
upper.limit = 1500 # maximum (theoretical) number of bills
bin.width = 50 # for plotting the posterio
# Defining and drawing from the prior distribution
N.bills = sample (x:1500, N.draw, replace=TRUE) # number of marked bills in each experiment
barplot(table(cut(N.bills, seq(x, upper.limit, bin.width)))/ length(N.bills), col = "gray")``````

``````# Defining the generative model
pick.bills <- function(N.bills) {
bills <- rep(0:1, c(N.bills - x, x)) # 0 for un-marked, 1 for marked in the inital capture
prob.pick <- ifelse(bills == 0, 1.0, u)
sum(sample(bills, y, prob = prob.pick)) # sampling y bills in the second capture
}
# Simulating the data
N.marked <- rep(NA, N.draw)
for(i in 1:N.draw) {
N.marked[i] <- pick.bills(N.bills[i])
}
# Filtering out parameter values N.marked that were not equal to w
post.bills <- N.bills[N.marked == w]
# Plotting the posterior distribution
length(post.bills)``````
``[1] 4277``
``barplot(table(cut(post.bills, seq(x,upper.limit,bin.width))) / length(post.bills), col = "blue")``

``````# Statistics
min(post.bills)``````
``[1] 933``
``mean(post.bills)``
``[1] 1135``
``median(post.bills)``
``[1] 1132``
``max(post.bills)``
``[1] 1427``

# Model 3 â€“ Listen to the Banker

``````N.draw = 500000 # number of replicates
x = 500 # number of bills marked in the initial capture
y = 300 # number of bills sampled in the second capture
w = 127 # number of marked bills in the sample
u = 0.9 # probability that marked bills will be retired
banker.mean = 1000
upper.limit = 1500 # maximum (theoretical) number of bills
bin.width = 50 # for plotting the posterio
# Defining and drawing from the prior distribution
N.bills = rnbinom(N.draw, mu = banker.mean - x, size = w) + x # number of marked bills in each experiment, using the banker's experience
barplot(table(cut(N.bills, seq(x, upper.limit, bin.width)))/ length(N.bills), col = "gray")``````

``````# Defining the generative model
pick.bills <- function(N.bills) {
bills <- rep(0:1, c(N.bills - x, x)) # 0 for un-marked, 1 for marked in the inital capture
prob.pick <- ifelse(bills == 0, 1.0, u)
sum(sample(bills, y, prob = prob.pick)) # sampling y bills in the second capture
}
# Simulating the data
N.marked <- rep(NA, N.draw)
for(i in 1:N.draw) {
N.marked[i] <- pick.bills(N.bills[i])
}
# Filtering out parameter values N.marked that were not equal to w
post.bills <- N.bills[N.marked == w]
# Plotting the posterior distribution
length(post.bills)``````
``[1] 5296``
``barplot(table(cut(post.bills, seq(x,upper.limit,bin.width))) / length(post.bills), col = "blue")``