@@ -337,8 +337,8 @@ posterior mean of the regression effects together with the 2.5\% and 97.5\%
337337quantiles of the posterior distribution in the following table.
338338
339339``` {r}
340- res_conj1 <- regression_conjugate(y, X, b0 = 0, B0 = 10, c0= 2.5, C0= 1.5)
341- post.sd.conj1 = sqrt(diag((res_conj1$CN / res_conj1$cN) * res_conj1$BN))
340+ res_conj1 <- regression_conjugate(y, X, b0 = 0, B0 = 10, c0 = 2.5, C0 = 1.5)
341+ post.sd.conj1 <- sqrt(diag((res_conj1$CN / res_conj1$cN) * res_conj1$BN))
342342
343343knitr::kable(round(cbind(
344344 qt(0.025, df = 2 * res_conj1$cN) * post.sd.conj1 + res_conj1$bN,
@@ -439,8 +439,11 @@ reg_semiconj <- function(y, X, b0 = 0, B0 = 10000, c0 = 2.5, C0 = 1.5,
439439 sigma2s <- rep(NA_real_, M)
440440
441441 # starting value for sigma2
442- sigma2 <-start.sigma2# var(y) / 2
443-
442+ if (missing(start.sigma2)) {
443+ start.sigma2 <- var(y) / 2
444+ }
445+ sigma2 <- start.sigma2
446+
444447 for (m in 1:(burnin + M)) {
445448 # sample beta from the full conditional
446449 BN <- solve(B0inv + XX / sigma2)
@@ -461,27 +464,31 @@ reg_semiconj <- function(y, X, b0 = 0, B0 = 10000, c0 = 2.5, C0 = 1.5,
461464}
462465```
463466### Example 6.5: Movie data
464- We run the sampler
467+
468+ We run the sampler for 1000 draws starting with a large value for the
469+ innovation variance.
470+
465471``` {r}
466472set.seed(1)
467473M <- 1000L # number of draws after burn-in
468474post.draws <- reg_semiconj(y, X, b0 = 0, B0 = 10000, c0 = 2.5, C0 = 1.5,
469- burnin = 0L, M = M,start.sigma2= 10^6)
475+ burnin = 0L, M = M, start.sigma2 = 10^6)
470476```
471477
472- We plot the draws
478+ We insepct the draws using trace plots.
479+
473480``` {r, echo = -c(1:2)}
474481if (pdfplots) {
475482 pdf("6-4_2a.pdf", width = 6, height = 6)
476483 par(mar = c(2.5, 1.5, 1.5, .1), mgp = c(1.6, .6, 0))
477484}
478- par(mfrow= c(2,2))
479- for (i in (1:d)) {
480- plot(post.draws$betas[,i], type= "l", xlab= "m", ylab= "",
481- main = colnames(post.draws$betas)[i])
485+ par(mfrow = c(2, 2))
486+ for (i in seq_len(ncol(post.draws$betas))) {
487+ plot(post.draws$betas[, i], type = "l", xlab = "m", ylab = "",
488+ main = colnames(post.draws$betas)[i])
482489}
483- plot(post.draws$sigma2s,type= "l",xlab= "m", ylab= "",
484- main = expression(sigma^2))
490+ plot(post.draws$sigma2s, type = "l", xlab = "m", ylab = "",
491+ main = expression(sigma^2))
485492```
486493
487494### Example 6.6: Movie data
@@ -514,7 +521,7 @@ Next, we define the prior parameters and run the sampler.
514521set.seed(1)
515522M <- 20000L # number of draws after burn-in
516523post.draws <- reg_semiconj(y, X, b0 = 0, B0 = 10000, c0 = 2.5, C0 = 1.5,
517- burnin = 1000L, M = M,start.sigma2=var(y)/2 )
524+ burnin = 1000L, M = M)
518525```
519526
520527To summarize the results nicely, we compute equal-tailed 95% credible
@@ -876,12 +883,12 @@ if (pdfplots) {
876883
877884# Some functions that return the log-marginal densities
878885## Log-marginal and marginal for Triple Gamma
879- logmarginal_TG <- function(x, a, c, kappa2){
886+ logmarginal_TG <- function(x, a, c, kappa2) {
880887 -0.5*log(4*pi/kappa2) - lbeta(a, c) + 0.5*(log(a) - log(c)) + lgamma(c + 0.5) +
881888 log(gsl::hyperg_U(c + 0.5, -a + 1.5, kappa2*a*x^2/(4*c), give = FALSE, strict = TRUE))
882889}
883890## Log-marginal and marginal for Double Gamma
884- logmarginal_DG <- function(x, a, kappa2){
891+ logmarginal_DG <- function(x, a, kappa2) {
885892 0.5*(a + 0.5)*log(a*kappa2) - 0.5*log(pi) - (a-0.5)*log(2) - lgamma(a) + (a - 0.5)*log(abs(x)) +
886893 log(besselK(sqrt(a*kappa2)*abs(x), a-0.5))
887894}
0 commit comments