# unbuffer the output

# set the working directory to the place You downloaded
# the raw data file in the family income case study

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

print( n <- length( y ) )

m <- 100000

monte.carlo.data.set <- matrix( NA, m, 4 )

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

set.seed( 1 )

alpha <- 0.05

print( theta.true <- mean( y ) )

system.time( {

  for ( i in 1:m ) {

    y.star <- sample( y, n, replace = T )
  
    monte.carlo.data.set[ i, 1 ] <- mean( y.star )

    clt.confidence.interval <- c( mean( y.star ) - 
      qnorm( 1 - alpha / 2 ) * sd( y.star ) / sqrt( n ), 
      mean( y.star ) + qnorm( 1 - alpha / 2 ) * 
      sd( y.star ) / sqrt( n ) )

    monte.carlo.data.set[ i, 2:3 ] <- clt.confidence.interval

    monte.carlo.data.set[ i, 4 ] <- ifelse ( 
      ( clt.confidence.interval[ 1 ] < theta.true ) & 
      ( clt.confidence.interval[ 2 ] > theta.true ), 1, 0 )

    if ( ( i %% 10000 ) == 0 ) print( i )

  }

} 

)

print( clt.95.percent.interval.coverage <- 
  mean( monte.carlo.data.set[ , 4 ] ) )

y.bar.star <- monte.carlo.data.set[ , 1 ]

qqnorm( y.bar.star, 
  main = 'Checking Normality of the Sample Mean With n = 842' )

abline( mean( y.bar.star ), sd( y.bar.star ), lwd = 2, 
  col = 'red' )

hist( y.bar.star, breaks = 100, prob = T, xlab = 'theta',
  main = 'Bootstrap Approximate Posterior Distribution For theta' )

print( theta.hat.bootstrap <- mean( y.bar.star ) )

print( se.theta.hat.bootstrap <- sd( y.bar.star ) )

print( theta.hat.bootstrap.95.percent.interval <-
  quantile( y.bar.star, probs = c( 0.025, 0.975 ) ) )