############
### DATA ###
############

MincerCrossdata <- matrix(c(
     96,   34,   13,    8,   11,
     54,   50,   45,   14,    6,
     11,   86,   83,   38,   10,
      9,   18,   74,  130,   58,
     14,   27,   11,   87,  375), 5, 5, byrow=T)

MincerPercentiles <- matrix(c(
   14.21, 14.64, 15.53, 16.76, 19.60,
   14.38, 15.02, 16.09, 17.26, 20.70,
   15.39, 15.75, 17.34, 18.40, 20.91,
   15.96, 16.64, 17.92, 19.43, 20.74,
   17.58, 18.82, 20.02, 21.52, 23.18), 5, 5, byrow=T)

OttowCountdata <- c(84,  27,   1,  32, 182)

OttowPercentiles <- matrix(c(
   12.05, 13.05, 13.84, 15.02, 19.15,
   13.68, 15.18, 15.94, 16.98, 17.88,
   17.77, 17.77, 17.77, 17.77, 17.77,
   16.13, 16.79, 17.28, 17.89, 24.84,
   17.46, 19.85, 21.56, 23.15, 24.98), 5, 5, byrow=T)

tableDataLucas <- c(0,  0,  0,  4, 25, 27, 39, 44, 50, 50, 50, 50, 50, 50, 50, 50, 50, 50, 50, 50)

getData <- function () 
{
  #returns the table of actual
  #observations
  # rows: teeth
  # columns: knees, 
  # 1: mature
  # 2: not mature
  # 3: "cannot tell"
  matrix(c(4176, 1735, 1364, 
           348, 1087, 237, 
           187, 83, 63), 
         3, 3)
}

getDataFemales <- function () 
  {
    #returns the table of actual
    #observations
    # rows: teeth
    # columns: knees, 
    # 1: mature
    # 2: not mature
    # 3: "cannot tell"
    matrix(c(190, 79, 40, 
             1, 12, 1, 
             6, 5, 3), 
           3, 3)
}

# Global parameter
alpha = 3
#alpha  = 10000

#########################
## Utility functions ####
#########################

rDirichlet <- 
function (alpha) 
{ # Assuming alpha is the parameter vector
  y <- rgamma(length(alpha), alpha, 1)
  y/sum(y)
}

makePriorAge <- function() 
{
  nAge <- 100
  
  # THE FOLLOWING IS THE ALTERNATIVE PRIOR: 
      mmin <- pnorm(15, 22.5, 4)
      mmax <- pnorm(30, 22.5, 4)
      steps <- mmin + (mmax - mmin)*(1:nAge - 0.5)/nAge
      agevals <- qnorm(steps, 22.5, 4)
  
  # THE FOLLOWING IS THE STANDARD PRIOR: 
  #mmax <- pgamma(15, 4, 1)
  #steps <- mmax *(1:nAge - 0.5)/nAge
  #agevals <- qgamma(steps, 4, 1) + 15
  
  cbind(agevals, 1/nAge)
}

getPriorValue <- function(x, prior) {
  sum(dnorm(x[1:2], prior[1:2], prior[3:4], log=TRUE))
}


getProbs <- function (param, ages) 
{
  result <- matrix(0, 3, length(ages))
  result[3,] <- param[3] + (ages-20)*param[4]
  result[1,] <- pnorm((ages - param[1])/param[2])
  result[2,] <- (1-result[1,])*(1-result[3,])
  result[1,] <- result[1,]*(1-result[3,])
  result
}

myDataSummary <- function(data) {
  # Assumes the data is of the form produced by the function reconstructing dat
  # from papers. Outputs means, sd, and percentiles. 
  for (i in 1:5) {
    print(i)
    d <- data[data[,2]==i,1]
    print(quantile(d, probs = c(0, 0.1, 0.25, 0.5, 0.75, 0.9, 1)))
    print(mean(d))
    print(sd(d))
  }
  
}


##################################################
## Functions for computations in paper        ####
##################################################

