Statistical Rethinking Chapter 3 Hard Question 1

Rethinking_ch3_H1.R

# Chapter three

## Hard questions


#' These data indicate the gender (male=1, female=0) of officially reported
#' first and second born children in 100 two-child families.

birth1 =c(1,0,0,0,1,1,0,1,0,1,0,0,1,1,0,1,1,0,0,0,1,0,0,0,1,0,
           0,0,0,1,1,1,0,1,0,1,1,1,0,1,0,1,1,0,1,0,0,1,1,0,1,0,0,0,0,0,0,0,
           1,1,0,1,0,0,1,0,0,0,1,0,0,1,1,1,1,0,1,0,1,1,1,1,1,0,0,1,0,1,1,0,
           1,0,1,1,1,0,1,1,1,1)
birth2 =c(0,1,0,1,0,1,1,1,0,0,1,1,1,1,1,0,0,1,1,1,0,0,1,1,1,0,
           1,1,1,0,1,1,1,0,1,0,0,1,1,1,1,0,0,1,0,1,1,1,1,1,1,1,1,1,1,1,1,1,
           1,1,1,0,1,1,0,1,1,0,1,1,1,0,0,0,0,0,0,1,0,0,0,1,1,0,0,1,0,0,1,1,
           0,0,0,1,1,1,0,0,0,0)

### Question 1

#' Using grid approximation, compute the posterior distribution for the probability
#' of a birth beinga boy. Assume a uniform prior probability. Which parameter value
#' maximises the posterior probability?

# Define parameter values
p = seq(from = 0, to = 1,
         length.out = 8128 #'arbitrary'
        )

# Define prior (uniform in this case)
prior = rep(0.6, #'arbitrary'
             length(p))

# Compute likelihoods

## sum number of boys (i.e. what we observe)
boys = sum(birth1) + sum(birth2)

## get likelihood (binom because binary e.g.)
likelihood = dbinom(boys,
                     size = 200,
                     prob = p)

# Get posterior
posterior = likelihood * prior
posterior = posterior / sum(posterior)

# Find the value that maximises posterior probability
p[which.max(posterior)]