#------------------------------------------------------------ 
  # R functions to do Metropolis-Hastings sampling 
  # for the NB10 data
  # 	
  # prior:       sigma2 ~ SI-chisq( nu.p, sigma2.p )
  # data model:  ( y_i | sigma2 ) ~IID N( mu, sigma2 ), i = 1, ..., n
  #
  # 
#------------------------------------------------------------ 
MH.normal.variance = function( y, mu, nu.p, sigma2.p, 
  sigma2.0, nu.star, n.burnin, n.monitor, n.thin, seed ) {

  # Main routine
 
  sigma2.old = sigma2.0
  R=0
  set.seed( seed )

  for ( i in 1:n.monitor) {
    sigma2.star = generate.CGD( nu.star, ( nu.star - 2 ) * sigma2.old / nu.star )
    u = runif( 1 )
    b = ( u <= alpha( sigma2.old, sigma2.star, y, mu, nu.p, sigma2.p, nu.star ) )
    sigma2.new = sigma2.star * b + sigma2.old * ( 1 - b )
    if ( i > n.burnin ) R = R + b
    if ( ( i > n.burnin ) & ( ( i - n.burnin ) %% n.thin == 0 ) ) 
    write( c( ( i - n.burnin ) / n.thin, signif(sigma2.new, digits = 5  )), 
           file="nb10.output", ncol = 2, append = T )
    }
  return( R / (n.monitor-n.burnin) )
}
#--------------------------------------------------------
generate.CGD = function( nu, sigma2 ) {

# candidate generating distribution
	
  return(               )

}
#---------------------------------------------------------
alpha = function( sigma2.old, sigma2.star, y, mu, nu.p, sigma2.p, nu.1 ) {
	
  # Acceptance probability calculation

  return(              )
}
#-----------------------------------------------------------
log.post = function( sigma2, y, mu, nu.p, sigma2.p ) {

  # log( posterior ) calculation

  return(                   )
}
#------------------------------------------------------------

log.lik = function( sigma2, y, mu ) {

  # log( likelihood ) calculation

  n = length( y )
  return(                         )
}
#------------------------------------------------------------
log.prior = function( sigma2 , nu.p, sigma2.p ) {

  # log( prior ) calculation

  return(                               )

}
#------------------------------------------------------------
log.CGD = function( sigma2.old, nu.1, sigma2 ) {

  # log( candidate generating density ) calculation

  return(                               )
}