rm( list = ls( ) ) # BNP analyses of the NB10 data # install DPpackage from CRAN library( DPpackage ) # unbuffer the output 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. ) ########## run 1 with PTlm nb10.sorted <- sort( y ) state.nb10.PTlm.1 <- NULL nburn.nb10.PTlm.1 <- 5000 nsave.nb10.PTlm.1 <- 10000 nskip.nb10.PTlm.1 <- 20 ndisplay.nb10.PTlm.1 <- 100 mcmc.nb10.PTlm.1 <- list( nburn = nburn.nb10.PTlm.1, nsave = nsave.nb10.PTlm.1, nskip = nskip.nb10.PTlm.1, ndisplay = ndisplay.nb10.PTlm.1 ) prior.nb10.PTlm.1 <- list( alpha = 1, beta0 = 0, Sbeta0 = diag( 1000, 1 ), tau1 = 0.01, tau2 = 0.01, M = 6 ) fit.nb10.PTlm.1 <- PTlm( formula = nb10.sorted ~ 1, prior = prior.nb10.PTlm.1, mcmc = mcmc.nb10.PTlm.1, state = state.nb10.PTlm.1, status = TRUE ) # 200,000 iterations (thinning factor 20) takes about 50 seconds summary( fit.nb10.PTlm.1 ) # Bayesian Semiparametric Regression Model # Call: # PTlm.default(formula = nb10.sorted ~ 1, prior = prior.nb10.PTlm.1, # mcmc = mcmc.nb10.PTlm.1, state = state.nb10.PTlm.1, status = TRUE) # Posterior Predictive Distributions (log): # Min. 1st Qu. Median Mean 3rd Qu. Max. # -13.030 -2.911 -2.552 -3.017 -2.392 -2.199 # Regression coefficients: # Mean Median Std. Dev. Naive Std.Error 95%HPD-Low # (Intercept) 4.042e+02 4.043e+02 5.835e-01 5.835e-03 4.030e+02 # mu 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 # sigma2 7.337e+01 7.134e+01 1.297e+01 1.297e-01 4.864e+01 # 95%HPD-Upp # (Intercept) 4.051e+02 # mu 0.000e+00 # sigma2 9.942e+01 # Acceptance Rate for Metropolis Step = 0.4921628 0 0.4420186 # Number of Observations: 100 plot( fit.nb10.PTlm.1, ask = TRUE ) colnames( fit.nb10.PTlm.1$save.state$thetasave ) # "(Intercept)" "mu" "sigma2" "alpha" theta.star.nb10.PTlm.1 <- fit.nb10.PTlm.1$save.state$thetasave[ , 1 ] sigma.2.star.nb10.PTlm.1 <- fit.nb10.PTlm.1$save.state$thetasave[ , 3 ] sigma.star.nb10.PTlm.1 <- sqrt( sigma.2.star.nb10.PTlm.1 ) par( mfrow = c( 1, 2 ) ) plot( density( theta.star.nb10.PTlm.1, adjust = 1 ), type = 'l', xlab = 'theta', main = '' ) plot( density( sigma.star.nb10.PTlm.1, adjust = 1 ), type = 'l', xlab = 'sigma', main = '' ) mean( theta.star.nb10.PTlm.1 ) # 404.1834 sd( theta.star.nb10.PTlm.1 ) # 0.5767217 print( rho.hat.theta.nb10.PTlm.1 <- acf( theta.star.nb10.PTlm.1, plot = F )$acf[ 2 ] ) # 0.7490921 print( mcmc.se.theta.nb10.PTlm.1 <- sqrt( ( 1 + rho.hat.theta.nb10.PTlm.1 ) / ( 1 - rho.hat.theta.nb10.PTlm.1 ) ) * sqrt( mean( sigma.star.nb10.PTlm.1 ) / nsave.nb10.PTlm.1 ) ) # 0.07703349 ########## run 2 with PTlm state.nb10.PTlm.2 <- NULL nburn.nb10.PTlm.2 <- 5000 nsave.nb10.PTlm.2 <- 60000 nskip.nb10.PTlm.2 <- 200 ndisplay.nb10.PTlm.2 <- 100 mcmc.nb10.PTlm.2 <- list( nburn = nburn.nb10.PTlm.2, nsave = nsave.nb10.PTlm.2, nskip = nskip.nb10.PTlm.2, ndisplay = ndisplay.nb10.PTlm.2 ) prior.nb10.PTlm.2 <- list( alpha = 1, beta0 = 0, Sbeta0 = diag( 1000, 1 ), tau1 = 0.01, tau2 = 0.01, M = 6 ) fit.nb10.PTlm.2 <- PTlm( formula = nb10.sorted ~ 1, prior = prior.nb10.PTlm.2, mcmc = mcmc.nb10.PTlm.2, state = state.nb10.PTlm.2, status = TRUE ) summary( fit.nb10.PTlm.2 ) # 12,000,000 iterations (thinning factor 200) takes about 50 minutes # Call: # PTlm.default(formula = nb10.sorted ~ 1, prior = prior.nb10.PTlm.2, # mcmc = mcmc.nb10.PTlm.2, state = state.nb10.PTlm.2, status = TRUE) # Posterior Predictive Distributions (log): # Min. 1st Qu. Median Mean 3rd Qu. Max. # -13.450 -2.909 -2.545 -3.024 -2.388 -2.195 # Regression coefficients: # Mean Median Std. Dev. Naive Std.Error 95%HPD-Low 95%HPD-Upp # (Intercept) 404.15258 404.25189 0.58051 0.00237 403.00593 405.00001 # mu 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 # sigma2 71.34025 70.48559 12.52918 0.05115 47.14353 97.60013 # Acceptance Rate for Metropolis Step = 0.4436739 0 0.4481326 # Number of Observations: 100 plot( fit.nb10.PTlm.2, ask = TRUE ) colnames( fit.nb10.PTlm.2$save.state$thetasave ) # "(Intercept)" "mu" "sigma2" "alpha" theta.star.nb10.PTlm.2 <- fit.nb10.PTlm.2$save.state$thetasave[ , 1 ] sigma.2.star.nb10.PTlm.2 <- fit.nb10.PTlm.2$save.state$thetasave[ , 3 ] sigma.star.nb10.PTlm.2 <- sqrt( sigma.2.star.nb10.PTlm.2 ) par( mfrow = c( 1, 2 ) ) plot( density( theta.star.nb10.PTlm.2, adjust = 1 ), type = 'l', xlab = 'theta', main = '' ) plot( density( sigma.star.nb10.PTlm.2, adjust = 1 ), type = 'l', xlab = 'sigma', main = '' ) mean( theta.star.nb10.PTlm.2 ) # 404.1526 sd( theta.star.nb10.PTlm.2 ) # 0.5805114 print( quantile( theta.star.nb10.PTlm.2, probs = c( 0.025, 0.975 ) ) ) print( rho.hat.theta.nb10.PTlm.2 <- acf( theta.star.nb10.PTlm.2, plot = F )$acf[ 2 ] ) # 0.05662473 print( mcmc.se.theta.nb10.PTlm.2 <- sqrt( ( 1 + rho.hat.theta.nb10.PTlm.2 ) / ( 1 - rho.hat.theta.nb10.PTlm.2 ) ) * sqrt( mean( sigma.star.nb10.PTlm.2 ) / nsave.nb10.PTlm.2 ) ) # 0.01253291 ###################### end NB10 PTlm analysis ########## run 1 with PTdensity nb10 <- 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. ) state <- NULL nburn <- 500 nsave <- 10000 nskip <- 20 ndisplay <- 500 mcmc.nb10.PTdensity.1 <- list( nburn = nburn, nsave = nsave, nskip = nskip, ndisplay = ndisplay, tune1 = 0.06 ) prior.nb10.PTdensity.1 <- list( a0 = NULL, alpha = 1, M = 6 ) fit.nb10.PTdensity.1 <- PTdensity( y = nb10, ngrid = 1000, prior = prior.nb10.PTdensity.1, mcmc = mcmc.nb10.PTdensity.1, state = state, status = TRUE ) # 200,000 iterations (thinning factor 20) took about 35 seconds summary( fit.nb10.PTdensity.1 ) Bayesian Density Estimation Using MPT # Call: # PTdensity.default(y = nb10, ngrid = 1000, prior = prior.nb10.PTdensity.1, # mcmc = mcmc.nb10.PTdensity.1, state = state, status = TRUE) # Posterior Predictive Distributions (log): # Min. 1st Qu. Median Mean 3rd Qu. Max. # -13.740 -2.933 -2.612 -3.035 -2.426 -2.207 # Baseline parameters: # Mean Median Std. Dev. Naive Std.Error 95%HPD-Low # mu:nb10 4.040e+02 4.039e+02 7.829e-01 7.829e-03 4.025e+02 # sigma:nb10 8.752e+00 8.598e+00 9.205e-01 9.205e-03 7.080e+00 # 95%HPD-Upp # mu:nb10 4.056e+02 # sigma:nb10 1.061e+01 # Acceptance Rate for Metropolis Step = 0.4514489 0.1685368 # Number of Observations: 100 # Number of Variables: 1 plot( fit.nb10.PTdensity.1, ask = TRUE ) points( nb10, rep( 0, length( nb10 ) ) ) plot( fit.nb10.PTdensity.1, ask = TRUE, output = "param" ) density.estimate.nb10.PTdensity.1$x1 <- cbind( fit.nb10.PTdensity.1$x1, fit.nb10.PTdensity.1$dens ) colnames( fit.nb10.PTdensity.1$save.state$thetasave ) [1] "mu:y" "sigma:y" "alpha" mu.star.nb10.PTdensity.1 <- fit.nb10.PTdensity.1$save.state$thetasave[ , 1 ] sigma.star.nb10.PTdensity.1 <- fit.nb10.PTdensity.1$save.state$thetasave[ , 2 ] alpha.star.nb10.PTdensity.1 <- fit.nb10.PTdensity.1$save.state$thetasave[ , 3 ] mean( mu.star.nb10.PTdensity.1 ) # 403.969 sd( mu.star.nb10.PTdensity.1 ) # 0.7829 rho.hat <- acf( mu.star.nb10.PTdensity.1, plot = F )$acf[ 2 ] print( mcmc.se <- sqrt( ( 1 + rho.hat ) / ( 1 - rho.hat ) ) * mean( sigma.star.nb10.PTdensity.1 ) / sqrt( nsave ) ) # 1.742979 ################ run 2 # not fully edited below this point nsave <- 5000 mcmc.nb10.PTdensity.2 <- list( nburn = 500, nsave = nsave, nskip = 500, ndisplay = 100, tune1 = 0.06, tune2 = 0.37 ) prior.nb10.PTdensity.2 <- list( a0 = 1, b0 = 1, M = 6 ) fit.nb10.PTdensity.2 <- PTdensity( y = nb10, ngrid = 100, prior = prior.nb10.PTdensity.2, mcmc = mcmc.nb10.PTdensity.2, state = NULL, status = TRUE ) fit.nb10.PTdensity.2 Bayesian Density Estimation Using MPT (7 min) Call: PTdensity.default(y = nb10, ngrid = 100, prior = prior.nb10.PTdensity.2, mcmc = mcmc.nb10.PTdensity.2, state = NULL, status = TRUE) Posterior Predictive Distributions (log): Min. 1st Qu. Median Mean 3rd Qu. Max. -16.8000 -2.2150 -1.8020 -2.3680 -1.2660 -0.9151 Posterior Inference of Parameters: mu:nb10 sigma:nb10 alpha 405.12102 6.23550 0.02075 Acceptance Rate for Metropolis Step = 0.3482985 0.2048569 0.5701457 Number of Observations: 100 Number of Variables: 1 plot( fit.nb10.PTdensity.2, ask = TRUE ) points( nb10, rep( 0, length( nb10 ) ) ) plot( fit.nb10.PTdensity.2, ask = TRUE, output = "param" ) density.estimate.nb10.PTdensity.2 <- cbind( fit.nb10.PTdensity.2$x1, fit.nb10.PTdensity.2$dens ) colnames( fit.nb10.PTdensity.2$save.state$thetasave ) # [1] "mu:nb10" "sigma:nb10" "alpha" mu.star.nb10.PTdensity.2 <- fit.nb10.PTdensity.2$save.state$thetasave[ , 1 ] sigma.star.nb10.PTdensity.2 <- fit.nb10.PTdensity.2$save.state$thetasave[ , 2 ] alpha.star.nb10.PTdensity.2 <- fit.nb10.PTdensity.2$save.state$thetasave[ , 3 ] mean( mu.star.nb10.PTdensity.2 ) # 405.121 sd( mu.star.nb10.PTdensity.2 ) # 0.1402037 print( rho.hat.nb10.PTdensity.2 <- acf( mu.star.nb10.PTdensity.2, plot = F )$acf[ 2 ] ) # 0.9503532 print( mcmc.se.nb10.PTdensity.2 <- sqrt( ( 1 + rho.hat.nb10.PTdensity.2 ) / ( 1 - rho.hat.nb10.PTdensity.2 ) ) * sqrt( mean( sigma.star.nb10.PTdensity.2 ) / nsave ) ) # 0.2213409 ################ run 3 nsave <- 15000 mcmc.nb10.PTdensity.3 <- list( nburn = 500, nsave = nsave, nskip = 5000, ndisplay = 100, tune1 = 0.06, tune2 = 0.2, tune3 = 0.9 ) prior.nb10.PTdensity.3 <- list( a0 = 1, b0 = 1, M = 6 ) fit.nb10.PTdensity.3 <- PTdensity( y = nb10, ngrid = 100, prior = prior.nb10.PTdensity.3, mcmc = mcmc.nb10.PTdensity.3, state = NULL, status = TRUE ) fit.nb10.PTdensity.3 Bayesian Density Estimation Using MPT (75 million iterations in 3.5 hours) Call: PTdensity.default(y = nb10, ngrid = 100, prior = prior.nb10.PTdensity.3, mcmc = mcmc.nb10.PTdensity.3, state = NULL, status = TRUE) Posterior Predictive Distributions (log): Min. 1st Qu. Median Mean 3rd Qu. Max. -19.5600 -2.1930 -1.7570 -2.3790 -1.2180 -0.8817 Posterior Inference of Parameters: mu:nb10 sigma:nb10 alpha 404.79958 5.98549 0.01919 Acceptance Rate for Metropolis Step = 0.3200288 0.3195898 0.6236641 Number of Observations: 100 Number of Variables: 1 plot( fit.nb10.PTdensity.3, ask = TRUE ) points( nb10, rep( 0, length( nb10 ) ) ) plot( fit.nb10.PTdensity.3, ask = TRUE, output = "param" ) density.estimate.nb10.PTdensity.3 <- cbind( fit.nb10.PTdensity.3$x1, fit.nb10.PTdensity.3$dens ) colnames( fit.nb10.PTdensity.3$save.state$thetasave ) # [1] "mu:nb10" "sigma:nb10" "alpha" mu.star.nb10.PTdensity.3 <- fit.nb10.PTdensity.3$save.state$thetasave[ , 1 ] sigma.star.nb10.PTdensity.3 <- fit.nb10.PTdensity.3$save.state$thetasave[ , 2 ] alpha.star.nb10.PTdensity.3 <- fit.nb10.PTdensity.3$save.state$thetasave[ , 3 ] mean( mu.star.nb10.PTdensity.3 ) # 404.7996 sd( mu.star.nb10.PTdensity.3 ) # 0.2587729 print( rho.hat.nb10.PTdensity.3 <- acf( mu.star.nb10.PTdensity.3, plot = F )$acf[ 2 ] ) # 0.9824521 print( mcmc.se.nb10.PTdensity.3 <- sqrt( ( 1 + rho.hat.nb10.PTdensity.3 ) / ( 1 - rho.hat.nb10.PTdensity.3 ) ) * sqrt( mean( sigma.star.nb10.PTdensity.3 ) / nsave ) ) # 0.2123212 ################ run 4 nsave <- 1000 mcmc.nb10.PTdensity.4 <- list( nburn = 500, nsave = nsave, nskip = 500, ndisplay = 100, tune1 = 0.05, tune2 = 0.16, tune3 = 1.3 ) prior.nb10.PTdensity.4 <- list( a0 = NULL, alpha = 10, M = 6, m0 = 404.6, S0 = 0.43, tau1 = 98.4, tau2 = 4159.2 ) fit.nb10.PTdensity.4 <- PTdensity( y = nb10, ngrid = 100, prior = prior.nb10.PTdensity.4, mcmc = mcmc.nb10.PTdensity.4, state = NULL, status = TRUE ) fit.nb10.PTdensity.4 Bayesian Density Estimation Using MPT (500,000 iterations in 1 minute) Call: PTdensity.default(y = nb10, ngrid = 100, prior = prior.nb10.PTdensity.4, mcmc = mcmc.nb10.PTdensity.4, state = NULL, status = TRUE) Posterior Predictive Distributions (log): Min. 1st Qu. Median Mean 3rd Qu. Max. -15.620 -3.237 -2.731 -3.235 -2.680 -2.620 Posterior Inference of Parameters: mu:nb10 sigma:nb10 alpha 404.610 6.762 10.000 Acceptance Rate for Metropolis Step = 0.7447617 0.796991 Number of Observations: 100 Number of Variables: 1 plot( fit.nb10.PTdensity.4, ask = TRUE ) points( nb10, rep( 0, length( nb10 ) ) ) plot( fit.nb10.PTdensity.4, ask = TRUE, output = "param" ) density.estimate.nb10.PTdensity.4 <- cbind( fit.nb10.PTdensity.4$x1, fit.nb10.PTdensity.4$dens ) colnames( fit.nb10.PTdensity.4$save.state$thetasave ) # [1] "mu:nb10" "sigma:nb10" "alpha" mu.star.nb10.PTdensity.4 <- fit.nb10.PTdensity.4$save.state$thetasave[ , 1 ] sigma.star.nb10.PTdensity.4 <- fit.nb10.PTdensity.4$save.state$thetasave[ , 2 ] alpha.star.nb10.PTdensity.4 <- fit.nb10.PTdensity.4$save.state$thetasave[ , 3 ] mean( mu.star.nb10.PTdensity.4 ) # 404.6103 sd( mu.star.nb10.PTdensity.4 ) # 0.6717502 print( rho.hat.nb10.PTdensity.4 <- acf( mu.star.nb10.PTdensity.4, plot = F )$acf[ 2 ] ) # 0.766459 print( mcmc.se.nb10.PTdensity.4 <- sqrt( ( 1 + rho.hat.nb10.PTdensity.4 ) / ( 1 - rho.hat.nb10.PTdensity.4 ) ) * sqrt( mean( sigma.star.nb10.PTdensity.4 ) / nsave ) ) # 0.2261615 ################ run 5 nb10 <- sort( nb10 ) nb10 # [1] 375 392 393 397 398 398 399 399 399 399 399 399 399 400 400 400 400 401 401 401 # [21] 401 401 401 401 401 401 401 401 401 402 402 402 402 402 402 402 402 403 403 403 # [41] 403 403 403 404 404 404 404 404 404 404 404 404 405 405 405 405 405 406 406 406 # [61] 406 406 406 406 406 406 406 406 406 407 407 407 407 407 407 407 407 408 408 408 # [81] 408 408 409 409 409 409 409 410 410 410 410 411 412 412 412 413 415 418 423 437 nsave <- 1000 mcmc.nb10.PTdensity.5 <- list( nburn = 500, nsave = nsave, nskip = 500, ndisplay = 100, tune1 = 0.1, tune2 = 0.3 ) prior.nb10.PTdensity.5 <- list( a0 = NULL, alpha = 5, M = 6, m0 = 404.6, S0 = 0.43, tau1 = 98.4, tau2 = 4159.2 ) fit.nb10.PTdensity.5 <- PTdensity( y = nb10, ngrid = 100, prior = prior.nb10.PTdensity.5, mcmc = mcmc.nb10.PTdensity.5, state = NULL, status = TRUE ) fit.nb10.PTdensity.5 Bayesian Density Estimation Using MPT (500,000 iterations in 1 minute) Call: PTdensity.default(y = nb10, ngrid = 100, prior = prior.nb10.PTdensity.5, mcmc = mcmc.nb10.PTdensity.5, state = NULL, status = TRUE) Posterior Predictive Distributions (log): Min. 1st Qu. Median Mean 3rd Qu. Max. -15.220 -3.230 -2.650 -3.193 -2.617 -2.554 Posterior Inference of Parameters: mu:nb10 sigma:nb10 alpha 404.656 6.931 5.000 Acceptance Rate for Metropolis Step = 0.5080758 0.6038644 Number of Observations: 100 Number of Variables: 1 plot( fit.nb10.PTdensity.5, ask = TRUE ) points( nb10, rep( 0, length( nb10 ) ) ) plot( fit.nb10.PTdensity.5, ask = TRUE, output = "param" ) density.estimate.nb10.PTdensity.5 <- cbind( fit.nb10.PTdensity.5$x1, fit.nb10.PTdensity.5$dens ) colnames( fit.nb10.PTdensity.5$save.state$thetasave ) # [1] "mu:nb10" "sigma:nb10" "alpha" mu.star.nb10.PTdensity.5 <- fit.nb10.PTdensity.5$save.state$thetasave[ , 1 ] sigma.star.nb10.PTdensity.5 <- fit.nb10.PTdensity.5$save.state$thetasave[ , 2 ] alpha.star.nb10.PTdensity.5 <- fit.nb10.PTdensity.5$save.state$thetasave[ , 3 ] mean( mu.star.nb10.PTdensity.5 ) # 404.6556 sd( mu.star.nb10.PTdensity.5 ) # 0.7420994 print( rho.hat.nb10.PTdensity.5 <- acf( mu.star.nb10.PTdensity.5, plot = F )$acf[ 2 ] ) # 0.756568 print( mcmc.se.nb10.PTdensity.5 <- sqrt( ( 1 + rho.hat.nb10.PTdensity.5 ) / ( 1 - rho.hat.nb10.PTdensity.5 ) ) * sqrt( mean( sigma.star.nb10.PTdensity.5 ) / nsave ) ) # 0.2236329 fit.nb10.PTdensity.5$cpo [1] 1.650863e-06 7.069208e-03 9.084301e-03 2.098114e-02 2.654275e-02 2.654275e-02 [7] 3.953982e-02 3.953982e-02 3.953982e-02 3.953982e-02 3.953982e-02 3.953982e-02 [13] 3.953982e-02 4.866045e-02 4.866045e-02 4.866045e-02 4.866045e-02 7.299144e-02 [19] 7.299144e-02 7.299144e-02 7.299144e-02 7.299144e-02 7.299144e-02 7.299144e-02 [25] 7.299144e-02 7.299144e-02 7.299144e-02 7.299144e-02 7.299144e-02 7.185608e-02 [31] 7.185608e-02 7.185608e-02 7.185608e-02 7.185608e-02 7.185608e-02 7.185608e-02 [37] 7.185608e-02 7.063841e-02 7.063841e-02 7.063841e-02 7.063841e-02 7.063841e-02 [43] 7.063841e-02 7.780352e-02 7.780352e-02 7.780352e-02 7.780352e-02 7.780352e-02 [49] 7.780352e-02 7.780352e-02 7.780352e-02 7.780352e-02 7.263345e-02 7.263345e-02 [55] 7.263345e-02 7.263345e-02 7.263345e-02 7.745945e-02 7.745945e-02 7.745945e-02 [61] 7.745945e-02 7.745945e-02 7.745945e-02 7.745945e-02 7.745945e-02 7.745945e-02 [67] 7.745945e-02 7.745945e-02 7.745945e-02 6.574401e-02 6.574401e-02 6.574401e-02 [73] 6.574401e-02 6.574401e-02 6.574401e-02 6.574401e-02 6.574401e-02 5.658003e-02 [79] 5.658003e-02 5.658003e-02 5.658003e-02 5.658003e-02 4.686425e-02 4.686425e-02 [85] 4.686425e-02 4.686425e-02 4.686425e-02 3.455696e-02 3.455696e-02 3.455696e-02 [91] 3.455696e-02 2.695459e-02 2.353519e-02 2.353519e-02 2.353519e-02 1.868259e-02 [97] 1.242316e-02 5.923828e-03 1.056741e-03 2.455426e-07 hist( nb10, nclass = 100, probability = T, main = '' ) lines( nb10, fit.nb10.PTdensity.5$cpo, lty = 2 ) hist( nb10, nclass = 50, probability = T, main = '' ) lines( nb10, fit.nb10.PTdensity.5$cpo, lty = 2 ) hist( nb10, nclass = 25, probability = T, main = '' ) lines( nb10, fit.nb10.PTdensity.5$cpo, lty = 2 ) colnames( fit.nb10.PTdensity.5$save.state$randsave ) # doesn't work dim( fit.nb10.PTdensity.5$save.state$randsave ) ################ run 6 nsave <- 1000 mcmc.nb10.PTdensity.6 <- list( nburn = 500, nsave = nsave, nskip = 500, ndisplay = 100, tune1 = 0.125, tune2 = 0.4 ) prior.nb10.PTdensity.6 <- list( a0 = NULL, alpha = 2.5, M = 6, m0 = 404.6, S0 = 0.43, tau1 = 98.4, tau2 = 4159.2 ) fit.nb10.PTdensity.6 <- PTdensity( y = nb10, ngrid = 100, prior = prior.nb10.PTdensity.6, mcmc = mcmc.nb10.PTdensity.6, state = NULL, status = TRUE ) fit.nb10.PTdensity.6 Bayesian Density Estimation Using MPT (500,000 iterations in 1.2 minutes) Call: PTdensity.default(y = nb10, ngrid = 100, prior = prior.nb10.PTdensity.6, mcmc = mcmc.nb10.PTdensity.6, state = NULL, status = TRUE) Posterior Predictive Distributions (log): Min. 1st Qu. Median Mean 3rd Qu. Max. -14.920 -3.097 -2.632 -3.132 -2.484 -2.399 Posterior Inference of Parameters: mu:nb10 sigma:nb10 alpha 404.358 7.164 2.500 Acceptance Rate for Metropolis Step = 0.3559342 0.4449212 Number of Observations: 100 Number of Variables: 1 plot( fit.nb10.PTdensity.6, ask = TRUE ) points( nb10, rep( 0, length( nb10 ) ) ) plot( fit.nb10.PTdensity.6, ask = TRUE, output = "param" ) density.estimate.nb10.PTdensity.6 <- cbind( fit.nb10.PTdensity.6$x1, fit.nb10.PTdensity.6$dens ) colnames( fit.nb10.PTdensity.6$save.state$thetasave ) # [1] "mu:nb10" "sigma:nb10" "alpha" mu.star.nb10.PTdensity.6 <- fit.nb10.PTdensity.6$save.state$thetasave[ , 1 ] sigma.star.nb10.PTdensity.6 <- fit.nb10.PTdensity.6$save.state$thetasave[ , 2 ] alpha.star.nb10.PTdensity.6 <- fit.nb10.PTdensity.6$save.state$thetasave[ , 3 ] mean( mu.star.nb10.PTdensity.6 ) # 404.3582 sd( mu.star.nb10.PTdensity.6 ) # 0.7076105 print( rho.hat.nb10.PTdensity.6 <- acf( mu.star.nb10.PTdensity.6, plot = F )$acf[ 2 ] ) # 0.8230109 print( mcmc.se.nb10.PTdensity.6 <- sqrt( ( 1 + rho.hat.nb10.PTdensity.6 ) / ( 1 - rho.hat.nb10.PTdensity.6 ) ) * sqrt( mean( sigma.star.nb10.PTdensity.6 ) / nsave ) ) # 0.2716344 hist( nb10, nclass = 100, probability = T, main = '' ) lines( nb10, fit.nb10.PTdensity.6$cpo, lty = 2 ) hist( nb10, nclass = 50, probability = T, main = '' ) lines( nb10, fit.nb10.PTdensity.6$cpo, lty = 2 ) hist( nb10, nclass = 25, probability = T, main = '' ) lines( nb10, fit.nb10.PTdensity.6$cpo, lty = 2 ) colnames( fit.nb10.PTdensity.6$save.state$randsave ) # doesn't work dim( fit.nb10.PTdensity.6$save.state$randsave ) # doesn't work ################ run 7 nsave <- 1000 mcmc.nb10.PTdensity.7 <- list( nburn = 500, nsave = nsave, nskip = 500, ndisplay = 100, tune1 = 0.11, tune2 = 0.4 ) prior.nb10.PTdensity.7 <- list( a0 = NULL, alpha = 1, M = 6, m0 = 404.6, S0 = 0.43, tau1 = 98.4, tau2 = 4159.2 ) fit.nb10.PTdensity.7 <- PTdensity( y = nb10, ngrid = 100, prior = prior.nb10.PTdensity.7, mcmc = mcmc.nb10.PTdensity.7, state = NULL, status = TRUE ) fit.nb10.PTdensity.7 Bayesian Density Estimation Using MPT (500,000 iterations in 1 minute) Call: PTdensity.default(y = nb10, ngrid = 100, prior = prior.nb10.PTdensity.7, mcmc = mcmc.nb10.PTdensity.7, state = NULL, status = TRUE) Posterior Predictive Distributions (log): Min. 1st Qu. Median Mean 3rd Qu. Max. -15.150 -2.904 -2.620 -3.047 -2.369 -2.115 Posterior Inference of Parameters: mu:nb10 sigma:nb10 alpha 404.161 7.302 1.000 Acceptance Rate for Metropolis Step = 0.2979123 0.3507119 Number of Observations: 100 Number of Variables: 1 plot( fit.nb10.PTdensity.7, ask = TRUE ) points( nb10, rep( 0, length( nb10 ) ) ) plot( fit.nb10.PTdensity.7, ask = TRUE, output = "param" ) density.estimate.nb10.PTdensity.7 <- cbind( fit.nb10.PTdensity.7$x1, fit.nb10.PTdensity.7$dens ) colnames( fit.nb10.PTdensity.7$save.state$thetasave ) # [1] "mu:nb10" "sigma:nb10" "alpha" mu.star.nb10.PTdensity.7 <- fit.nb10.PTdensity.7$save.state$thetasave[ , 1 ] sigma.star.nb10.PTdensity.7 <- fit.nb10.PTdensity.7$save.state$thetasave[ , 2 ] alpha.star.nb10.PTdensity.7 <- fit.nb10.PTdensity.7$save.state$thetasave[ , 3 ] mean( mu.star.nb10.PTdensity.7 ) # 404.1607 sd( mu.star.nb10.PTdensity.7 ) # 0.611553 print( rho.hat.nb10.PTdensity.7 <- acf( mu.star.nb10.PTdensity.7, plot = F )$acf[ 2 ] ) # 0.9296005 print( mcmc.se.nb10.PTdensity.7 <- sqrt( ( 1 + rho.hat.nb10.PTdensity.7 ) / ( 1 - rho.hat.nb10.PTdensity.7 ) ) * sqrt( mean( sigma.star.nb10.PTdensity.7 ) / nsave ) ) # 0.4473693 hist( nb10, nclass = 100, probability = T, main = '' ) lines( nb10, fit.nb10.PTdensity.7$cpo, lty = 2 ) hist( nb10, nclass = 50, probability = T, main = '' ) lines( nb10, fit.nb10.PTdensity.7$cpo, lty = 2 ) hist( nb10, nclass = 25, probability = T, main = '' ) lines( nb10, fit.nb10.PTdensity.7$cpo, lty = 2 ) colnames( fit.nb10.PTdensity.7$save.state$randsave ) # doesn't work dim( fit.nb10.PTdensity.7$save.state$randsave ) # doesn't work ################ run 8 nsave <- 1000 mcmc.nb10.PTdensity.8 <- list( nburn = 500, nsave = nsave, nskip = 500, ndisplay = 100, tune1 = 0.08, tune2 = 0.35 ) prior.nb10.PTdensity.8 <- list( a0 = NULL, alpha = 0.5, M = 6, m0 = 404.6, S0 = 0.43, tau1 = 98.4, tau2 = 4159.2 ) fit.nb10.PTdensity.8 <- PTdensity( y = nb10, ngrid = 100, prior = prior.nb10.PTdensity.8, mcmc = mcmc.nb10.PTdensity.8, state = NULL, status = TRUE ) fit.nb10.PTdensity.8 Bayesian Density Estimation Using MPT (500,000 iterations in 1 minute) Call: PTdensity.default(y = nb10, ngrid = 100, prior = prior.nb10.PTdensity.8, mcmc = mcmc.nb10.PTdensity.8, state = NULL, status = TRUE) Posterior Predictive Distributions (log): Min. 1st Qu. Median Mean 3rd Qu. Max. -14.510 -2.936 -2.545 -2.949 -2.243 -1.880 Posterior Inference of Parameters: mu:nb10 sigma:nb10 alpha 404.667 7.386 0.500 Acceptance Rate for Metropolis Step = 0.3522433 0.3430449 Number of Observations: 100 Number of Variables: 1 plot( fit.nb10.PTdensity.8, ask = TRUE ) points( nb10, rep( 0, length( nb10 ) ) ) plot( fit.nb10.PTdensity.8, ask = TRUE, output = "param" ) density.estimate.nb10.PTdensity.8 <- cbind( fit.nb10.PTdensity.8$x1, fit.nb10.PTdensity.8$dens ) colnames( fit.nb10.PTdensity.8$save.state$thetasave ) # [1] "mu:nb10" "sigma:nb10" "alpha" mu.star.nb10.PTdensity.8 <- fit.nb10.PTdensity.8$save.state$thetasave[ , 1 ] sigma.star.nb10.PTdensity.8 <- fit.nb10.PTdensity.8$save.state$thetasave[ , 2 ] alpha.star.nb10.PTdensity.8 <- fit.nb10.PTdensity.8$save.state$thetasave[ , 3 ] mean( mu.star.nb10.PTdensity.8 ) # sd( mu.star.nb10.PTdensity.8 ) # print( rho.hat.nb10.PTdensity.8 <- acf( mu.star.nb10.PTdensity.8, plot = F )$acf[ 2 ] ) # print( mcmc.se.nb10.PTdensity.8 <- sqrt( ( 1 + rho.hat.nb10.PTdensity.8 ) / ( 1 - rho.hat.nb10.PTdensity.8 ) ) * sqrt( mean( sigma.star.nb10.PTdensity.8 ) / nsave ) ) # hist( nb10, nclass = 100, probability = T, main = '' ) lines( nb10, fit.nb10.PTdensity.8$cpo, lty = 2 ) hist( nb10, nclass = 50, probability = T, main = '' ) lines( nb10, fit.nb10.PTdensity.8$cpo, lty = 2 ) hist( nb10, nclass = 25, probability = T, main = '' ) lines( nb10, fit.nb10.PTdensity.8$cpo, lty = 2 ) colnames( fit.nb10.PTdensity.8$save.state$randsave ) # doesn't work dim( fit.nb10.PTdensity.8$save.state$randsave ) # doesn't work ################ run 9 nsave <- 10000 mcmc.nb10.PTdensity.9 <- list( nburn = 500, nsave = nsave, nskip = 5000, ndisplay = 100, tune1 = 0.08, tune2 = 0.35 ) prior.nb10.PTdensity.9 <- list( a0 = NULL, alpha = 1, M = 6, m0 = 404.6, S0 = 0.43, tau1 = 98.4, tau2 = 4159.2 ) fit.nb10.PTdensity.9 <- PTdensity( y = nb10, ngrid = 100, prior = prior.nb10.PTdensity.9, mcmc = mcmc.nb10.PTdensity.9, state = NULL, status = TRUE ) fit.nb10.PTdensity.9 Bayesian Density Estimation Using MPT (5,000,000 iterations in 1.7 hours) Call: PTdensity.default(y = nb10, ngrid = 100, prior = prior.nb10.PTdensity.9, mcmc = mcmc.nb10.PTdensity.9, state = NULL, status = TRUE) Posterior Predictive Distributions (log): Min. 1st Qu. Median Mean 3rd Qu. Max. -14.790 -2.932 -2.613 -3.046 -2.373 -2.138 Posterior Inference of Parameters: mu:nb10 sigma:nb10 alpha 404.285 7.346 1.000 Acceptance Rate for Metropolis Step = 0.3838885 0.3873463 Number of Observations: 100 Number of Variables: 1 plot( fit.nb10.PTdensity.9 ) points( nb10, rep( 0, length( nb10 ) ) ) plot( fit.nb10.PTdensity.9, ask = TRUE, output = "param" ) density.estimate.nb10.PTdensity.9 <- cbind( fit.nb10.PTdensity.9$x1, fit.nb10.PTdensity.9$dens ) hist( nb10, nclass = 25, probability = T, main = '' ) lines( density.estimate.nb10.PTdensity.9[ , 1 ], density.estimate.nb10.PTdensity.9[ , 2 ], lty = 2 ) colnames( fit.nb10.PTdensity.9$save.state$thetasave ) # [1] "mu:nb10" "sigma:nb10" "alpha" mu.star.nb10.PTdensity.9 <- fit.nb10.PTdensity.9$save.state$thetasave[ , 1 ] sigma.star.nb10.PTdensity.9 <- fit.nb10.PTdensity.9$save.state$thetasave[ , 2 ] alpha.star.nb10.PTdensity.9 <- fit.nb10.PTdensity.9$save.state$thetasave[ , 3 ] mean( mu.star.nb10.PTdensity.9 ) # 404.2855 sd( mu.star.nb10.PTdensity.9 ) # 0.7107078 print( rho.hat.nb10.PTdensity.9 <- acf( mu.star.nb10.PTdensity.9, plot = F )$acf[ 2 ] ) # 0.7100531 print( mcmc.se.nb10.PTdensity.9 <- sqrt( ( 1 + rho.hat.nb10.PTdensity.9 ) / ( 1 - rho.hat.nb10.PTdensity.9 ) ) * sqrt( mean( sigma.star.nb10.PTdensity.9 ) / nsave ) ) # 0.06582072 plot( density( mu.star.nb10.PTdensity.9 ), type = 'l', main = '', xlab = 'mu' ) plot( density( sigma.star.nb10.PTdensity.9 ), type = 'l', main = '', xlab = 'sigma' ) hist( nb10, nclass = 100, probability = T, main = '' ) lines( nb10, fit.nb10.PTdensity.9$cpo, lty = 2 ) hist( nb10, nclass = 50, probability = T, main = '' ) lines( nb10, fit.nb10.PTdensity.9$cpo, lty = 2 ) hist( nb10, nclass = 25, probability = T, main = '' ) lines( nb10, fit.nb10.PTdensity.9$cpo, lty = 2 ) write( t( density.estimate.nb10.PTdensity.9 ), 'density.estimate.nb10.PTdensity.9.txt', ncol = 2 ) write( t( mu.star.nb10.PTdensity.9 ), 'mu.star.nb10.PTdensity.9.txt', ncol = 1 ) write( t( sigma.star.nb10.PTdensity.9 ), 'sigma.star.nb10.PTdensity.9.txt', ncol = 1 ) ################ run 10 nsave <- 50000 mcmc.nb10.PTdensity.10 <- list( nburn = 500, nsave = nsave, nskip = 5000, ndisplay = 100, tune1 = 0.08, tune2 = 0.35 ) prior.nb10.PTdensity.10 <- list( a0 = NULL, alpha = 1, M = 6, m0 = 404.6, S0 = 0.43, tau1 = 98.4, tau2 = 4159.2 ) fit.nb10.PTdensity.10 <- PTdensity( y = nb10, ngrid = 100, prior = prior.nb10.PTdensity.10, mcmc = mcmc.nb10.PTdensity.10, state = NULL, status = TRUE ) fit.nb10.PTdensity.10 Bayesian Density Estimation Using MPT (250,000,000 iterations in 8.5 hours) Call: PTdensity.default(y = nb10, ngrid = 100, prior = prior.nb10.PTdensity.10, mcmc = mcmc.nb10.PTdensity.10, state = NULL, status = TRUE) Posterior Predictive Distributions (log): Min. 1st Qu. Median Mean 3rd Qu. Max. -14.750 -2.931 -2.605 -3.045 -2.376 -2.135 Posterior Inference of Parameters: mu:nb10 sigma:nb10 alpha 404.302 7.343 1.000 Acceptance Rate for Metropolis Step = 0.3853494 0.388329 Number of Observations: 100 Number of Variables: 1 plot( fit.nb10.PTdensity.10 ) points( nb10, rep( 0, length( nb10 ) ) ) plot( fit.nb10.PTdensity.10, ask = TRUE, output = "param" ) density.estimate.nb10.PTdensity.10 <- cbind( fit.nb10.PTdensity.10$x1, fit.nb10.PTdensity.10$dens ) hist( nb10, nclass = 25, probability = T, main = '' ) lines( density.estimate.nb10.PTdensity.10[ , 1 ], density.estimate.nb10.PTdensity.10[ , 2 ], lty = 2 ) colnames( fit.nb10.PTdensity.10$save.state$thetasave ) # [1] "mu:nb10" "sigma:nb10" "alpha" mu.star.nb10.PTdensity.10 <- fit.nb10.PTdensity.10$save.state$thetasave[ , 1 ] sigma.star.nb10.PTdensity.10 <- fit.nb10.PTdensity.10$save.state$thetasave[ , 2 ] alpha.star.nb10.PTdensity.10 <- fit.nb10.PTdensity.10$save.state$thetasave[ , 3 ] mean( mu.star.nb10.PTdensity.10 ) # 404.302 sd( mu.star.nb10.PTdensity.10 ) # 0.7225608 print( rho.hat.nb10.PTdensity.10 <- acf( mu.star.nb10.PTdensity.10, plot = F )$acf[ 2 ] ) # 0.7122207 print( mcmc.se.nb10.PTdensity.10 <- sqrt( ( 1 + rho.hat.nb10.PTdensity.10 ) / ( 1 - rho.hat.nb10.PTdensity.10 ) ) * sqrt( mean( sigma.star.nb10.PTdensity.10 ) / nsave ) ) # 0.02956072 plot( density( mu.star.nb10.PTdensity.10 ), type = 'l', main = '', xlab = 'mu' ) plot( density( sigma.star.nb10.PTdensity.10 ), type = 'l', main = '', xlab = 'sigma' ) hist( nb10, nclass = 100, probability = T, main = '' ) lines( nb10, fit.nb10.PTdensity.10$cpo, lty = 2 ) hist( nb10, nclass = 50, probability = T, main = '' ) lines( nb10, fit.nb10.PTdensity.10$cpo, lty = 2 ) hist( nb10, nclass = 25, probability = T, main = '' ) lines( nb10, fit.nb10.PTdensity.10$cpo, lty = 2 ) write( t( density.estimate.nb10.PTdensity.10 ), 'density.estimate.nb10.PTdensity.10.txt', ncol = 2 ) write( t( mu.star.nb10.PTdensity.10 ), 'mu.star.nb10.PTdensity.10.txt', ncol = 1 ) write( t( sigma.star.nb10.PTdensity.10 ), 'sigma.star.nb10.PTdensity.10.txt', ncol = 1 ) ###################### end NB10 PTdensity analysis ###################### begin simulated NB10 PTdensity analysis n.nb10 <- length( nb10 ) nb10.simulated <- mean( nb10 ) + sd( nb10 ) * qnorm( ( 1:n.nb10 - 0.5 ) / n.nb10 ) ################ run 1 nsave <- 1000 mcmc.nb10.simulated.PTdensity.1 <- list( nburn = 500, nsave = nsave, nskip = 500, ndisplay = 100, tune1 = 0.08, tune2 = 0.35 ) prior.nb10.simulated.PTdensity.1 <- list( a0 = NULL, alpha = 1, M = 6, m0 = 404.6, S0 = 0.43, tau1 = 98.4, tau2 = 4159.2 ) fit.nb10.simulated.PTdensity.1 <- PTdensity( y = nb10.simulated, ngrid = 100, prior = prior.nb10.simulated.PTdensity.1, mcmc = mcmc.nb10.simulated.PTdensity.1, state = NULL, status = TRUE ) fit.nb10.simulated.PTdensity.1 Bayesian Density Estimation Using MPT (500,000 iterations in 5.6 minutes) Call: PTdensity.default(y = nb10.simulated, ngrid = 100, prior = prior.nb10.simulated.PTdensity.1, mcmc = mcmc.nb10.simulated.PTdensity.1, state = NULL, status = TRUE) Posterior Predictive Distributions (log): Min. 1st Qu. Median Mean 3rd Qu. Max. -6.238 -3.540 -3.102 -3.383 -2.940 -2.881 Posterior Inference of Parameters: mu:nb10.simulated sigma:nb10.simulated alpha 404.562 6.541 1.000 Acceptance Rate for Metropolis Step = 0.6692463 0.6408435 Number of Observations: 100 Number of Variables: 1 ################ run 2 mcmc.nb10.simulated.PTdensity.2 <- list( nburn = 500, nsave = 1000, nskip = 500, ndisplay = 100, tune1 = 0.16, tune2 = 0.70 ) prior.nb10.simulated.PTdensity.2 <- list( a0 = NULL, alpha = 1, M = 6, m0 = 404.6, S0 = 0.43, tau1 = 98.4, tau2 = 4159.2 ) fit.nb10.simulated.PTdensity.2 <- PTdensity( y = nb10.simulated, ngrid = 100, prior = prior.nb10.simulated.PTdensity.2, mcmc = mcmc.nb10.simulated.PTdensity.2, state = NULL, status = TRUE ) fit.nb10.simulated.PTdensity.2 Bayesian Density Estimation Using MPT (500,000 iterations in 5.6 minutes) Call: PTdensity.default(y = nb10.simulated, ngrid = 100, prior = prior.nb10.simulated.PTdensity.2, mcmc = mcmc.nb10.simulated.PTdensity.2, state = NULL, status = TRUE) Posterior Predictive Distributions (log): Min. 1st Qu. Median Mean 3rd Qu. Max. -6.239 -3.539 -3.098 -3.384 -2.943 -2.859 Posterior Inference of Parameters: mu:nb10.simulated sigma:nb10.simulated alpha 404.599 6.531 1.000 Acceptance Rate for Metropolis Step = 0.6132124 0.5617188 Number of Observations: 100 Number of Variables: 1 ################ run 3 nsave <- 1000 mcmc.nb10.simulated.PTdensity.3 <- list( nburn = 500, nsave = nsave, nskip = 500, ndisplay = 100, tune1 = 1.0, tune2 = 1.2 ) prior.nb10.simulated.PTdensity.3 <- list( a0 = NULL, alpha = 1, M = 6, m0 = 404.6, S0 = 0.43, tau1 = 98.4, tau2 = 4159.2 ) fit.nb10.simulated.PTdensity.3 <- PTdensity( y = nb10.simulated, ngrid = 100, prior = prior.nb10.simulated.PTdensity.3, mcmc = mcmc.nb10.simulated.PTdensity.3, state = NULL, status = TRUE ) fit.nb10.simulated.PTdensity.3 Bayesian Density Estimation Using MPT (500,000 iterations in 5.6 minutes) Call: PTdensity.default(y = nb10.simulated, ngrid = 100, prior = prior.nb10.simulated.PTdensity.3, mcmc = mcmc.nb10.simulated.PTdensity.3, state = NULL, status = TRUE) Posterior Predictive Distributions (log): Min. 1st Qu. Median Mean 3rd Qu. Max. -6.260 -3.554 -3.100 -3.381 -2.934 -2.859 Posterior Inference of Parameters: mu:nb10.simulated sigma:nb10.simulated alpha 404.588 6.501 1.000 Acceptance Rate for Metropolis Step = 0.4317906 0.4787677 Number of Observations: 100 Number of Variables: 1 plot( fit.nb10.simulated.PTdensity.3 ) points( nb10, rep( 0, length( nb10 ) ) ) plot( fit.nb10.simulated.PTdensity.3, ask = TRUE, output = "param" ) density.estimate.nb10.simulated.PTdensity.3 <- cbind( fit.nb10.simulated.PTdensity.3$x1, fit.nb10.simulated.PTdensity.3$dens ) hist( nb10, nclass = 25, probability = T, main = '' ) lines( density.estimate.nb10.simulated.PTdensity.3[ , 1 ], density.estimate.nb10.simulated.PTdensity.3[ , 2 ], lty = 2 ) colnames( fit.nb10.simulated.PTdensity.3$save.state$thetasave ) # [1] "mu:nb10" "sigma:nb10" "alpha" mu.star.nb10.simulated.PTdensity.3 <- fit.nb10.simulated.PTdensity.3$save.state$thetasave[ , 1 ] sigma.star.nb10.simulated.PTdensity.3 <- fit.nb10.simulated.PTdensity.3$save.state$thetasave[ , 2 ] alpha.star.nb10.simulated.PTdensity.3 <- fit.nb10.simulated.PTdensity.3$save.state$thetasave[ , 3 ] mean( mu.star.nb10.simulated.PTdensity.3 ) # 404.5878 sd( mu.star.nb10.simulated.PTdensity.3 ) # 0.6155323 print( rho.hat.nb10.simulated.PTdensity.3 <- acf( mu.star.nb10.simulated.PTdensity.3, plot = F )$acf[ 2 ] ) # -0.0519037 print( mcmc.se.nb10.simulated.PTdensity.3 <- sqrt( ( 1 + rho.hat.nb10.simulated.PTdensity.3 ) / ( 1 - rho.hat.nb10.simulated.PTdensity.3 ) ) * sqrt( mean( sigma.star.nb10.simulated.PTdensity.3 ) / nsave ) ) # 0.0765482 hist( nb10.simulated, nclass = 100, probability = T, main = '' ) lines( nb10.simulated, fit.nb10.simulated.PTdensity.3$cpo, lty = 2 ) hist( nb10.simulated, nclass = 50, probability = T, main = '' ) lines( nb10.simulated, fit.nb10.simulated.PTdensity.3$cpo, lty = 2 ) hist( nb10.simulated, nclass = 25, probability = T, main = '' ) lines( nb10.simulated, fit.nb10.simulated.PTdensity.3$cpo, lty = 2 ) plot( density( mu.star.nb10.simulated.PTdensity.3 ), main = '', xlab = 'mu' ) plot( density( sigma.star.nb10.simulated.PTdensity.3 ), main = '', xlab = 'sigma' ) acf( mu.star.nb10.simulated.PTdensity.3, main = 'mu' ) ################ run 4 nsave <- 10000 mcmc.nb10.simulated.PTdensity.4 <- list( nburn = 500, nsave = nsave, nskip = 500, ndisplay = 100, tune1 = 1.0, tune2 = 1.2 ) prior.nb10.simulated.PTdensity.4 <- list( a0 = NULL, alpha = 1, M = 6, m0 = 404.6, S0 = 0.43, tau1 = 98.4, tau2 = 4159.2 ) fit.nb10.simulated.PTdensity.4 <- PTdensity( y = nb10.simulated, ngrid = 100, prior = prior.nb10.simulated.PTdensity.4, mcmc = mcmc.nb10.simulated.PTdensity.4, state = NULL, status = TRUE ) fit.nb10.simulated.PTdensity.4 Bayesian Density Estimation Using MPT (5,000,000 iterations in 56 minutes) Call: PTdensity.default(y = nb10.simulated, ngrid = 100, prior = prior.nb10.simulated.PTdensity.4, mcmc = mcmc.nb10.simulated.PTdensity.4, state = NULL, status = TRUE) Posterior Predictive Distributions (log): Min. 1st Qu. Median Mean 3rd Qu. Max. -6.250 -3.544 -3.103 -3.383 -2.939 -2.876 Posterior Inference of Parameters: mu:nb10.simulated sigma:nb10.simulated alpha 404.597 6.515 1.000 Acceptance Rate for Metropolis Step = 0.4318256 0.4782411 Number of Observations: 100 Number of Variables: 1 plot( fit.nb10.simulated.PTdensity.4 ) points( nb10.simulated, rep( 0, length( nb10 ) ) ) plot( fit.nb10.simulated.PTdensity.4, ask = TRUE, output = "param" ) density.estimate.nb10.simulated.PTdensity.4 <- cbind( fit.nb10.simulated.PTdensity.4$x1, fit.nb10.simulated.PTdensity.4$dens ) hist( nb10, nclass = 25, probability = T, main = '' ) lines( density.estimate.nb10.simulated.PTdensity.4[ , 1 ], density.estimate.nb10.simulated.PTdensity.4[ , 2 ], lty = 2 ) colnames( fit.nb10.simulated.PTdensity.4$save.state$thetasave ) # [1] "mu:nb10" "sigma:nb10" "alpha" mu.star.nb10.simulated.PTdensity.4 <- fit.nb10.simulated.PTdensity.4$save.state$thetasave[ , 1 ] sigma.star.nb10.simulated.PTdensity.4 <- fit.nb10.simulated.PTdensity.4$save.state$thetasave[ , 2 ] alpha.star.nb10.simulated.PTdensity.4 <- fit.nb10.simulated.PTdensity.4$save.state$thetasave[ , 3 ] mean( mu.star.nb10.simulated.PTdensity.4 ) # 404.5967 sd( mu.star.nb10.simulated.PTdensity.4 ) # 0.6282579 print( rho.hat.nb10.simulated.PTdensity.4 <- acf( mu.star.nb10.simulated.PTdensity.4, plot = F )$acf[ 2 ] ) # -0.005189169 print( mcmc.se.nb10.simulated.PTdensity.4 <- sqrt( ( 1 + rho.hat.nb10.simulated.PTdensity.4 ) / ( 1 - rho.hat.nb10.simulated.PTdensity.4 ) ) * sqrt( mean( sigma.star.nb10.simulated.PTdensity.4 ) / nsave ) ) # 0.02539234 hist( nb10.simulated, nclass = 100, probability = T, main = '' ) lines( nb10.simulated, fit.nb10.simulated.PTdensity.4$cpo, lty = 2 ) hist( nb10.simulated, nclass = 50, probability = T, main = '' ) lines( nb10.simulated, fit.nb10.simulated.PTdensity.4$cpo, lty = 2 ) hist( nb10.simulated, nclass = 25, probability = T, main = '' ) lines( nb10.simulated, fit.nb10.simulated.PTdensity.4$cpo, lty = 2 ) plot( density( mu.star.nb10.simulated.PTdensity.4 ), main = '', xlab = 'mu' ) plot( density( sigma.star.nb10.simulated.PTdensity.4 ), main = '', xlab = 'sigma' ) acf( mu.star.nb10.simulated.PTdensity.4, main = 'mu' ) ################ run 5 nsave <- 1000 mcmc.nb10.simulated.PTdensity.5 <- list( nburn = 500, nsave = nsave, nskip = 500, ndisplay = 100, tune1 = 1.0, tune2 = 1.2 ) prior.nb10.simulated.PTdensity.5 <- list( a0 = 1, b0 = 1, M = 6 ) fit.nb10.simulated.PTdensity.5 <- PTdensity( y = nb10.simulated, ngrid = 100, prior = prior.nb10.simulated.PTdensity.5, mcmc = mcmc.nb10.simulated.PTdensity.5, state = NULL, status = TRUE ) fit.nb10.simulated.PTdensity.5 Bayesian Density Estimation Using MPT (500,000 iterations in 7.0 minutes) Call: PTdensity.default(y = nb10.simulated, ngrid = 100, prior = prior.nb10.simulated.PTdensity.5, mcmc = mcmc.nb10.simulated.PTdensity.5, state = NULL, status = TRUE) Posterior Predictive Distributions (log): Min. 1st Qu. Median Mean 3rd Qu. Max. -6.286 -3.527 -3.076 -3.348 -2.901 -2.822 Posterior Inference of Parameters: mu:nb10.simulated sigma:nb10.simulated alpha 404.542 6.670 3.369 Acceptance Rate for Metropolis Step = 0.5893659 0.6336451 0.6456211 Number of Observations: 100 Number of Variables: 1 ################ run 6 nsave <- 1000 mcmc.nb10.simulated.PTdensity.6 <- list( nburn = 500, nsave = nsave, nskip = 500, ndisplay = 100, tune1 = 1.5, tune2 = 1.8, tune3 = 1.9 ) prior.nb10.simulated.PTdensity.6 <- list( a0 = 1, b0 = 1, M = 6 ) fit.nb10.simulated.PTdensity.6 <- PTdensity( y = nb10.simulated, ngrid = 100, prior = prior.nb10.simulated.PTdensity.6, mcmc = mcmc.nb10.simulated.PTdensity.6, state = NULL, status = TRUE ) fit.nb10.simulated.PTdensity.6 Bayesian Density Estimation Using MPT (500,000 iterations in 7.0 minutes) Call: PTdensity.default(y = nb10.simulated, ngrid = 100, prior = prior.nb10.simulated.PTdensity.6, mcmc = mcmc.nb10.simulated.PTdensity.6, state = NULL, status = TRUE) Posterior Predictive Distributions (log): Min. 1st Qu. Median Mean 3rd Qu. Max. -6.258 -3.515 -3.075 -3.346 -2.903 -2.812 Posterior Inference of Parameters: mu:nb10.simulated sigma:nb10.simulated alpha 404.605 6.674 3.406 Acceptance Rate for Metropolis Step = 0.4982433 0.5494975 0.4744148 Number of Observations: 100 Number of Variables: 1 ################ run 7 nsave <- 1000 mcmc.nb10.simulated.PTdensity.7 <- list( nburn = 500, nsave = nsave, nskip = 500, ndisplay = 100, tune1 = 1.75, tune2 = 2.4, tune3 = 2.0 ) prior.nb10.simulated.PTdensity.7 <- list( a0 = 1, b0 = 1, M = 6 ) fit.nb10.simulated.PTdensity.7 <- PTdensity( y = nb10.simulated, ngrid = 100, prior = prior.nb10.simulated.PTdensity.7, mcmc = mcmc.nb10.simulated.PTdensity.7, state = NULL, status = TRUE ) fit.nb10.simulated.PTdensity.7 Bayesian Density Estimation Using MPT (500,000 iterations in 7.0 minutes) Call: PTdensity.default(y = nb10.simulated, ngrid = 100, prior = prior.nb10.simulated.PTdensity.7, mcmc = mcmc.nb10.simulated.PTdensity.7, state = NULL, status = TRUE) Posterior Predictive Distributions (log): Min. 1st Qu. Median Mean 3rd Qu. Max. -6.266 -3.518 -3.078 -3.348 -2.903 -2.812 Posterior Inference of Parameters: mu:nb10.simulated sigma:nb10.simulated alpha 404.603 6.685 3.375 Acceptance Rate for Metropolis Step = 0.4594776 0.4781535 0.4575992 Number of Observations: 100 Number of Variables: 1 plot( fit.nb10.simulated.PTdensity.7 ) points( nb10.simulated, rep( 0, length( nb10 ) ) ) plot( fit.nb10.simulated.PTdensity.7, ask = TRUE, output = "param" ) density.estimate.nb10.simulated.PTdensity.7 <- cbind( fit.nb10.simulated.PTdensity.7$x1, fit.nb10.simulated.PTdensity.7$dens ) hist( nb10, nclass = 25, probability = T, main = '' ) lines( density.estimate.nb10.simulated.PTdensity.7[ , 1 ], density.estimate.nb10.simulated.PTdensity.7[ , 2 ], lty = 2 ) colnames( fit.nb10.simulated.PTdensity.7$save.state$thetasave ) # [1] "mu:nb10" "sigma:nb10" "alpha" mu.star.nb10.simulated.PTdensity.7 <- fit.nb10.simulated.PTdensity.7$save.state$thetasave[ , 1 ] sigma.star.nb10.simulated.PTdensity.7 <- fit.nb10.simulated.PTdensity.7$save.state$thetasave[ , 2 ] alpha.star.nb10.simulated.PTdensity.7 <- fit.nb10.simulated.PTdensity.7$save.state$thetasave[ , 3 ] mean( mu.star.nb10.simulated.PTdensity.7 ) # 404.6032 sd( mu.star.nb10.simulated.PTdensity.7 ) # 1.677472 print( rho.hat.nb10.simulated.PTdensity.7 <- acf( mu.star.nb10.simulated.PTdensity.7, plot = F )$acf[ 2 ] ) # -0.03938924 print( mcmc.se.nb10.simulated.PTdensity.7 <- sqrt( ( 1 + rho.hat.nb10.simulated.PTdensity.7 ) / ( 1 - rho.hat.nb10.simulated.PTdensity.7 ) ) * sqrt( mean( sigma.star.nb10.simulated.PTdensity.7 ) / nsave ) ) # 0.07860199 hist( nb10.simulated, nclass = 100, probability = T, main = '' ) lines( nb10.simulated, fit.nb10.simulated.PTdensity.7$cpo, lty = 2 ) hist( nb10.simulated, nclass = 50, probability = T, main = '' ) lines( nb10.simulated, fit.nb10.simulated.PTdensity.7$cpo, lty = 2 ) hist( nb10.simulated, nclass = 25, probability = T, main = '' ) lines( nb10.simulated, fit.nb10.simulated.PTdensity.7$cpo, lty = 2 ) plot( density( mu.star.nb10.simulated.PTdensity.7 ), main = '', xlab = 'mu' ) plot( density( sigma.star.nb10.simulated.PTdensity.7 ), main = '', xlab = 'sigma' ) acf( mu.star.nb10.simulated.PTdensity.7, main = 'mu' ) ################ run 8 nsave <- 10000 mcmc.nb10.simulated.PTdensity.8 <- list( nburn = 500, nsave = nsave, nskip = 5000, ndisplay = 100, tune1 = 1.75, tune2 = 2.5, tune3 = 2.0 ) prior.nb10.simulated.PTdensity.8 <- list( a0 = 1, b0 = 1, M = 6 ) fit.nb10.simulated.PTdensity.8 <- PTdensity( y = nb10.simulated, ngrid = 100, prior = prior.nb10.simulated.PTdensity.8, mcmc = mcmc.nb10.simulated.PTdensity.8, state = NULL, status = TRUE ) fit.nb10.simulated.PTdensity.8 Bayesian Density Estimation Using MPT (50,000,000 iterations in 11.6 hours) Call: PTdensity.default(y = nb10.simulated, ngrid = 100, prior = prior.nb10.simulated.PTdensity.8, mcmc = mcmc.nb10.simulated.PTdensity.8, state = NULL, status = TRUE) Posterior Predictive Distributions (log): Min. 1st Qu. Median Mean 3rd Qu. Max. -6.270 -3.524 -3.074 -3.347 -2.896 -2.819 Posterior Inference of Parameters: mu:nb10.simulated sigma:nb10.simulated alpha 404.617 6.669 3.386 Acceptance Rate for Metropolis Step = 0.4584765 0.4672574 0.4577286 Number of Observations: 100 Number of Variables: 1 plot( fit.nb10.simulated.PTdensity.8 ) points( nb10.simulated, rep( 0, length( nb10 ) ) ) plot( fit.nb10.simulated.PTdensity.8, ask = TRUE, output = "param" ) density.estimate.nb10.simulated.PTdensity.8 <- cbind( fit.nb10.simulated.PTdensity.8$x1, fit.nb10.simulated.PTdensity.8$dens ) hist( nb10, nclass = 25, probability = T, main = '' ) lines( density.estimate.nb10.simulated.PTdensity.8[ , 1 ], density.estimate.nb10.simulated.PTdensity.8[ , 2 ], lty = 2 ) colnames( fit.nb10.simulated.PTdensity.8$save.state$thetasave ) # [1] "mu:nb10" "sigma:nb10" "alpha" mu.star.nb10.simulated.PTdensity.8 <- fit.nb10.simulated.PTdensity.8$save.state$thetasave[ , 1 ] sigma.star.nb10.simulated.PTdensity.8 <- fit.nb10.simulated.PTdensity.8$save.state$thetasave[ , 2 ] alpha.star.nb10.simulated.PTdensity.8 <- fit.nb10.simulated.PTdensity.8$save.state$thetasave[ , 3 ] mean( mu.star.nb10.simulated.PTdensity.8 ) # 404.617 sd( mu.star.nb10.simulated.PTdensity.8 ) # 1.659243 print( rho.hat.nb10.simulated.PTdensity.8 <- acf( mu.star.nb10.simulated.PTdensity.8, plot = F )$acf[ 2 ] ) # -0.008240254 print( mcmc.se.nb10.simulated.PTdensity.8 <- sqrt( ( 1 + rho.hat.nb10.simulated.PTdensity.8 ) / ( 1 - rho.hat.nb10.simulated.PTdensity.8 ) ) * sqrt( mean( sigma.star.nb10.simulated.PTdensity.8 ) / nsave ) ) # 0.02561224 hist( nb10.simulated, nclass = 100, probability = T, main = '' ) lines( nb10.simulated, fit.nb10.simulated.PTdensity.8$cpo, lty = 2 ) hist( nb10.simulated, nclass = 50, probability = T, main = '' ) lines( nb10.simulated, fit.nb10.simulated.PTdensity.8$cpo, lty = 2 ) hist( nb10.simulated, nclass = 25, probability = T, main = '' ) lines( nb10.simulated, fit.nb10.simulated.PTdensity.8$cpo, lty = 2 ) plot( density( mu.star.nb10.simulated.PTdensity.8 ), main = '', xlab = 'mu' ) plot( density( sigma.star.nb10.simulated.PTdensity.8 ), main = '', xlab = 'sigma' ) acf( mu.star.nb10.simulated.PTdensity.8, main = 'mu' ) print( c( mean( alpha.star.nb10.simulated.PTdensity.8 ), sd( alpha.star.nb10.simulated.PTdensity.8 ) ) ) # 3.386029 1.541352 print( quantile( alpha.star.nb10.simulated.PTdensity.8 ) ) # 0% 25% 50% 75% 100% # 0.2576925 2.2811824 3.1045414 4.2101597 12.5203937 print( mean( alpha.star.nb10.simulated.PTdensity.8 < 1 ) ) # 0.0113 ###################### end simulated NB10 PTdensity analysis ###################### begin simulated NB10 DPdensity analysis nb10 <- 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. ) n.nb10 <- length( nb10 ) nb10.simulated <- mean( nb10 ) + sd( nb10 ) * qnorm( ( 1:n.nb10 - 0.5 ) / n.nb10 ) nsave <- 10000 mcmc.nb10.simulated.DPdensity.1 <- list( nburn = 1000, nsave = nsave, nskip = 10, ndisplay = 100 ) prior.nb10.simulated.DPdensity.1 <- list( alpha = 1, m2 = rep( 0, 1 ), s2 = diag( 100000, 1 ), psiinv2 = solve( diag( 0.5, 1 ) ), nu1 = 4, nu2 = 4, tau1 = 1, tau2 = 100 ) fit.nb10.simulated.DPdensity.1 <- DPdensity( y = nb10.simulated, prior = prior.nb10.simulated.DPdensity.1, mcmc = mcmc.nb10.simulated.DPdensity.1, state = NULL, status = TRUE ) fit.nb10.simulated.DPdensity.1 DPM of normals model for density estimation (100,000 iterations in 2.5 minutes) Call: DPdensity.default(y = nb10.simulated, prior = prior.nb10.simulated.DPdensity.1, mcmc = mcmc.nb10.simulated.DPdensity.1, state = NULL, status = TRUE) Posterior Predictive Distributions (log): Min. 1st Qu. Median Mean 3rd Qu. Max. -6.392 -3.485 -3.026 -3.311 -2.838 -2.784 Posterior Inference of Parameters: m1-nb10.simulated k0 psi1-nb10.simulated 367.24596 0.01724 0.32215 ncluster 2.10980 Number of Observations: 100 Number of Variables: 1 plot( fit.nb10.simulated.DPdensity.1, ask = FALSE ) density.nb10.simulated.DPdensity.1 <- cbind( fit.nb10.simulated.DPdensity.1$x1, fit.nb10.simulated.DPdensity.1$dens ) hist( nb10.simulated, breaks = 30, probability = T, main = '' ) lines( density.nb10.simulated.DPdensity.1[ , 1 ], density.nb10.simulated.DPdensity.1[ , 2 ], lwd = 2 ) plot( fit.nb10.simulated.DPdensity.1, ask = T, output = "param" ) DPrandom( fit.nb10.simulated.DPdensity.1 ) plot( DPrandom( fit.nb10.simulated.DPdensity.1, predictive = TRUE ), ask = T ) plot( DPrandom( fit.nb10.simulated.DPdensity.1, predictive = TRUE ), ask = T, hpd = FALSE ) plot( DPrandom( fit.nb10.simulated.DPdensity.1 ), ask = T, hpd = FALSE ) colnames( fit.nb10.simulated.DPdensity.1$save.state$thetasave ) # "m1-nb10.simulated" "k0" "psi1-nb10.simulated" "ncluster" "alpha" mu.star.nb10.simulated.DPdensity.1 <- fit.nb10.simulated.DPdensity.1$save.state$thetasave[ , 1 ] print( mean( mu.star.nb10.simulated.DPdensity.1 ) ) # print( sd( mu.star.nb10.simulated.DPdensity.1 ) ) # ###################### end simulated NB10 DPdensity analysis ###################### begin simulated NB10 WinBUGS PT analysis ########## round 1, with the code that doesn't monitor \int y G ( dy ) ################################################## # MPT code, uses "ones" trick # c ~ gamma(5,1), mu ~ dnorm(0,0.001) # sigma ~ dgamma(0.01,0.01) ################################################## model{ for(i in 1:n){ ones[i] <- 1 z[i] <- (y[i]-mu)/sigma ind[i] <- trunc(sets*phi(z[i])+1) like[i] <- coi*exp(-0.5*z[i]*z[i])*p[ind[i]]/sigma # E{cpo.inv[i]|data}=1/CPO[i] cpo.inv[i] <- 1/(like[i]*sets) ones[i] ~ dbern(like[i]) } x[1,1] <- 0.5; x[1,2] <- 0.5 for(k in 1:j.int[Jt]){ m[1,k] <- log(x[1,1]) m[1,j.int[Jt]+k] <- log(x[1,2]) } for(i in 2:Jt){ par.lev[i] <- c*i*i for(j in 1:j.int[i]){ x[i,2*j-1] ~ dbeta(par.lev[i],par.lev[i]); x[i,2*j] <- 1-x[i,2*j-1] for(k in 1:j.int[Jt-i+1]){ m[i,(2*j-2)*j.int[Jt-i+1]+k] <- log(x[i,2*j-1]) m[i,(2*j-1)*j.int[Jt-i+1]+k] <- log(x[i,2*j]) }}} for(i in 1:sets){ p[i] <- exp(sum(m[1:Jt,i])) } for(i in 1:200){ grid[i]<-left+(right-left)*i/201 z.f[i] <- (grid[i]-mu)/sigma ind.f[i] <- trunc(sets*phi(z.f[i])+1) # f[i] is predictive density over grid[i] f[i] <- sets*coi*exp(-0.5*z.f[i]*z.f[i])*p[ind.f[i]]/sigma } mu ~ dnorm(0,0.001) sigma ~ dgamma(0.01,0.01) coi <- 0.3989422804 c ~ dgamma( 5, 1 ) left <- mean( y[ ] ) - 3.25 * sd( y[ ] ) right <- mean( y[ ] ) + 3.25 * sd( y[ ] ) } ################################################## # starting values ################################################## list( mu = 404.59, sigma = 6.458053, c = 2 ) ################################################## # data; n=100 is sample size; Jt=5 is level of finite Polya tree # j.int is number of independent beta distributions at each level # left and right are interval endpoints to evaluate predictive density ################################################## list(y=c(387.9325, 390.5564, 391.9152, 392.8727, 393.6261, 394.2547, 394.7985, 395.2808, 395.7162, 396.1147, 396.4834, 396.8275, 397.1509, 397.4567, 397.7473, 398.0247, 398.2906, 398.5462, 398.7926, 399.0310, 399.2620, 399.4864, 399.7048, 399.9178, 400.1259, 400.3294, 400.5288, 400.7244, 400.9165, 401.1054, 401.2914, 401.4747, 401.6556, 401.8342, 402.0107, 402.1853, 402.3581, 402.5294, 402.6993, 402.8678, 403.0352, 403.2016, 403.3670, 403.5316, 403.6956, 403.8590, 404.0219, 404.1845, 404.3468, 404.5089, 404.6711, 404.8332, 404.9955, 405.1581, 405.3210, 405.4844, 405.6484, 405.8130, 405.9784, 406.1448, 406.3122, 406.4807, 406.6506, 406.8219, 406.9947, 407.1693, 407.3458, 407.5244, 407.7053, 407.8886, 408.0746, 408.2635, 408.4556, 408.6512, 408.8506, 409.0541, 409.2622, 409.4752, 409.6936, 409.9180, 410.1490, 410.3874, 410.6338, 410.8894, 411.1553, 411.4327, 411.7233, 412.0291, 412.3525, 412.6966, 413.0653, 413.4638, 413.8992, 414.3815, 414.9253, 415.5539, 416.3073, 417.2648, 418.6236, 421.2475), n=100, Jt=5,sets=32, j.int=c(1,2,4,8,16)) ########## R post-processing nb10.simulated.pt.predictive.1 <- matrix( scan( "nb10-simulated-pt-predictive-1.txt" ), 2000000, 2, byrow = T )[ , 2 ] nb10.simulated.pt.predictive.1 <- matrix( nb10.simulated.pt.predictive.1, 10000, 200, byrow = F ) apply( nb10.simulated.pt.predictive.1, 2, mean ) left <- mean( nb10.simulated ) - 3.25 * sd( nb10.simulated ) right <- mean( nb10.simulated ) + 3.25 * sd( nb10.simulated ) y.grid <- matrix( left + ( right - left ) * ( 1:200 ) / 201, 200, 1 ) sum( nb10.simulated.pt.predictive.1 * as.vector( y.grid ) / sum( nb10.simulated.pt.predictive.1 ) ) sum( nb10.simulated.pt.predictive.1[ 2, ] ) predictive.mean.star <- as.vector( nb10.simulated.pt.predictive.1 %*% y.grid ) * ( y.grid[ 2 ] - y.grid[ 1 ] ) # mean( predictive.mean.star ) # sd( predictive.mean.star ) ########## round 2, with the code that does monitor \int y G ( dy ) ######### MPT WinBUGS code starts below this line ################################################## # MPT code, uses "ones" trick # c ~ gamma(5,1), mu ~ dnorm(0,0.001) # sigma ~ dgamma(0.01,0.01) ################################################## model{ for(i in 1:n){ ones[i] <- 1 z[i] <- (y[i]-mu)/sigma ind[i] <- trunc(sets*phi(z[i])+1) like[i] <- coi*exp(-0.5*z[i]*z[i])*p[ind[i]]/sigma # E{cpo.inv[i]|data}=1/CPO[i] cpo.inv[i] <- 1/(like[i]*sets) ones[i]~dbern(like[i]) } # CHANGED BELOW TO NOT BE MEDIAN-MU ANYMORE AND UPDATE X[1,1] # x[1,1] <- 0.5; x[1,2] <- 0.5 # for(k in 1:j.int[Jt]){ # m[1,k] <- log(x[1,1]) # m[1,j.int[Jt]+k] <- log(x[1,2]) # } for(i in 1:Jt){ par.lev[i] <- c*i*i for(j in 1:j.int[i]){ x[i,2*j-1] ~ dbeta(par.lev[i],par.lev[i]); x[i,2*j] <- 1-x[i,2*j-1] for(k in 1:j.int[Jt-i+1]){ m[i,(2*j-2)*j.int[Jt-i+1]+k] <- log(x[i,2*j-1]) m[i,(2*j-1)*j.int[Jt-i+1]+k] <- log(x[i,2*j]) }}} for(i in 1:sets){ p[i] <- exp(sum(m[1:Jt,i])) g.prod[i]<-p[i]*m.g[i] } G.mean<-sum(g.prod[1:sets])*sigma+mu for(i in 1:200){ grid[i]<-left+(right-left)*i/201 z.f[i] <- (grid[i]-mu)/sigma ind.f[i] <- trunc(sets*phi(z.f[i])+1) # f[i] is predictive density over grid[i] f[i] <- sets*coi*exp(-0.5*z.f[i]*z.f[i])*p[ind.f[i]]/sigma } mu ~ dnorm(0,0.001) sigma ~ dgamma(0.01,0.01) coi <- 0.3989422804 c ~ dgamma(5,1) } ################################################## # starting values ################################################## list(mu=21,sigma=6,c=2) ################################################## # data; n=82 is sample size; Jt=5 is level of finite Polya tree # j.int is number of independent beta distributions at each level # left and right are interval endpoints to evaluate predictive density ################################################## list(y=c(9.172,9.35,9.483,9.558,9.775,10.227,10.406,16.084, 16.17,18.419,18.552,18.6,18.927,19.052,19.07,19.33,19.343, 19.349,19.44,19.473,19.529,19.541,19.547,19.663,19.846,19.856, 19.863,19.914,19.918,19.973,19.989,20.166,20.175,20.179,20.196,20.215,20.221,20.415, 20.629,20.795,20.821,20.846,20.875,20.986,21.137,21.492,21.701,21.814,21.921,21.96, 22.185,22.209,22.242,22.249,22.314,22.374,22.495,22.746,22.747,22.888,22.914,23.206, 23.241,23.263,23.484,23.538,23.542,23.666,23.706,23.711,24.129,24.285,24.289,24.366, 24.717,24.99,25.633,26.96,26.995,32.065,32.789,34.279), n=82, Jt=5,sets=32, j.int=c(1,2,4,8,16),left=0,right=40, m.g=c(-0.0703,-0.0526,-0.0443,-0.0384,-0.0336,-0.0296,-0.0259, -0.0226,-0.0195,-0.0166,-0.0139,-0.0112,-0.0086,-0.0061, -0.0036,-0.0012,0.0012,0.0036,0.0061,0.0086,0.0112,0.0139, 0.0166,0.0195,0.0226,0.0259,0.0296,0.0336,0.0384,0.0443,0.0526,0.0703)) ### post-processing of a run with c = 1 f.raw <- scan( "f.txt" ) f.raw <- f.raw[ 2 * ( 1:22000000 ) ] f <- matrix( f.raw, 110000, 200 ) f.mean <- apply( f, 2, mean ) left <- mean( nb10.simulated ) - 3.25 * sd( nb10.simulated ) right <- mean( nb10.simulated ) + 3.25 * sd( nb10.simulated ) nb10.simulated.grid <- seq( left, right, length.out = 200 ) plot( nb10.simulated.grid, f.mean, type = 'l' ) c G.mean mu sigma largest monitoring mean SD mean SD mean SD mcse run 1 404.6 0.8755 404.6 0.9001 6.606 0.7509 0.01105 107000 try again use 50,000 monitoring after a burn-in of 1,000 (MCSEs on the order of 0.05) prior on mu is N ( \mu_0, precision \tau ) prior on sigma is U ( 0, A ), A chosen not to truncate likelihood; if necessary tune A based on initial monitoring run of 5,000 Gk = \Gamma( k, 1 ) c = 1 mcse (mean) = 0.15 1* = with new simulated data (mean 404.18, sd 6.535489) PN = parametric normal model G.mean mu sigma c c mu_0 \tau A mean SD mean SD mean SD mean SD G1 0.0 1.0E-6 20 404.6 1.693 404.6 1.738 6.747 0.7153 3.262 1.534 1 0.0 1.0E-6 20 404.5 3.128 404.5 3.22 7.114 1.11 -- -- 1* 0.0 1.0E-6 20 403.1 2.211 403.0 2.275 6.915 0.9255 -- -- ---------------- G5 0.0 1.0E-6 20 404.6 1.315 404.6 1.346 6.658 0.5869 6.374 2.539 G10,2 0.0 1.0E-6 20 404.6 1.626 404.6 1.39 6.663 0.599 5.736 1.626 ---------------- G10 0.0 1.0E-6 20 404.6 1.086 404.6 1.108 6.607 0.5482 10.93 3.208 10 0.0 1.0E-6 20 404.6 1.107 404.6 1.129 6.619 0.5479 -- -- ---------------- G20 0.0 1.0E-6 20 404.6 0.9574 404.6 0.9733 6.589 0.5148 20.63 4.531 ---------------- G50 0.0 1.0E-6 20 404.6 0.8153 404.6 0.8244 6.562 0.4905 50.24 7.076 ---------------- G100 0.0 1.0E-6 20 404.6 0.7479 404.6 0.7533 6.554 0.4838 100.1 9.986 100 0.0 1.0E-6 20 404.6 0.7429 404.6 0.7482 6.553 0.4817 -- -- ---------------- G200 0.0 1.0E-6 20 404.6 0.7084 404.6 0.7114 6.544 0.4736 200.1 14.22 ---------------- G500 0.0 1.0E-6 20 404.6 0.6783 404.6 0.6795 6.542 0.4754 499.8 22.33 ---------------- PN 0.0 1.0E-6 20 -- -- 404.6 0.6528 6.538 0.4714 ---------------- U( 0, 1000) 0.0 1.0E-6 20 404.6 0.6939 404.6 0.6961 6.549 0.4752 528.1 277.7 mcse for mean of c is 6.997 posterior for c is close to prior nb10.simulated.data.2 <- sort( rnorm( 100, mean( y ), sd( y ) ) ) dput( nb10.simulated.data.2, "nb10-simulated-data-2.txt" ) ###################### end simulated NB10 WinBUGS PT analysis