Bayesian HW5

Read in the data

lizards <- read.csv("~/Desktop/hornedLizards.csv")
y <- lizards

2a)

NormalNormal <- function(y,S=10000,mu0=0,tau2=100,a=1,b=1){
  n = length(y)
  
  # starting values 
  mu <- 0
  sig2inv <- 1
  
  mu_keep <- rep(NA,S)
  sig2inv_keep <- rep(NA,S)
  
  # MCMC 
  for(s in 1:S){
    # sample mu
    v <- 1/(n*sig2inv + 1/tau2)
    m <- v * (n*sig2inv*mean(y) + mu0/tau2)
    mu <- rnorm(1, m,sqrt(v))
    
    # sample sig2inv
    sig2inv <- rgamma(1, a+n/2 , b+0.5*sum((y-mu)^2))
    
    # save values
    mu_keep[s] <- mu
    sig2inv_keep[s] <- sig2inv
    
  }
  out <- list( mu = mu_keep, sigma=1/sqrt(sig2inv_keep))
  
  ## remove burnin
  mu_keep <- mu_keep[-c(1:round(S/2))]
  sigma_keep <- 1/sqrt(sig2inv_keep[-c(1:round(S/2))])
  
  ## make table
  results <- data.frame(
    mean = c(mean(mu_keep),mean(sigma_keep)),
    sd = c(sd(mu_keep),sd(sigma_keep)),
    lower = c(quantile(mu_keep,0.025),quantile(sigma_keep,0.025)),
    upper = c(quantile(mu_keep,0.975),quantile(sigma_keep,0.975))
  )
  row.names(results) <- c("mu", "sigma")
  out$results <- results
  return(out)
}

2b)

fit <- NormalNormal(y=lizards$squamosalHornLength, S=10000, mu0=0,tau2=100)
fit$results

##            mean        sd     lower     upper
## mu    23.894954 0.2058958 23.487111 24.295188
## sigma  2.768433 0.1445616  2.507819  3.077187

3a)

living <- subset(lizards, lizards$Survival=="living")
killed <- subset(lizards, lizards$Survival=="killed")

3b)

fit2 <- NormalNormal(y=living$squamosalHornLength , S=10000, mu0=0,tau2=100)
fit2$results

##            mean        sd     lower     upper
## mu    24.268440 0.2113857 23.843518 24.676152
## sigma  2.625905 0.1489444  2.352563  2.944811

fit3 <- NormalNormal(y=killed$squamosalHornLength, S=10000, mu0=0,tau2=100)
fit3$results

##            mean        sd    lower     upper
## mu    21.918450 0.4947286 20.96443 22.904673
## sigma  2.690542 0.3486966  2.11127  3.475563

3c)

mean(killed$squamosalHornLength < living$squamosalHornLength)

## Warning in killed$squamosalHornLength < living$squamosalHornLength: longer
## object length is not a multiple of shorter object length

## [1] 0.7012987