ParameterEstimationFromPapers <- function(MincerCd = MincerCrossdata,
                                          MincerPerc = MincerPercentiles, 
                                          OttowCountd = OttowCountdata, 
                                          OttowPerc = OttowPercentiles, 
                                          tableDataL = tableDataLucas) {
  resultparams <- matrix(0, 4, 2)
  likelihood <- function (theta, data) 
  {
    # This LOGlikelihoodfunction 
    # is for numerical optimization
    x <- pnorm(data[,1], theta[1], exp(theta[2]))
    x[data[,2]<5] <- 1 - x[data[,2]<5]
    -sum(log(x))
  }
  # 
  # LUCAS: 
  # 
  dataLucas <- matrix(1, 1000, 2)
  for (i in 1:20) dataLucas[((i-1)*50+1):((i-1)*50+50),1] <- 16 + (i-1)*0.5 + 0.25
  for (i in 4:20) dataLucas[((i-1)*50+1):((i-1)*50 + tableDataL[i]),2] <- 5
  res <- nlm(likelihood, c(19, log(1)), dataLucas)
  resultparams[1,1] <- round(res$estimate[1],1)
  resultparams[1,2] <- round(exp(res$estimate[2]),1)
  #
  # Construct Mincer Rawdata
  # 
  N <- 271
  grsize <- apply(MincerCd, 2, sum)
  grsize <- grsize/sum(grsize)*0.999 
  gr <- round(N*grsize, 0)
  start <- cumsum(c(1, gr))
  result <- matrix(0, N, 2)
  perp <- qnorm(c(0.1, 0.25, 0.5, 0.75, 0.9))
  linExtend <- function (x0, x1, x2, y1, y2) 
  {
    #A linear extrapolation
    (x0-x1)/(x2-x1)*(y2-y1) + y1
  }
  for (i in 1:5) {
    result[start[i]:(start[i+1]-1),2] <- i
    xout <- qnorm(((1:gr[i])-0.5)/gr[i])
    ymin <- linExtend(xout[1], perp[1], perp[2], 
                      MincerPercentiles[i,1], MincerPerc[i,2])
    ymax <- linExtend(xout[length(xout)], perp[4], perp[5], 
                      MincerPercentiles[i,4], MincerPerc[i,5])
    ymin <- max(ymin, 14.1)
    ymax <- min(ymax, 24.5)
    pp <- c(xout[1], perp, xout[length(xout)])
    ppy <- c(ymin, MincerPercentiles[i,], ymax)
    result[start[i]:(start[i+1]-1),1] <- 
      approx(x=pp, y=ppy, xout = xout)$y 
  }
  # Do a numerical optimization: 
  res <- nlm(likelihood, c(19, log(2)), result)
  resultparams[2,1] <- round(res$estimate[1],1)
  resultparams[2,2] <- round(exp(res$estimate[2]),1)
  # 
  # OTTOW: 
  # 
  gr <- OttowCountd
  start <- cumsum(c(1, gr))
  N <- sum(gr)
  result <- matrix(0, N, 2)
  perp <- qnorm(c(0.25, 0.5, 0.75))
  for (i in 1:5) {
    result[start[i]:(start[i+1]-1),2] <- i
    xout <- qnorm(((1:gr[i])-0.5)/gr[i])
    pp <- c(xout[1], perp, xout[length(xout)])
    ppy <- OttowPerc[i,]
    result[start[i]:(start[i+1]-1),1] <- 
      approx(x=pp, y=ppy, xout = xout)$y 
  }
  # Do a numerical optimization: 
  res <- nlm(likelihood, c(19, log(2)), result)
  resultparams[3,1] <- round(res$estimate[1],1)
  resultparams[3,2] <- round(exp(res$estimate[2]),1)
  #
  # Adjusted Ottow results: 
  #
  # There are 32 persons assigned category 4. 
  # In this run, we re-assign 16 of them to category 5. 
  ind <- sample(start[4]:(start[5]-1), 16)
  result[ind,2] <- 5
  
  # Do a numerical optimization: 
  res <- nlm(likelihood, c(19, log(2)), result)
  resultparams[4,1] <- round(res$estimate[1],1)
  resultparams[4,2] <- round(exp(res$estimate[2]),1)
  
  resultparams
}

FinalMakeFigure1 <- function(prior1, prior2) {
  age <- makePriorAge()
  x <- age[,1]
  y <- pnorm((x - prior1[1])/prior1[2])
  par(mfrow=c(1,2))
  plot(x, y, xlab="Age", ylab="Probability", type="l")
  for (i in 1:10) {
    theta1 <- rnorm(1, prior1[1], prior1[3])
    theta2 <- -1
    while (theta2<=0) theta2 <- rnorm(1, prior1[2], prior1[4])
    y <- pnorm((x - theta1)/theta2)
    lines(x, y, lty=2)
  }
  y <- pnorm((x - prior2[1])/prior2[2])
  plot(x, y, xlab="Age", ylab="Probability", type="l")
  for (i in 1:10) {
    theta1 <- rnorm(1, prior2[1], prior2[3])
    theta2 <- -1
    while (theta2<=0) theta2 <- rnorm(1, prior2[2], prior2[4])
    y <- pnorm((x - theta1)/theta2)
    lines(x, y, lty=2)
  }
}

