###################### begin NB10 PTlm analysis

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

nb10 <- y

state <- NULL

nburn <- 5000

nsave <- 10000

nskip <- 20

ndisplay <- 100

mcmc <- list( nburn = nburn, nsave = nsave, nskip = nskip,
  ndisplay = ndisplay )

prior.1 <- list( alpha = 1, beta0 = 0, Sbeta0 = diag( 1000, 1 ),
  tau1 = 0.01, tau2 = 0.01, M = 6 )

fit.nb10.1 <- PTlm( formula = y ~ 1, prior = prior.1, mcmc = mcmc,
  state = state, status = TRUE )

summary( fit.nb10.1 )

Bayesian Semiparametric Regression Model

Call:
PTlm.default(formula = y ~ 1, prior = prior.1, mcmc = mcmc, state = state, 
    status = TRUE )

Posterior Predictive Distributions (log):
   Min.  1st Qu.   Median     Mean  3rd Qu.     Max.  
-13.100   -2.915   -2.551   -3.018   -2.393   -2.196  

Regression coefficients:
             Mean       Median     Std. Dev.  Naive Std.Error  95%HPD-Low
(Intercept)  404.18094  404.29339    0.57301    0.00573        403.00049 
mu             0.00000    0.00000    0.00000    0.00000          0.00000 
sigma2        73.41452   71.47945   12.70997    0.12710         48.35440 
             95%HPD-Upp
(Intercept)  405.01635 
mu             0.00000 
sigma2        99.53967 
 
Acceptance Rate for Metropolis Step =  0.4896698 0 0.447507 

Number of Observations: 100

summary( fit.nb10.1, hpd = FALSE )

Bayesian Semiparametric Regression Model

Call:
PTlm.default(formula = y ~ 1, prior = prior.1, mcmc = mcmc, state = state, 
    status = TRUE )

Posterior Predictive Distributions (log):
   Min.  1st Qu.   Median     Mean  3rd Qu.     Max.  
-13.100   -2.915   -2.551   -3.018   -2.393   -2.196  

Regression coefficients:
             Mean       Median     Std. Dev.  Naive Std.Error  95%CI-Low
(Intercept)  404.18094  404.29339    0.57301    0.00573        403.14151
mu             0.00000    0.00000    0.00000    0.00000          0.00000
sigma2        73.41452   71.47945   12.70997    0.12710         50.85340
             95%CI-Upp
(Intercept)  405.47888
mu             0.00000
sigma2       103.20073
 
Acceptance Rate for Metropolis Step =  0.4896698 0 0.447507 

Number of Observations: 100

plot( fit.nb10.1, ask = TRUE )

# Table of Pseudo Contour Probabilities

anova( fit.nb10.1 ) # didn't work
 
colnames( fit.nb10.1$save.state$thetasave )

[1] "(Intercept)" "mu"          "sigma2"      "alpha"  

theta.star <- fit.nb10.1$save.state$thetasave[ , 1 ]

sigma.2.star <- fit.nb10.1$save.state$thetasave[ , 3 ]

par( mfrow = c( 1, 1 ) )

plot( density( theta.star, adjust = 4 ), type = 'l', xlab = 'theta',
  main = '' )

plot( density( sigma.2.star, adjust = 2 ), type = 'l', xlab = 'sigma^2',
  main = '' )

prior.2 <- list( alpha = 1, beta0 = 404.59, Sbeta0 = diag( 1, 1 ),
  tau1 = 1, tau2 = 1, M = 6, mub = 404.59, Sb = 1, frstlprob = FALSE )

fit.nb10.2 <- PTlm( formula = y ~ 1, prior = prior.2, mcmc = mcmc,
  state = state, status = TRUE )

prior.2 <- list( alpha = 1, beta0 = 0, Sbeta0 = diag( 1000, 1 ),
  tau1 = 0.01, tau2 = 0.01, M = 6, mub = 0, Sb = 1000, frstlprob = FALSE )

###################### end NB10 PTlm analysis

###################### begin NB10 PTdensity 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. )

state <- NULL

nburn <- 500

nsave <- 100000

nskip <- 50

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 )

fit.nb10.PTdensity.1

Bayesian Density Estimation Using MPT (took 59.5 min)

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.  
-12.200   -3.039   -2.583   -3.035   -2.400   -2.184  

Posterior Inference of Parameters:
   mu:nb10  sigma:nb10       alpha  
   405.384       9.012       1.000  

Acceptance Rate for Metropolis Step =  0.4484321 0.1586613 

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 )

sd( mu.star.nb10.PTdensity.1 )

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 ) / sqrt( nsave ) )

################ run 2

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 (37.5 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 18.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
  5.6 minutes)

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 )

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
  5.5 minutes)

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

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
  6.6 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
  5.5 minutes)

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
  5.5 minutes)

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
  9.0 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

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 <- 1000

mcmc.nb10.PTdensity.10 <- list( nburn = 500, nsave = nsave, 
  nskip = 500, 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 (500,000 iterations in
  x.x minutes)

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.9$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 )

# 

sd( mu.star.nb10.PTdensity.9 )

# 

print( rho.hat.nb10.PTdensity.9 <- 
  acf( mu.star.nb10.PTdensity.9, plot = F )$acf[ 2 ] )

# 

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 ) )

# 

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

z <- list( call = cl, y = y, varnames = varnames, modelname = model.name,
  cpo = foo$cpo, prior = prior, mcmc = mcmc, state = state, 
  save.state = save.state, nrec = foo$nrec, nvar = nvar, crand = crand,
  coefficients = coeff, f = f, grid1 = grid1, grid2 = grid2, fun1 = fun1,
  fun2 = fun2, x1 = x1, x2 = x2, dens = dens, acrate = acrate ) 

return( z )