qsample.RdThis 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 = NULL,
selection = NULL,
verbose = inla.getOption("verbose"),
.debug = FALSE
)Number of samples. Only used if missing(sample)
The precision matrix or a filename containing it.
The linear term
The mu term
A matrix of optional samples where each column is a sample. If set, then evaluate the log-density for each sample only.
Optional linear constraints; see ?INLA::f and argument
extraconstr
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.
Control the RNG. If seed=0L then GMRFLib will set the
seed intelligently/at 'random', and this is and should be the default
behaviour. 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.
PLEASE NOTE1: If seed!=0 then the computations will run in serial
mode, over-riding whatever is set in num.threads (a warning might be
issued).
PLEASE NOTE2: If the PARDISO sparse matrix library is used, continuity of the samples with respect to small changes in the precision matrix, can be expected but is not guaranteed. If this feature is required, please use the TAUCS sparse matrix library.
If TRUE, compute also the log-density of each sample.
Note that the output format then change.
If TRUE, compute also the (constrained) mean.
Note that the output format then change.
Maximum number of threads the inla-program will
use, or as 'A:B' defining the number threads in the outer (A) and inner (B)
layer for nested parallelism. seed!=0 requires serial comptuations.
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.
Logical. Run in verbose mode or not.
Logical. Internal debug-mode.
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.
g = system.file("demodata/germany.graph", package="INLA")
Q = inla.graph2matrix(g)
diag(Q) = dim(Q)[1L]
x = inla.qsample(10, Q)
if (FALSE) matplot(x) # \dontrun{}
x = inla.qsample(10, Q, logdens=TRUE)
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)
xx = inla.qsample(Q=Q, sample = x$sample)
print(x$logdens - xx$logdens)
#> logdens1 logdens2 logdens3
#> 0 0 0
## 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)
#> sample:1 sample:2 sample:3
#> [1,] 4.059321e-09 -3.429512e-09 -1.543883e-08
## control the RNG (require serial mode)
x = inla.qsample(n, Q=Q, seed = 123, num.threads="1:1")
## restart from same seed, only sample 1
xx = inla.qsample(n=1, Q=Q, seed = 123, num.threads="1:1")
## continue from the save state, sample the remaining 2
xxx = inla.qsample(n=n-1, Q=Q, seed = -1, num.threads="1:1")
## should be 0
print(x - cbind(xx, xxx))
#> sample:1 sample:2 sample:3
#> x:1 0 0 0
#> x:2 0 0 0
#> x:3 0 0 0