## Probability of a girl birth given placenta previa (BDA3 p. 37).

Calculate the posterior distribution on a discrete grid of points by multiplying the likelihood and a non-conjugate prior at each point, and normalizing over the points. Simulate samples from the resulting non-standard posterior distribution using inverse cdf using the discrete grid.

ggplot2 and gridExtra are used for plotting, tidyr for manipulating data frames

library(ggplot2)
theme_set(theme_minimal())
library(gridExtra)
library(tidyr)

### Evaluating posterior with non-conjugate prior in grid

Posterior with observations (437,543) and uniform prior (Beta(1,1))

a <- 437
b <- 543

Evaluate densities at evenly spaced points between 0.1 and 1

df1 <- data.frame(theta = seq(0.1, 1, 0.001))
df1$con <- dbeta(df1$theta, a, b)

Compute the density of non-conjugate prior in discrete points, i.e. in a grid this non-conjugate prior is the same as in figure 2.4 in the book

pp <- rep(1, nrow(df1))
pi <- sapply(c(0.388, 0.488, 0.588), function(pi) which(df1$theta == pi)) pm <- 11 pp[pi[1]:pi[2]] <- seq(1, pm, length.out = length(pi[1]:pi[2])) pp[pi[3]:pi[2]] <- seq(1, pm, length.out = length(pi[3]:pi[2])) normalize the prior df1$nc_p <- pp / sum(pp)

compute the un-normalized non-conjugate posterior in a grid

po <- dbeta(df1$theta, a, b) * pp normalize the posterior df1$nc_po <- po / sum(po)

Plot posterior with uniform prior, non-conjugate prior and the corresponding non-conjugate posterior

# gather the data frame into key-value pairs
# and change variable names for plotting
df2 <- gather(df1, grp, p, -theta, factor_key = T) #%>%
levels(df2$grp) <- c('Posterior with uniform prior', 'Non-conjugate prior', 'Non-conjugate posterior') ggplot(data = df2) + geom_line(aes(theta, p)) + facet_wrap(~grp, ncol = 1, scales = 'free_y') + coord_cartesian(xlim = c(0.35,0.6)) + scale_y_continuous(breaks=NULL) + labs(x = '', y = '') ### Inverse cdf sampling compute the cumulative density in a grid df1$cs_po <- cumsum(df1$nc_po) Sample from uniform distribution U(0,1) # set.seed(seed) is used to set seed for the randon number generator set.seed(2601) # runif(k) returns k uniform random numbers from interval [0,1] r <- runif(10000) Inverse-cdf sampling # function to find the value smallest value theta at which the cumulative # sum of the posterior densities is greater than r. invcdf <- function(r, df) df$theta[sum(df\$cs_po < r) + 1]
# sapply function for each sample r. The returned values s are now
# random draws from the distribution.
s <- sapply(r, invcdf, df1)

Create three plots: p1 is the posterior, p2 is the cdf of the posterior and p3 is the histogram of posterior samples (drawn using inv-cdf)

p1 <- ggplot(data = df1) +
geom_line(aes(theta, nc_po)) +
coord_cartesian(xlim = c(0.35, 0.6)) +
labs(title = 'Non-conjugate posterior', x = '', y = '') +
scale_y_continuous(breaks = NULL)
p2 <- ggplot(data = df1) +
geom_line(aes(theta, cs_po)) +
coord_cartesian(xlim = c(0.35, 0.6)) +
labs(title = 'Posterior-cdf', x = '', y = '') +
scale_y_continuous(breaks = NULL)
p3 <- ggplot() +
geom_histogram(aes(s), binwidth = 0.003) +
coord_cartesian(xlim = c(0.35, 0.6)) +
labs(title = 'Histogram of posterior samples', x = '', y = '') +
scale_y_continuous(breaks = NULL)
# combine the plots
grid.arrange(p1, p2, p3)

