# This does Metropolis sampling in the one-sample Gaussian problem 
# (known mean, unknown variance) with the NB10 data

metropolis.example <- function( y, mu, sigma.star, sigma2.0, M, B ) {

  n <- length( y )

  sigma2.hat <- sum( ( y - mu )^2 ) / n

  lambda.hat <- log( sigma2.hat )

  mcmc.data.set <- matrix( NA, M + B + 1, 2 )

  mcmc.data.set[ 1, ] <- c( log( sigma2.0 ), sigma2.0 )

  acceptance.count <- 0

  for ( i in 2:( M + B + 1 ) ) {

    lambda.current <- mcmc.data.set[ i - 1, 1 ]

    lambda.star <- rnorm( 1, lambda.current, sigma.star )

    u <- runif( 1 )

    if ( u <= acceptance.probability( lambda.star, lambda.current,
      lambda.hat, n ) ) {

        mcmc.data.set[ i, 1 ] <- lambda.star

        mcmc.data.set[ i, 2 ] <- exp( lambda.star )

        acceptance.count <- acceptance.count + 1

    }

    else {

      mcmc.data.set[ i, 1 ] <- lambda.current

      mcmc.data.set[ i, 2 ] <- exp( lambda.current )

    }

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

  }

  print( acceptance.rate <- acceptance.count / ( M + B ) )

  return( mcmc.data.set )

}

acceptance.probability <- function( lambda.star, lambda.current, 
  lambda.hat, n ) {

  return( exp( log.posterior( lambda.star, lambda.hat, n ) -
    log.posterior( lambda.current, lambda.hat, n ) ) )

}

log.posterior <- function( lambda, lambda.hat, n ) {

  return( log.prior( lambda ) + log.likelihood( lambda, lambda.hat, 
    n ) )

}

log.likelihood <- function( lambda, lambda.hat, n ) {

  return( - lambda * n / 2 - n * exp( lambda.hat ) / 
    ( 2 * exp( lambda ) ) )

}

log.prior <- function( lambda ) {

  return( 0 )

}

y <- c( 409., 400., 406., 399., 402., 406., 401., 403., 401., 403., 398.,
  403., 407., 402., 401., 399., 400., 401., 405., 402., 408., 399., 399.,
  402., 399., 397., 407., 401., 399., 401., 403., 400., 410., 401., 407.,
  423., 406., 406., 402., 405., 405., 409., 399., 402., 407., 406., 413.,
  409., 404., 402., 404., 406., 407., 405., 411., 410., 410., 410., 401.,
  402., 404., 405., 392., 407., 406., 404., 403., 408., 404., 407., 412.,
  406., 409., 400., 408., 404., 401., 404., 408., 406., 408., 406., 401.,
  412., 393., 437., 418., 415., 404., 401., 401., 407., 412., 375., 409.,
  406., 398., 406., 403., 404. )

print( mu <- mean( y ) )

# 404.59

sigma.star <- 1

print( sigma2.0 <- var( y ) )

# 41.8201

M <- 100000

B <- 500

mcmc.data.set.1 <- metropolis.example( y, mu, sigma.star, sigma2.0, M, B )

# 0.1747861

mean( mcmc.data.set.1[ , 2 ] )

# 42.19285

sd( mcmc.data.set.1[ , 2 ] )

# 6.038515

plot( 1:100, mcmc.data.set.1[ 1:100, 2 ], type = 'l',
  xlab = 'Iteration Number', ylab = 'Draws From p ( sigma^2 | y )' )

abline( h = 42.19285, lty = 2 )

rsichi2 <- function( n, nu, s2 ) {

  return( nu * s2 / rchisq( n, nu ) )

}

n <- 100

print( sigma2.hat <- ( n - 1 ) * var( y ) / n )

# 41.4019

direct.simulation <- rsichi2( 100000, 100, 41.4019 )

n * sigma2.hat / ( n - 2 )

# 42.24684

qqplot( mcmc.data.set.1[ , 2 ], direct.simulation, xlab = 'MCMC Draws',
  ylab = 'Scaled Inverse Chi Square Draws' )

abline( 0, 1 )

sigma.star <- 0.05

mcmc.data.set.2 <- metropolis.example( y, mu, sigma.star, sigma2.0, M, B )

# 0.8885174

plot( 1:100, mcmc.data.set.2[ 1:100, 2 ], type = 'l',
  xlab = 'Iteration Number', ylab = 'Draws From p ( sigma^2 | y )' )

2.4 * sd( mcmc.data.set.1[ , 1 ] )

# 0.3415815

sigma.star <- 0.34

mcmc.data.set.3 <- metropolis.example( y, mu, sigma.star, sigma2.0, M, B )

# 0.4420597

