Posterior predictive checking of normal model for light data

ggplot2 is used for plotting, tidyr for manipulating data frames

library(ggplot2)
theme_set(theme_minimal())
library(tidyr)
library(latex2exp)
library(rprojroot)
root<-has_file(".BDA_R_demos_root")$make_fix_file()

Data

y <- read.table(root("demos_ch6","light.txt"))$V1

Sufficient statistics

n <- length(y)
s <- sd(y)
my <- mean(y)

Create 9 random replicate data sets from the posterior predictive density. Each set has same number of virtual observations as the original data set.

sampt <- replicate(9, rt(n, n-1)*sqrt(1+1/n)*s+my) %>%
  as.data.frame()

Replace one of the replicates with observed data. If you can spot which one has been replaced, it means that the replicates do not resemble the original data and thus the model has a defect

ind <- sample(9, 1)
sampt_y <- replace(sampt, ind, y) %>% setNames(1:9) %>% gather()
ggplot(data = sampt_y) +
  geom_histogram(aes(x = value), fill = 'steelblue',
                 color = 'black', binwidth = 4) +
  facet_wrap(~key, nrow = 3) +
  coord_cartesian(xlim = c(-55, 65)) +
  labs(x = '', y = '') +
  scale_y_continuous(breaks=NULL) +
  theme(strip.background = element_blank())

Generate 1000 replicate data sets and compute test statistic. The test statistic here is the smallest observation in the data or in the replicated data.

sampt1000 <- replicate(1000, rt(n, n-1)*sqrt(1+1/n)*s+my) %>%
  as.data.frame()
minvals <- data.frame(x = sapply(sampt1000, min))

Plot test statistic for the data and the replicated data sets

title1 <- 'Smallest observation in the replicated
data (hist.) vs in the original data (vertical line)'
ggplot(data = minvals) +
  geom_histogram(aes(x = x), fill = 'steelblue',
                 color = 'black', binwidth = 2) +
  geom_vline(aes(xintercept = min(x)), data = data.frame(x = y),
             color = 'red') +
  coord_cartesian(xlim = c(-50, 20)) +
  labs(x = TeX('Minimum of \\mathit{y} and \\mathit{y}^{\\mathrm{rep}}'),
       y = '', title = title1) +
  scale_y_continuous(breaks=NULL)