LS0tCnRpdGxlOiAiQmF5ZXNpYW4gZGF0YSBhbmFseXNpcyBkZW1vIDIuNCIKYXV0aG9yOiAiQWtpIFZlaHRhcmksIE1hcmt1cyBQYWFzaW5pZW1pIgpkYXRlOiAiYHIgZm9ybWF0KFN5cy5EYXRlKCkpYCIKb3V0cHV0OgogIGh0bWxfZG9jdW1lbnQ6CiAgICB0aGVtZTogcmVhZGFibGUKICAgIGNvZGVfZG93bmxvYWQ6IHRydWUKLS0tCiMjIFByb2JhYmlsaXR5IG9mIGEgZ2lybCBiaXJ0aCBnaXZlbiBwbGFjZW50YSBwcmV2aWEgKEJEQTMgcC4gMzcpLgoKQ2FsY3VsYXRlIHRoZSBwb3N0ZXJpb3IgZGlzdHJpYnV0aW9uIG9uIGEgZGlzY3JldGUgZ3JpZCBvZiBwb2ludHMgYnkKbXVsdGlwbHlpbmcgdGhlIGxpa2VsaWhvb2QgYW5kIGEgbm9uLWNvbmp1Z2F0ZSBwcmlvciBhdCBlYWNoIHBvaW50LAphbmQgbm9ybWFsaXppbmcgb3ZlciB0aGUgcG9pbnRzLiBTaW11bGF0ZSBzYW1wbGVzIGZyb20gdGhlIHJlc3VsdGluZwpub24tc3RhbmRhcmQgcG9zdGVyaW9yIGRpc3RyaWJ1dGlvbiB1c2luZyBpbnZlcnNlIGNkZiB1c2luZyB0aGUKZGlzY3JldGUgZ3JpZC4KCmdncGxvdDIgYW5kIGdyaWRFeHRyYSBhcmUgdXNlZCBmb3IgcGxvdHRpbmcsIHRpZHlyIGZvciBtYW5pcHVsYXRpbmcgZGF0YSBmcmFtZXMKCmBgYHtyIHNldHVwLCBtZXNzYWdlPUZBTFNFLCBlcnJvcj1GQUxTRSwgd2FybmluZz1GQUxTRX0KbGlicmFyeShnZ3Bsb3QyKQp0aGVtZV9zZXQodGhlbWVfbWluaW1hbCgpKQpsaWJyYXJ5KGdyaWRFeHRyYSkKbGlicmFyeSh0aWR5cikKYGBgCgojIyMgRXZhbHVhdGluZyBwb3N0ZXJpb3Igd2l0aCBub24tY29uanVnYXRlIHByaW9yIGluIGdyaWQKClBvc3RlcmlvciB3aXRoIG9ic2VydmF0aW9ucyAoNDM3LDU0MykgYW5kIHVuaWZvcm0gcHJpb3IgKEJldGEoMSwxKSkKCmBgYHtyIH0KYSA8LSA0MzcKYiA8LSA1NDMKYGBgCgpFdmFsdWF0ZSBkZW5zaXRpZXMgYXQgZXZlbmx5IHNwYWNlZCBwb2ludHMgYmV0d2VlbiAwLjEgYW5kIDEKCmBgYHtyIH0KZGYxIDwtIGRhdGEuZnJhbWUodGhldGEgPSBzZXEoMC4xLCAxLCAwLjAwMSkpCmRmMSRjb24gPC0gZGJldGEoZGYxJHRoZXRhLCBhLCBiKQpgYGAKCkNvbXB1dGUgdGhlIGRlbnNpdHkgb2Ygbm9uLWNvbmp1Z2F0ZSBwcmlvciBpbiBkaXNjcmV0ZSBwb2ludHMsIGkuZS4gaW4gYSBncmlkCnRoaXMgbm9uLWNvbmp1Z2F0ZSBwcmlvciBpcyB0aGUgc2FtZSBhcyBpbiBmaWd1cmUgMi40IGluIHRoZSBib29rCgpgYGB7ciB9CnBwIDwtIHJlcCgxLCBucm93KGRmMSkpCnBpIDwtIHNhcHBseShjKDAuMzg4LCAwLjQ4OCwgMC41ODgpLCBmdW5jdGlvbihwaSkgd2hpY2goZGYxJHRoZXRhID09IHBpKSkKcG0gPC0gMTEKcHBbcGlbMV06cGlbMl1dIDwtIHNlcSgxLCBwbSwgbGVuZ3RoLm91dCA9IGxlbmd0aChwaVsxXTpwaVsyXSkpCnBwW3BpWzNdOnBpWzJdXSA8LSBzZXEoMSwgcG0sIGxlbmd0aC5vdXQgPSBsZW5ndGgocGlbM106cGlbMl0pKQpgYGAKCm5vcm1hbGl6ZSB0aGUgcHJpb3IKCmBgYHtyIH0KZGYxJG5jX3AgPC0gcHAgLyBzdW0ocHApCmBgYAoKY29tcHV0ZSB0aGUgdW4tbm9ybWFsaXplZCBub24tY29uanVnYXRlIHBvc3RlcmlvciBpbiBhIGdyaWQKCmBgYHtyIH0KcG8gPC0gZGJldGEoZGYxJHRoZXRhLCBhLCBiKSAqIHBwCmBgYAoKbm9ybWFsaXplIHRoZSBwb3N0ZXJpb3IKCmBgYHtyIH0KZGYxJG5jX3BvIDwtIHBvIC8gc3VtKHBvKQpgYGAKClBsb3QgcG9zdGVyaW9yIHdpdGggdW5pZm9ybSBwcmlvciwgbm9uLWNvbmp1Z2F0ZQpwcmlvciBhbmQgdGhlIGNvcnJlc3BvbmRpbmcgbm9uLWNvbmp1Z2F0ZSBwb3N0ZXJpb3IKCmBgYHtyIH0KIyBnYXRoZXIgdGhlIGRhdGEgZnJhbWUgaW50byBrZXktdmFsdWUgcGFpcnMKIyBhbmQgY2hhbmdlIHZhcmlhYmxlIG5hbWVzIGZvciBwbG90dGluZwpkZjIgPC0gZ2F0aGVyKGRmMSwgZ3JwLCBwLCAtdGhldGEsIGZhY3Rvcl9rZXkgPSBUKSAjJT4lCmxldmVscyhkZjIkZ3JwKSA8LSBjKCdQb3N0ZXJpb3Igd2l0aCB1bmlmb3JtIHByaW9yJywKICAgICAgICAgICAgICAgICAgICAgJ05vbi1jb25qdWdhdGUgcHJpb3InLCAnTm9uLWNvbmp1Z2F0ZSBwb3N0ZXJpb3InKQpnZ3Bsb3QoZGF0YSA9IGRmMikgKwogIGdlb21fbGluZShhZXModGhldGEsIHApKSArCiAgZmFjZXRfd3JhcCh+Z3JwLCBuY29sID0gMSwgc2NhbGVzID0gJ2ZyZWVfeScpICsKICBjb29yZF9jYXJ0ZXNpYW4oeGxpbSA9IGMoMC4zNSwwLjYpKSArCiAgc2NhbGVfeV9jb250aW51b3VzKGJyZWFrcz1OVUxMKSArCiAgbGFicyh4ID0gJycsIHkgPSAnJykKYGBgCgojIyMgSW52ZXJzZSBjZGYgc2FtcGxpbmcKCmNvbXB1dGUgdGhlIGN1bXVsYXRpdmUgZGVuc2l0eSBpbiBhIGdyaWQKCmBgYHtyIH0KZGYxJGNzX3BvIDwtIGN1bXN1bShkZjEkbmNfcG8pCmBgYAoKU2FtcGxlIGZyb20gdW5pZm9ybSBkaXN0cmlidXRpb24gVSgwLDEpCgpgYGB7ciB9CiMgc2V0LnNlZWQoc2VlZCkgaXMgdXNlZCB0byBzZXQgc2VlZCBmb3IgdGhlIHJhbmRvbiBudW1iZXIgZ2VuZXJhdG9yCnNldC5zZWVkKDI2MDEpCiMgcnVuaWYoaykgcmV0dXJucyBrIHVuaWZvcm0gcmFuZG9tIG51bWJlcnMgZnJvbSBpbnRlcnZhbCBbMCwxXQpyIDwtIHJ1bmlmKDEwMDAwKQpgYGAKCkludmVyc2UtY2RmIHNhbXBsaW5nCgpgYGB7ciB9CiMgZnVuY3Rpb24gdG8gZmluZCB0aGUgdmFsdWUgc21hbGxlc3QgdmFsdWUgdGhldGEgYXQgd2hpY2ggdGhlIGN1bXVsYXRpdmUKIyBzdW0gb2YgdGhlIHBvc3RlcmlvciBkZW5zaXRpZXMgaXMgZ3JlYXRlciB0aGFuIHIuCmludmNkZiA8LSBmdW5jdGlvbihyLCBkZikgZGYkdGhldGFbc3VtKGRmJGNzX3BvIDwgcikgKyAxXQojIHNhcHBseSBmdW5jdGlvbiBmb3IgZWFjaCBzYW1wbGUgci4gVGhlIHJldHVybmVkIHZhbHVlcyBzIGFyZSBub3cKIyByYW5kb20gZHJhd3MgZnJvbSB0aGUgZGlzdHJpYnV0aW9uLgpzIDwtIHNhcHBseShyLCBpbnZjZGYsIGRmMSkKYGBgCgpDcmVhdGUgdGhyZWUgcGxvdHM6IHAxIGlzIHRoZSBwb3N0ZXJpb3IsIHAyIGlzIHRoZSBjZGYgb2YgdGhlIHBvc3RlcmlvcgphbmQgcDMgaXMgdGhlIGhpc3RvZ3JhbSBvZiBwb3N0ZXJpb3Igc2FtcGxlcyAoZHJhd24gdXNpbmcgaW52LWNkZikKCmBgYHtyIH0KcDEgPC0gZ2dwbG90KGRhdGEgPSBkZjEpICsKICBnZW9tX2xpbmUoYWVzKHRoZXRhLCBuY19wbykpICsKICBjb29yZF9jYXJ0ZXNpYW4oeGxpbSA9IGMoMC4zNSwgMC42KSkgKwogIGxhYnModGl0bGUgPSAnTm9uLWNvbmp1Z2F0ZSBwb3N0ZXJpb3InLCB4ID0gJycsIHkgPSAnJykgKwogIHNjYWxlX3lfY29udGludW91cyhicmVha3MgPSBOVUxMKQpwMiA8LSBnZ3Bsb3QoZGF0YSA9IGRmMSkgKwogIGdlb21fbGluZShhZXModGhldGEsIGNzX3BvKSkgKwogIGNvb3JkX2NhcnRlc2lhbih4bGltID0gYygwLjM1LCAwLjYpKSArCiAgbGFicyh0aXRsZSA9ICdQb3N0ZXJpb3ItY2RmJywgeCA9ICcnLCB5ID0gJycpICsKICBzY2FsZV95X2NvbnRpbnVvdXMoYnJlYWtzID0gTlVMTCkKcDMgPC0gZ2dwbG90KCkgKwogIGdlb21faGlzdG9ncmFtKGFlcyhzKSwgYmlud2lkdGggPSAwLjAwMykgKwogIGNvb3JkX2NhcnRlc2lhbih4bGltID0gYygwLjM1LCAwLjYpKSArCiAgbGFicyh0aXRsZSA9ICdIaXN0b2dyYW0gb2YgcG9zdGVyaW9yIHNhbXBsZXMnLCB4ID0gJycsIHkgPSAnJykgKwogIHNjYWxlX3lfY29udGludW91cyhicmVha3MgPSBOVUxMKQojIGNvbWJpbmUgdGhlIHBsb3RzCmdyaWQuYXJyYW5nZShwMSwgcDIsIHAzKQpgYGAKCg==