FinalMakeTableSuppl <- function(teeth0, knees0) {

  age    <- makePriorAge()
  data   <- getData()
  loglik <- function(theta) {
    # computes the loglikelihood of data given the remaining parameters
    teeth    <- c(teeth0, theta[1:2])
    knees    <- c(knees0, theta[3:4])
    probsT   <- getProbs(teeth, age[,1])
    probsK   <- getProbs(knees, age[,1])
    result   <- 0
    for (i in 1:3)
      for (j in 1:3) 
        result <- result - data[i,j]*
      log(sum(probsT[i,]*probsK[j,]*age[,2]))
    result 
  }
  res <- nlm(loglik,c(0.2, 0, 0.05, 0))$estimate
  print("Maximum Likelihood estimates for missing parameters:")
  print(res)
  teeth <- c(teeth0, res[1:2])
  knees <- c(knees0, res[3:4])
  probsT   <- getProbs(teeth, age[,1])
  probsK   <- getProbs(knees, age[,1])
  myprobs  <- matrix(0, 3, 3)
  for (i in 1:3)
    for (j in 1:3) 
      myprobs[i,j] <- sum(probsT[i,]*probsK[j,]*age[,2])
  print(chisq.test(as.vector(data), p = as.vector(myprobs)))
  predtable <- myprobs*sum(data)
  predtable <- round(predtable)
  print("Data:")
  print(data)
  print("Predicted data:")
  print(predtable)
  print(apply(predtable, 1, sum))
  print(apply(predtable, 2, sum))
  print(sum(predtable))
  pvals <- matrix(0, 3, 3)
  for (i in 1:3)
    for (j in 1:3) 
      pvals[i,j] <- binom.test(data[i,j], sum(data), myprobs[i,j])$p.value
  print(pvals)
}

FinalGibbsSampler <- function(priorT, priorK) 
{
#  alpha <- 3
  
  # Set up the age structure; draw an initial age distribution, and 
  # compute (as age[,2]) the probabilities at the discretized ages 
  # using this age profile: 
  age <- makePriorAge()
  Nage <- dim(age)[1]
  
  #set up starting variables for Gibbs sampling
  teeth    <- c(19, 1.5, 0.2, 0) # was c(19, 1.5, 0.2, 0)
  knees    <- c(18, 1.5, 0.05, 0)  # was c(21, 2, 0.05, 0)
  counts   <- array(0, c(3, 3, Nage))
  data     <- getData()
  
  # get value for prior density at this parameter
  dpriorT  <- getPriorValue(teeth, priorT)
  # get value for prior density at this parameter
  dpriorK  <- getPriorValue(knees, priorK)
  
  probsT   <- getProbs(teeth, age[,1]) # returns 3xnAge matrix
  probsK   <- getProbs(knees, age[,1]) # returns 3xnAge matrix
  
  #Gibbs sampling
  maxstep = 1000000
  result <- list(teeth = matrix(0, maxstep, 4),
                 knees = matrix(0, maxstep, 4),
                 popul = matrix(0, as.integer(maxstep/100), Nage),
                 cumprobs = matrix(0, 4, Nage),
                 tableprobs = matrix(0, maxstep, 9))
  
  result$cumprobs[1:2,] <- probsT[2:3,]
  result$cumprobs[3:4,] <- probsK[2:3,]
  counter <- 1
  
  countT <- 0
  countK <- 0
  for (step in 1:maxstep) {
    # simulate counts given obs model and population model
    for (i in 1:3)
      for (j in 1:3)
        counts[i,j,] <- rmultinom(
          1, data[i,j], 
          prob = probsT[i,]*probsK[j,]*age[,2])
      
      # monitoring convergence: 
      if (step/1000 == round(step/1000)) {
        mytab <- matrix(0, 3, 3)
        for (i in 1:3)
          for (j in 1:3)
            mytab[i,j] <- sum(probsT[i,]*probsK[j,]*age[,2])
          mytab <- mytab/sum(mytab)*sum(data)
          print(round(mytab))
          print(paste("====", step))   
      }
      
      # Simulate a new teeth model: 
      newteeth  <- teeth + rnorm(4, 0, c(0.05, 0.05, 0.01, 0.001))
      newProbsT <- getProbs(newteeth, age[,1])
      if (all(newProbsT > 0)) { # reject directly nonsensical proposals
        newpriorT <- getPriorValue(newteeth, priorT)
        A  <- apply(counts, c(1,3), sum)
        B  <- sum((log(newProbsT) - log(probsT)) * A, na.rm=T) 
        if (runif(1) < exp(B + newpriorT - dpriorT)) {
          countT <- countT + 1
          teeth <- newteeth
          probsT <- newProbsT
          dpriorT <- newpriorT                
        }
      }
      
      # Simulate a new knee model: 
      newknees  <- knees + rnorm(4, 0, c(0.05, 0.05, 0.01, 0.001))
      newProbsK <- getProbs(newknees, age[,1])
      if (all(newProbsK > 0)) { # reject directly nonsensical proposals
        newpriorK <- getPriorValue(newknees, priorK)
        A  <- apply(counts, c(2,3), sum)
        B  <- sum((log(newProbsK) - log(probsK)) * A, na.rm=T) 
        if (runif(1) < exp(B + newpriorK - dpriorK)) {
          countK <- countK + 1
          knees <- newknees
          probsK <- newProbsK
          dpriorK <- newpriorK
        }
      }
      
      # simulate population given counts
      ct <- apply(counts, 3, sum)
      age[,2] <- rDirichlet(ct + alpha/Nage)
      
      
      # Collect results
      result$teeth[step,] <- teeth
      result$knees[step,] <- knees
      result$cumprobs <- result$cumprobs + 
        rbind(probsT[-1,], probsK[-1,])
      result$tableprobs[step,] <- as.vector(
        apply(counts[,,age[,1]<=18], 1:2, sum))
      A <- step/100
      if (A == as.integer(A)) {
        result$popul[counter,] <- cumsum(age[,2])
        counter <- counter + 1
      }
      
  } 
      print(paste("Accep. rate teeth", countT/maxstep))
      print(paste("Accep. rate knees", countK/maxstep))
  result
}