LS0tCnRpdGxlOiAiQmF5ZXNpYW4gZGF0YSBhbmFseXNpcyBkZW1vIDYuMSIKYXV0aG9yOiAiQWtpIFZlaHRhcmksIE1hcmt1cyBQYWFzaW5pZW1pIgpkYXRlOiAiYHIgZm9ybWF0KFN5cy5EYXRlKCkpYCIKb3V0cHV0OgogIGh0bWxfZG9jdW1lbnQ6CiAgICB0aGVtZTogcmVhZGFibGUKICAgIGNvZGVfZG93bmxvYWQ6IHRydWUKLS0tCiMjIFBvc3RlcmlvciBwcmVkaWN0aXZlIGNoZWNraW5nIG9mIG5vcm1hbCBtb2RlbCBmb3IgbGlnaHQgZGF0YQoKZ2dwbG90MiBpcyB1c2VkIGZvciBwbG90dGluZywgdGlkeXIgZm9yIG1hbmlwdWxhdGluZyBkYXRhIGZyYW1lcwoKYGBge3Igc2V0dXAsIG1lc3NhZ2U9RkFMU0UsIGVycm9yPUZBTFNFLCB3YXJuaW5nPUZBTFNFfQpsaWJyYXJ5KGdncGxvdDIpCnRoZW1lX3NldCh0aGVtZV9taW5pbWFsKCkpCmxpYnJhcnkodGlkeXIpCmxpYnJhcnkobGF0ZXgyZXhwKQpsaWJyYXJ5KHJwcm9qcm9vdCkKcm9vdDwtaGFzX2ZpbGUoIi5CREFfUl9kZW1vc19yb290IikkbWFrZV9maXhfZmlsZSgpCmBgYAoKRGF0YQoKYGBge3IgfQp5IDwtIHJlYWQudGFibGUocm9vdCgiZGVtb3NfY2g2IiwibGlnaHQudHh0IikpJFYxCmBgYAoKU3VmZmljaWVudCBzdGF0aXN0aWNzCgpgYGB7ciB9Cm4gPC0gbGVuZ3RoKHkpCnMgPC0gc2QoeSkKbXkgPC0gbWVhbih5KQpgYGAKCkNyZWF0ZSA5IHJhbmRvbSByZXBsaWNhdGUgZGF0YSBzZXRzIGZyb20gdGhlIHBvc3RlcmlvcgpwcmVkaWN0aXZlIGRlbnNpdHkuCkVhY2ggc2V0IGhhcyBzYW1lIG51bWJlciBvZiB2aXJ0dWFsIG9ic2VydmF0aW9ucyBhcyB0aGUKb3JpZ2luYWwgZGF0YSBzZXQuCgpgYGB7ciB9CnNhbXB0IDwtIHJlcGxpY2F0ZSg5LCBydChuLCBuLTEpKnNxcnQoMSsxL24pKnMrbXkpICU+JQogIGFzLmRhdGEuZnJhbWUoKQpgYGAKClJlcGxhY2Ugb25lIG9mIHRoZSByZXBsaWNhdGVzIHdpdGggb2JzZXJ2ZWQgZGF0YS4KSWYgeW91IGNhbiBzcG90IHdoaWNoIG9uZSBoYXMgYmVlbiByZXBsYWNlZCwgaXQgbWVhbnMKdGhhdCB0aGUgcmVwbGljYXRlcyBkbyBub3QgcmVzZW1ibGUgdGhlIG9yaWdpbmFsIGRhdGEKYW5kIHRodXMgdGhlIG1vZGVsIGhhcyBhIGRlZmVjdAoKYGBge3IgfQppbmQgPC0gc2FtcGxlKDksIDEpCnNhbXB0X3kgPC0gcmVwbGFjZShzYW1wdCwgaW5kLCB5KSAlPiUgc2V0TmFtZXMoMTo5KSAlPiUgZ2F0aGVyKCkKZ2dwbG90KGRhdGEgPSBzYW1wdF95KSArCiAgZ2VvbV9oaXN0b2dyYW0oYWVzKHggPSB2YWx1ZSksIGZpbGwgPSAnc3RlZWxibHVlJywKICAgICAgICAgICAgICAgICBjb2xvciA9ICdibGFjaycsIGJpbndpZHRoID0gNCkgKwogIGZhY2V0X3dyYXAofmtleSwgbnJvdyA9IDMpICsKICBjb29yZF9jYXJ0ZXNpYW4oeGxpbSA9IGMoLTU1LCA2NSkpICsKICBsYWJzKHggPSAnJywgeSA9ICcnKSArCiAgc2NhbGVfeV9jb250aW51b3VzKGJyZWFrcz1OVUxMKSArCiAgdGhlbWUoc3RyaXAuYmFja2dyb3VuZCA9IGVsZW1lbnRfYmxhbmsoKSkKYGBgCgpHZW5lcmF0ZSAxMDAwIHJlcGxpY2F0ZSBkYXRhIHNldHMgYW5kIGNvbXB1dGUgdGVzdCBzdGF0aXN0aWMuIFRoZQp0ZXN0IHN0YXRpc3RpYyBoZXJlIGlzIHRoZSBzbWFsbGVzdCBvYnNlcnZhdGlvbiBpbiB0aGUgZGF0YSBvciBpbgp0aGUgcmVwbGljYXRlZCBkYXRhLgoKYGBge3IgfQpzYW1wdDEwMDAgPC0gcmVwbGljYXRlKDEwMDAsIHJ0KG4sIG4tMSkqc3FydCgxKzEvbikqcytteSkgJT4lCiAgYXMuZGF0YS5mcmFtZSgpCm1pbnZhbHMgPC0gZGF0YS5mcmFtZSh4ID0gc2FwcGx5KHNhbXB0MTAwMCwgbWluKSkKYGBgCgpQbG90IHRlc3Qgc3RhdGlzdGljIGZvciB0aGUgZGF0YSBhbmQgdGhlIHJlcGxpY2F0ZWQgZGF0YSBzZXRzCgpgYGB7ciB9CnRpdGxlMSA8LSAnU21hbGxlc3Qgb2JzZXJ2YXRpb24gaW4gdGhlIHJlcGxpY2F0ZWQKZGF0YSAoaGlzdC4pIHZzIGluIHRoZSBvcmlnaW5hbCBkYXRhICh2ZXJ0aWNhbCBsaW5lKScKZ2dwbG90KGRhdGEgPSBtaW52YWxzKSArCiAgZ2VvbV9oaXN0b2dyYW0oYWVzKHggPSB4KSwgZmlsbCA9ICdzdGVlbGJsdWUnLAogICAgICAgICAgICAgICAgIGNvbG9yID0gJ2JsYWNrJywgYmlud2lkdGggPSAyKSArCiAgZ2VvbV92bGluZShhZXMoeGludGVyY2VwdCA9IG1pbih4KSksIGRhdGEgPSBkYXRhLmZyYW1lKHggPSB5KSwKICAgICAgICAgICAgIGNvbG9yID0gJ3JlZCcpICsKICBjb29yZF9jYXJ0ZXNpYW4oeGxpbSA9IGMoLTUwLCAyMCkpICsKICBsYWJzKHggPSBUZVgoJ01pbmltdW0gb2YgXFxtYXRoaXR7eX0gYW5kIFxcbWF0aGl0e3l9XntcXG1hdGhybXtyZXB9fScpLAogICAgICAgeSA9ICcnLCB0aXRsZSA9IHRpdGxlMSkgKwogIHNjYWxlX3lfY29udGludW91cyhicmVha3M9TlVMTCkKYGBgCgo=