@@ -371,11 +371,99 @@ probabilities (under equal prior model probabilities).
371371
372372``` r
373373logBF <- logmarglikM1 - logmarglikM3
374- PrM2 <- 0.5 / (0.5 + 0.5 * exp(logBF ))
375- PrM1 <- 1 - PrM2
376- knitr :: kable(cbind(BF = exp(logBF ), logBF = logBF , PrM1 , PrM2 ))
374+ PrM3 <- 0.5 / (0.5 + 0.5 * exp(logBF ))
375+ PrM1 <- 1 - PrM3
376+ knitr :: kable(cbind(logML1 = logmarglikM1 , logML3 = logmarglikM3 ,
377+ BF = exp(logBF ), logBF = logBF , Pr1 = PrM1 , Pr3 = PrM3 ))
378+ ```
379+
380+ | logML1 | logML3 | BF | logBF | Pr1 | Pr3 |
381+ | ----------:| ----------:| ----:| ----------:| ----:| ----:|
382+ | -3456.384 | -3305.694 | 0 | -150.6898 | 0 | 1 |
383+
384+ ### Example 10.10: Testing homogeneity against unobserved heterogeneity
385+
386+ We begin by loading the data and specifying the hyperparameters.
387+
388+ ``` r
389+ data(eyetracking , package = " BayesianLearningCode" )
390+ y <- eyetracking $ anomalies
391+ N <- length(y )
392+
393+ a0_tmp <- c(0.1 , 0.5 , 1 , 2 )
394+ m0_tmp <- c(1 , mean(y ), 5 , 10 , 20 )
395+ grid <- expand.grid(a0 = a0_tmp , m0 = m0_tmp )
396+ a0 <- grid $ a0
397+ b0 <- grid $ a0 / grid $ m0
398+ ```
399+
400+ Now we can compute and print the log marginal likelihoods.
401+
402+ ``` r
403+ logmarglikM1 <- a0 * log(b0 ) + lgamma(a0 + sum(y )) -
404+ lgamma(a0 ) - (a0 + sum(y )) * log(b0 + N ) -
405+ sum(lgamma(y + 1 ))
406+
407+ logmarglikM2 <- rep(NA_real_ , length(a0 ))
408+ for (i in seq_along(a0 )) {
409+ logmarglikM2 [i ] <- N * a0 [i ] * log(b0 [i ]) + sum(lgamma(a0 [i ] + y )) -
410+ N * lgamma(a0 [i ]) - sum(a0 [i ] + y ) * log(b0 [i ] + 1 ) -
411+ sum(lgamma(y + 1 ))
412+ }
413+
414+ knitr :: kable(matrix (logmarglikM1 , nrow = length(a0_tmp ), ncol = length(m0_tmp ),
415+ dimnames = list (a0 = a0_tmp , m0 = round(m0_tmp , 2 ))), digits = 2 )
416+ ```
417+
418+ | | 1 | 3.52 | 5 | 10 | 20 |
419+ | :----| --------:| --------:| --------:| --------:| --------:|
420+ | 0.1 | -472.91 | -472.78 | -472.78 | -472.82 | -472.87 |
421+ | 0.5 | -472.25 | -471.62 | -471.64 | -471.81 | -472.07 |
422+ | 1 | -472.45 | -471.20 | -471.25 | -471.59 | -472.11 |
423+ | 2 | -473.31 | -470.81 | -470.92 | -471.60 | -472.63 |
424+
425+ ``` r
426+
427+ knitr :: kable(matrix (logmarglikM2 , nrow = length(a0_tmp ), ncol = length(m0_tmp ),
428+ dimnames = list (a0 = a0_tmp , m0 = round(m0_tmp , 2 ))), digits = 2 )
429+ ```
430+
431+ | | 1 | 3.52 | 5 | 10 | 20 |
432+ | :----| --------:| --------:| --------:| --------:| --------:|
433+ | 0.1 | -249.61 | -237.68 | -238.22 | -241.61 | -246.80 |
434+ | 0.5 | -271.04 | -223.77 | -226.24 | -242.34 | -267.54 |
435+ | 1 | -316.77 | -241.38 | -245.87 | -276.12 | -324.87 |
436+ | 2 | -381.56 | -273.80 | -281.39 | -335.39 | -426.86 |
437+
438+ ### Example 10.11: Testing Poisson vs. negative binomial
439+
440+ In the negative binomial case, if $a_ {0}$ is fixed and $b_ {0}$ follows a
441+ beta prime prior, we have a closed-form expression for the marginal
442+ likelihood.
443+
444+ ``` r
445+ alpha_b0 <- rep(1 , length(a0 ))
446+ beta_b0 <- 1 + alpha_b0 / b0
447+ sdb0 <- sqrt(alpha_b0 * (alpha_b0 + beta_b0 - 1 ) /
448+ ((beta_b0 - 2 ) * (beta_b0 - 1 )^ 2 ))
449+ # > Warning in sqrt(alpha_b0 * (alpha_b0 + beta_b0 - 1)/((beta_b0 - 2) * (beta_b0 -
450+ # > : NaNs produced
451+ alpha_bN <- alpha_b0 + N * a0
452+ beta_bN <- beta_b0 + sum(y )
453+
454+ logmarglikM3 <- rep(NA_real_ , length(a0 ))
455+ for (i in seq_along(a0 )) {
456+ logmarglikM3 [i ] <- lbeta(alpha_bN [i ], beta_bN [i ]) - lbeta(alpha_b0 [i ], beta_b0 [i ]) -
457+ N * lgamma(a0 [i ]) + sum(lgamma(a0 [i ] + y )) - sum(lgamma(y + 1 ))
458+ }
459+
460+ knitr :: kable(matrix (logmarglikM3 , nrow = length(a0_tmp ), ncol = length(m0_tmp ),
461+ dimnames = list (a0 = a0_tmp , m0 = round(m0_tmp , 2 ))), digits = 2 )
377462```
378463
379- | BF | logBF | PrM1 | PrM2 |
380- | ----:| ----------:| -----:| -----:|
381- | 0 | -150.6898 | 0 | 1 |
464+ | | 1 | 3.52 | 5 | 10 | 20 |
465+ | :----| --------:| --------:| --------:| --------:| --------:|
466+ | 0.1 | -239.42 | -238.96 | -239.02 | -239.61 | -241.09 |
467+ | 0.5 | -226.13 | -225.82 | -225.89 | -226.55 | -228.38 |
468+ | 1 | -243.96 | -243.78 | -243.86 | -244.49 | -246.27 |
469+ | 2 | -276.60 | -276.55 | -276.65 | -277.22 | -278.83 |
0 commit comments