printTable2 <- function(result, ind) {
  print(round(apply(result$teeth[ind,1:2], 2, mean),1))
  print(round(apply(result$knees[ind,1:2], 2, mean),1))
  print(round(quantile(result$teeth[ind,1], c(0.025, 0.975)),1))
  print(round(quantile(result$teeth[ind,2], c(0.025, 0.975)),1))
  print(round(quantile(result$knees[ind,1], c(0.025, 0.975)),1))
  print(round(quantile(result$knees[ind,2], c(0.025, 0.975)),1))
}

FinalMakeTable2 <- function(resultM2L = resM2L, 
                            resultM2M = resM2M, 
                            resultM3  = resM3, 
                            resultM4L = resM4L, 
                            resultM4M = resM4M, 
                            resultSIX = NULL, 
                            NoOfInputs = 5) {
  
  burnin <- 20000
  ind <- (burnin + 1):(dim(resultM2L$teeth)[1])
  printTable2(resultM2L, ind)
  printTable2(resultM2M, ind)
  printTable2(resultM3, ind)
  printTable2(resultM4L, ind)
  if (NoOfInputs>4)
    printTable2(resultM4M, ind)
  if (NoOfInputs>5)
    printTable2(resultSIX, ind)
}

plotFigure2 <- function(result, burnin, age, title) {
  teeth <- apply(result$teeth[-(1:burnin),], 2, mean)
  knees <- apply(result$knees[-(1:burnin),], 2, mean)
  prbT <- getProbs(teeth, age[,1])
  prbK <- getProbs(knees, age[,1])
  plot(age[,1], c(0, rep(1, dim(age)[1]-1)), type="n", 
       xlab="Age", ylab="Probability", main=title)
  lines(age[,1], 1-prbT[2,])
  lines(age[,1], prbT[3,])
  lines(age[,1], 1-prbK[2,], lty=2)
  lines(age[,1], prbK[3,], lty=2)
}


