Generate samples from a GMRF using the GMRFLib implementation

This function generate samples from a GMRF using the GMRFLib implementation

inla.qsample(
        n = 1L,
        Q,
        b,
        mu, 
        sample,
        constr,
        reordering = INLA::inla.reorderings(),
        seed = 0L,
        logdens = ifelse(missing(sample), FALSE, TRUE),
        compute.mean = ifelse(missing(sample), FALSE, TRUE),
        num.threads = if (seed == 0L) 1L else NULL,
        selection = NULL, verbose = FALSE)

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 for TAUCS; 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.
compute.mean: If TRUE, compute also the (constrained) mean. Note that the output format then change.
num.threads: The number of threads that can be used. num.threads>1L requires seed = 0L.
selection: A vector of indices of each sample to return. NULL means return the whole sample. (Note that the log-density retured, is for the whole sample.) The use of selection cannot be combined with the use of sample.
verbose: Logical. Run in verbose mode or not.

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 or compute.mean is TRUE, then a list with names sample, logdens and mean 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 as the vector logdens. The mean (include corrections for the constraints, if any) is store in the vector mean.

Author

Havard Rue hrue@r-inla.org

Examples

 g = system.file("demodata/germany.graph", package="INLA")
 Q = inla.graph2matrix(g)
 diag(Q) = dim(Q)[1L]
 x = inla.qsample(10, Q)
#> Warning: error in running command
#> Warning: cannot open file '/tmp/RtmpB7Z7ym/file9c6e3cd9600c/file9c6e1be642f9': No such file or directory
#> Error in file(rng.file, "rb"): cannot open the connection
 if (FALSE) matplot(x) # \dontrun{}
 x = inla.qsample(10, Q, logdens=TRUE)
#> Warning: error in running command
#> Warning: cannot open file '/tmp/RtmpB7Z7ym/file9c6e3c63cfdd/file9c6ed85c95': No such file or directory
#> Error in file(rng.file, "rb"): cannot open the connection
 if (FALSE) matplot(x$sample) # \dontrun{}

 n = 3
 Q = diag(n)
 ns = 2
 
 ## sample and evaluate a sample
 x = inla.qsample(n, Q=Q, logdens=TRUE)
#> Warning: error in running command
#> Warning: cannot open file '/tmp/RtmpB7Z7ym/file9c6e29101545/file9c6e5046e884': No such file or directory
#> Error in file(rng.file, "rb"): cannot open the connection
 xx = inla.qsample(Q=Q,  sample = x$sample)
#> Error: object 'x' not found
 print(x$logdens - xx$logdens)
#> Error in h(simpleError(msg, call)): error in evaluating the argument 'x' in selecting a method for function 'print': object 'x' not found
 
 ## the use of a constraint
 constr = list(A = matrix(rep(1, n), 1, n), e = 0)
 x = inla.qsample(n, Q=Q, constr=constr)
#> Warning: error in running command
#> Warning: cannot open file '/tmp/RtmpB7Z7ym/file9c6e58a75c8c/file9c6ec3fbbcc': No such file or directory
#> Error in file(rng.file, "rb"): cannot open the connection
 print(constr$A %*% x)
#> Error in h(simpleError(msg, call)): error in evaluating the argument 'x' in selecting a method for function 'print': object 'x' not found
 
 ## control the RNG
 x = inla.qsample(n, Q=Q, seed = 123)
#> Warning: error in running command
#> Warning: cannot open file '/tmp/RtmpB7Z7ym/file9c6e20bf7d66/file9c6e37b0d397': No such file or directory
#> Error in file(rng.file, "rb"): cannot open the connection
 ## restart from same seed,  only sample 1
 xx = inla.qsample(n=1, Q=Q, seed = 123)
#> Warning: error in running command
#> Warning: cannot open file '/tmp/RtmpB7Z7ym/file9c6e5e875fbf/file9c6e15a01cda': No such file or directory
#> Error in file(rng.file, "rb"): cannot open the connection
 ## continue from the save state, sample the remaining 2
 xxx = inla.qsample(n=n-1, Q=Q, seed = -1)
#> Warning: error in running command
#> Warning: cannot open file '/tmp/RtmpB7Z7ym/file9c6e4c056740/file9c6e1eebced6': No such file or directory
#> Error in file(rng.file, "rb"): cannot open the connection
 ## should be 0
 print(x - cbind(xx, xxx))
#> Error in h(simpleError(msg, call)): error in evaluating the argument 'x' in selecting a method for function 'print': object 'x' not found