A Markov chain Monte Carlo (MCMC) sampler for location-scale regression
models from the lmls() function. The sampler uses Gibbs updates for the
location coefficients and the Riemann manifold Metropolis-adjusted Langevin
algorithm (MMALA) from Girolami and Calderhead (2011) with the Fisher-Rao
metric tensor for the scale coefficients. The priors for the regression
coefficients are assumed to be flat.
To find the optimal step size for the MMALA updates, the dual averaging algorithm from Nesterov (2009) is used during a warm-up phase.
Arguments
- m
A location-scale regression model from the
lmls()function.- num_samples
The number of MCMC samples after the warm-up. Defaults to 1000.
- num_warmup
The number of MCMC samples for the warm-up. Defaults to 1000.
- target_accept
The target acceptance rate for the dual averaging algorithm used for the warm-up. Defaults to 0.8.
Value
An lmls S3 object, see lmls(). The entry mcmc with the matrices
of MCMC samples is added to the object as a list with the names location
and scale.
References
Girolami, M. and Calderhead, B. (2011), Riemann manifold Langevin and Hamiltonian Monte Carlo methods. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 73: 123-214. doi:10.1111/j.1467-9868.2010.00765.x
Nesterov, Y. (2009), Primal-dual subgradient methods for convex problems. Mathematical Programming, 120: 221–259. doi:10.1007/s10107-007-0149-x
Examples
library(lmls)
m <- lmls(y ~ poly(x, 2), ~ x, data = abdom, light = FALSE)
m <- mcmc(m)
summary(m, type = "mcmc")
#>
#> Call:
#> lmls(location = y ~ poly(x, 2), scale = ~x, data = abdom, light = FALSE)
#>
#> Deviance residuals:
#> Min. 1st Qu. Median Mean 3rd Qu. Max.
#> -3.48600 -0.76150 -0.05559 -0.01225 0.62090 4.20700
#>
#> Location coefficients (identity link):
#> Mean 2.5% 50% 97.5%
#> (Intercept) 226.74 225.64 226.76 227.79
#> poly(x, 2)1 2160.52 2130.65 2160.82 2189.83
#> poly(x, 2)2 -98.48 -122.70 -98.78 -73.89
#>
#> Scale coefficients (log link):
#> Mean 2.5% 50% 97.5%
#> (Intercept) 1.35344 1.17007 1.35456 1.538
#> x 0.04255 0.03633 0.04243 0.049
#>
#> Residual degrees of freedom: 605
#> Log-likelihood: -2397.38
#> AIC: 4804.76
#> BIC: 4826.83
#>
plot(m$mcmc$scale[, 2], type = "l")