FinalMakeFigure2 <- function(resultM2L = resM2L, 
                             resultM2M = resM2M, 
                             resultM3  = resM3, 
                             resultM4L = resM4L, 
                             resultM4M = resM4M, 
                             resultSIX = NULL, 
                             headings) {
  NoOfInputs <- length(headings)
  if (NoOfInputs == 4)
    par(mfrow=c(2,2))
  else
    par(mfrow=c(3,2))
  age <- makePriorAge()
  # Plot prior 
  if (NoOfInputs == 5) {
  prbT <- getProbs(c(18.62, 0.74, 0.19, 0.023), age[,1])
  prbK <- getProbs(c(18.49, 1.52, 0.03, -0.003), age[,1])
  plot(age[,1], c(0, rep(1, dim(age)[1]-1)), type="n", 
       xlab="Age", ylab="Probability", main = "Fixed parameters")
  lines(age[,1], 1-prbT[2,])
  lines(age[,1], prbT[3,])
  lines(age[,1], 1-prbK[2,], lty=2)
  lines(age[,1], prbK[3,], lty=2)
  }
  
  plotFigure2(resultM2L, 20000, age, headings[1])
  plotFigure2(resultM2M, 20000, age, headings[2])
  plotFigure2(resultM3,  20000, age, headings[3])
  plotFigure2(resultM4L, 20000, age, headings[4])
  if (NoOfInputs > 4)
    plotFigure2(resultM4M, 20000, age, headings[5])
  if (NoOfInputs > 5)
    plotFigure2(resultSIX, 20000, age, headings[6])
}

plotFigure3 <- function(result, age, title) {
  y <- apply(result$popul, 2, median)
  plot(age[,1], y, type="l", xlab="Age", ylab="Probaility", main=title)
  y <- apply(result$popul, 2, quantile, probs = 0.25, names=FALSE)
  lines(age[,1], y, lty=2)
  y <- apply(result$popul, 2, quantile, probs = 0.75, names=FALSE)
  lines(age[,1], y, lty=2)
  y <- apply(result$popul, 2, quantile, probs = 0.025, names=FALSE)
  lines(age[,1], y, lty=3)
  y <- apply(result$popul, 2, quantile, probs = 0.975, names=FALSE)
  lines(age[,1], y, lty=3)
}

FinalMakeFigure3 <- function(resultM2L = resM2L, 
                             resultM2M = resM2M, 
                             resultM3  = resM3, 
                             resultM4L = resM4L, 
                             resultM4M = resM4M, 
                             resultSIX = NULL, 
                             headings) {

  result <- resultM2L
  result$popul <- result$popul[-(1:200),]
  age <- makePriorAge()
  nAge <- dim(age)[1]
  if (length(headings)==4)
    par(mfrow=c(2,2))
  else 
    par(mfrow=c(3,2))
  if (length(headings)==5) {
  # Plot prior: 
  plot(age[,1], (1:nAge)/nAge, type="l", lty=1, xlab="Age", ylab="Probability", main="Prior")
  lines(age[,1], qbeta(0.25, 1:nAge*alpha/nAge, (nAge-1):0*alpha/nAge), lty=2)
  lines(age[,1], qbeta(0.75, 1:nAge*alpha/nAge, (nAge-1):0*alpha/nAge), lty=2)
  lines(age[,1], qbeta(0.025, 1:nAge*alpha/nAge, (nAge-1):0*alpha/nAge), lty=3)
  lines(age[,1], qbeta(0.975, 1:nAge*alpha/nAge, (nAge-1):0*alpha/nAge), lty=3)
  }
  # Plot posterior: 
  # 
  plotFigure3(resultM2L, age, headings[1])
  plotFigure3(resultM2M, age, headings[2])
  plotFigure3(resultM3,  age, headings[3])
  plotFigure3(resultM4L, age, headings[4])
  if (length(headings)>4)
    plotFigure3(resultM4M, age, headings[5])
  if (length(headings)> 5)
    plotFigure3(resultSIX, age, headings[6])
}





