dta <- read.csv("Lect13-2013.csv")
ycorn <- dta[,6]

like <- function(b) {
        ll <- log(dgamma(ycorn[1],b[1],1/b[2]))
        for (j in 2:nrow(dta)) {
          ll <- ll+log(dgamma(ycorn[j],b[1],1/b[2])) }
        return(-ll)}

li <- 0.0
while (li <= 0.0) {
  rdraw <- rchisq(2,1,ncp=0)
  astar <- 96.60*rdraw[1]
  bstar <- 1.41*rdraw[2]
  li <- exp(-like(cbind(astar,bstar))) }

drw <- cbind(astar,bstar,li,astar*li,bstar*li)

for (j in 2:5000) {
  li <- 0
  while (li <= 0.0) {
    rdraw <- rchisq(2,1,ncp=0)
    astar <- 96.60*rdraw[1]
    bstar <- 1.41*rdraw[2]
    li <- exp(-like(cbind(astar,bstar))) }
  drw <- rbind(drw,cbind(astar,bstar,li,astar*li,bstar*li)) }
print(drw[1:15,])

pdrw <- sum(drw[,3])
adrw <- sum(drw[,4])/pdrw
bdrw <- sum(drw[,5])/pdrw

print(cbind(pdrw,adrw,bdrw,adrw*bdrw,adrw*bdrw^2))

ml <- optim(cbind(92.0,1.25), like, method = "L-BFGS-B", lower = 0.0001)

print(ml)