Estimating the speed of light using normal model BDA3 p. 66

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

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

Data

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

Sufficient statistics

n <- length(y)
s2 <- var(y)
my <- mean(y)

Positive values only

y_pos <- y[y > 0]

Sufficient statistics

n_pos <- length(y_pos)
s2_pos <- var(y_pos)
my_pos <- mean(y_pos)

Compute the density of mu in these points

tl1 <- c(18, 34)
df1 <- data.frame(t1 = seq(tl1[1], tl1[2], length.out = 1000))

Compute the exact marginal density for mu

# multiplication by 1./sqrt(s2/n) is due to the transformation of variable
# z=(x-mean(y))/sqrt(s2/n), see BDA3 p. 21
df1$pm_mu <- dt((df1$t1 - my) / sqrt(s2/n), n-1) / sqrt(s2/n)

Compute the exact marginal density for mu for the positive data

df1$pm_mu_pos = dt((df1$t1 - my_pos) / sqrt(s2_pos/n_pos), n_pos-1) / sqrt(s2_pos/n_pos)

Create a histogram of the measurements

p1 <- ggplot() +
  geom_histogram(aes(y), binwidth = 2, fill = 'steelblue', color = 'black') +
  coord_cartesian(xlim = c(-42, 42)) +
  scale_y_continuous(breaks = NULL) +
  labs(title = 'Newcomb\'s measurements', x = 'y')

Create a plot of the normal model

# gather the data points into key-value pairs
df2 <- gather(df1, grp, p, -t1)
# legend labels
labs2 <- c('Posterior of mu', 'Posterior of mu given y > 0', 'Modern estimate')
p2 <- ggplot(data = df2) +
  geom_line(aes(t1, p, color = grp)) +
  geom_vline(aes(xintercept = 33, color = 'q'),
             linetype = 'dashed', show.legend = F) +
  coord_cartesian(xlim = c(-42, 42)) +
  scale_y_continuous(breaks = NULL) +
  labs(title = 'Normal model', x = 'mu', y = '') +
  scale_color_manual(values = c('blue', 'darkgreen', 'black')) +
  guides(color=FALSE) +
  annotate("text", x=24, y=0.26, label=labs2[1], hjust="right", size=5) +
  annotate("text", x=26, y=0.58, label=labs2[2], hjust="right", size=5) +
  annotate("text", x=32, y=0.7, label=labs2[3], hjust="right", size=5)

Combine the plots

grid.arrange(p1, p2)

