inla.sample {INLA} | R Documentation |
This function generate samples, and functions of those, from an approximated posterior of a fitted model (an inla-object)
inla.posterior.sample(n = 1L, result, selection = list(), intern = FALSE, use.improved.mean = TRUE, skew.corr = TRUE, add.names = TRUE, seed = 0L, num.threads = NULL, parallel.configs = TRUE, verbose=FALSE) inla.posterior.sample.eval(fun, samples, return.matrix = TRUE, ...)
n |
Number of samples. |
result |
The inla-object, ie the output from an |
selection |
Select what part of the sample to return. By default, the whole sample
is returned. |
intern |
Logical. If |
use.improved.mean |
Logical. If |
skew.corr |
Logical. If |
add.names |
Logical. If |
seed |
See the same argument in |
num.threads |
The number of threads to use in the format 'A:B' defining the number threads in the
outer (A) and inner (B) layer for nested parallelism. A '0' will be replaced
intelligently.
|
parallel.configs |
Logical. If |
verbose |
Logical. Run in verbose mode or not. |
fun |
The function to evaluate for each sample. Upon entry, the variable names
defined in the model are defined as the value of the sample.
The list of names are defined in |
samples |
|
return.matrix |
Logical. If |
... |
Additional arguments to |
The hyperparameters are sampled from the configurations used to do the
numerical integration, hence if you want a higher resolution, you need to
to change the int.stratey
variable and friends. The latent field is
sampled from the Gaussian approximation conditioned on the hyperparameters,
but with a correction for the mean (default),
and optional (and by default) corrected for the estimated skewness.
The log.density report is only correct when there is no constraints. With constraints, it correct the Gaussian part of the sample for the constraints.
After the sample is (optional) skewness corrected, the log.density is is not exact for correcting for constraints, but the error is very small in most cases.
inla.posterior.sample
returns a list of the samples,
where each sample is a list with
names hyperpar
and latent
, and with their marginal
densities in logdens$hyperpar
and logdens$latent
and the joint density is in logdens$joint
.
inla.posterior.sample.eval
return a list or a matrix of
fun
applied to each sample.
Havard Rue hrue@r-inla.org and Cristian Chiuchiolo cristian.chiuchiolo@kaust.edu.sa
r = inla(y ~ 1 ,data = data.frame(y=rnorm(1)), control.compute = list(config=TRUE)) samples = inla.posterior.sample(2,r) ## reproducible results: inla.seed = as.integer(runif(1)*.Machine$integer.max) set.seed(12345) x = inla.posterior.sample(10, r, seed = inla.seed, num.threads="1:1") set.seed(12345) xx = inla.posterior.sample(10, r, seed = inla.seed, num.threads="1.1") all.equal(x, xx) set.seed(1234) n = 25 xx = rnorm(n) yy = rev(xx) z = runif(n) y = rnorm(n) r = inla(y ~ 1 + z + f(xx) + f(yy, copy="xx"), data = data.frame(y, z, xx, yy), control.compute = list(config=TRUE), family = "gaussian") r.samples = inla.posterior.sample(10, r) fun = function(...) { mean(xx) - mean(yy) } f1 = inla.posterior.sample.eval(fun, r.samples) fun = function(...) { c(exp(Intercept), exp(Intercept + z)) } f2 = inla.posterior.sample.eval(fun, r.samples) fun = function(...) { return (theta[1]/(theta[1] + theta[2])) } f3 = inla.posterior.sample.eval(fun, r.samples) ## Predicting nz new observations, and ## comparing the estimated one with the true one set.seed(1234) n = 100 alpha = beta = s = 1 z = rnorm(n) y = alpha + beta * z + rnorm(n, sd = s) r = inla(y ~ 1 + z, data = data.frame(y, z), control.compute = list(config=TRUE), family = "gaussian") r.samples = inla.posterior.sample(10^3, r) ## just return samples of the intercept intercepts = inla.posterior.sample.eval("Intercept", r.samples) nz = 3 znew = rnorm(nz) fun = function(zz = NA) { ## theta[1] is the precision return (Intercept + z * zz + rnorm(length(zz), sd = sqrt(1/theta[1]))) } par(mfrow=c(1, nz)) f1 = inla.posterior.sample.eval(fun, r.samples, zz = znew) for(i in 1:nz) { hist(f1[i, ], n = 100, prob = TRUE) m = alpha + beta * znew[i] xx = seq(m-4*s, m+4*s, by = s/100) lines(xx, dnorm(xx, mean=m, sd = s), lwd=2) } ## ## Be aware that using non-clean variable names might be a little tricky ## n <- 100 X <- matrix(rnorm(n^2), n, 2) x <- X[, 1] xx <- X[, 2] xxx <- x*xx y <- 1 + 2*x + 3*xx + 4*xxx + rnorm(n, sd = 0.01) r <- inla(y ~ X[, 1]*X[, 2], data = list(y = y, X = X), control.compute = list(config = TRUE)) print(round(dig = 4, r$summary.fixed[,"mean"])) sam <- inla.posterior.sample(100, r) sam.extract <- inla.posterior.sample.eval( (function(...) { beta.1 <- get("X[, 1]") beta.2 <- get("X[, 2]") beta.12 <- get("X[, 1]:X[, 2]") return(c(Intercept, beta.1, beta.2, beta.12)) }), sam) print(round(dig = 4, rowMeans(sam.extract))) ## a simpler form can also be used here, and in the examples below sam.extract <- inla.posterior.sample.eval( c("Intercept", "X[, 1]", "X[, 2]", "X[, 1]:X[, 2]"), sam) print(round(dig = 4, rowMeans(sam.extract))) r <- inla(y ~ x + xx + xxx, data = list(y = y, x = x, xx = xx, xxx = xxx), control.compute = list(config = TRUE)) sam <- inla.posterior.sample(100, r) sam.extract <- inla.posterior.sample.eval( (function(...) { return(c(Intercept, x, xx, xxx)) }), sam) print(round(dig = 4, rowMeans(sam.extract))) sam.extract <- inla.posterior.sample.eval(c("Intercept", "x", "xx", "xxx"), sam) print(round(dig = 4, rowMeans(sam.extract))) r <- inla(y ~ x*xx, data = list(y = y, x = x, xx = xx), control.compute = list(config = TRUE)) sam <- inla.posterior.sample(100, r) sam.extract <- inla.posterior.sample.eval( (function(...) { return(c(Intercept, x, xx, get("x:xx"))) }), sam) print(round(dig = 4, rowMeans(sam.extract))) sam.extract <- inla.posterior.sample.eval(c("Intercept", "x", "xx", "x:xx"), sam) print(round(dig = 4, rowMeans(sam.extract)))