rm( list = ls( ) ) # change directory to where the code and data files are source( "LNDPMM-code-R.txt" ) IB <- 500 B <- 2000 y <- scan( "sim-data.txt" ) a.alpha <- 2.0 b.alpha <- 1.8 a.mu <- 0.0 b.mu <- 4.0 a.tau2 <- 2.0 b.tau2 <- 4.0 a.phi <- 2.0 b.phi <- 4.0 # unbuffer the output mcmc( IB, B, y, a.alpha, b.alpha, a.mu, b.mu, a.tau2, b.tau2, a.phi, b.phi ) y.0.post <- scan( "post_pred_y.txt" ) y.0.prior <- scan( "prior_pred_y.txt" ) hist( y, probability = T, xlim = c( -8, 7 ), ylim = c( 0, 0.2 ), main = "", ylab = "Density", xlab = 'prior dashed, posterior solid' ) lines( density( y.0.prior ), lty = 5, xlab = "p(y_0) dashed, p(y_0|data) solid", col = 'blue' ) lines( density( y.0.post ), lty = 1, col = 'red', lwd = 2 )