demo3_5.Rmd

---
title: "Bayesian data analysis demo 3.5"
author: "Aki Vehtari, Markus Paasiniemi"
date: "`r format(Sys.Date())`"
output:
  html_document:
    theme: readable
    code_download: true
---
## Estimating the speed of light using normal model BDA3 p. 66

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

```{r setup, message=FALSE, error=FALSE, warning=FALSE}
library(ggplot2)
theme_set(theme_minimal())
library(gridExtra)
library(tidyr)
library(rprojroot)
root<-has_file(".BDA_R_demos_root")$make_fix_file()
```

Data

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

Sufficient statistics

```{r }
n <- length(y)
s2 <- var(y)
my <- mean(y)
```

Positive values only

```{r }
y_pos <- y[y > 0]
```

Sufficient statistics

```{r }
n_pos <- length(y_pos)
s2_pos <- var(y_pos)
my_pos <- mean(y_pos)
```

Compute the density of mu in these points

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

Compute the exact marginal density for mu

```{r }
# 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

```{r }
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

```{r }
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

```{r }
# 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

```{r }
grid.arrange(p1, p2)
```