# Code to generate the contour plots on page 7 of lecture 11 overheads

rats <- read.table("rats_data.txt", header=T)
y <- rats$y
n <- rats$n


u <- (-220:-140)/100
v <- (100:500)/100
lenu <- length(u)
lenv <- length(v)


logpostden2020 <- matrix(0,ncol=lenv, nrow=lenu)
for(i in 1:lenu) {
  for(j in 1:lenv) {
    beta <- exp(v[j])/(1+exp(u[i]))
    alpha <- exp(u[i]) * beta
    logpostden2020[i,j] <- log(alpha) + log(beta) + sum(lbeta(alpha+y,beta+n-y)-lbeta(alpha,beta))
  }
}

postden2020 <- exp(logpostden2020-max(logpostden2020))
for(i in 1:lenu) {
  for(j in 1:lenv) {
    beta <- exp(v[j])/(1+exp(u[i]))
    alpha <- exp(u[i]) * beta
    postden2020[i,j] <- postden2020[i,j] * (alpha < 20) * (beta < 20)
  }
}

alpha <- (1:100)/10
beta <- (1:80)/2
lena <- length(alpha)
lenb <- length(beta)

logpostden2020a <- matrix(0, ncol=lenb, nrow=lena)
for(i in 1:lena) {
  for(j in 1:lenb) {
    logpostden2020a[i,j] <- sum(lbeta(alpha[i]+y,beta[j]+n-y)-lbeta(alpha[i],beta[j]))
  }
}

postden2020a <- exp(logpostden2020a-max(logpostden2020a))
for(i in 1:lena) {
  for(j in 1:lenb) {
    postden2020a[i,j] <- postden2020a[i,j] * (alpha[i] < 20) * (beta[j] < 20)
  }
}


postscript("../ratpost2020.eps",horiz=F,width=8,height=4)
par(mfrow=c(1,2), pty="m", mar=c(4,4,2,1)+0.1)

contour(u,v,postden2020,levels=seq(0.05,0.95,0.1),drawlabels=F,
  main=expression(paste(alpha," ~ U(0,20), ", beta," ~ U(0,20)")),
  xlab=expression(log(alpha/beta)), ylab=expression(log(alpha+beta)))
contour(alpha,beta,postden2020a,levels=seq(0.05,0.95,0.1),drawlabels=F,
  main=expression(paste(alpha," ~ U(0,20), ", beta," ~ U(0,20)")),
  xlab=expression(alpha), ylab=expression(beta))
dev.off()