Marginal posterior predictive checking

Tail area probabilities of marginal predictive distributions, aka probability integral transformation (PIT).

Normal model for light speed 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)

Tail area probabilities of marginal predictive distributions, aka probability integral transformation (PIT)

Ty <- data.frame(x = pt((y - my)/(sqrt(1+1/n)*s), n-1))

#* Plot histogram of PIT values. Ideally histogram should be close to uniform.
title1 <- 'Light speed example
distribution of marginal posterior tail-values'
ggplot(data = Ty) +
  geom_histogram(aes(x = x), fill = 'steelblue',
                 color = 'black', binwidth = 0.05) +
  coord_cartesian(xlim = c(0, 1)) +
  labs(x = TeX('\\mathit{p}(\\mathit{y}^{\\mathrm{rep}}_{\\mathit{i}} < \\mathit{y_i} | \\mathit{y})'),
       y = '', title = title1) +
  scale_y_continuous(breaks=NULL)

Repeat the PIT checking after removing two "outliers"

y <- y[y>0]

Sufficient statistics

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

Tail area probabilities of marginal predictive distributions, aka probability integral transformation (PIT)

Ty <- data.frame(x = pt((y - my)/(sqrt(1+1/n)*s), n-1))

Plot histogram of PIT values. Ideally histogram should be close to uniform.

title1 <- 'Light speed example
distribution of marginal posterior tail-values'
ggplot(data = Ty) +
  geom_histogram(aes(x = x), fill = 'steelblue',
                 color = 'black', binwidth = 0.05) +
  coord_cartesian(xlim = c(0, 1)) +
  labs(x = TeX('\\mathit{p}(\\mathit{y}^{\\mathrm{rep}}_{\\mathit{i}} < \\mathit{y_i} | \\mathit{y})'),
       y = '', title = title1) +
  scale_y_continuous(breaks=NULL)

