# unbuffer the output

# set the working directory to the place You downloaded
# the raw data file in the family income case study,
# and ensure that that place also contains the rjags
# model files for the lognormal and gamma models

library( rjags ) 

y <- scan( 'U.S.-family-income-2009.txt' )

print( n <- length( y ) )

alpha <- 0.05

theta.true <- mean( y )

income.lognormal.initial.values <- 
  list( mu = 0.0, log.tau = 0.0, y.new = mean( y ) )

income.gamma.initial.values <- 
  list( log.alpha = 0.0, log.beta = 0.0, y.new = mean( y ) )

lognormal.number.of.hits.so.far <- 0

cumulative.lognormal.hit.rate <- 0

gamma.number.of.hits.so.far <- 0

cumulative.gamma.hit.rate <- 0

# the next four commands are up to you

# use Your own personal random number seed instead of 1
# in the set.seed command below

set.seed( 9470233 )

# use a bigger value of m (e.g., 200 or 500 or 1000) 
# in the next command, if you have time

m <- 100

# use a bigger value of n.burnin (e.g., 1000) 
# in the next command, if you have time

n.burnin <- 100

# use a bigger value of n.monitor (e.g., 20000 or 50000 or 100000) 
# in the next command, if you have time

n.monitor <- 10000

system.time( {

  for ( i in 1:m ) {

    y.star <- sample( y, n, replace = T )

    income.bootstrapped.data <- list( y = y.star, 
      n = length( y.star ) )
  
    income.lognormal.calibration.run <- jags.model( file = 
      'ams-206-income-lognormal-analysis-rjags-model.txt', 
      data = income.bootstrapped.data, 
      inits = income.lognormal.initial.values )

    update( income.lognormal.calibration.run, n.iter = n.burnin )

    income.lognormal.calibration.run.results <- 
      jags.samples( income.lognormal.calibration.run, 
      variable.names = c( 'theta' ), n.iter = n.monitor )

    lognormal.theta.star <- 
      income.lognormal.calibration.run.results$theta

    lognormal.interval.endpoints <- quantile( lognormal.theta.star, 
      probs = c( alpha / 2, 1 - alpha / 2 ) ) 

    lognormal.hit <- ( lognormal.interval.endpoints[ 1 ] < 
      theta.true ) & ( theta.true < lognormal.interval.endpoints[ 2 ] )

    lognormal.number.of.hits.so.far <- 
      lognormal.number.of.hits.so.far + lognormal.hit

    cumulative.lognormal.hit.rate <- 
      signif( as.numeric( lognormal.number.of.hits.so.far / i ), 4 )

    income.gamma.calibration.run <- jags.model( file = 
      'ams-206-income-gamma-analysis-rjags-model.txt', 
      data = income.bootstrapped.data, 
      inits = income.gamma.initial.values )

    update( income.gamma.calibration.run, n.iter = n.burnin )

    income.gamma.calibration.run.results <- 
      jags.samples( income.gamma.calibration.run, 
      variable.names = c( 'theta' ), n.iter = n.monitor )

    gamma.theta.star <- 
      income.gamma.calibration.run.results$theta

    gamma.interval.endpoints <- quantile( gamma.theta.star, 
      probs = c( alpha / 2, 1 - alpha / 2 ) ) 

    gamma.hit <- ( gamma.interval.endpoints[ 1 ] < 
      theta.true ) & ( theta.true < gamma.interval.endpoints[ 2 ] )

    gamma.number.of.hits.so.far <- 
      gamma.number.of.hits.so.far + gamma.hit

    cumulative.gamma.hit.rate <- 
      signif( as.numeric( gamma.number.of.hits.so.far / i ), 4 )

    cat( '\n' )

    print( c( 'iteration number', i, 'cumulative lognormal hit rate',
      cumulative.lognormal.hit.rate, 'cumulative gamma hit rate',
      cumulative.gamma.hit.rate ) )

    cat( '\n' )

  }

}

)

# on my desktop at home, with a single-threaded cpu at 3.4 GHz, 
# m = 10 replications took about 503 seconds with n.burnin 
# and n.monitor set to the unusually low values above, so if
# your processor has about the same GHz as mine you can figure
# on about 50 seconds per replication with those same values
# of n.burnin and n.monitor; this may help you decide
# on your chosen values of m, n.burnin and n.monitor

# a small future software engineering project would be to use
# foreach and %dopar% (or snow, or ...) to parallelize the
# code above, to make use of all available threads on your
# machine simultaneously (that should not be hard to do);
# if anybody wants to do this and send me the code
# for extra-extra credit, that would be terrific