This is a set of S+, BUGS, Maple, and C programs to perform the calculations in Chapter 2 of Draper D (1998), Bayesian Hierarchical Modeling. You will need to edit the text below to extract the programs you want. A long line of asterisks is used to separate each program. Below in order you will find the following: * S+ functions to do Hastings sampling in the model (2.7) * S+ function (preCODA) that prepares the output of all the S+ samplers supplied here for CODA diagnostic processing * S+ functions to do Metropolis sampling in the model (2.7) * S+ functions to do Gibbs sampling in the model (2.31, 2.32) * S+ functions to do Metropolis sampling in the model (2.46, 2.48) * BUGS files for Gibbs sampling in the t example of Section 2.6 * Maple code to find the posterior mode and compute its approximate covariance matrix (based on the Hessian at the mode), in the NB10 t model of Section 2.6 * C functions to do Hastings sampling in the model (2.7) Please notify David Draper (d.draper@maths.bath.ac.uk) of any problems you encounter with any of these programs. 27 Jan 98 **************************************************************************** # S+ functions to do Hastings sampling in the model (2.7): # # sigma2 ~ SI-chisq( nu.p, sigma2.p ) # ( y_i | sigma2 ) ~IID N( mu, sigma2 ), i = 1, ..., n # # (written by DD, with trick to avoid S+ memory- # management problems due to Brian Ripley) # # Inputs: # # y = data vector, of length n = sample size # mu = known mean in Gaussian likelihood # nu.p = prior effective sample size # sigma2.p = prior estimate of sigma2 # sigma2.0 = initial value for sigma2 in Hastings iterations # nu.star = scaling factor for Hastings proposal distribution # (affects the acceptance rate R; to increase R, increase # nu.star). nu.star must be > 2 in this implementation; # values near 20-25 lead to good mixing # n.burnin = length of burn-in period # n.monitor = length of monitoring period # n.thin = thinning constant (only every n.thin-th iteration # in the monitoring period will be written to disk) # seed = random number seed (for generating repeatable # sequences of Hastings iterations); must be an integer # from 0 to 1000 # output.file.prefix = character string naming where you want # the MCMC data set to go; for example, output.file.prefix # = "NB10" would write the MCMC data set to the file # "NB10.d" # # Outputs: # # Acceptance rate R returned when iterations are finished # A file called paste( output.file.prefix, ".d", sep = "" ) is # written (in the same directory where S+ has been called) # containing one row for each monitored iteration and three # columns: the monitored iteration number (from 1 to # n.monitor / n.thin), the simulated draw from the posterior # for sigma2 for that iteration, and the corresponding # simulated draw from the predictive distribution for a new # y.star. If the output file exists before the function is # invoked, it will be over-written hastings.gaussian.variance <- function( y, mu, nu.p, sigma2.p, sigma2.0, nu.star, n.burnin, n.monitor, n.thin, seed, output.file.prefix ) { # Main routine # # Includes trick due to Ripley to overcome S+ memory-management # problems sigma2.old <- sigma2.0 R <- 0 write( c( sigma2.old, R ), "loop.result", append = F ) set.seed( seed ) for ( i in 1:( n.burnin + n.monitor ) ) { null <- loop( y, mu, nu.p, sigma2.p, nu.star, output.file.prefix, i ) } loop.result <- scan( "loop.result" ) R <- loop.result[2] return( R / n.monitor ) } loop <- function( y, mu, nu.p, sigma2.p, nu.star, output.file.prefix, i ) { # Ripley idea: put everything inside an explicit loop into a # function that returns nothing, reading from and writing to # disk as needed to maintain communication. This will keep S+ # from accumulating dynamic memory as it goes around the loop loop.result <- scan( "loop.result" ) sigma2.old <- loop.result[1] R <- loop.result[2] sigma2.star <- PD.sim( 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 ) y.new <- rnorm( 1, mu, sqrt( sigma2.new ) ) 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( c( sigma2.new, y.new ), digits = 5 ) ), paste( output.file.prefix, ".d", sep = "" ), ncol = 3, append = ( i > n.burnin + n.thin ) ) sigma2.old <- sigma2.new write( c( sigma2.old, R ), "loop.result", append = F ) return( NULL ) } PD.sim <- function( nu, sigma2 ) { # Proposal distribution simulation return( nu * sigma2 / rchisq( 1, nu ) ) } alpha <- function( sigma2.old, sigma2.new, y, mu, nu.p, sigma2.p, nu.star ) { # Acceptance probability calculation return( min( 1, exp( log.post( sigma2.new, y, mu, nu.p, sigma2.p ) + log.PD( sigma2.old, nu.star, sigma2.new ) - log.PD( sigma2.new, nu.star, sigma2.old ) - log.post( sigma2.old, y, mu, nu.p, sigma2.p ) ) ) ) } log.post <- function( sigma2, y, mu, nu.p, sigma2.p ) { # log( posterior ) calculation return( log.lik( sigma2, y, mu ) + log.prior( sigma2, nu.p, sigma2.p ) ) } log.lik <- function( sigma2, y, mu ) { # log( likelihood ) calculation n <- length( y ) return( ( - n / 2 ) * log( sigma2 ) - sum( ( y - mu )^2 ) / ( 2 * sigma2 ) ) } log.prior <- function( sigma2 , nu.p, sigma2.p ) { # log( prior ) calculation return( ( -1 - nu.p / 2 ) * log( sigma2 ) - nu.p * sigma2.p / ( 2 * sigma2 ) ) } log.PD <- function( sigma2, nu.star, sigma2.star ) { # log( proposal distribution ) calculation return( ( nu.star / 2 ) * log( sigma2.star ) - ( 1 + nu.star / 2 ) * log( sigma2 ) - ( nu.star - 2 ) * sigma2.star / ( 2 * sigma2 ) ) } **************************************************************************** # S+ function to prepare the output of all the S+ and C samplers # supplied here for CODA diagnostic processing # # (written by DD) # # Inputs: # # n.monitor = length of monitoring period used to produce # input file # n.thin = thinning constant used to produce input file # p = number of variables monitored in input file (this is # one less than the number of columns in that file, since # the first column contains the iteration number) # input.file.prefix = character string naming the file with # the MCMC data set in it; for example, input.file.prefix # = "NB10" assumes that the MCMC data set is in a file # called "NB10.d" in the directory where S+ is running # var.names = character vector of length p supplying (in order) # the names of the variables monitored in MCMC data set # output.file.prefix = character string naming where you want # the output of preCODA to go; for example, # output.file.prefix = "NB10" will write the MCMC data set # out into a file called "NB10.out" in a format suitable for # reading by CODA, and will also create a file called # "NB10.ind" that CODA uses to figure out the format of # "NB10.out" # # Outputs: # # Files called paste( output.file.prefix, ".out", sep = "" ) # and paste( output.file.prefix, ".ind", sep = "" ) # are written (in the same directory where S+ has been # called); see under output.file.prefix above for a # description of the two files. If these files exist before # preCODA is invoked, they will be over-written preCODA <- function( n.monitor, n.thin, p, input.file.prefix, var.names, output.file.prefix ) { # Prepares MCMC output for reading by CODA MCMC.data <- matrix( scan( paste( input.file.prefix, ".d", sep = "" ) ), n.monitor / n.thin, p + 1, byrow = T ) for ( i in 1:p ) { write( t( cbind( MCMC.data[, 1], MCMC.data[, i + 1] ) ), paste( output.file.prefix, ".out", sep = "" ), ncol = 2, append = ( i > 1 ) ) write( c( var.names[i], as.character( 1 + ( i - 1 ) * n.monitor / n.thin ), as.character( i * n.monitor / n.thin ) ), paste( output.file.prefix, ".ind", sep = "" ), ncol = 3, append = ( i > 1 ) ) } return( paste( output.file.prefix, ".out", " written", sep = "" ) ) } **************************************************************************** # # S+ functions to do Metropolis sampling in the model (2.7): # # sigma2 ~ SI-chisq( nu.p, sigma2.p ) # ( y_i | sigma2 ) ~IID N( mu, sigma2 ), i = 1, ..., n # # Inputs: # # y = data vector, of length n = sample size # mu = known mean in Gaussian likelihood # nu.p = prior effective sample size # sigma2.p = prior estimate of sigma2 # sigma2.0 = initial value for sigma2 in Metropolis iterations # kappa = scaling factor for Metropolis proposal distribution # (affects the acceptance rate R; to increase R, decrease # kappa) # n.burnin = length of burn-in period # n.monitor = length of monitoring period # n.thin = thinning constant (only every n.thin-th iteration # in the monitoring period will be written to disk) # seed = random number seed (for generating repeatable # sequences of Hastings iterations); must be an integer # from 0 to 1000 # output.file.prefix = character string naming where you want # the MCMC data set to go; for example, output.file.prefix # = "NB10" would write the MCMC data set to the file # "NB10.d" # # Outputs: # # Acceptance rate R returned when iterations are finished # A file called paste( output.file.prefix, ".d", sep = "" ) # is written (in the same directory where S+ has been # called) containing one row for each monitored iteration # and four columns: the monitored iteration number (from 1 # to n.monitor / n.thin), the simulated draws from the # posterior for lambda = log(sigma2) and sigma2 for that # iteration, and the corresponding simulated draw from the # predictive distribution for a new y.star. If the output # file exists before the function is invoked, it will be # over-written metropolis.gaussian.variance <- function( y, mu, nu.p, sigma2.p, sigma2.0, kappa, n.burnin, n.monitor, n.thin, seed, output.file.prefix ) { # Main routine lambda.old <- log( sigma2.0 ) R <- 0 write( c( lambda.old, R ), "loop.result", append = F ) set.seed( seed ) for ( i in 1:( n.burnin + n.monitor ) ) { null <- loop( y, mu, nu.p, sigma2.p, kappa, output.file.prefix, i ) } loop.result <- scan( "loop.result" ) R <- loop.result[2] return( R / n.monitor ) } loop <- function( y, mu, nu.p, sigma2.p, kappa, output.file.prefix, i ) { loop.result <- scan( "loop.result" ) lambda.old <- loop.result[1] R <- loop.result[2] n <- length( y ) lambda.star <- PD.sim( lambda.old, kappa, n ) u <- runif( 1 ) b <- ( u <= alpha( lambda.old, lambda.star, y, mu, nu.p, sigma2.p ) ) lambda.new <- lambda.star * b + lambda.old * ( 1 - b ) y.new <- rnorm( 1, mu, sqrt( exp( lambda.new ) ) ) 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( c( lambda.new, exp( lambda.new ), y.new ), digits = 5 ) ), paste( output.file.prefix, ".d", sep = "" ), ncol = 4, append = ( i > n.burnin + n.thin ) ) lambda.old <- lambda.new write( c( lambda.old, R ), "loop.result", append = F ) return( NULL ) } PD.sim <- function( lambda, kappa, n ) { # Proposal distribution simulation return( rnorm( 1, lambda, sqrt( 2.0 * kappa / n ) ) ) } alpha <- function( lambda.old, lambda.new, y, mu, nu.p, sigma2.p ) { # Acceptance probability calculation return( min( 1, exp( log.post( lambda.new, y, mu, nu.p, sigma2.p ) - log.post( lambda.old, y, mu, nu.p, sigma2.p ) ) ) ) } log.post <- function( lambda, y, mu, nu.p, sigma2.p ) { # log( posterior ) calculation return( log.lik( lambda, y, mu ) + log.prior( lambda, nu.p, sigma2.p ) ) } log.lik <- function( lambda, y, mu ) { # log( likelihood ) calculation n <- length( y ) sigma2 <- exp( lambda ) return( ( - n / 2 ) * log( sigma2 ) - sum( ( y - mu )^2 ) / ( 2 * sigma2 ) ) } log.prior <- function( lambda , nu.p, sigma2.p ) { # log( prior ) calculation (including Jacobian) return( ( - nu.p / 2 ) * lambda - nu.p * sigma2.p / ( 2 * exp( lambda ) ) ) } **************************************************************************** # S+ functions to do Gibbs sampling in the model (2.31, 2.32): # # mu ~ U( L, H ) # ( sigma | mu ) ~ U( 0, min( mu - L, H - mu ) ) # ( y_i | mu, sigma ) ~IID U( mu - sigma, mu + sigma ), # i = 1, ..., n # # (written by DD) # # Inputs: # # y = data vector, of length n = sample size # L = known lower bound for true range # H = known upper bound for true range # sigma.0 = initial value for sigma in Gibbs iterations # n.burnin = length of burn-in period # n.monitor = length of monitoring period # n.thin = thinning constant (only every n.thin-th iteration # in the monitoring period will be written to disk) # seed = random number seed (for generating repeatable # sequences of Gibbs iterations); must be an integer from 0 # to 1000 # output.file.prefix = character string naming where you want # the MCMC data set to go; for example, output.file.prefix # = "ammonite" would write the MCMC data set to the file # "ammonite.d" # # Outputs: # # A file called paste( output.file.prefix, ".d", sep = "" ) is # written (in the same directory where S+ has been called) # containing one row for each monitored iteration and five # columns: the monitored iteration number (from 1 to # n.monitor / n.thin), the simulated draws from the marginal # posteriors for mu and sigma for that iteration, and the # corresponding simulated draws from the marginal posteriors # for A = mu - sigma and B = mu + sigma. If the output file # exists before the function is invoked, it will be # over-written gibbs.uniform <- function( y, L, H, sigma.0, n.burnin, n.monitor, n.thin, seed, output.file.prefix ) { # Main routine sigma.previous <- sigma.0 for ( i in 1:( n.burnin + n.monitor ) ) { mu.hat <- sample.mu( y, sigma.previous, L, H ) sigma.hat <- sample.sigma( y, mu.hat, L, H ) A.hat <- mu.hat - sigma.hat B.hat <- mu.hat + sigma.hat if ( ( i > n.burnin ) & ( ( i - n.burnin ) %% n.thin == 0 ) ) write( c( ( i - n.burnin ) / n.thin, signif( c( mu.hat, sigma.hat, A.hat, B.hat ), digits = 6 ) ), paste( output.file.prefix, ".d", sep = "" ), ncol = 5, append = ( i > n.burnin + n.thin ) ) sigma.previous <- sigma.hat } return( paste( n.monitor / n.thin, " iterations written to ", output.file.prefix, ".d", sep = "" ) ) } sample.mu <- function( y, sigma, L, H ) { # Sampling from the full conditional for mu c.1 <- max( L + sigma, max( y ) - sigma ) c.2 <- min( y ) + sigma c.3 <- min( H - sigma, c.2 ) sigma.star <- ( H - L ) / 2.0 c.4 <- 1.0 / log( sigma.star^2 / ( ( c.1 - L ) * ( H - c.3 ) ) ) mu.star <- ( H + L ) / 2.0 U <- runif( 1 ) if ( sigma < mu.star - min( y ) ) { return( L + ( c.1 - L ) * ( ( c.2 - L ) / ( c.1 - L ) )^U ) } else if ( U < c.4 * log( sigma.star / ( c.1 - L ) ) ) { return( L + ( c.1 - L ) * exp( U / c.4 ) ) } else { return( H - sigma.star / exp( U / c.4 - log( sigma.star / ( c.1 - L ) ) ) ) } } sample.sigma <- function( y, mu, L, H ) { # Sampling from the full conditional for sigma n <- length( y ) c.5 <- max( mu - min( y ), max( y ) - mu ) c.6 <- min( mu - L, H - mu ) U <- runif( 1 ) return( ( ( 1.0 - U ) * c.5^( 1.0 - n ) + U * c.6^( 1.0 - n ) )^ ( 1.0 / ( 1.0 - n ) ) ) } **************************************************************************** # # S+ functions to do Metropolis sampling in the model (2.46, # 2.48): # # theta = ( mu, eta, lambda ) # ( theta ) ~ N( mu.0, sigma.mu2 ) * N( eta.0, sigma.eta2 ) * # N( lambda.0, sigma.lambda2 ) # ( y_i | theta ) ~IID t.exp( lambda ) ( mu, exp( 2 * eta ) ), # i = 1, ..., n # # Inputs: # # y = data vector, of length n = sample size # mu.0 = prior mean for mu # sigma.mu = prior SD for mu # eta.0 = prior mean for eta # sigma.eta = prior SD for eta # lambda.0 = prior mean for lambda # sigma.lambda = prior SD for lambda # kappa = scaling factor for Metropolis proposal distribution # (affects the acceptance rate R; to increase R, decrease # kappa) # Sigma = proposal distribution covariance matrix # theta.0 = initial value for ( mu, eta, lambda ) in # Metropolis iterations # n.burnin = length of burn-in period # n.monitor = length of monitoring period # n.thin = thinning constant (only every n.thin-th iteration # in the monitoring period will be written to disk) # seed = random number seed (for generating repeatable # sequences of Hastings iterations); must be an integer # from 0 to 1000 # output.file.prefix = character string naming where you want # the MCMC data set to go; for example, output.file.prefix # = "NB10" would write the MCMC data set to the file # "NB10.d" # # Outputs: # # Acceptance rate R returned when iterations are finished # A file called paste( output.file.prefix, ".d", sep = "" ) is # written (in the same directory where S+ has been called) # containing one row for each monitored iteration and six # columns: the monitored iteration number (from 1 to # n.monitor/n.thin), the simulated draws from the posterior # for theta = ( mu, eta, lambda ) for that iteration, and # the corresponding simulated draws from the posterior for # sigma = exp( eta ) and nu = exp( lambda ). If the output # file exists before the function is invoked, it will be # over-written metropolis.t <- function( y, mu.0, sigma.mu, eta.0, sigma.eta, lambda.0, sigma.lambda, kappa, Sigma, theta.0, n.burnin, n.monitor, n.thin, seed, output.file.prefix ) { # Main routine theta.old <- theta.0 p <- length( theta.old ) R <- 0 write( c( theta.old, R ), "loop.result", append = F ) set.seed( seed ) L <- t( chol( Sigma ) ) L.kappa <- sqrt( kappa ) * L for ( i in 1:( n.burnin + n.monitor ) ) { null <- loop( p, L.kappa, y, mu.0, sigma.mu, eta.0, sigma.eta, lambda.0, sigma.lambda, output.file.prefix, i ) } loop.result <- scan( "loop.result" ) R <- loop.result[p + 1] return( R / n.monitor ) } loop <- function( p, L.kappa, y, mu.0, sigma.mu, eta.0, sigma.eta, lambda.0, sigma.lambda, output.file.prefix, i ) { loop.result <- scan( "loop.result" ) theta.old <- loop.result[1:p] R <- loop.result[p+1] n <- length( y ) theta.star <- PD.sim( theta.old, p, L.kappa ) u <- runif( 1 ) b <- ( u <= alpha( theta.old, theta.star, y, mu.0, sigma.mu, eta.0, sigma.eta, lambda.0, sigma.lambda ) ) theta.new <- theta.star * b + theta.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( c( theta.new, exp( theta.new[c( 2, 3 )] ) ), digits = 5 ) ), paste( output.file.prefix, ".d", sep = "" ), ncol = p + 3, append = ( i > n.burnin + n.thin ) ) theta.old <- theta.new write( c( theta.old, R ), "loop.result", append = F ) return( NULL ) } PD.sim <- function( theta, p, L.kappa ) { # Proposal distribution simulation Z <- matrix( rnorm( p ), p, 1 ) Mu <- matrix( theta, p, 1 ) theta.star <- c( Mu + ( L.kappa %*% Z ) ) return( theta.star ) } alpha <- function( theta.old, theta.new, y, mu.0, sigma.mu, eta.0, sigma.eta, lambda.0, sigma.lambda ) { # Acceptance probability calculation return( min( 1, exp( log.post( theta.new, y, mu.0, sigma.mu, eta.0, sigma.eta, lambda.0, sigma.lambda ) - log.post( theta.old, y, mu.0, sigma.mu, eta.0, sigma.eta, lambda.0, sigma.lambda ) ) ) ) } log.post <- function( theta, y, mu.0, sigma.mu, eta.0, sigma.eta, lambda.0, sigma.lambda ) { # log( posterior ) calculation return( log.prior( theta, mu.0, sigma.mu, eta.0, sigma.eta, lambda.0, sigma.lambda ) + log.lik( theta, y ) ) } log.prior <- function( theta, mu.0, sigma.mu, eta.0, sigma.eta, lambda.0, sigma.lambda ) { # log( prior ) calculation (including Jacobian) mu <- theta[1] eta <- theta[2] lambda <- theta[3] return( -0.5 * ( ( mu - mu.0 ) / sigma.mu )^2 - 0.5 * ( ( eta - eta.0 ) / sigma.eta )^2 - 0.5 * ( ( lambda - lambda.0 ) / sigma.lambda )^2 ) } log.lik <- function( theta, y ) { # log( likelihood ) calculation mu <- theta[1] eta <- theta[2] lambda <- theta[3] n <- length( y ) return( n * lgamma( 0.5 * ( exp( lambda ) + 1.0 ) ) - n * eta - n * lgamma( 0.5 * exp( lambda ) ) - 0.5 * n * lambda - 0.5 * ( exp( lambda ) + 1.0 ) * sum( log( 1.0 + exp( - ( lambda + 2.0 * eta ) ) * ( y - mu )^2 ) ) ) } **************************************************************************** Next comes five files needed by BUGS to do the Gibbs sampling in the t example of Section 2.6. The files are delimited by [[ ]] with a blank line before and after [[ and ]]; in each case below you should strip away the line naming the file (containing [[), the blank line immediately following that line, the line consisting just of ]] at the end of the file, and the blank line immediately preceding that line. the file nb10.3.bug [[ model nb10.3; # Naming the model. const n = 100, g = 101; # Defining the constants. var mu, tau, u, grid[g], # Specifying the variables. nu, y[n], sigma; data y in "nb10-y.dat", # Reading in the data grid in "nb10-grid.dat"; # and initial values. inits in "nb10.3.in"; # { mu ~ dnorm( 0.0, 4.0E-6 ); # Specifying the priors tau ~ dgamma( 0.25, 0.12 ); # for mu, sigma, and nu u ~ dcat( grid[] ); # (see text). nu <- 1.0 + u / 7.0; # for ( i in 1:n ) { # y[i] ~ dt( mu, tau, nu ); # Specifying the likelihood. } # sigma <- 1.0 / sqrt( tau ); # Defining a derived quantity. } ]] the file NB10-y.dat [[ 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 ]] the file NB10-grid.dat [[ 0.00105225 0.00169802 0.00250893 0.00346489 0.00453795 0.00569633 0.00690769 0.00814134 0.00936973 0.0105692 0.0117201 0.0128071 0.0138185 0.0147461 0.0155848 0.0163318 0.0169867 0.0175504 0.0180256 0.0184155 0.0187245 0.0189572 0.0191188 0.0192144 0.0192493 0.0192288 0.0191579 0.0190417 0.0188848 0.0186918 0.0184668 0.0182138 0.0179365 0.0176383 0.0173223 0.0169915 0.0166484 0.0162954 0.0159347 0.0155684 0.0151981 0.0148254 0.0144518 0.0140785 0.0137067 0.0133373 0.0129712 0.0126091 0.0122518 0.0118997 0.0115534 0.0112132 0.0108796 0.0105527 0.0102328 0.00992013 0.00961471 0.00931669 0.0090261 0.00874299 0.00846733 0.00819911 0.00793827 0.00768476 0.00743848 0.00719935 0.00696724 0.00674205 0.00652365 0.00631191 0.00610668 0.00590783 0.00571521 0.00552866 0.00534804 0.0051732 0.00500399 0.00484024 0.00468182 0.00452857 0.00438034 0.00423698 0.00409836 0.00396432 0.00383474 0.00370947 0.00358837 0.00347132 0.00335819 0.00324885 0.00314317 0.00304105 0.00294236 0.002847 0.00275484 0.00266579 0.00257974 0.00249658 0.00241623 0.00233859 0.000977389 ]] the file NB10.3.cmd [[ compile( "nb10.5.bug" ) update( 1000 ) monitor( mu, 14 ) monitor( sigma, 14 ) monitor( nu, 14 ) update( 70000 ) q( ) ]] the file NB10.3.in [[ list( mu = 405.03, tau = 0.18233, u = 30 ) ]] **************************************************************************** # # Maple code to find the posterior mode and compute its # approximate covariance matrix (based on the Hessian at the # mode), in the NB10 $t$ model of Section 2.6. # # (written by DD, with some help from Riccardo Gatto) # # Input: # # Reads from a file called "nb10.dat" which contains the NB10 # measurements, 1 per line for 100 lines # # Outputs: # # Obtains the MAP (maximum a posteriori) equations by # differentiating the log posterior function symbolically, # solves them numerically to get the mode, calculates the # Hessian symbolically, and evaluates it numerically at the mode # n := 100; readlib( readdata ); y := readdata( `nb10.dat`, float, 1 ); mu.0 := 0.0; sigma.mu := 500.0; eta.0 := 3.80; sigma.eta := 0.77; lambda.0 := 1.84; sigma.lambda := 0.59; log.prior := -0.5 * ( ( mu - mu.0 ) / sigma.mu )^2 - 0.5 * \ ( ( eta - eta.0 ) / sigma.eta )^2 - 0.5 * ( ( lambda - \ lambda.0 ) / sigma.lambda )^2; log.likelihood := n * lnGAMMA( 0.5 * ( exp( lambda ) + \ 1.0 ) ) - n * eta - n * lnGAMMA( 0.5 * exp( lambda ) ) - \ 0.5 * n * lambda - 0.5 * ( exp( lambda ) + 1.0 ) * sum( \ log( 1.0 + exp( - ( lambda + 2.0 * eta ) ) * ( y[i] - \ mu )^2 ), i = 1 .. n ); log.posterior := log.prior + log.likelihood; map.eq1 := diff( log.posterior, mu ); map.eq2 := diff( log.posterior, eta ); map.eq3 := diff( log.posterior, lambda ); fsolve( { map.eq1, map.eq2, map.eq3 }, { mu, eta, lambda }, \ { mu = 403 .. 406, eta = 1 .. 3, lambda = 0.5 .. 3.0 } ); # At this point Maple takes about 2.5 seconds at 333Mhz to # iteratively solve the MAP equations, obtaining the values of # mu, eta, and lambda listed below. with( linalg ); H := hessian( log.posterior, [mu, eta, lambda] ); mu := 404.2956374; eta := 1.346258072; lambda := 1.259790967; He := matrix( 3, 3, ( i,j ) -> 0 ): for i from 1 to 3 do for j from 1 to 3 do He[i,j] := evalf( H[i,j] ): od: od: Sigma := inverse( - He ); # At this point Maple returns the covariance matrix: # # [ .2159336246 .002989379323 .008083662383] # [ ] # Sigma := [.002989379323 .01193790730 .01486075432 ] # [ ] # [.008083662383 .01486075432 .07490518257 ] **************************************************************************** /* * C functions to do Hastings sampling in the model (2.7): * * sigma2 ~ SI-chisq( nu.p, sigma2.p ) * ( y_i | sigma2 ) ~IID N( mu, sigma2 ), i = 1, ..., n * * (random number generators written by William Browne; Hastings * code written by DD and Dimitris Fouskakis and edited by WB) * * Inputs are contained in a file called * "hastings.gaussian.variance.in", in the directory where the * program is to be run, with the contents listed below (one * input per line, in the order listed, from mu to n; then the * data vector y is input, with one element on each line) * * mu = known mean in Gaussian likelihood * nu.p = prior effective sample size * sigma2.p = prior estimate of sigma2 * sigma2.0 = initial value for sigma2 in Hastings iterations * nu.star = scaling factor for Hastings proposal distribution * (affects the acceptance rate R; to increase R, increase * nu.star). nu.star must be > 2 in this implementation; * values near 20-25 lead to good mixing * n.burnin = length of burn-in period * n.monitor = length of monitoring period * n.thin = thinning constant (only every n.thin-th iteration * in the monitoring period will be written to disk) * seed = random number seed (for generating repeatable * sequences of Hastings iterations); must be an integer * n = sample size = length of data vector * y = data vector * * Outputs: * * Acceptance rate R printed when iterations are finished * A file called "hastings.gaussian.variance.d" is written (in * the same directory where the program has been run) * containing one row for each monitored iteration and three * columns: the monitored iteration number (from 1 to * n.monitor / n.thin), the simulated draw from the posterior * for sigma2 for that iteration, and the corresponding * simulated draw from the predictive distribution for a new * y.star. If "hastings.gaussian.variance.d" exists before * the function is invoked, it will be over-written */ #include #include /* * Defined constants */ #define SEED1 13 #define SEED2 4 #define SEED3 1972 #define PI 3.1415927 #define E 2.71828182 long int seed1; /* Seeds declared externally to avoid passing */ long int seed2; /* each time wichmann is called */ long int seed3; double wichmann() /* * Random number generator for U(0,1) distribution. */ { extern long int seed1, seed2, seed3; double random; seed1 = (171 * seed1)%30269; seed2 = (172 * seed2)%30307; seed3 = (170 * seed3)%30323; random = fmod(seed1/30269.0 + seed2/30307.0 + seed3/30323.0, 1.0); return random; } double rexp(double lambda) /* * Generates from an exponential distribution */ { double random, uniform; uniform = wichmann(); random = - (1/lambda) * log(uniform); return random; } double rgamma1(double alpha) /* * Generates from a gamma distribution with alpha < 1 */ { double uniform0, uniform1; double random, x; int done = 0; uniform0 = wichmann(); uniform1 = wichmann(); if (uniform0 > E/(alpha + E)) { random = -log((alpha + E)*(1-uniform0)/(alpha*E)); if ( uniform1 > pow(random,alpha - 1)) return -1; else return random; } else { x = (alpha + E) * uniform0 / E; random = pow(x,1/alpha); if ( uniform1 > exp(-random)) return -1; else return random; } } double rgamma2(double alpha) /* * Generates from a gamma distribution with alpha > 1 */ { double uniform1,uniform2; double c1,c2,c3,c4,c5,w; double random; int done = 1; c1 = alpha - 1; c2 = (alpha - 1/(6 * alpha))/c1; c3 = 2 / c1; c4 = c3 + 2; c5 = 1 / sqrt(alpha); do { uniform1 = wichmann(); uniform2 = wichmann(); if (alpha > 2.5) { uniform1 = uniform2 + c5 * (1 - 1.86 * uniform1); } } while ((uniform1 >= 1) || (uniform1 <= 0)); w = c2 * uniform2 / uniform1; if ((c3 * uniform1 + w + 1/w) > c4) { if ((c3 * log(uniform1) - log(w) + w) >= 1) { done = 0; } } if (done == 0) return -1; random = c1 * w; return random; } double rgamma(double alpha, double beta) /* * Generates from a general gamma(alpha,beta) distribution */ { double random; if (alpha < 1) do { random = rgamma1(alpha)/beta; } while (random < 0 ); if (alpha == 1) random = rexp(1)/beta; if (alpha > 1) do { random = rgamma2(alpha)/beta; } while (random < 0); return random; } double rstd_normal() /* * Generates from a standard normal(0,1) distribution */ { double uniform1,uniform2; double theta,r; double random; uniform1 = wichmann(); uniform2 = wichmann(); theta = 2 * PI * uniform1; r = sqrt(2 * ( - log(uniform2))); random = r * cos(theta); return random; } double rnormal(double mean, double sd) /* * Generates from a general normal(mu,sigma2) distribution */ { double random; random = mean + sd * rstd_normal(); return random; } /* * Hastings code begins here */ double PD_sim( nu, sigma2 ) double nu, sigma2; /* * Proposal distribution simulation */ { double result = nu * sigma2 / rgamma( nu / 2.0, 0.5 ); return result; } double log_lik( sigma2, y, mu, n) double sigma2, *y, mu; long int n; /* * log( likelihood ) calculation */ { int i; double result = ( - n / 2.0 ) * log( sigma2 ); for ( i = 0; i < n; i++ ) { result = result - pow( y[i] - mu, 2.0 ) / ( 2.0 * sigma2 ); } return result; } double log_prior( sigma2 , nu_p, sigma2_p ) double sigma2 , nu_p, sigma2_p; /* * log( prior ) calculation */ { double result = ( -1.0 - nu_p / 2.0 ) * log (sigma2 ) - nu_p * sigma2_p / ( 2.0 * sigma2 ); return result; } double log_post( sigma2, y, mu, nu_p, sigma2_p, n ) double sigma2, *y, mu, nu_p, sigma2_p; long int n; /* * log( posterior ) calculation */ { double result = log_lik( sigma2, y, mu, n ) + log_prior( sigma2, nu_p, sigma2_p ); return result; } double log_PD( sigma2, nu_star, sigma2_star ) double sigma2, nu_star, sigma2_star; /* * log( proposal distribution ) calculation */ { double result = ( nu_star / 2.0 ) * log( sigma2_star ) - ( 1.0 + nu_star / 2.0 ) * log( sigma2 ) - (nu_star - 2.0) * sigma2_star / ( 2.0 * sigma2 ); return result; } double alpha( sigma2_old, sigma2_new, y, mu, nu_p, sigma2_p, nu_star, n ) double sigma2_old, sigma2_new, *y, mu, nu_p, sigma2_p, nu_star; long int n; /* * acceptance probability calculation */ { double ratio = exp( log_post( sigma2_new, y, mu, nu_p, sigma2_p, n ) + log_PD( sigma2_old, nu_star, sigma2_new ) - log_PD( sigma2_new, nu_star, sigma2_old ) - log_post( sigma2_old, y, mu, nu_p, sigma2_p, n ) ); if ( ratio > 1.0 ) return 1.0; else return ratio; } void main( ) /* * Main routine */ { long int burnin, monitor, thin, n_run, i, n; double accept = 0.0, *y, sigma2_old, sigma2_star, nu_star, u, mu, nu_p, sigma2_p, sigma2_new, y_new, arg2, b; FILE *fp,*fpout; fp = fopen( "hastings.gaussian.variance.in", "r" ); fscanf(fp,"%lf",&mu); fscanf(fp,"%lf",&nu_p); fscanf(fp,"%lf",&sigma2_p); fscanf(fp,"%lf",&sigma2_old); fscanf(fp,"%lf",&nu_star); fscanf(fp,"%ld",&burnin); fscanf(fp,"%ld",&monitor); fscanf(fp,"%ld",&thin); fscanf(fp,"%ld",&seed1); fscanf(fp,"%ld",&n); y = (double *)calloc(n,sizeof(double)); for(i=0;i burnin ) accept = accept + b; if ( ( i > burnin ) && ( (i - burnin) % thin == 0 ) ) fprintf( fpout, "%d %10.4f %10.4f\n", ( i - burnin ) / thin, sigma2_new, y_new ); sigma2_old = sigma2_new; } printf( "%f\n", accept / monitor ); } ****************************************************************************