plot( 1:100, mcmc.data.set.3[ 1:100, 2 ], type = 'l',
  xlab = 'Iteration Number', ylab = 'Draws From p ( sigma^2 | y )' )

par( mfrow = c( 2, 2 ) )

plot( 1:( M + B + 1 ), mcmc.data.set.1[ , 2 ], type = 'l',
  xlab = 'Iteration Number', ylab = 'Draws from p ( sigma^2 | y )',
  main = 'sigma.star = 1.0' )

plot( density( mcmc.data.set.1[ , 2 ], adjust = 2 ), xlab = 'sigma^2',
  ylab = 'Density', main = '' )

acf( mcmc.data.set.1[ , 2 ], main = '' )

pacf( mcmc.data.set.1[ , 2 ], main = '' )

par( mfrow = c( 2, 2 ) )

plot( 1:( M + B + 1 ), mcmc.data.set.2[ , 2 ], type = 'l',
  xlab = 'Iteration Number', ylab = 'Draws from p ( sigma^2 | y )',
  main = 'sigma.star = 0.05' )

plot( density( mcmc.data.set.2[ , 2 ], adjust = 2 ), xlab = 'sigma^2',
  ylab = 'Density', main = '' )

acf( mcmc.data.set.2[ , 2 ], main = '' )

pacf( mcmc.data.set.2[ , 2 ], main = '' )

par( mfrow = c( 2, 2 ) )

plot( 1:( M + B + 1 ), mcmc.data.set.3[ , 2 ], type = 'l',
  xlab = 'Iteration Number', ylab = 'Draws from p ( sigma^2 | y )',
  main = 'sigma.star = 0.34' )

plot( density( mcmc.data.set.3[ , 2 ], adjust = 2 ), xlab = 'sigma^2',
  ylab = 'Density', main = '' )

acf( mcmc.data.set.3[ , 2 ], main = '' )

pacf( mcmc.data.set.3[ , 2 ], main = '' )

print( rho.1 <- acf( mcmc.data.set.1[ , 2 ], plot = F )$acf[ 2 ] )

# 0.7852402

print( rho.2 <- acf( mcmc.data.set.2[ , 2 ], plot = F )$acf[ 2 ] )

# 0.6560033

print( rho.3 <- acf( mcmc.data.set.3[ , 2 ], plot = F )$acf[ 2 ] )

# 0.9526438

print( rho.4 <- acf( mcmc.data.set.4[ , 2 ], plot = F )$acf[ 2 ] )

# 0.6365326

print( posterior.mean.1 <- mean( mcmc.data.set.1[ , 2 ] ) )

# 42.19285

print( posterior.sd.1 <- sd( mcmc.data.set.1[ , 2 ] ) )

# 6.038515

print( MCSE.1 <- ( posterior.sd.1 / sqrt( M ) ) * 
  sqrt( ( 1 + rho.1 ) / ( 1 - rho.1 ) ) )

# 0.05505566

print( posterior.mean.2 <- mean( mcmc.data.set.2[ , 2 ] ) )

# 42.30085

print( posterior.sd.2 <- sd( mcmc.data.set.2[ , 2 ] ) )

# 6.081415

print( MCSE.2 <- ( posterior.sd.2 / sqrt( M ) ) * 
  sqrt( ( 1 + rho.2 ) / ( 1 - rho.2 ) ) )

# 0.04219472

print( posterior.mean.3 <- mean( mcmc.data.set.3[ , 2 ] ) )

# 42.32991

print( posterior.sd.3 <- sd( mcmc.data.set.3[ , 2 ] ) )

# 6.09088

print( MCSE.3 <- ( posterior.sd.3 / sqrt( M ) ) * 
  sqrt( ( 1 + rho.3 ) / ( 1 - rho.3 ) ) )

# 0.1236811

print( posterior.mean.4 <- mean( mcmc.data.set.4[ , 2 ] ) )

# 42.19168

print( posterior.sd.4 <- sd( mcmc.data.set.4[ , 2 ] ) )

# 6.100921

print( MCSE.4 <- ( posterior.sd.4 / sqrt( M ) ) * 
  sqrt( ( 1 + rho.4 ) / ( 1 - rho.4 ) ) )

# 0.04093785

#       acceptance  posterior     MCSE of     posterior
# PDSD     rate       mean    posterior mean     SD

# 0.05      0.89      42.3        0.124         6.09
# 0.34      0.44      42.2        0.041         6.10
# 1.00      0.18      42.2        0.055         6.04

print( correct.posterior.mean <- sigma2.hat * n / ( n - 2 ) )

# 42.24684

print( correct.posterior.sd <- sigma2.hat * sqrt( 2 * n^2 /
  ( ( n - 2 )^2 * ( n - 4 ) ) ) )

# 6.097806