Bayesian Data Analysis Tutorial
(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] 4699barplot(table(cut(post.bills, seq(x,upper.limit,bin.width))) / length(post.bills), col = "blue")
# Statistics
min(post.bills)[1] 996mean(post.bills)[1] 1194median(post.bills)[1] 1191max(post.bills)[1] 1450Model 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] 4277barplot(table(cut(post.bills, seq(x,upper.limit,bin.width))) / length(post.bills), col = "blue")
# Statistics
min(post.bills)[1] 933mean(post.bills)[1] 1135median(post.bills)[1] 1132max(post.bills)[1] 1427Model 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] 5296barplot(table(cut(post.bills, seq(x,upper.limit,bin.width))) / length(post.bills), col = "blue")
# Statistics
min(post.bills)[1] 948mean(post.bills)[1] 1058median(post.bills)[1] 1057max(post.bills)[1] 1217LS0tCnRpdGxlOiAiQmF5ZXNpYW4gRGF0YSBBbmFseXNpcyBUdXRvcmlhbCIKb3V0cHV0OiBodG1sX25vdGVib29rCi0tLQooYmFzZWQgb24gUi4gQmFhdGgncyBVc2VSISAyMDE1IHR1dG9yaWFsKQojIE1vZGVsIDEgLS0gU2ltcGxlCmBgYHtyfQpOLmRyYXcgPSA1MDAwMDAgIyBudW1iZXIgb2YgcmVwbGljYXRlcwp4ID0gNTAwICMgbnVtYmVyIG9mIGJpbGxzIG1hcmtlZCBpbiB0aGUgaW5pdGlhbCBjYXB0dXJlCnkgPSAzMDAgIyBudW1iZXIgb2YgYmlsbHMgc2FtcGxlZCBpbiB0aGUgc2Vjb25kIGNhcHR1cmUKdyA9IDEyNyAjIG51bWJlciBvZiBtYXJrZWQgYmlsbHMgaW4gdGhlIHNhbXBsZQp1cHBlci5saW1pdCA9IDE1MDAgIyBtYXhpbXVtICh0aGVvcmV0aWNhbCkgbnVtYmVyIG9mIGJpbGxzCmJpbi53aWR0aCA9IDUwICMgZm9yIHBsb3R0aW5nIHRoZSBwb3N0ZXJpbwoKIyBEZWZpbmluZyBhbmQgZHJhd2luZyBmcm9tIHRoZSBwcmlvciBkaXN0cmlidXRpb24KTi5iaWxscyA9IHNhbXBsZSAoeDoxNTAwLCBOLmRyYXcsIHJlcGxhY2U9VFJVRSkgIyBudW1iZXIgb2YgbWFya2VkIGJpbGxzIGluIGVhY2ggZXhwZXJpbWVudApiYXJwbG90KHRhYmxlKGN1dChOLmJpbGxzLCBzZXEoeCwgdXBwZXIubGltaXQsIGJpbi53aWR0aCkpKS8gbGVuZ3RoKE4uYmlsbHMpLCBjb2wgPSAiZ3JheSIpCgojIERlZmluaW5nIHRoZSBnZW5lcmF0aXZlIG1vZGVsCnBpY2suYmlsbHMgPC0gZnVuY3Rpb24oTi5iaWxscykgewogIGJpbGxzIDwtIHJlcCgwOjEsIGMoTi5iaWxscyAtIHgsIHgpKSAjIDAgZm9yIHVuLW1hcmtlZCwgMSBmb3IgbWFya2VkIGluIHRoZSBpbml0YWwgY2FwdHVyZQogIHN1bShzYW1wbGUoYmlsbHMsIHkpKSAjIHNhbXBsaW5nIHkgYmlsbHMgaW4gdGhlIHNlY29uZCBjYXB0dXJlICAgICAgIAp9CgojIFNpbXVsYXRpbmcgdGhlIGRhdGEKTi5tYXJrZWQgPC0gcmVwKE5BLCBOLmRyYXcpIAogIGZvcihpIGluIDE6Ti5kcmF3KSB7CiAgICBOLm1hcmtlZFtpXSA8LSBwaWNrLmJpbGxzKE4uYmlsbHNbaV0pCiAgfQoKIyBGaWx0ZXJpbmcgb3V0IHBhcmFtZXRlciB2YWx1ZXMgTi5tYXJrZWQgdGhhdCB3ZXJlIG5vdCBlcXVhbCB0byB3CnBvc3QuYmlsbHMgPC0gTi5iaWxsc1tOLm1hcmtlZCA9PSB3XQoKIyBQbG90dGluZyB0aGUgcG9zdGVyaW9yIGRpc3RyaWJ1dGlvbgpsZW5ndGgocG9zdC5iaWxscykKYmFycGxvdCh0YWJsZShjdXQocG9zdC5iaWxscywgc2VxKHgsdXBwZXIubGltaXQsYmluLndpZHRoKSkpIC8gbGVuZ3RoKHBvc3QuYmlsbHMpLCBjb2wgPSAiYmx1ZSIpCgojIFN0YXRpc3RpY3MKbWluKHBvc3QuYmlsbHMpCm1lYW4ocG9zdC5iaWxscykKbWVkaWFuKHBvc3QuYmlsbHMpCm1heChwb3N0LmJpbGxzKQpgYGAKIyBNb2RlbCAyIC0tIE1hcmtlZCBCaWxscyBBcmUgQnJpdHRsZQpgYGB7cn0KTi5kcmF3ID0gNTAwMDAwICMgbnVtYmVyIG9mIHJlcGxpY2F0ZXMKeCA9IDUwMCAjIG51bWJlciBvZiBiaWxscyBtYXJrZWQgaW4gdGhlIGluaXRpYWwgY2FwdHVyZQp5ID0gMzAwICMgbnVtYmVyIG9mIGJpbGxzIHNhbXBsZWQgaW4gdGhlIHNlY29uZCBjYXB0dXJlCncgPSAxMjcgIyBudW1iZXIgb2YgbWFya2VkIGJpbGxzIGluIHRoZSBzYW1wbGUKdSA9IDAuOSAjIHByb2JhYmlsaXR5IHRoYXQgbWFya2VkIGJpbGxzIHdpbGwgYmUgcmV0aXJlZAp1cHBlci5saW1pdCA9IDE1MDAgIyBtYXhpbXVtICh0aGVvcmV0aWNhbCkgbnVtYmVyIG9mIGJpbGxzCmJpbi53aWR0aCA9IDUwICMgZm9yIHBsb3R0aW5nIHRoZSBwb3N0ZXJpbwoKIyBEZWZpbmluZyBhbmQgZHJhd2luZyBmcm9tIHRoZSBwcmlvciBkaXN0cmlidXRpb24KTi5iaWxscyA9IHNhbXBsZSAoeDoxNTAwLCBOLmRyYXcsIHJlcGxhY2U9VFJVRSkgIyBudW1iZXIgb2YgbWFya2VkIGJpbGxzIGluIGVhY2ggZXhwZXJpbWVudApiYXJwbG90KHRhYmxlKGN1dChOLmJpbGxzLCBzZXEoeCwgdXBwZXIubGltaXQsIGJpbi53aWR0aCkpKS8gbGVuZ3RoKE4uYmlsbHMpLCBjb2wgPSAiZ3JheSIpCgojIERlZmluaW5nIHRoZSBnZW5lcmF0aXZlIG1vZGVsCnBpY2suYmlsbHMgPC0gZnVuY3Rpb24oTi5iaWxscykgewogIGJpbGxzIDwtIHJlcCgwOjEsIGMoTi5iaWxscyAtIHgsIHgpKSAjIDAgZm9yIHVuLW1hcmtlZCwgMSBmb3IgbWFya2VkIGluIHRoZSBpbml0YWwgY2FwdHVyZQogIHByb2IucGljayA8LSBpZmVsc2UoYmlsbHMgPT0gMCwgMS4wLCB1KQogIHN1bShzYW1wbGUoYmlsbHMsIHksIHByb2IgPSBwcm9iLnBpY2spKSAjIHNhbXBsaW5nIHkgYmlsbHMgaW4gdGhlIHNlY29uZCBjYXB0dXJlICAgICAgIAp9CgojIFNpbXVsYXRpbmcgdGhlIGRhdGEKTi5tYXJrZWQgPC0gcmVwKE5BLCBOLmRyYXcpIAogIGZvcihpIGluIDE6Ti5kcmF3KSB7CiAgICBOLm1hcmtlZFtpXSA8LSBwaWNrLmJpbGxzKE4uYmlsbHNbaV0pCiAgfQoKIyBGaWx0ZXJpbmcgb3V0IHBhcmFtZXRlciB2YWx1ZXMgTi5tYXJrZWQgdGhhdCB3ZXJlIG5vdCBlcXVhbCB0byB3CnBvc3QuYmlsbHMgPC0gTi5iaWxsc1tOLm1hcmtlZCA9PSB3XQoKIyBQbG90dGluZyB0aGUgcG9zdGVyaW9yIGRpc3RyaWJ1dGlvbgpsZW5ndGgocG9zdC5iaWxscykKYmFycGxvdCh0YWJsZShjdXQocG9zdC5iaWxscywgc2VxKHgsdXBwZXIubGltaXQsYmluLndpZHRoKSkpIC8gbGVuZ3RoKHBvc3QuYmlsbHMpLCBjb2wgPSAiYmx1ZSIpCgojIFN0YXRpc3RpY3MKbWluKHBvc3QuYmlsbHMpCm1lYW4ocG9zdC5iaWxscykKbWVkaWFuKHBvc3QuYmlsbHMpCm1heChwb3N0LmJpbGxzKQpgYGAKCiMgTW9kZWwgMyAtLSBMaXN0ZW4gdG8gdGhlIEJhbmtlcgpgYGB7cn0KTi5kcmF3ID0gNTAwMDAwICMgbnVtYmVyIG9mIHJlcGxpY2F0ZXMKeCA9IDUwMCAjIG51bWJlciBvZiBiaWxscyBtYXJrZWQgaW4gdGhlIGluaXRpYWwgY2FwdHVyZQp5ID0gMzAwICMgbnVtYmVyIG9mIGJpbGxzIHNhbXBsZWQgaW4gdGhlIHNlY29uZCBjYXB0dXJlCncgPSAxMjcgIyBudW1iZXIgb2YgbWFya2VkIGJpbGxzIGluIHRoZSBzYW1wbGUKdSA9IDAuOSAjIHByb2JhYmlsaXR5IHRoYXQgbWFya2VkIGJpbGxzIHdpbGwgYmUgcmV0aXJlZApiYW5rZXIubWVhbiA9IDEwMDAKdXBwZXIubGltaXQgPSAxNTAwICMgbWF4aW11bSAodGhlb3JldGljYWwpIG51bWJlciBvZiBiaWxscwpiaW4ud2lkdGggPSA1MCAjIGZvciBwbG90dGluZyB0aGUgcG9zdGVyaW8KCiMgRGVmaW5pbmcgYW5kIGRyYXdpbmcgZnJvbSB0aGUgcHJpb3IgZGlzdHJpYnV0aW9uCk4uYmlsbHMgPSBybmJpbm9tKE4uZHJhdywgbXUgPSBiYW5rZXIubWVhbiAtIHgsIHNpemUgPSB3KSArIHggIyBudW1iZXIgb2YgbWFya2VkIGJpbGxzIGluIGVhY2ggZXhwZXJpbWVudCwgdXNpbmcgdGhlIGJhbmtlcidzIGV4cGVyaWVuY2UKYmFycGxvdCh0YWJsZShjdXQoTi5iaWxscywgc2VxKHgsIHVwcGVyLmxpbWl0LCBiaW4ud2lkdGgpKSkvIGxlbmd0aChOLmJpbGxzKSwgY29sID0gImdyYXkiKQoKIyBEZWZpbmluZyB0aGUgZ2VuZXJhdGl2ZSBtb2RlbApwaWNrLmJpbGxzIDwtIGZ1bmN0aW9uKE4uYmlsbHMpIHsKICBiaWxscyA8LSByZXAoMDoxLCBjKE4uYmlsbHMgLSB4LCB4KSkgIyAwIGZvciB1bi1tYXJrZWQsIDEgZm9yIG1hcmtlZCBpbiB0aGUgaW5pdGFsIGNhcHR1cmUKICBwcm9iLnBpY2sgPC0gaWZlbHNlKGJpbGxzID09IDAsIDEuMCwgdSkKICBzdW0oc2FtcGxlKGJpbGxzLCB5LCBwcm9iID0gcHJvYi5waWNrKSkgIyBzYW1wbGluZyB5IGJpbGxzIGluIHRoZSBzZWNvbmQgY2FwdHVyZSAgICAgICAKfQoKIyBTaW11bGF0aW5nIHRoZSBkYXRhCk4ubWFya2VkIDwtIHJlcChOQSwgTi5kcmF3KSAKICBmb3IoaSBpbiAxOk4uZHJhdykgewogICAgTi5tYXJrZWRbaV0gPC0gcGljay5iaWxscyhOLmJpbGxzW2ldKQogIH0KCiMgRmlsdGVyaW5nIG91dCBwYXJhbWV0ZXIgdmFsdWVzIE4ubWFya2VkIHRoYXQgd2VyZSBub3QgZXF1YWwgdG8gdwpwb3N0LmJpbGxzIDwtIE4uYmlsbHNbTi5tYXJrZWQgPT0gd10KCiMgUGxvdHRpbmcgdGhlIHBvc3RlcmlvciBkaXN0cmlidXRpb24KbGVuZ3RoKHBvc3QuYmlsbHMpCmJhcnBsb3QodGFibGUoY3V0KHBvc3QuYmlsbHMsIHNlcSh4LHVwcGVyLmxpbWl0LGJpbi53aWR0aCkpKSAvIGxlbmd0aChwb3N0LmJpbGxzKSwgY29sID0gImJsdWUiKQoKIyBTdGF0aXN0aWNzCm1pbihwb3N0LmJpbGxzKQptZWFuKHBvc3QuYmlsbHMpCm1lZGlhbihwb3N0LmJpbGxzKQptYXgocG9zdC5iaWxscykKYGBgCg==