FinalMakeTable3 <- function(resultM2L = resM2L, 
                            resultM2M = resM2M, 
                            resultM3  = resM3,
                            resultM4L = resM4L,
                            resultM4M = resM4M) {
  data <- getData()
  expectedRate <- matrix(0, 8, 5)
  lowerRate    <- expectedRate
  upperRate    <- expectedRate
  #
  info <- resultM2L$tableprobs
  A <- matrix(apply(info, 2, mean), 3, 3)
  B <- round(100*A/data)
  expectedRate[1:3,1] <- B[1:3,1]
  expectedRate[4:5,1] <- B[1,2:3]
  expectedRate[6:7,1] <- 100 - B[2:3,2]
  expectedRate[8  ,1] <- 100 - B[2,3]
  A <- matrix(apply(info, 2, quantile, probs=0.025, names=FALSE), 3, 3)
  B <- round(100*A/data)
  lowerRate[1:3,1]    <- B[1:3,1]
  lowerRate[4:5,1]    <- B[1,2:3]
  upperRate[6:7,1]    <- 100 - B[2:3,2]
  upperRate[8  ,1]    <- 100 - B[2,3]
  A <- matrix(apply(info, 2, quantile, probs=0.975, names=FALSE), 3, 3)
  B <- round(100*A/data)
  upperRate[1:3,1]    <- B[1:3,1]
  upperRate[4:5,1]    <- B[1,2:3]
  lowerRate[6:7,1]    <- 100 - B[2:3,2]
  lowerRate[8  ,1]    <- 100 - B[2,3]
  #
  info <- resultM2M$tableprobs
  A <- matrix(apply(info, 2, mean), 3, 3)
  B <- round(100*A/data)
  expectedRate[1:3,2] <- B[1:3,1]
  expectedRate[4:5,2] <- B[1,2:3]
  expectedRate[6:7,2] <- 100 - B[2:3,2]
  expectedRate[8  ,2] <- 100 - B[2,3]
  A <- matrix(apply(info, 2, quantile, probs=0.025, names=FALSE), 3, 3)
  B <- round(100*A/data)
  lowerRate[1:3,2]    <- B[1:3,1]
  lowerRate[4:5,2]    <- B[1,2:3]
  upperRate[6:7,2]    <- 100 - B[2:3,2]
  upperRate[8  ,2]    <- 100 - B[2,3]
  A <- matrix(apply(info, 2, quantile, probs=0.975, names=FALSE), 3, 3)
  B <- round(100*A/data)
  upperRate[1:3,2]    <- B[1:3,1]
  upperRate[4:5,2]    <- B[1,2:3]
  lowerRate[6:7,2]    <- 100 - B[2:3,2]
  lowerRate[8  ,2]    <- 100 - B[2,3]
  #
  info <- resultM3$tableprobs
  A <- matrix(apply(info, 2, mean), 3, 3)
  B <- round(100*A/data)
  expectedRate[1:3,3] <- B[1:3,1]
  expectedRate[4:5,3] <- B[1,2:3]
  expectedRate[6:7,3] <- 100 - B[2:3,2]
  expectedRate[8  ,3] <- 100 - B[2,3]
  A <- matrix(apply(info, 2, quantile, probs=0.025, names=FALSE), 3, 3)
  B <- round(100*A/data)
  lowerRate[1:3,3]    <- B[1:3,1]
  lowerRate[4:5,3]    <- B[1,2:3]
  upperRate[6:7,3]    <- 100 - B[2:3,2]
  upperRate[8  ,3]    <- 100 - B[2,3]
  A <- matrix(apply(info, 2, quantile, probs=0.975, names=FALSE), 3, 3)
  B <- round(100*A/data)
  upperRate[1:3,3]    <- B[1:3,1]
  upperRate[4:5,3]    <- B[1,2:3]
  lowerRate[6:7,3]    <- 100 - B[2:3,2]
  lowerRate[8  ,3]    <- 100 - B[2,3]
  #
  info <- resultM4L$tableprobs
  A <- matrix(apply(info, 2, mean), 3, 3)
  B <- round(100*A/data)
  expectedRate[1:3,4] <- B[1:3,1]
  expectedRate[4:5,4] <- B[1,2:3]
  expectedRate[6:7,4] <- 100 - B[2:3,2]
  expectedRate[8  ,4] <- 100 - B[2,3]
  A <- matrix(apply(info, 2, quantile, probs=0.025, names=FALSE), 3, 3)
  B <- round(100*A/data)
  lowerRate[1:3,4]    <- B[1:3,1]
  lowerRate[4:5,4]    <- B[1,2:3]
  upperRate[6:7,4]    <- 100 - B[2:3,2]
  upperRate[8  ,4]    <- 100 - B[2,3]
  A <- matrix(apply(info, 2, quantile, probs=0.975, names=FALSE), 3, 3)
  B <- round(100*A/data)
  upperRate[1:3,4]    <- B[1:3,1]
  upperRate[4:5,4]    <- B[1,2:3]
  lowerRate[6:7,4]    <- 100 - B[2:3,2]
  lowerRate[8  ,4]    <- 100 - B[2,3]
  #
  info <- resultM4M$tableprobs
  A <- matrix(apply(info, 2, mean), 3, 3)
  B <- round(100*A/data)
  expectedRate[1:3,5] <- B[1:3,1]
  expectedRate[4:5,5] <- B[1,2:3]
  expectedRate[6:7,5] <- 100 - B[2:3,2]
  expectedRate[8  ,5] <- 100 - B[2,3]
  A <- matrix(apply(info, 2, quantile, probs=0.025, names=FALSE), 3, 3)
  B <- round(100*A/data)
  lowerRate[1:3,5]    <- B[1:3,1]
  lowerRate[4:5,5]    <- B[1,2:3]
  upperRate[6:7,5]    <- 100 - B[2:3,2]
  upperRate[8  ,5]    <- 100 - B[2,3]
  A <- matrix(apply(info, 2, quantile, probs=0.975, names=FALSE), 3, 3)
  B <- round(100*A/data)
  upperRate[1:3,5]    <- B[1:3,1]
  upperRate[4:5,5]    <- B[1,2:3]
  lowerRate[6:7,5]    <- 100 - B[2:3,2]
  lowerRate[8  ,5]    <- 100 - B[2,3]
  
  
  print(expectedRate)
  print(lowerRate)
  print(upperRate)
}