LS0tCnRpdGxlOiAiQmF5ZXNpYW4gZGF0YSBhbmFseXNpcyBkZW1vIDYuNCIKYXV0aG9yOiAiQWtpIFZlaHRhcmksIE1hcmt1cyBQYWFzaW5pZW1pIgpkYXRlOiAiYHIgZm9ybWF0KFN5cy5EYXRlKCkpYCIKb3V0cHV0OgogIGh0bWxfZG9jdW1lbnQ6CiAgICB0aGVtZTogcmVhZGFibGUKICAgIGNvZGVfZG93bmxvYWQ6IHRydWUKLS0tCiMjIE1hcmdpbmFsIHBvc3RlcmlvciBwcmVkaWN0aXZlIGNoZWNraW5nCgpUYWlsIGFyZWEgcHJvYmFiaWxpdGllcyBvZiBtYXJnaW5hbCBwcmVkaWN0aXZlIGRpc3RyaWJ1dGlvbnMsCmFrYSBwcm9iYWJpbGl0eSBpbnRlZ3JhbCB0cmFuc2Zvcm1hdGlvbiAoUElUKS4KCk5vcm1hbCBtb2RlbCBmb3IgbGlnaHQgc3BlZWQgZGF0YS4KCmdncGxvdDIgaXMgdXNlZCBmb3IgcGxvdHRpbmcsIHRpZHlyIGZvciBtYW5pcHVsYXRpbmcgZGF0YSBmcmFtZXMKCmBgYHtyIHNldHVwLCBtZXNzYWdlPUZBTFNFLCBlcnJvcj1GQUxTRSwgd2FybmluZz1GQUxTRX0KbGlicmFyeShnZ3Bsb3QyKQp0aGVtZV9zZXQodGhlbWVfbWluaW1hbCgpKQpsaWJyYXJ5KHRpZHlyKQpsaWJyYXJ5KGxhdGV4MmV4cCkKbGlicmFyeShycHJvanJvb3QpCnJvb3Q8LWhhc19maWxlKCIuQkRBX1JfZGVtb3Nfcm9vdCIpJG1ha2VfZml4X2ZpbGUoKQpgYGAKCkRhdGEKCmBgYHtyIH0KeSA8LSByZWFkLnRhYmxlKHJvb3QoImRlbW9zX2NoNiIsImxpZ2h0LnR4dCIpKSRWMQpgYGAKClN1ZmZpY2llbnQgc3RhdGlzdGljcwoKYGBge3IgfQpuIDwtIGxlbmd0aCh5KQpzIDwtIHNkKHkpCm15IDwtIG1lYW4oeSkKYGBgCgpUYWlsIGFyZWEgcHJvYmFiaWxpdGllcyBvZiBtYXJnaW5hbCBwcmVkaWN0aXZlIGRpc3RyaWJ1dGlvbnMsCmFrYSBwcm9iYWJpbGl0eSBpbnRlZ3JhbCB0cmFuc2Zvcm1hdGlvbiAoUElUKQoKYGBge3IgfQpUeSA8LSBkYXRhLmZyYW1lKHggPSBwdCgoeSAtIG15KS8oc3FydCgxKzEvbikqcyksIG4tMSkpCgojKiBQbG90IGhpc3RvZ3JhbSBvZiBQSVQgdmFsdWVzLiBJZGVhbGx5IGhpc3RvZ3JhbSBzaG91bGQgYmUgY2xvc2UgdG8gdW5pZm9ybS4KdGl0bGUxIDwtICdMaWdodCBzcGVlZCBleGFtcGxlCmRpc3RyaWJ1dGlvbiBvZiBtYXJnaW5hbCBwb3N0ZXJpb3IgdGFpbC12YWx1ZXMnCmdncGxvdChkYXRhID0gVHkpICsKICBnZW9tX2hpc3RvZ3JhbShhZXMoeCA9IHgpLCBmaWxsID0gJ3N0ZWVsYmx1ZScsCiAgICAgICAgICAgICAgICAgY29sb3IgPSAnYmxhY2snLCBiaW53aWR0aCA9IDAuMDUpICsKICBjb29yZF9jYXJ0ZXNpYW4oeGxpbSA9IGMoMCwgMSkpICsKICBsYWJzKHggPSBUZVgoJ1xcbWF0aGl0e3B9KFxcbWF0aGl0e3l9XntcXG1hdGhybXtyZXB9fV97XFxtYXRoaXR7aX19IDwgXFxtYXRoaXR7eV9pfSB8IFxcbWF0aGl0e3l9KScpLAogICAgICAgeSA9ICcnLCB0aXRsZSA9IHRpdGxlMSkgKwogIHNjYWxlX3lfY29udGludW91cyhicmVha3M9TlVMTCkKYGBgCgpSZXBlYXQgdGhlIFBJVCBjaGVja2luZyBhZnRlciByZW1vdmluZyB0d28gIm91dGxpZXJzIgoKYGBge3IgfQp5IDwtIHlbeT4wXQpgYGAKClN1ZmZpY2llbnQgc3RhdGlzdGljcwoKYGBge3IgfQpuIDwtIGxlbmd0aCh5KQpzIDwtIHNkKHkpCm15IDwtIG1lYW4oeSkKYGBgCgpUYWlsIGFyZWEgcHJvYmFiaWxpdGllcyBvZiBtYXJnaW5hbCBwcmVkaWN0aXZlIGRpc3RyaWJ1dGlvbnMsCmFrYSBwcm9iYWJpbGl0eSBpbnRlZ3JhbCB0cmFuc2Zvcm1hdGlvbiAoUElUKQoKYGBge3IgfQpUeSA8LSBkYXRhLmZyYW1lKHggPSBwdCgoeSAtIG15KS8oc3FydCgxKzEvbikqcyksIG4tMSkpCmBgYAoKUGxvdCBoaXN0b2dyYW0gb2YgUElUIHZhbHVlcy4gSWRlYWxseSBoaXN0b2dyYW0gc2hvdWxkIGJlIGNsb3NlIHRvIHVuaWZvcm0uCgpgYGB7ciB9CnRpdGxlMSA8LSAnTGlnaHQgc3BlZWQgZXhhbXBsZQpkaXN0cmlidXRpb24gb2YgbWFyZ2luYWwgcG9zdGVyaW9yIHRhaWwtdmFsdWVzJwpnZ3Bsb3QoZGF0YSA9IFR5KSArCiAgZ2VvbV9oaXN0b2dyYW0oYWVzKHggPSB4KSwgZmlsbCA9ICdzdGVlbGJsdWUnLAogICAgICAgICAgICAgICAgIGNvbG9yID0gJ2JsYWNrJywgYmlud2lkdGggPSAwLjA1KSArCiAgY29vcmRfY2FydGVzaWFuKHhsaW0gPSBjKDAsIDEpKSArCiAgbGFicyh4ID0gVGVYKCdcXG1hdGhpdHtwfShcXG1hdGhpdHt5fV57XFxtYXRocm17cmVwfX1fe1xcbWF0aGl0e2l9fSA8IFxcbWF0aGl0e3lfaX0gfCBcXG1hdGhpdHt5fSknKSwKICAgICAgIHkgPSAnJywgdGl0bGUgPSB0aXRsZTEpICsKICBzY2FsZV95X2NvbnRpbnVvdXMoYnJlYWtzPU5VTEwpCmBgYAoK