LS0tCnRpdGxlOiAiQmF5ZXNpYW4gZGF0YSBhbmFseXNpcyBkZW1vIDMuNSIKYXV0aG9yOiAiQWtpIFZlaHRhcmksIE1hcmt1cyBQYWFzaW5pZW1pIgpkYXRlOiAiYHIgZm9ybWF0KFN5cy5EYXRlKCkpYCIKb3V0cHV0OgogIGh0bWxfZG9jdW1lbnQ6CiAgICB0aGVtZTogcmVhZGFibGUKICAgIGNvZGVfZG93bmxvYWQ6IHRydWUKLS0tCiMjIEVzdGltYXRpbmcgdGhlIHNwZWVkIG9mIGxpZ2h0IHVzaW5nIG5vcm1hbCBtb2RlbCBCREEzIHAuIDY2CgpnZ3Bsb3QyLCBhbmQgZ3JpZEV4dHJhIGFyZSB1c2VkIGZvciBwbG90dGluZywgdGlkeXIgZm9yCm1hbmlwdWxhdGluZyBkYXRhIGZyYW1lcwoKYGBge3Igc2V0dXAsIG1lc3NhZ2U9RkFMU0UsIGVycm9yPUZBTFNFLCB3YXJuaW5nPUZBTFNFfQpsaWJyYXJ5KGdncGxvdDIpCnRoZW1lX3NldCh0aGVtZV9taW5pbWFsKCkpCmxpYnJhcnkoZ3JpZEV4dHJhKQpsaWJyYXJ5KHRpZHlyKQpsaWJyYXJ5KHJwcm9qcm9vdCkKcm9vdDwtaGFzX2ZpbGUoIi5CREFfUl9kZW1vc19yb290IikkbWFrZV9maXhfZmlsZSgpCmBgYAoKRGF0YQoKYGBge3IgfQp5IDwtIHJlYWQudGFibGUocm9vdCgiZGVtb3NfY2gzIiwibGlnaHQudHh0IikpJFYxCmBgYAoKU3VmZmljaWVudCBzdGF0aXN0aWNzCgpgYGB7ciB9Cm4gPC0gbGVuZ3RoKHkpCnMyIDwtIHZhcih5KQpteSA8LSBtZWFuKHkpCmBgYAoKUG9zaXRpdmUgdmFsdWVzIG9ubHkKCmBgYHtyIH0KeV9wb3MgPC0geVt5ID4gMF0KYGBgCgpTdWZmaWNpZW50IHN0YXRpc3RpY3MKCmBgYHtyIH0Kbl9wb3MgPC0gbGVuZ3RoKHlfcG9zKQpzMl9wb3MgPC0gdmFyKHlfcG9zKQpteV9wb3MgPC0gbWVhbih5X3BvcykKYGBgCgpDb21wdXRlIHRoZSBkZW5zaXR5IG9mIG11IGluIHRoZXNlIHBvaW50cwoKYGBge3IgfQp0bDEgPC0gYygxOCwgMzQpCmRmMSA8LSBkYXRhLmZyYW1lKHQxID0gc2VxKHRsMVsxXSwgdGwxWzJdLCBsZW5ndGgub3V0ID0gMTAwMCkpCmBgYAoKQ29tcHV0ZSB0aGUgZXhhY3QgbWFyZ2luYWwgZGVuc2l0eSBmb3IgbXUKCmBgYHtyIH0KIyBtdWx0aXBsaWNhdGlvbiBieSAxLi9zcXJ0KHMyL24pIGlzIGR1ZSB0byB0aGUgdHJhbnNmb3JtYXRpb24gb2YgdmFyaWFibGUKIyB6PSh4LW1lYW4oeSkpL3NxcnQoczIvbiksIHNlZSBCREEzIHAuIDIxCmRmMSRwbV9tdSA8LSBkdCgoZGYxJHQxIC0gbXkpIC8gc3FydChzMi9uKSwgbi0xKSAvIHNxcnQoczIvbikKYGBgCgpDb21wdXRlIHRoZSBleGFjdCBtYXJnaW5hbCBkZW5zaXR5IGZvciBtdSBmb3IgdGhlIHBvc2l0aXZlIGRhdGEKCmBgYHtyIH0KZGYxJHBtX211X3BvcyA9IGR0KChkZjEkdDEgLSBteV9wb3MpIC8gc3FydChzMl9wb3Mvbl9wb3MpLCBuX3Bvcy0xKSAvIHNxcnQoczJfcG9zL25fcG9zKQpgYGAKCkNyZWF0ZSBhIGhpc3RvZ3JhbSBvZiB0aGUgbWVhc3VyZW1lbnRzCgpgYGB7ciB9CnAxIDwtIGdncGxvdCgpICsKICBnZW9tX2hpc3RvZ3JhbShhZXMoeSksIGJpbndpZHRoID0gMiwgZmlsbCA9ICdzdGVlbGJsdWUnLCBjb2xvciA9ICdibGFjaycpICsKICBjb29yZF9jYXJ0ZXNpYW4oeGxpbSA9IGMoLTQyLCA0MikpICsKICBzY2FsZV95X2NvbnRpbnVvdXMoYnJlYWtzID0gTlVMTCkgKwogIGxhYnModGl0bGUgPSAnTmV3Y29tYlwncyBtZWFzdXJlbWVudHMnLCB4ID0gJ3knKQpgYGAKCkNyZWF0ZSBhIHBsb3Qgb2YgdGhlIG5vcm1hbCBtb2RlbAoKYGBge3IgfQojIGdhdGhlciB0aGUgZGF0YSBwb2ludHMgaW50byBrZXktdmFsdWUgcGFpcnMKZGYyIDwtIGdhdGhlcihkZjEsIGdycCwgcCwgLXQxKQojIGxlZ2VuZCBsYWJlbHMKbGFiczIgPC0gYygnUG9zdGVyaW9yIG9mIG11JywgJ1Bvc3RlcmlvciBvZiBtdSBnaXZlbiB5ID4gMCcsICdNb2Rlcm4gZXN0aW1hdGUnKQpwMiA8LSBnZ3Bsb3QoZGF0YSA9IGRmMikgKwogIGdlb21fbGluZShhZXModDEsIHAsIGNvbG9yID0gZ3JwKSkgKwogIGdlb21fdmxpbmUoYWVzKHhpbnRlcmNlcHQgPSAzMywgY29sb3IgPSAncScpLAogICAgICAgICAgICAgbGluZXR5cGUgPSAnZGFzaGVkJywgc2hvdy5sZWdlbmQgPSBGKSArCiAgY29vcmRfY2FydGVzaWFuKHhsaW0gPSBjKC00MiwgNDIpKSArCiAgc2NhbGVfeV9jb250aW51b3VzKGJyZWFrcyA9IE5VTEwpICsKICBsYWJzKHRpdGxlID0gJ05vcm1hbCBtb2RlbCcsIHggPSAnbXUnLCB5ID0gJycpICsKICBzY2FsZV9jb2xvcl9tYW51YWwodmFsdWVzID0gYygnYmx1ZScsICdkYXJrZ3JlZW4nLCAnYmxhY2snKSkgKwogIGd1aWRlcyhjb2xvcj1GQUxTRSkgKwogIGFubm90YXRlKCJ0ZXh0IiwgeD0yNCwgeT0wLjI2LCBsYWJlbD1sYWJzMlsxXSwgaGp1c3Q9InJpZ2h0Iiwgc2l6ZT01KSArCiAgYW5ub3RhdGUoInRleHQiLCB4PTI2LCB5PTAuNTgsIGxhYmVsPWxhYnMyWzJdLCBoanVzdD0icmlnaHQiLCBzaXplPTUpICsKICBhbm5vdGF0ZSgidGV4dCIsIHg9MzIsIHk9MC43LCBsYWJlbD1sYWJzMlszXSwgaGp1c3Q9InJpZ2h0Iiwgc2l6ZT01KQpgYGAKCkNvbWJpbmUgdGhlIHBsb3RzCgpgYGB7ciB9CmdyaWQuYXJyYW5nZShwMSwgcDIpCmBgYAoK