Now we consider densities that are bounded but have sufficiently
heavy tails that one or more of \(x_i ^ {r d +
1} f(x) ^ r, i = 1, \ldots, d\) are unbounded if \(r = 1/2\) is used.
Cauchy density
Suppose that \(X\) has a standard
Cauchy density, that is \(f(x) \propto (1 +
x^2)^{-1}, -\infty < x < \infty\). The bounding box cannot
be constructed if \(r < 1\). For
\(r = 1\) the bounding box parameters
\(b_1^-(r) = -1\) and \(b_1^+(r) = 1\) are attained in the limits
as \(x\) decreases/increases to
infinity respectively. This is fine in theory but using \(r > 1\) avoids this problem and the
largest probability of acceptance is obtained for \(r \approx 1.26\). Below are the estimated
bounding boxes and probabilities of acceptance for \(r = 1\) (left) and \(r = 1.26\) (right).
ptr_c <- create_xptr("logcauchy")
cauchy1 <- ru_rcpp(logf = ptr_c, log = TRUE, init = 0, r = 1, n = 1000)
cauchy126 <- ru_rcpp(logf = ptr_c, log = TRUE, init = 0, r = 1.26, n = 1000)
cbind(cauchy1$box, cauchy126$box)
#> box vals1 conv box vals1 conv
#> a 1.0000000 0.000 0 1.0000000 0.000000 0
#> b1minus -0.9999999 -2509.398 0 -0.8310221 -2.948272 0
#> b1plus 0.9999999 2509.398 0 0.8310221 2.948272 0
c(cauchy1$pa, cauchy126$pa)
#> [1] 0.7987220 0.8488964
For \(r = 1\) the algorithm finds
essentially the correct values of \(b_1^-\) and \(b_1^+\) even though these do not occur at a
turning point of the objective functions. The estimated probability of
acceptance is indeed larger for \(r =
1.26\).
Hierarchical 1-way ANOVA with few groups
Consider the hierarchical 1-way ANOVA model \[\begin{equation*}
Y_{ij} = \mu + \alpha_i + \epsilon_{ij}, \quad\mbox{for} \,\, i = 1,
\ldots, I, \,
j = 1, \ldots, n_i, \label{eqn:1way}
\end{equation*}\] where \(\alpha_i \sim
N(0, \sigma_\alpha^2)\) and \(\epsilon_{ij} \sim N(0, \sigma^2)\) and all
random variables are independent.
If the number \(I\) of groups is
small (in the example we consider below \(I =
4\)) then there is little information about the parameter \(\sigma_\alpha\) that summarizes the
variability in the mean of the response variable \(Y\) across groups. Unless strong prior
information about \(\sigma_\alpha\) is
provided the posterior distribution for \(\sigma_\alpha\) will tend to have a heavy
right tail (Gelman 2006). Indeed, for a
commonly-used `noninformative’ prior \(\pi(\mu, \sigma_\alpha, \log\sigma) \propto
1\) (Gelman et al. 2014) a proper
posterior density is not obtained if \(I <
3\) (Gelman 2006).
To illustrate the problem that a heavy-tailed density can cause a
naive implementation of the generalized ratio-of-uniforms algorithm we
use the blood coagulation time data presented in Table 11.2 of Gelman et al. (2014). These data are available
as the dataframe coagulation in the bang
package (Northrop and Hall 2017): the
first column holds the response values (the coagulation times of blood
drawn from 24 animals), the second column a factor indicating to which
of 4 diets the animal is allocated.
The bang function hanova1 samples from
the marginal posterior distribution of \((\mu,
\sigma_\alpha, \sigma)\) given data based on a user-supplied
prior distribution. The default prior is \(\pi(\mu, \sigma_\alpha, \log\sigma) \propto
1\). Under this prior, or indeed any prior in which \(\mu\) is normally distributed and
independent of \((\sigma_\alpha,
\sigma)\) a priori, the generalized ratio-of-uniforms
method can be used to sample from the marginal posterior distribution of
\((\sigma_\alpha, \sigma)\). By default
(argument param = "trans") then hanova1
parameterizes this marginal posterior in terms of \((\log \sigma_\alpha, \log \sigma)\). If
instead we use param = "original", so that this posterior
is parameterized in terms of \((\mu,
\sigma_\alpha, \sigma)\), then, with the default \(r = 1/2\), we find that the bounding box
cannot be found because the right tail of the posterior for \(\sigma_\alpha\) is heavy enough to prevent
this. However, if we use \(r = 1\) then
the bounding box can be found. The two successful approaches
(reparameterization or use of \(r=1\))
are illustrated below.
library(bang)
coag1 <- hanova1(resp = coagulation[, 1], fac = coagulation[, 2], n = 1000)
coag2 <- hanova1(resp = coagulation[, 1], fac = coagulation[, 2], n = 1000,
param = "original", r = 1)
probs <- c(2.5, 25, 50, 75, 97.5) / 100
all1 <- cbind(coag1$theta_sim_vals, coag1$sim_vals)
all2 <- cbind(coag2$theta_sim_vals, coag2$sim_vals)
round(t(apply(all1, 2, quantile, probs = probs)), 1)
#> 2.5% 25% 50% 75% 97.5%
#> theta[1] 58.8 60.4 61.2 62.0 63.7
#> theta[2] 63.9 65.2 65.8 66.5 67.8
#> theta[3] 65.7 67.1 67.8 68.4 69.8
#> theta[4] 59.5 60.6 61.2 61.7 62.9
#> mu 54.0 62.2 63.9 65.8 72.8
#> sigma[alpha] 2.0 3.6 5.3 8.0 26.3
#> sigma 1.8 2.2 2.4 2.7 3.4
round(t(apply(all2, 2, quantile, probs = probs)), 1)
#> 2.5% 25% 50% 75% 97.5%
#> theta[1] 58.7 60.4 61.3 62.1 63.8
#> theta[2] 63.9 65.3 65.9 66.6 67.9
#> theta[3] 65.6 67.1 67.8 68.4 69.7
#> theta[4] 59.6 60.6 61.1 61.7 62.8
#> mu 54.3 62.3 64.0 65.7 72.2
#> sigma[alpha] 1.9 3.5 5.0 7.6 26.6
#> sigma 1.8 2.2 2.4 2.7 3.5
These posterior summaries are similar to those presented in Table
11.3 of Gelman et al. (2014) (where \(\sigma_\alpha\) is denoted \(\tau\)), which were obtained using Gibbs
sampling.
The reparameterization strategy has the higher estimated probability
of acceptance.
coag1$pa
#> [1] 0.5089059
coag2$pa
#> [1] 0.3059039
