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

437 girls and 543 boys have been observed. Calculate and plot the posterior distribution of the proportion of girls $$\theta$$, using uniform prior on $$\theta$$.

ggplot2 is used for plotting to install new packages, type e.g. install.packages(‘ggplot2’)

library(ggplot2)
theme_set(theme_minimal())

Posterior is Beta(438,544)

# seq creates evenly spaced values
df1 <- data.frame(theta = seq(0.375, 0.525, 0.001))
a <- 438
b <- 544
# dbeta computes the posterior density
df1$p <- dbeta(df1$theta, a, b)

compute also 95% central interval

# seq creates evenly spaced values from 2.5% quantile
# to 97.5% quantile (i.e., 95% central interval)
# qbeta computes the value for a given quantile given parameters a and b
df2 <- data.frame(theta = seq(qbeta(0.025, a, b), qbeta(0.975, a, b), length.out = 100))
# compute the posterior density
df2$p <- dbeta(df2$theta, a, b)

Plot posterior (Beta(438,544)) and 48.8% line for population average

ggplot(mapping = aes(theta, p)) +
geom_line(data = df1) +
# Add a layer of colorized 95% posterior interval
geom_area(data = df2, aes(fill='1')) +
# Add the proportion of girl babies in general population
geom_vline(xintercept = 0.488, linetype='dotted') +
# Decorate the plot a little
labs(title='Uniform prior -> Posterior is Beta(438,544)', y = '') +
scale_y_continuous(expand = c(0, 0.1), breaks = NULL) +
scale_fill_manual(values = 'lightblue', labels = '95% posterior interval') +
theme(legend.position = 'bottom', legend.title = element_blank())