####################################################################
### RUNNING COMPUTATIONS                                         ###
####################################################################


# Fix these values, for now: 
#resultparams <- ParameterEstimationFromPapers()

# The calculations above may vary because of a random sampling inthe adjusted
# Ottow data. Thus we present below the numbers used in the paper: 
resultparams  <- matrix(c(
  18.6, 0.7, 
  20.0, 3.2, 
  18.5, 1.5, 
  17.8, 1.7), 4, 2, byrow=T)

#print("Parameters estimated from Lucas:")
#print(resultparams[1,])
#print("Parameters estimated from Mincer:")
#print(resultparams[2,])
#print("Parameters estimated from Ottow:")
#print(resultparams[3,])
#print("Parameters estimated from adjusted Ottow data:")
#print(resultparams[4,])

#FinalMakeTableSuppl(resultparams[1,], resultparams[3,])
#FinalMakeTableSuppl(resultparams[2,], resultparams[3,])

# Set the priors based on estimated parameters
TL  <- c(resultparams[1,], 0.001, 0.001)
TM  <- c(resultparams[2,], 0.001, 0.001)
TLM <- c(apply(resultparams[1:2,], 2, mean), 0.001, 0.001)
T1  <- c(resultparams[1,], 0.2, 0.2)
T2  <- c(resultparams[2,], 0.2, 0.2)
T3  <- c(apply(resultparams[1:2,], 2, mean), 0.8, 0.8)
KO  <- c(resultparams[3,], 0.001, 0.001)
KOA <- c(resultparams[4,], 0.001, 0.001)
K1  <- c(resultparams[3,], 0.2, 0.2)
K2  <- c(resultparams[3,], 0.8, 0.8)
K3  <- c(resultparams[4,], 0.2, 0.2)
Ksoc1 <- c(18.6, 0.7, 0.2, 0.2)
Ksoc2 <- c(19.6, 0.7, 0.2, 0.2)
Ksoc3 <- c(20.0, 3.2, 0.2, 0.2)
Ksoc4 <- c(21.0, 3.2, 0.2, 0.2)
prior <- list(TL=TL, TM=TM, TLM=TLM, T1=T1, T2=T2, T3=T3, 
              KO=KO, KOA=KOA, K1=K1, K2=K2, K3=K3, 
              Ksoc1=Ksoc1, Ksoc2=Ksoc2, Ksoc3=Ksoc3, Ksoc4=Ksoc4)
#print(prior)

#FinalMakeFigure1(prior$K1, prior$K2)
#x11()

#print(system.time(resML   <- FinalGibbsSampler(prior$TL, prior$KO)))
#print(system.time(resMM   <- FinalGibbsSampler(prior$TM, prior$KO)))

#print(system.time(resM2L <- FinalGibbsSampler(prior$T1, prior$K1)))
#print(system.time(resM2M <- FinalGibbsSampler(prior$T2, prior$K1)))
#print(system.time(resM3  <- FinalGibbsSampler(prior$T3, prior$K2)))
#print(system.time(resM4L <- FinalGibbsSampler(prior$T1, prior$K3)))
#print(system.time(resM4M <- FinalGibbsSampler(prior$T2, prior$K3)))

#FinalMakeTable2()
#FinalMakeFigure2(headings = c("Model 1", "Model 2", "Model 3",
#                              "Model 4", "Model 5"))
#x11()

#FinalMakeFigure3(headings = c("Model 1", "Model 2", "Model 3", 
#                              "Model 4", "Model 5"))
#FinalMakeTable3()

################################################
### Code for supplementary material 
################################################

