Loading [MathJax]/jax/output/HTML-CSS/jax.js

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 θ, using uniform prior on θ.

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==