LS0tCnRpdGxlOiAiQmF5ZXNpYW4gZGF0YSBhbmFseXNpcyBkZW1vIDIuMSIKYXV0aG9yOiAiQWtpIFZlaHRhcmksIE1hcmt1cyBQYWFzaW5pZW1pIgpkYXRlOiAiYHIgZm9ybWF0KFN5cy5EYXRlKCkpYCIKb3V0cHV0OgogIGh0bWxfZG9jdW1lbnQ6CiAgICB0aGVtZTogcmVhZGFibGUKICAgIGNvZGVfZG93bmxvYWQ6IHRydWUKLS0tCiMjIFByb2JhYmlsaXR5IG9mIGEgZ2lybCBiaXJ0aCBnaXZlbiBwbGFjZW50YSBwcmV2aWEgKEJEQTMgcC4gMzcpLgoKNDM3IGdpcmxzIGFuZCA1NDMgYm95cyBoYXZlIGJlZW4gb2JzZXJ2ZWQuIENhbGN1bGF0ZSBhbmQgcGxvdCB0aGUKcG9zdGVyaW9yIGRpc3RyaWJ1dGlvbiBvZiB0aGUgcHJvcG9ydGlvbiBvZiBnaXJscyAkXHRoZXRhJCwgdXNpbmcKdW5pZm9ybSBwcmlvciBvbiAkXHRoZXRhJC4KCmdncGxvdDIgaXMgdXNlZCBmb3IgcGxvdHRpbmcKdG8gaW5zdGFsbCBuZXcgcGFja2FnZXMsIHR5cGUgZS5nLiBpbnN0YWxsLnBhY2thZ2VzKCdnZ3Bsb3QyJykKCmBgYHtyIHNldHVwLCBtZXNzYWdlPUZBTFNFLCBlcnJvcj1GQUxTRSwgd2FybmluZz1GQUxTRX0KbGlicmFyeShnZ3Bsb3QyKQp0aGVtZV9zZXQodGhlbWVfbWluaW1hbCgpKQpgYGAKClBvc3RlcmlvciBpcyBCZXRhKDQzOCw1NDQpCgpgYGB7ciB9CiMgc2VxIGNyZWF0ZXMgZXZlbmx5IHNwYWNlZCB2YWx1ZXMKZGYxIDwtIGRhdGEuZnJhbWUodGhldGEgPSBzZXEoMC4zNzUsIDAuNTI1LCAwLjAwMSkpIAphIDwtIDQzOApiIDwtIDU0NAojIGRiZXRhIGNvbXB1dGVzIHRoZSBwb3N0ZXJpb3IgZGVuc2l0eQpkZjEkcCA8LSBkYmV0YShkZjEkdGhldGEsIGEsIGIpCmBgYAoKY29tcHV0ZSBhbHNvIDk1JSBjZW50cmFsIGludGVydmFsCgpgYGB7ciB9CiMgc2VxIGNyZWF0ZXMgZXZlbmx5IHNwYWNlZCB2YWx1ZXMgZnJvbSAyLjUlIHF1YW50aWxlCiMgdG8gOTcuNSUgcXVhbnRpbGUgKGkuZS4sIDk1JSBjZW50cmFsIGludGVydmFsKQojIHFiZXRhIGNvbXB1dGVzIHRoZSB2YWx1ZSBmb3IgYSBnaXZlbiBxdWFudGlsZSBnaXZlbiBwYXJhbWV0ZXJzIGEgYW5kIGIKZGYyIDwtIGRhdGEuZnJhbWUodGhldGEgPSBzZXEocWJldGEoMC4wMjUsIGEsIGIpLCBxYmV0YSgwLjk3NSwgYSwgYiksIGxlbmd0aC5vdXQgPSAxMDApKQojIGNvbXB1dGUgdGhlIHBvc3RlcmlvciBkZW5zaXR5CmRmMiRwIDwtIGRiZXRhKGRmMiR0aGV0YSwgYSwgYikKYGBgCgpQbG90IHBvc3RlcmlvciAoQmV0YSg0MzgsNTQ0KSkKYW5kIDQ4LjglIGxpbmUgZm9yIHBvcHVsYXRpb24gYXZlcmFnZQoKYGBge3IgfQpnZ3Bsb3QobWFwcGluZyA9IGFlcyh0aGV0YSwgcCkpICsKICBnZW9tX2xpbmUoZGF0YSA9IGRmMSkgKwogICMgQWRkIGEgbGF5ZXIgb2YgY29sb3JpemVkIDk1JSBwb3N0ZXJpb3IgaW50ZXJ2YWwKICBnZW9tX2FyZWEoZGF0YSA9IGRmMiwgYWVzKGZpbGw9JzEnKSkgKwogICMgQWRkIHRoZSBwcm9wb3J0aW9uIG9mIGdpcmwgYmFiaWVzIGluIGdlbmVyYWwgcG9wdWxhdGlvbgogIGdlb21fdmxpbmUoeGludGVyY2VwdCA9IDAuNDg4LCBsaW5ldHlwZT0nZG90dGVkJykgKwogICMgRGVjb3JhdGUgdGhlIHBsb3QgYSBsaXR0bGUKICBsYWJzKHRpdGxlPSdVbmlmb3JtIHByaW9yIC0+IFBvc3RlcmlvciBpcyBCZXRhKDQzOCw1NDQpJywgeSA9ICcnKSArCiAgc2NhbGVfeV9jb250aW51b3VzKGV4cGFuZCA9IGMoMCwgMC4xKSwgYnJlYWtzID0gTlVMTCkgKwogIHNjYWxlX2ZpbGxfbWFudWFsKHZhbHVlcyA9ICdsaWdodGJsdWUnLCBsYWJlbHMgPSAnOTUlIHBvc3RlcmlvciBpbnRlcnZhbCcpICsKICB0aGVtZShsZWdlbmQucG9zaXRpb24gPSAnYm90dG9tJywgbGVnZW5kLnRpdGxlID0gZWxlbWVudF9ibGFuaygpKQpgYGAKCg==