#alpha <- 100000 #fix the population profile
#print(system.time(supp1  <- FinalGibbsSampler(prior$T1, prior$K1)))
#print(system.time(supp2  <- FinalGibbsSampler(prior$T1, prior$K2)))
#print(system.time(supp3  <- FinalGibbsSampler(prior$T1, prior$K3)))
#print(system.time(supp4  <- FinalGibbsSampler(prior$T2, prior$K1)))
#print(system.time(supp5  <- FinalGibbsSampler(prior$T2, prior$K2)))
#print(system.time(supp6  <- FinalGibbsSampler(prior$T2, prior$K3)))
#print(system.time(supp7  <- FinalGibbsSampler(prior$T3, prior$K1)))
#print(system.time(supp8  <- FinalGibbsSampler(prior$T3, prior$K2)))
#print(system.time(supp9  <- FinalGibbsSampler(prior$T3, prior$K3)))

#FinalMakeTable2(supp1, supp2, supp3, supp4, supp5)
#FinalMakeTable2(supp6, supp7, supp8, supp9, supp9)

#alpha <- 3 # reset to hierarchical prior
#print(system.time(supp10  <- FinalGibbsSampler(prior$TL, prior$KO)))
#print(system.time(supp11  <- FinalGibbsSampler(prior$TM, prior$KO)))
#print(system.time(supp12  <- FinalGibbsSampler(prior$TLM, prior$KO)))
#print(system.time(supp13  <- FinalGibbsSampler(prior$TL, prior$KOA)))
#print(system.time(supp14  <- FinalGibbsSampler(prior$TM, prior$KOA)))
#print(system.time(supp15  <- FinalGibbsSampler(prior$TLM, prior$KOA)))
#FinalMakeFigure3(supp10, supp11, supp12, supp13, supp14, supp15, 
#                 headings = c("Model 10", "Model 11", "Model 12", 
#                              "Model 13", "Model 14", "Model 15"))

#print(system.time(supp16  <- FinalGibbsSampler(prior$T1, prior$K2)))
#print(system.time(supp17  <- FinalGibbsSampler(prior$T2, prior$K2)))
#print(system.time(supp18  <- FinalGibbsSampler(prior$T3, prior$K1)))
#print(system.time(supp19  <- FinalGibbsSampler(prior$T3, prior$K3)))
#FinalMakeFigure3(supp16, supp17, supp18, supp19, headings=c("Model 16", 
#                                                            "Model 17", 
#                                                            "Model 18", 
#                                                            "Model 19"))
#x11()
#FinalMakeFigure2(supp16, supp17, supp18, supp19, headings=c("Model 16", 
#                                                            "Model 17", 
#                                                            "Model 18", 
#                                                            "Model 19"))
#FinalMakeTable2(supp16, supp17, supp18, supp19, NoOfInputs = 4)
#FinalMakeTable3(supp16, supp17, supp18, supp19)

#print(system.time(supp20  <- FinalGibbsSampler(prior$T1, prior$Ksoc1)))
#print(system.time(supp21  <- FinalGibbsSampler(prior$T1, prior$Ksoc2)))
#print(system.time(supp22  <- FinalGibbsSampler(prior$T2, prior$Ksoc3)))
#print(system.time(supp23  <- FinalGibbsSampler(prior$T2, prior$Ksoc4)))
#FinalMakeTable2(supp20, supp21, supp22, supp23, NoOfInputs = 4)
#FinalMakeFigure3(supp20, supp21, supp22, supp23, headings = c("Model 20", 
#                                                              "Model 21", 
#                                                              "Model 22", 
#                                                              "Model 23"))
#x11()
#FinalMakeFigure2(supp20, supp21, supp22, supp23, headings = c("Model 20", 
#                                                              "Model 21", 
#                                                              "Model 22", 
#                                                              "Model 23"))
#FinalMakeTable3(supp20, supp21, supp22, supp23)


# RESET THE STARTING POINT FOR THE PRIOR FOR THE AGE DISTRIBUTION
# BY EDITING makePriorAge

# Recompute paper results: 
# print(system.time(resM2L <- FinalGibbsSampler(prior$T1, prior$K1)))
# print(system.time(resM2M <- FinalGibbsSampler(prior$T2, prior$K1)))
# print(system.time(resM3  <- FinalGibbsSampler(prior$T3, prior$K2)))
# print(system.time(resM4L <- FinalGibbsSampler(prior$T1, prior$K3)))
# print(system.time(resM4M <- FinalGibbsSampler(prior$T2, prior$K3)))

# FinalMakeTable2()
# FinalMakeFigure2(headings = c("Model 1", "Model 2", "Model 3",
#                               "Model 4", "Model 5"))
# x11()
# FinalMakeFigure3(headings = c("Model 1", "Model 2", "Model 3", 
#                               "Model 4", "Model 5"))
# FinalMakeTable3()