Sample, transform and evaluate from a joint marginal approximation

Sample, transform and evaluate from from a joint marginal approximation as returned using argument selection in inla.

inla.rjmarginal(n, jmarginal, constr)

inla.rjmarginal.eval(fun, samples, ...)

# S3 method for class 'inla.jmarginal'
print(x, ...)

# S3 method for class 'inla.jmarginal'
summary(object, ...)

# S3 method for class 'summary.inla.jmarginal'
print(x, ...)

inla.tjmarginal(jmarginal, A)

inla.1djmarginal(jmarginal)

Arguments

n

The number of samples

jmarginal

A marginal object given either by a inla object or result$selection

constr

Optional linear constraints; see ?INLA::f and argument extraconstr

fun

A function which is evaluated for each sample, similar to inla.posterior.sample.eval: please see the documentation for this functions for details.

samples

The samples, as in the form of the output from inla.rjmarginal

...

Arguments passed on to other methods (printing and summarising)

x

Object to be printed

object

Object to be summarised

A

A matrix used for the linear combination

Value

THESE FUNCTIONS ARE EXPERIMENTAL FOR THE MOMENT (JULY 2020)

inla.rjmarginal returns a list with the samples in samples (matrix) and the corresponding log-densities in log.density (vector). Each column in samples contains one sample.

inla.rjmarginal.eval returns a matrix, where each row is the (vector) function evaluated at each sample.

inla.tjmarginal returns a inla.jmarginal-object of the linear combination defined by the matrix A.

inla.1djmarginal return the marginal densities from a joint approximation.

See also

inla()

Author

Cristian Chiuchiolo and Havard Rue hrue@r-inla.org

Examples

 n = 10
 x = 1+rnorm(n)
 xx = 3 + rnorm(n)
 y = 1 + x + xx + rnorm(n)
 selection = list(xx=1, x=1)
 r = inla(y ~ 1 + x + xx,
          data = data.frame(y, x, xx),
          selection = selection)
 ns = 100
 xx = inla.rjmarginal(ns, r)

 print(cbind(mean = r$selection$mean, sample.mean = rowMeans(xx$samples)))
#>           mean sample.mean
#> xx:1 1.5895074   1.6096240
#> x:1  0.8244295   0.8222172
 print("cov matrix")
#> [1] "cov matrix"
 print(round(r$selection$cov.matrix, dig=3))
#>       [,1]  [,2]
#> [1,] 0.078 0.012
#> [2,] 0.012 0.045
 print("sample cov matrix")
#> [1] "sample cov matrix"
 print(round(cov(t(xx$samples)), dig=3))
#>       xx:1   x:1
#> xx:1 0.080 0.014
#> x:1  0.014 0.041

 skew = function(z) mean((z-mean(z))^3)/var(z)^1.5
 print(round(cbind(skew = r$selection$skewness,
                   sample.skew = apply(xx$samples, 1, skew)), digits = 3))
#>        skew sample.skew
#> xx:1 -0.001      -0.315
#> x:1   0.000      -0.251

 ## illustrating the eval function
 n = 10
 x = rnorm(n)
 eta = 1 + x
 y = eta + rnorm(n, sd=0.1)
 selection = list(x = 1, '(Intercept)' = 1)
 r = inla(y ~ 1 + x,
          data = data.frame(y, x),
          selection = selection)
 xx = inla.rjmarginal(100,  r)
 xx.eval = inla.rjmarginal.eval(function() c(x, Intercept),  xx)
#> Warning: Function 'inla.rjmarginal.eval()' is experimental.
 print(cbind(xx$samples[, 1]))
#>                    [,1]
#> x:1           1.0179514
#> (Intercept):1 0.9856987
 print(cbind(xx.eval[, 1]))
#>             [,1]
#> fun[1] 1.0179514
#> fun[2] 0.9856987

 constr <- list(A = matrix(1, ncol = nrow(xx$samples), nrow = 1), e = 1)
 x <- inla.rjmarginal(10, r, constr = constr)

 A <- matrix(rnorm(nrow(xx$samples)^2), nrow(xx$samples), nrow(xx$samples))
 b <- inla.tjmarginal(r, A)
 b.marg <- inla.1djmarginal(b)