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

# pivot the data points into key-value pairs
df2 <- df1 %>% pivot_longer(cols = !t1, names_to="grp", values_to="p")
# 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="none") +
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)

LS0tCnRpdGxlOiAiQmF5ZXNpYW4gZGF0YSBhbmFseXNpcyBkZW1vIDMuNSIKYXV0aG9yOiAiQWtpIFZlaHRhcmksIE1hcmt1cyBQYWFzaW5pZW1pIgpkYXRlOiAiYHIgZm9ybWF0KFN5cy5EYXRlKCkpYCIKb3V0cHV0OgogIGh0bWxfZG9jdW1lbnQ6CiAgICB0aGVtZTogcmVhZGFibGUKICAgIGNvZGVfZG93bmxvYWQ6IHRydWUKLS0tCiMjIEVzdGltYXRpbmcgdGhlIHNwZWVkIG9mIGxpZ2h0IHVzaW5nIG5vcm1hbCBtb2RlbCBCREEzIHAuIDY2CgpnZ3Bsb3QyLCBhbmQgZ3JpZEV4dHJhIGFyZSB1c2VkIGZvciBwbG90dGluZywgdGlkeXIgZm9yCm1hbmlwdWxhdGluZyBkYXRhIGZyYW1lcwoKYGBge3Igc2V0dXAsIG1lc3NhZ2U9RkFMU0UsIGVycm9yPUZBTFNFLCB3YXJuaW5nPUZBTFNFfQpsaWJyYXJ5KGdncGxvdDIpCnRoZW1lX3NldCh0aGVtZV9taW5pbWFsKCkpCmxpYnJhcnkoZ3JpZEV4dHJhKQpsaWJyYXJ5KHRpZHlyKQpsaWJyYXJ5KHJwcm9qcm9vdCkKcm9vdDwtaGFzX2ZpbGUoIi5CREFfUl9kZW1vc19yb290IikkbWFrZV9maXhfZmlsZSgpCmBgYAoKRGF0YQoKYGBge3IgfQp5IDwtIHJlYWQudGFibGUocm9vdCgiZGVtb3NfY2gzIiwibGlnaHQudHh0IikpJFYxCmBgYAoKU3VmZmljaWVudCBzdGF0aXN0aWNzCgpgYGB7ciB9Cm4gPC0gbGVuZ3RoKHkpCnMyIDwtIHZhcih5KQpteSA8LSBtZWFuKHkpCmBgYAoKUG9zaXRpdmUgdmFsdWVzIG9ubHkKCmBgYHtyIH0KeV9wb3MgPC0geVt5ID4gMF0KYGBgCgpTdWZmaWNpZW50IHN0YXRpc3RpY3MKCmBgYHtyIH0Kbl9wb3MgPC0gbGVuZ3RoKHlfcG9zKQpzMl9wb3MgPC0gdmFyKHlfcG9zKQpteV9wb3MgPC0gbWVhbih5X3BvcykKYGBgCgpDb21wdXRlIHRoZSBkZW5zaXR5IG9mIG11IGluIHRoZXNlIHBvaW50cwoKYGBge3IgfQp0bDEgPC0gYygxOCwgMzQpCmRmMSA8LSBkYXRhLmZyYW1lKHQxID0gc2VxKHRsMVsxXSwgdGwxWzJdLCBsZW5ndGgub3V0ID0gMTAwMCkpCmBgYAoKQ29tcHV0ZSB0aGUgZXhhY3QgbWFyZ2luYWwgZGVuc2l0eSBmb3IgbXUKCmBgYHtyIH0KIyBtdWx0aXBsaWNhdGlvbiBieSAxLi9zcXJ0KHMyL24pIGlzIGR1ZSB0byB0aGUgdHJhbnNmb3JtYXRpb24gb2YgdmFyaWFibGUKIyB6PSh4LW1lYW4oeSkpL3NxcnQoczIvbiksIHNlZSBCREEzIHAuIDIxCmRmMSRwbV9tdSA8LSBkdCgoZGYxJHQxIC0gbXkpIC8gc3FydChzMi9uKSwgbi0xKSAvIHNxcnQoczIvbikKYGBgCgpDb21wdXRlIHRoZSBleGFjdCBtYXJnaW5hbCBkZW5zaXR5IGZvciBtdSBmb3IgdGhlIHBvc2l0aXZlIGRhdGEKCmBgYHtyIH0KZGYxJHBtX211X3BvcyA9IGR0KChkZjEkdDEgLSBteV9wb3MpIC8gc3FydChzMl9wb3Mvbl9wb3MpLCBuX3Bvcy0xKSAvIHNxcnQoczJfcG9zL25fcG9zKQpgYGAKCkNyZWF0ZSBhIGhpc3RvZ3JhbSBvZiB0aGUgbWVhc3VyZW1lbnRzCgpgYGB7ciB9CnAxIDwtIGdncGxvdCgpICsKICBnZW9tX2hpc3RvZ3JhbShhZXMoeSksIGJpbndpZHRoID0gMiwgZmlsbCA9ICdzdGVlbGJsdWUnLCBjb2xvciA9ICdibGFjaycpICsKICBjb29yZF9jYXJ0ZXNpYW4oeGxpbSA9IGMoLTQyLCA0MikpICsKICBzY2FsZV95X2NvbnRpbnVvdXMoYnJlYWtzID0gTlVMTCkgKwogIGxhYnModGl0bGUgPSAnTmV3Y29tYlwncyBtZWFzdXJlbWVudHMnLCB4ID0gJ3knKQpgYGAKCkNyZWF0ZSBhIHBsb3Qgb2YgdGhlIG5vcm1hbCBtb2RlbAoKYGBge3IgfQojIHBpdm90IHRoZSBkYXRhIHBvaW50cyBpbnRvIGtleS12YWx1ZSBwYWlycwpkZjIgPC0gZGYxICU+JSBwaXZvdF9sb25nZXIoY29scyA9ICF0MSwgbmFtZXNfdG89ImdycCIsIHZhbHVlc190bz0icCIpCiMgbGVnZW5kIGxhYmVscwpsYWJzMiA8LSBjKCdQb3N0ZXJpb3Igb2YgbXUnLCAnUG9zdGVyaW9yIG9mIG11IGdpdmVuIHkgPiAwJywgJ01vZGVybiBlc3RpbWF0ZScpCnAyIDwtIGdncGxvdChkYXRhID0gZGYyKSArCiAgZ2VvbV9saW5lKGFlcyh0MSwgcCwgY29sb3IgPSBncnApKSArCiAgZ2VvbV92bGluZShhZXMoeGludGVyY2VwdCA9IDMzLCBjb2xvciA9ICdxJyksCiAgICAgICAgICAgICBsaW5ldHlwZSA9ICdkYXNoZWQnLCBzaG93LmxlZ2VuZCA9IEYpICsKICBjb29yZF9jYXJ0ZXNpYW4oeGxpbSA9IGMoLTQyLCA0MikpICsKICBzY2FsZV95X2NvbnRpbnVvdXMoYnJlYWtzID0gTlVMTCkgKwogIGxhYnModGl0bGUgPSAnTm9ybWFsIG1vZGVsJywgeCA9ICdtdScsIHkgPSAnJykgKwogIHNjYWxlX2NvbG9yX21hbnVhbCh2YWx1ZXMgPSBjKCdibHVlJywgJ2RhcmtncmVlbicsICdibGFjaycpKSArCiAgZ3VpZGVzKGNvbG9yPSJub25lIikgKwogIGFubm90YXRlKCJ0ZXh0IiwgeD0yNCwgeT0wLjI2LCBsYWJlbD1sYWJzMlsxXSwgaGp1c3Q9InJpZ2h0Iiwgc2l6ZT01KSArCiAgYW5ub3RhdGUoInRleHQiLCB4PTI2LCB5PTAuNTgsIGxhYmVsPWxhYnMyWzJdLCBoanVzdD0icmlnaHQiLCBzaXplPTUpICsKICBhbm5vdGF0ZSgidGV4dCIsIHg9MzIsIHk9MC43LCBsYWJlbD1sYWJzMlszXSwgaGp1c3Q9InJpZ2h0Iiwgc2l6ZT01KQpgYGAKCkNvbWJpbmUgdGhlIHBsb3RzCgpgYGB7ciB9CmdyaWQuYXJyYW5nZShwMSwgcDIpCmBgYAoK