qsample: Generate samples from a GMRF using the GMRFLib implementation

In andrewzm/INLA: Functions which allow to perform full Bayesian analysis of latent Gaussian models using Integrated Nested Laplace Approximaxion

Description

This function generate samples from a GMRF using the GMRFLib implementation

Usage

 inla.qsample(
        n = 1L,
        Q,
        b,
        mu, 
        sample,
        constr,
        reordering = inla.reorderings(),
        seed = 0L,
        logdens = ifelse(missing(sample), FALSE, TRUE))
 

Arguments

n	Number of samples. Only used if missing(sample)
Q	The precision matrix or a filename containing it.
b	The linear term
mu	The mu term
sample	A matrix of optional samples where each column is a sample. If set, then evaluate the log-density for each sample only.
constr	Optional linear constraints; see ?INLA::f and argument extraconstr
reordering	The type of reordering algorithm to be used; either one of the names listed in inla.reorderings() or the output from inla.qreordering(Q). The default is "auto" which try several reordering algorithm and use the best one for this particular matrix.
seed	Control the RNG. If seed=0L then GMRFLib will set the seed intelligently/at 'random'. If seed < 0L then the saved state of the RNG will be reused if possible, otherwise, GMRFLib will set the seed intelligently/at 'random'. If seed > 0L then this value is used as the seed for the RNG.
logdens	If TRUE, compute also the log-density of each sample. Note that the output format then change.

Value

The log-density has form -1/2(x-mu)^T Q (x-mu) + b^T x

If logdens is FALSE, then inla.qsample returns the samples in a matrix, where each column is a sample. If logdens is TRUE, then a list with names sample and logdens is returned. The samples are stored in the matrix sample where each column is a sample, and the log densities of each sample are stored the vector logdens.

Author(s)

Havard Rue [e-mail]

Examples

 g = system.file("demodata/germany.graph", package="INLA")
 Q = inla.graph2matrix(g)
 diag(Q) = dim(Q)[1L]
 x = inla.qsample(10, Q)
 ## Not run: matplot(x)
 x = inla.qsample(10, Q, logdens=TRUE)
 ## Not run: matplot(x$sample)

 n = 3
 Q = diag(n)
 ns = 2
 
 ## sample and evaluate a sample
 x = inla.qsample(n, Q=Q, logdens=TRUE)
 xx = inla.qsample(Q=Q,  sample = x$sample)
 print(x$logdens - xx$logdens)
 
 ## the use of a constraint
 constr = list(A = matrix(rep(1, n), 1, n), e = 0)
 x = inla.qsample(n, Q=Q, constr=constr)
 print(constr$A %*% x)
 
 ## control the RNG
 x = inla.qsample(n, Q=Q, seed = 123)
 ## restart from same seed,  only sample 1
 xx = inla.qsample(n=1, Q=Q, seed = 123)
 ## continue from the save state, sample the remaining 2
 xxx = inla.qsample(n=n-1, Q=Q, seed = -1)
 ## should be 0
 print(x - cbind(xx, xxx))

andrewzm/INLA documentation built on May 10, 2019, 11:12 a.m.