R/marginal.R
In andrewzm/INLA: Functions which allow to perform full Bayesian analysis of latent Gaussian models using Integrated Nested Laplace Approximaxion
## Export: inla.dmarginal inla.pmarginal inla.qmarginal inla.rmarginal inla.hpdmarginal
## Export: inla.smarginal inla.emarginal inla.tmarginal inla.mmarginal inla.zmarginal
## NOT-YET-IN-USE-Export: summary!inla.marginal plot!inla.marginal
##! \name{marginal}
##! \alias{inla.marginal}
##! \alias{marginal}
##! \alias{pmarginal}
##! \alias{inla.pmarginal}
##! \alias{qmarginal}
##! \alias{inla.qmarginal}
##! \alias{dmarginal}
##! \alias{inla.dmarginal}
##! \alias{rmarginal}
##! \alias{inla.rmarginal}
##! \alias{inla.hpdmarginal}
##! \alias{hpdmarginal}
##! \alias{inla.expectation}
##! \alias{inla.emarginal}
##! \alias{emarginal}
##! \alias{inla.marginal.expectation}
##! \alias{marginal.expectation}
##! \alias{inla.spline}
##! \alias{inla.smarginal}
##! \alias{smarginal}
##! \alias{inla.marginal.transform}
##! \alias{marginal.transform}
##! \alias{inla.tmarginal}
##! \alias{inla.mmarginal}
##! \alias{mmarginal}
##! \alias{inla.zmarginal}
##! \alias{zmarginal}
##!
##! \title{Functions which operates on marginals}
##!
##! \description{Density, distribution function, quantile function, random
##! generation, hpd-interval, interpolation, expectations, mode and transformations of
##! marginals obtained by \code{inla} or \code{inla.hyperpar()}.
##! These functions computes the density (inla.dmarginal),
##! the distribution function (inla.pmarginal),
##! the quantile function (inla.qmarginal),
##! random generation (inla.rmarginal),
##! spline smoothing (inla.smarginal),
##! computes expected values (inla.emarginal),
##! computes the mode (inla.mmarginal),
##! transforms the marginal (inla.tmarginal), and provide summary statistics (inla.zmarginal).
##! }
##!
##! \usage{
##! inla.dmarginal(x, marginal, log = FALSE)
##! inla.pmarginal(q, marginal, normalize = TRUE, len = 1024)
##! inla.qmarginal(p, marginal, len = 1024)
##! inla.rmarginal(n, marginal)
##! inla.hpdmarginal(p, marginal, len = 1024)
##! inla.smarginal(marginal, log = FALSE, extrapolate = 0.0, keep.type = FALSE, factor=10L)
##! inla.emarginal(fun, marginal, ...)
##! inla.mmarginal(marginal)
##! inla.tmarginal(fun, marginal, n=1024, h.diff = .Machine$double.eps^(1/3),
##! method = c("quantile", "linear"), ...)
##! inla.zmarginal(marginal, silent = FALSE)
##! }
##! \arguments{
##!
##! \item{marginal}{A marginal object from either \code{inla} or
##! \code{inla.hyperpar()}, which is either \code{list(x=c(), y=c())}
##! with density values \code{y} at locations \code{x}, or a
##! \code{matrix(,n,2)} for which the density values are the second
##! column and the locations in the first column.
##! The\code{inla.hpdmarginal()}-function
##! assumes a unimodal density.}
##!
##! \item{fun}{A (vectorised) function like \code{function(x) exp(x)} to
##! compute the expectation against, or which define the transformation
##! new = fun(old)}
##!
##! \item{x}{Evaluation points}
##!
##! \item{q}{Quantiles}
##!
##! \item{p}{Probabilities}
##!
##! \item{n}{The number of observations. If \code{length(n) > 1}, the
##! length is taken to be the number required. For
##! \code{inla.marginal.transform}, its the number of points to use
##! in the new density.}
##!
##! \item{h.diff}{The step-length for the numerical differeniation inside \code{inla.marginal.transform}}
##!
##! \item{...}{Further arguments to be passed to function which
##! expectation is to be computed.}
##!
##! \item{log}{Return density or interpolated density in log-scale?}
##!
##! \item{normalize}{Renormalise the density after interpolation?}
##! \item{len}{Number of locations used to interpolate the distribution
##! function.}
##!
##! \item{keep.type}{If \code{FALSE} then return a \code{list(x=, y=)}, otherwise if \code{TRUE},
##! then return a matrix if the input is a matrix}
##!
##! \item{extrapolate}{How much to extrapolate on each side when computing the
##! interpolation. In fraction of the range.}
##!
##! \item{factor}{The number of points after interpolation is \code{factor} times the original number of points;
##! which is argument \code{n} in \code{spline}}
##!
##! \item{method}{Which method should be used to layout points for where the transformation is computed.}
##!
##! \item{silent}{Output the result visually (TRUE) or just through the call.}
##! }
##!
##! \value{%%
##! \code{inla.smarginal} returns \code{list=c(x=c(), y=c())} of
##! interpolated values do extrapolation using the factor given,
##! and the remaining function returns what they say they should do. }
##! %%
##!
##! \author{Havard Rue \email{hrue@math.ntnu.no}}
##!
##! \seealso{\code{\link{inla}}, \code{\link{inla.hyperpar}}}
##!
##! \examples{
##! ## a simple linear regression example
##! n = 10
##! x = rnorm(n)
##! sd = 0.1
##! y = 1+x + rnorm(n,sd=sd)
##! res = inla(y ~ 1 + x, data = data.frame(x,y),
##! control.family=list(initial = log(1/sd^2),fixed=TRUE))
##!
##! ## chose a marginal and compare the with the results computed by the
##! ## inla-program
##! r = res$summary.fixed["x",]
##! m = res$marginals.fixed$x
##!
##! ## compute the 95% HPD interval
##! inla.hpdmarginal(0.95, m)
##!
##! x = seq(-6, 6, len = 1000)
##! y = dnorm(x)
##! inla.hpdmarginal(0.95, list(x=x, y=y))
##!
##! ## compute the the density for exp(r), version 1
##! r.exp = inla.tmarginal(exp, m)
##! ## or version 2
##! r.exp = inla.tmarginal(function(x) exp(x), m)
##!
##! ## to plot the marginal, we use the inla.smarginal, which interpolates (in
##! ## log-scale). Compare with some samples.
##! plot(inla.smarginal(m), type="l")
##! s = inla.rmarginal(1000, m)
##! hist(inla.rmarginal(1000, m), add=TRUE, prob=TRUE)
##! lines(density(s), lty=2)
##!
##! m1 = inla.emarginal(function(x) x^1, m)
##! m2 = inla.emarginal(function(x) x^2, m)
##! stdev = sqrt(m2 - m1^2)
##! q = inla.qmarginal(c(0.025,0.975), m)
##!
##! ## inla-program results
##! print(r)
##!
##! ## inla.marginal-results (they shouldn't be perfect!)
##! print(c(mean=m1, sd=stdev, "0.025quant" = q[1], "0.975quant" = q[2]))
##! ## using the buildt-in function
##! inla.zmarginal(m)
##! }
##!
### functions to work with the marginal, either defined as a matrix
### x[2, n], or a list(x=, y=). NOTE:: there are NO EXTRAPOLATION, so
### the marginal is defined FROM...TO!
`inla.marginal.fix` = function(marginal)
{
## just remove points where the density is <= 0 or where the
## density is to small compared to the maximum density. (othewise
## we can get trouble with the spline interpolation). same with
## 'x'? No...
eps = sqrt(.Machine$double.eps)
if (is.matrix(marginal)) {
i = (marginal[, 2] > 0) & (abs(marginal[, 2]/max(marginal[, 2])) > eps)
m = list(x=marginal[i, 1], y=marginal[i, 2])
## i = c(diff(marginal[, 1]) > eps, TRUE)
## m = list(x=marginal[i, 1], y=marginal[i, 2])
} else {
i = (marginal$y > 0) & (abs(marginal$y/max(marginal$y)) > eps)
m = list(x = marginal$x[i], y = marginal$y[i])
## i = c(diff(marginal$x) > eps, TRUE)
## m = list(x = marginal$x[i], y = marginal$y[i])
}
return (m)
}
`inla.spline` = function(marginal, log = FALSE, extrapolate = 0.0, factor = 10L) {
return (inla.smarginal(marginal, log, extrapolate, factor))
}
`inla.smarginal` = function(marginal, log = FALSE, extrapolate = 0.0, keep.type = FALSE, factor=10L)
{
## for marginal in matrix MARGINAL, which is a marginal density,
## return the nice interpolated (x, y) where the interpolation is
## done in log(y)
is.mat = is.matrix(marginal)
m = inla.marginal.fix(marginal)
r = diff(range(m$x))
ans = spline(m$x, log(m$y), xmin = min(m$x) - extrapolate * r, xmax = max(m$x) + extrapolate * r,
n = factor*length(m$x), method = "natural")
if (!log) {
ans$y = exp(ans$y)
}
if (is.mat && keep.type) {
return (cbind(ans$x, ans$y))
} else {
return (ans)
}
}
`inla.sfmarginal` = function(marginal)
{
## for marginal in matrix MARGINAL, which is a marginal density,
## return the spline-function which returns the log-density
m = inla.marginal.fix(marginal)
r = range(m$x)
return (list(range = r, fun = splinefun(m$x, log(m$y))))
}
`inla.expectation` = function(fun, marginal, ...) {
return (inla.emarginal(fun, marginal, ...))
}
`inla.emarginal` = function(fun, marginal, ...)
{
## compute E(FUN(x)), where the marginal of x is given in
## `marginal'; see inla.smarginal(). Also work for FUN(x)
## returning a vector.
xx = inla.smarginal(marginal)
n = length(xx$x)
if (n%%2 == 0)
n = n -1
## use Simpsons integration rule
i.0 = c(1, n)
i.4 = seq(2, n-1, by=2)
i.2 = seq(3, n-2, by=2)
fun = match.fun(fun)
ff = fun(xx$x[1:n], ...) * xx$y[1:n]
nf = length(ff) %/% n
e = numeric(nf)
off = 0L
for(i in 1:nf) {
e[i] = sum(sum(ff[i.0 + off]) + 4*sum(ff[i.4 + off]) + 2*sum(ff[i.2 + off]))
off = off + n
}
## normalise, so that E(1) = 1
ff = 1 * xx$y[1:n]
e.1 = sum(sum(ff[i.0]) + 4*sum(ff[i.4]) + 2*sum(ff[i.2]))
return (e/e.1)
}
`inla.dmarginal` = function(x, marginal, log = FALSE)
{
## return the density of the marginal. if LOG, in log-scale. this
## density is not renormalized (but should be already normalized
## in some scale).
## return list(range=, fun=..)
f = inla.sfmarginal(inla.smarginal(marginal))
n = length(x)
d = numeric(n)
for(i in 1:n) {
if (x[i] >= f$range[1] && x[i] <= f$range[2]) {
if (log)
d[i] = f$fun(x[i])
else
d[i] = exp(f$fun(x[i]))
} else {
if (log)
d[i] = -Inf
else
d[i] = 0.0
}
}
return (d)
}
`inla.pmarginal` = function(q, marginal, normalize = TRUE, len = 1024)
{
f = inla.sfmarginal(inla.smarginal(marginal))
xx = seq(f$range[1], f$range[2], length = len)
d = cumsum(exp(f$fun(xx)))
d = d/d[length(d)]
## just spline-interpolate the mapping
fq = splinefun(xx, d, method = "monoH.FC")
## just make sure the p's are in [0, 1]
n = length(q)
xx = pmin(pmax(q, rep(f$range[1], n)), rep(f$range[2], n))
return (fq(xx))
}
`inla.qmarginal` = function(p, marginal, len = 1024)
{
f = inla.sfmarginal(inla.smarginal(marginal))
xx = seq(f$range[1], f$range[2], length = len)
d = cumsum(exp(f$fun(xx)))
d = d/d[length(d)]
## for the moment, we remove only duplicated zero's and one's.
eps = .Machine$double.eps * 1000.0
for(val in c(0.0, 1.0)) {
is.val = which(abs(d - val) <= eps)
if (length(is.val) > 1) {
## keep the first, ie remove it from the list to the removed
is.val = is.val[-1]
## and remove the rest
d = d[-is.val]
xx = xx[-is.val]
}
}
## just spline-interpolate the inverse mapping
fq = splinefun(d, xx, method = "monoH.FC")
## just make sure the p's are in [0, 1]
n = length(p)
pp = pmin(pmax(p, rep(0, n)), rep(1, n))
return (fq(pp))
}
`inla.hpdmarginal` = function(p, marginal, len = 2048L)
{
sm = inla.smarginal(marginal, keep.type = FALSE)
f = inla.sfmarginal(sm)
x.range = f$range
xx = seq(x.range[1], x.range[2], length = len)
## simply add 0 density outside the range.
dx = diff(xx)[1]
ef = exp(f$fun(xx))
xx = c(min(xx) - dx, xx, max(xx) + dx)
d = c(0, ef, 0)
d = cumsum(d)
d = d/d[length(d)]
## for the moment, we remove only duplicated zero's and one's.
eps = .Machine$double.eps * 1000.0
for(val in c(0.0, 1.0)) {
is.val = which(abs(d - val) <= eps)
if (length(is.val) > 1) {
## keep the first, ie remove it from the list to the removed
is.val = is.val[-1]
## and remove the rest
d = d[-is.val]
xx = xx[-is.val]
}
}
## just spline-interpolate the inverse mapping
fq = splinefun(d, xx, method = "monoH.FC")
## just make sure the p's are in [0, 1]
np = length(p)
pp = 1-pmin(pmax(p, rep(0, np)), rep(1, np))
## parts of this code below is taken from TeachingDemo::hpd
f = function(x, posterior.icdf, conf) {
return (posterior.icdf(1 - conf + x) - posterior.icdf(x))
}
tol= sqrt(.Machine$double.eps)
result = matrix(NA, np, 2)
for(i in 1:np) {
out = optimize(f, c(0, pp[i]), posterior.icdf = fq, conf = pp[i], tol = tol)
result[i, ] = c(fq(out$minimum), fq(1 - pp[i] + out$minimum))
## make sure its in the range of the marginal itself as we expanded the domain
result[i, ] = c(min(x.range[2], max(result[i, 1], x.range[1])),
min(x.range[2], max(result[i, 2], x.range[1])))
}
colnames(result) = c("low", "high")
rownames(result) = paste("level:", format(1-pp, digits=6, justify="left", trim=TRUE), sep="")
return (result)
}
`inla.rmarginal` = function(n, marginal)
{
return (inla.qmarginal(runif(n), marginal))
}
`inla.marginal.transform` = function(fun, marginal, n=1024, h.diff = .Machine$double.eps^(1/3),
method = c("quantile", "linear"), ...)
{
return (inla.tmarginal(fun, marginal, n, h.diff, method = method, ...))
}
`inla.tmarginal` = function(fun, marginal, n=1024, h.diff = .Machine$double.eps^(1/3),
method = c("quantile", "linear"), ...)
{
f = match.fun(fun)
ff = function(x) f(x, ...)
is.mat = is.matrix(marginal)
m = inla.smarginal(marginal)
r = range(m$x)
method = match.arg(method)
if (inla.strcasecmp(method, "quantile")) {
x = inla.qmarginal((1:n)/(n+1), marginal)
} else if (inla.strcasecmp(method, "linear")) {
x = seq(r[1], r[2], length = n)
} else {
stop("unknown method")
}
xx = ff(x)
## use the numDeriv library if present
present = inla.require("numDeriv")
if (present) {
## using numDeriv
dif = numeric(length(x))
for(i in 1:length(x)) {
dif[i] = numDeriv::grad(ff, x[i])
}
log.dens = inla.dmarginal(x, marginal, log=FALSE)/ abs(dif)
} else {
## use a simple algorithm
log.dens = inla.dmarginal(x, marginal, log=FALSE)/
abs((ff(x + h.diff) - ff(x-h.diff))/(2*h.diff))
}
## if decreasing function, then reverse
if (xx[1] > xx[n]) {
xx[1:n] = xx[n:1]
log.dens[1:n] = log.dens[n:1]
}
if (is.mat) {
ret = cbind(x = xx, y = log.dens)
} else {
ret = list(x = xx, y = log.dens)
}
if (FALSE) {
class(ret) = "inla.marginal"
attr(ret, "inla.tag") = paste(attr(marginal, "inla.tag"), "transformed")
}
return (ret)
}
`inla.mmarginal` = function(marginal)
{
p = seq(0.01, 0.99, by=0.01)
n = length(p)
x = inla.qmarginal(p, marginal)
d = inla.dmarginal(x, marginal, log=TRUE)
idx = which.max(d)
res = optimize(
f = inla.dmarginal,
interval = c(x[max(1L, idx-1L)], x[min(n, idx+1L)]),
maximum = TRUE,
## arguments to inla.dmarginal
marginal = marginal, log=TRUE)
return(res$maximum)
}
`inla.zmarginal` = function(marginal, silent = FALSE)
{
stopifnot(inla.is.marginal(marginal))
m = inla.emarginal(function(xx) c(xx, xx^2), marginal)
q = inla.qmarginal(c(0.025, 0.25, 0.5, 0.75, 0.975), marginal)
s = sqrt(max(0, m[2]-m[1]^2))
if (!silent) {
##cat(as.character(match.call()[-1]), "\n")
##cat("+--------+-----+------------\n")
cat("Mean ", format(m[1], digits=6), "\n")
cat("Stdev ", format(s, digits=6), "\n")
cat("Quantile 0.025", format(q[1], digits=6), "\n")
cat("Quantile 0.25 ", format(q[2], digits=6), "\n")
cat("Quantile 0.5 ", format(q[3], digits=6), "\n")
cat("Quantile 0.75 ", format(q[4], digits=6), "\n")
cat("Quantile 0.975", format(q[5], digits=6), "\n")
}
return (invisible(list(mean = m[1], sd = s, "quant0.025" = q[1],
"quant0.25" = q[2], "quant0.5" = q[3], "quant0.75"=q[4], "quant0.975" = q[5])))
}
`inla.is.marginal` = function(m)
{
return ((is.matrix(m) && ncol(m) == 2L && nrow(m) > 2L) ||
(is.list(m) && all(names(m) == c("x", "y"))))
}
## 'plot' and 'summary'-methods for marginals
plot.inla.marginal = function(x, y, ...)
{
m = inla.emarginal(function(xx) c(xx, xx^2), x)
xlab = paste("x (mean", format(m[1], digits=4),
", stdev", format(sqrt(max(0, m[2]-m[1]^2)), digits=4), ")")
plot(inla.smarginal(x, extrapolate=0), type = "l", xlab = xlab, ylab = "marginal density", ...)
return (invisible())
}
summary.inla.marginal = function(object, ...)
{
## option
silent=FALSE
inla.eval.dots(...)
m = inla.emarginal(function(xx) c(xx, xx^2), object)
q = inla.qmarginal(c(0.025, 0.25, 0.5, 0.75, 0.975), object)
s = sqrt(max(0, m[2]-m[1]^2))
if (!silent) {
cat("Properties of: ", "FIXME", "\n")
cat("+--------------+------------------\n")
cat("Mean ", format(m[1], digits=6), "\n")
cat("Stdev ", format(s, digits=6), "\n")
cat("Quantile 0.025", format(q[1], digits=6), "\n")
cat("Quantile 0.25 ", format(q[2], digits=6), "\n")
cat("Quantile 0.5 ", format(q[3], digits=6), "\n")
cat("Quantile 0.75 ", format(q[4], digits=6), "\n")
cat("Quantile 0.975", format(q[5], digits=6), "\n")
}
return (invisible(list(mean = m[1], sd = s, "quant0.025" = q[1],
"quant0.25" = q[2], "quant0.5" = q[3], "quant0.75"=q[4], "quant0.975" = q[5])))
}
`inla.plot.inla.marginals` = function(x, ...)
{
stop("NOT IN USE FOR THE MOMENT")
## input here is a list of marginals
n = length(x)
if (n == 1) {
## length is 1, so here we plot the marginal itself
plot(x[[1]], ...)
} else {
xx = sapply(x, summary, silent=TRUE)
stopifnot(dim(xx)[1] == 7)
stopifnot(dim(xx)[2] == n)
xx.min = min(as.numeric(xx))
xx.max = max(as.numeric(xx))
ylim = range(pretty(c(xx.min, xx.max)))
lab = c("mean", "quant0.5", "quant0.025", "quant0.975", "quant0.25", "quant0.75")
xval = 1:n
plot(xval, rep(ylim[2]*1000, n),
ylim = ylim, main = "FIXME",
xlab = "x", ylab = inla.paste(lab, sep=", "), ...)
polygon(c(xval, rev(xval)), c(xx[lab[3], ], rev(xx[lab[4], ])), col="lightgray", border=NA, ...)
polygon(c(xval, rev(xval)), c(xx[lab[5], ], rev(xx[lab[6], ])), col="gray", border=NA, ...)
## mean
j=1; lines(xval, xx[lab[j], ], lwd=1, lty=j, ...)
j=1; points(xval, xx[lab[j], ], pch=20, cex=0.4)
## median
j=2; lines(xval, xx[lab[j], ], lwd=1, lty=j, ...)
return (invisible())
}
}
andrewzm/INLA documentation built on May 10, 2019, 11:12 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
## Export: inla.dmarginal inla.pmarginal inla.qmarginal inla.rmarginal inla.hpdmarginal
## Export: inla.smarginal inla.emarginal inla.tmarginal inla.mmarginal inla.zmarginal
## NOT-YET-IN-USE-Export: summary!inla.marginal plot!inla.marginal
##! \name{marginal}
##! \alias{inla.marginal}
##! \alias{marginal}
##! \alias{pmarginal}
##! \alias{inla.pmarginal}
##! \alias{qmarginal}
##! \alias{inla.qmarginal}
##! \alias{dmarginal}
##! \alias{inla.dmarginal}
##! \alias{rmarginal}
##! \alias{inla.rmarginal}
##! \alias{inla.hpdmarginal}
##! \alias{hpdmarginal}
##! \alias{inla.expectation}
##! \alias{inla.emarginal}
##! \alias{emarginal}
##! \alias{inla.marginal.expectation}
##! \alias{marginal.expectation}
##! \alias{inla.spline}
##! \alias{inla.smarginal}
##! \alias{smarginal}
##! \alias{inla.marginal.transform}
##! \alias{marginal.transform}
##! \alias{inla.tmarginal}
##! \alias{inla.mmarginal}
##! \alias{mmarginal}
##! \alias{inla.zmarginal}
##! \alias{zmarginal}
##!
##! \title{Functions which operates on marginals}
##!
##! \description{Density, distribution function, quantile function, random
##! generation, hpd-interval, interpolation, expectations, mode and transformations of
##! marginals obtained by \code{inla} or \code{inla.hyperpar()}.
##! These functions computes the density (inla.dmarginal),
##! the distribution function (inla.pmarginal),
##! the quantile function (inla.qmarginal),
##! random generation (inla.rmarginal),
##! spline smoothing (inla.smarginal),
##! computes expected values (inla.emarginal),
##! computes the mode (inla.mmarginal),
##! transforms the marginal (inla.tmarginal), and provide summary statistics (inla.zmarginal).
##! }
##!
##! \usage{
##! inla.dmarginal(x, marginal, log = FALSE)
##! inla.pmarginal(q, marginal, normalize = TRUE, len = 1024)
##! inla.qmarginal(p, marginal, len = 1024)
##! inla.rmarginal(n, marginal)
##! inla.hpdmarginal(p, marginal, len = 1024)
##! inla.smarginal(marginal, log = FALSE, extrapolate = 0.0, keep.type = FALSE, factor=10L)
##! inla.emarginal(fun, marginal, ...)
##! inla.mmarginal(marginal)
##! inla.tmarginal(fun, marginal, n=1024, h.diff = .Machine$double.eps^(1/3),
##! method = c("quantile", "linear"), ...)
##! inla.zmarginal(marginal, silent = FALSE)
##! }
##! \arguments{
##!
##! \item{marginal}{A marginal object from either \code{inla} or
##! \code{inla.hyperpar()}, which is either \code{list(x=c(), y=c())}
##! with density values \code{y} at locations \code{x}, or a
##! \code{matrix(,n,2)} for which the density values are the second
##! column and the locations in the first column.
##! The\code{inla.hpdmarginal()}-function
##! assumes a unimodal density.}
##!
##! \item{fun}{A (vectorised) function like \code{function(x) exp(x)} to
##! compute the expectation against, or which define the transformation
##! new = fun(old)}
##!
##! \item{x}{Evaluation points}
##!
##! \item{q}{Quantiles}
##!
##! \item{p}{Probabilities}
##!
##! \item{n}{The number of observations. If \code{length(n) > 1}, the
##! length is taken to be the number required. For
##! \code{inla.marginal.transform}, its the number of points to use
##! in the new density.}
##!
##! \item{h.diff}{The step-length for the numerical differeniation inside \code{inla.marginal.transform}}
##!
##! \item{...}{Further arguments to be passed to function which
##! expectation is to be computed.}
##!
##! \item{log}{Return density or interpolated density in log-scale?}
##!
##! \item{normalize}{Renormalise the density after interpolation?}
##! \item{len}{Number of locations used to interpolate the distribution
##! function.}
##!
##! \item{keep.type}{If \code{FALSE} then return a \code{list(x=, y=)}, otherwise if \code{TRUE},
##! then return a matrix if the input is a matrix}
##!
##! \item{extrapolate}{How much to extrapolate on each side when computing the
##! interpolation. In fraction of the range.}
##!
##! \item{factor}{The number of points after interpolation is \code{factor} times the original number of points;
##! which is argument \code{n} in \code{spline}}
##!
##! \item{method}{Which method should be used to layout points for where the transformation is computed.}
##!
##! \item{silent}{Output the result visually (TRUE) or just through the call.}
##! }
##!
##! \value{%%
##! \code{inla.smarginal} returns \code{list=c(x=c(), y=c())} of
##! interpolated values do extrapolation using the factor given,
##! and the remaining function returns what they say they should do. }
##! %%
##!
##! \author{Havard Rue \email{hrue@math.ntnu.no}}
##!
##! \seealso{\code{\link{inla}}, \code{\link{inla.hyperpar}}}
##!
##! \examples{
##! ## a simple linear regression example
##! n = 10
##! x = rnorm(n)
##! sd = 0.1
##! y = 1+x + rnorm(n,sd=sd)
##! res = inla(y ~ 1 + x, data = data.frame(x,y),
##! control.family=list(initial = log(1/sd^2),fixed=TRUE))
##!
##! ## chose a marginal and compare the with the results computed by the
##! ## inla-program
##! r = res$summary.fixed["x",]
##! m = res$marginals.fixed$x
##!
##! ## compute the 95% HPD interval
##! inla.hpdmarginal(0.95, m)
##!
##! x = seq(-6, 6, len = 1000)
##! y = dnorm(x)
##! inla.hpdmarginal(0.95, list(x=x, y=y))
##!
##! ## compute the the density for exp(r), version 1
##! r.exp = inla.tmarginal(exp, m)
##! ## or version 2
##! r.exp = inla.tmarginal(function(x) exp(x), m)
##!
##! ## to plot the marginal, we use the inla.smarginal, which interpolates (in
##! ## log-scale). Compare with some samples.
##! plot(inla.smarginal(m), type="l")
##! s = inla.rmarginal(1000, m)
##! hist(inla.rmarginal(1000, m), add=TRUE, prob=TRUE)
##! lines(density(s), lty=2)
##!
##! m1 = inla.emarginal(function(x) x^1, m)
##! m2 = inla.emarginal(function(x) x^2, m)
##! stdev = sqrt(m2 - m1^2)
##! q = inla.qmarginal(c(0.025,0.975), m)
##!
##! ## inla-program results
##! print(r)
##!
##! ## inla.marginal-results (they shouldn't be perfect!)
##! print(c(mean=m1, sd=stdev, "0.025quant" = q[1], "0.975quant" = q[2]))
##! ## using the buildt-in function
##! inla.zmarginal(m)
##! }
##!
### functions to work with the marginal, either defined as a matrix
### x[2, n], or a list(x=, y=). NOTE:: there are NO EXTRAPOLATION, so
### the marginal is defined FROM...TO!
`inla.marginal.fix` = function(marginal)
{
## just remove points where the density is <= 0 or where the
## density is to small compared to the maximum density. (othewise
## we can get trouble with the spline interpolation). same with
## 'x'? No...
eps = sqrt(.Machine$double.eps)
if (is.matrix(marginal)) {
i = (marginal[, 2] > 0) & (abs(marginal[, 2]/max(marginal[, 2])) > eps)
m = list(x=marginal[i, 1], y=marginal[i, 2])
## i = c(diff(marginal[, 1]) > eps, TRUE)
## m = list(x=marginal[i, 1], y=marginal[i, 2])
} else {
i = (marginal$y > 0) & (abs(marginal$y/max(marginal$y)) > eps)
m = list(x = marginal$x[i], y = marginal$y[i])
## i = c(diff(marginal$x) > eps, TRUE)
## m = list(x = marginal$x[i], y = marginal$y[i])
}
return (m)
}
`inla.spline` = function(marginal, log = FALSE, extrapolate = 0.0, factor = 10L) {
return (inla.smarginal(marginal, log, extrapolate, factor))
}
`inla.smarginal` = function(marginal, log = FALSE, extrapolate = 0.0, keep.type = FALSE, factor=10L)
{
## for marginal in matrix MARGINAL, which is a marginal density,
## return the nice interpolated (x, y) where the interpolation is
## done in log(y)
is.mat = is.matrix(marginal)
m = inla.marginal.fix(marginal)
r = diff(range(m$x))
ans = spline(m$x, log(m$y), xmin = min(m$x) - extrapolate * r, xmax = max(m$x) + extrapolate * r,
n = factor*length(m$x), method = "natural")
if (!log) {
ans$y = exp(ans$y)
}
if (is.mat && keep.type) {
return (cbind(ans$x, ans$y))
} else {
return (ans)
}
}
`inla.sfmarginal` = function(marginal)
{
## for marginal in matrix MARGINAL, which is a marginal density,
## return the spline-function which returns the log-density
m = inla.marginal.fix(marginal)
r = range(m$x)
return (list(range = r, fun = splinefun(m$x, log(m$y))))
}
`inla.expectation` = function(fun, marginal, ...) {
return (inla.emarginal(fun, marginal, ...))
}
`inla.emarginal` = function(fun, marginal, ...)
{
## compute E(FUN(x)), where the marginal of x is given in
## `marginal'; see inla.smarginal(). Also work for FUN(x)
## returning a vector.
xx = inla.smarginal(marginal)
n = length(xx$x)
if (n%%2 == 0)
n = n -1
## use Simpsons integration rule
i.0 = c(1, n)
i.4 = seq(2, n-1, by=2)
i.2 = seq(3, n-2, by=2)
fun = match.fun(fun)
ff = fun(xx$x[1:n], ...) * xx$y[1:n]
nf = length(ff) %/% n
e = numeric(nf)
off = 0L
for(i in 1:nf) {
e[i] = sum(sum(ff[i.0 + off]) + 4*sum(ff[i.4 + off]) + 2*sum(ff[i.2 + off]))
off = off + n
}
## normalise, so that E(1) = 1
ff = 1 * xx$y[1:n]
e.1 = sum(sum(ff[i.0]) + 4*sum(ff[i.4]) + 2*sum(ff[i.2]))
return (e/e.1)
}
`inla.dmarginal` = function(x, marginal, log = FALSE)
{
## return the density of the marginal. if LOG, in log-scale. this
## density is not renormalized (but should be already normalized
## in some scale).
## return list(range=, fun=..)
f = inla.sfmarginal(inla.smarginal(marginal))
n = length(x)
d = numeric(n)
for(i in 1:n) {
if (x[i] >= f$range[1] && x[i] <= f$range[2]) {
if (log)
d[i] = f$fun(x[i])
else
d[i] = exp(f$fun(x[i]))
} else {
if (log)
d[i] = -Inf
else
d[i] = 0.0
}
}
return (d)
}
`inla.pmarginal` = function(q, marginal, normalize = TRUE, len = 1024)
{
f = inla.sfmarginal(inla.smarginal(marginal))
xx = seq(f$range[1], f$range[2], length = len)
d = cumsum(exp(f$fun(xx)))
d = d/d[length(d)]
## just spline-interpolate the mapping
fq = splinefun(xx, d, method = "monoH.FC")
## just make sure the p's are in [0, 1]
n = length(q)
xx = pmin(pmax(q, rep(f$range[1], n)), rep(f$range[2], n))
return (fq(xx))
}
`inla.qmarginal` = function(p, marginal, len = 1024)
{
f = inla.sfmarginal(inla.smarginal(marginal))
xx = seq(f$range[1], f$range[2], length = len)
d = cumsum(exp(f$fun(xx)))
d = d/d[length(d)]
## for the moment, we remove only duplicated zero's and one's.
eps = .Machine$double.eps * 1000.0
for(val in c(0.0, 1.0)) {
is.val = which(abs(d - val) <= eps)
if (length(is.val) > 1) {
## keep the first, ie remove it from the list to the removed
is.val = is.val[-1]
## and remove the rest
d = d[-is.val]
xx = xx[-is.val]
}
}
## just spline-interpolate the inverse mapping
fq = splinefun(d, xx, method = "monoH.FC")
## just make sure the p's are in [0, 1]
n = length(p)
pp = pmin(pmax(p, rep(0, n)), rep(1, n))
return (fq(pp))
}
`inla.hpdmarginal` = function(p, marginal, len = 2048L)
{
sm = inla.smarginal(marginal, keep.type = FALSE)
f = inla.sfmarginal(sm)
x.range = f$range
xx = seq(x.range[1], x.range[2], length = len)
## simply add 0 density outside the range.
dx = diff(xx)[1]
ef = exp(f$fun(xx))
xx = c(min(xx) - dx, xx, max(xx) + dx)
d = c(0, ef, 0)
d = cumsum(d)
d = d/d[length(d)]
## for the moment, we remove only duplicated zero's and one's.
eps = .Machine$double.eps * 1000.0
for(val in c(0.0, 1.0)) {
is.val = which(abs(d - val) <= eps)
if (length(is.val) > 1) {
## keep the first, ie remove it from the list to the removed
is.val = is.val[-1]
## and remove the rest
d = d[-is.val]
xx = xx[-is.val]
}
}
## just spline-interpolate the inverse mapping
fq = splinefun(d, xx, method = "monoH.FC")
## just make sure the p's are in [0, 1]
np = length(p)
pp = 1-pmin(pmax(p, rep(0, np)), rep(1, np))
## parts of this code below is taken from TeachingDemo::hpd
f = function(x, posterior.icdf, conf) {
return (posterior.icdf(1 - conf + x) - posterior.icdf(x))
}
tol= sqrt(.Machine$double.eps)
result = matrix(NA, np, 2)
for(i in 1:np) {
out = optimize(f, c(0, pp[i]), posterior.icdf = fq, conf = pp[i], tol = tol)
result[i, ] = c(fq(out$minimum), fq(1 - pp[i] + out$minimum))
## make sure its in the range of the marginal itself as we expanded the domain
result[i, ] = c(min(x.range[2], max(result[i, 1], x.range[1])),
min(x.range[2], max(result[i, 2], x.range[1])))
}
colnames(result) = c("low", "high")
rownames(result) = paste("level:", format(1-pp, digits=6, justify="left", trim=TRUE), sep="")
return (result)
}
`inla.rmarginal` = function(n, marginal)
{
return (inla.qmarginal(runif(n), marginal))
}
`inla.marginal.transform` = function(fun, marginal, n=1024, h.diff = .Machine$double.eps^(1/3),
method = c("quantile", "linear"), ...)
{
return (inla.tmarginal(fun, marginal, n, h.diff, method = method, ...))
}
`inla.tmarginal` = function(fun, marginal, n=1024, h.diff = .Machine$double.eps^(1/3),
method = c("quantile", "linear"), ...)
{
f = match.fun(fun)
ff = function(x) f(x, ...)
is.mat = is.matrix(marginal)
m = inla.smarginal(marginal)
r = range(m$x)
method = match.arg(method)
if (inla.strcasecmp(method, "quantile")) {
x = inla.qmarginal((1:n)/(n+1), marginal)
} else if (inla.strcasecmp(method, "linear")) {
x = seq(r[1], r[2], length = n)
} else {
stop("unknown method")
}
xx = ff(x)
## use the numDeriv library if present
present = inla.require("numDeriv")
if (present) {
## using numDeriv
dif = numeric(length(x))
for(i in 1:length(x)) {
dif[i] = numDeriv::grad(ff, x[i])
}
log.dens = inla.dmarginal(x, marginal, log=FALSE)/ abs(dif)
} else {
## use a simple algorithm
log.dens = inla.dmarginal(x, marginal, log=FALSE)/
abs((ff(x + h.diff) - ff(x-h.diff))/(2*h.diff))
}
## if decreasing function, then reverse
if (xx[1] > xx[n]) {
xx[1:n] = xx[n:1]
log.dens[1:n] = log.dens[n:1]
}
if (is.mat) {
ret = cbind(x = xx, y = log.dens)
} else {
ret = list(x = xx, y = log.dens)
}
if (FALSE) {
class(ret) = "inla.marginal"
attr(ret, "inla.tag") = paste(attr(marginal, "inla.tag"), "transformed")
}
return (ret)
}
`inla.mmarginal` = function(marginal)
{
p = seq(0.01, 0.99, by=0.01)
n = length(p)
x = inla.qmarginal(p, marginal)
d = inla.dmarginal(x, marginal, log=TRUE)
idx = which.max(d)
res = optimize(
f = inla.dmarginal,
interval = c(x[max(1L, idx-1L)], x[min(n, idx+1L)]),
maximum = TRUE,
## arguments to inla.dmarginal
marginal = marginal, log=TRUE)
return(res$maximum)
}
`inla.zmarginal` = function(marginal, silent = FALSE)
{
stopifnot(inla.is.marginal(marginal))
m = inla.emarginal(function(xx) c(xx, xx^2), marginal)
q = inla.qmarginal(c(0.025, 0.25, 0.5, 0.75, 0.975), marginal)
s = sqrt(max(0, m[2]-m[1]^2))
if (!silent) {
##cat(as.character(match.call()[-1]), "\n")
##cat("+--------+-----+------------\n")
cat("Mean ", format(m[1], digits=6), "\n")
cat("Stdev ", format(s, digits=6), "\n")
cat("Quantile 0.025", format(q[1], digits=6), "\n")
cat("Quantile 0.25 ", format(q[2], digits=6), "\n")
cat("Quantile 0.5 ", format(q[3], digits=6), "\n")
cat("Quantile 0.75 ", format(q[4], digits=6), "\n")
cat("Quantile 0.975", format(q[5], digits=6), "\n")
}
return (invisible(list(mean = m[1], sd = s, "quant0.025" = q[1],
"quant0.25" = q[2], "quant0.5" = q[3], "quant0.75"=q[4], "quant0.975" = q[5])))
}
`inla.is.marginal` = function(m)
{
return ((is.matrix(m) && ncol(m) == 2L && nrow(m) > 2L) ||
(is.list(m) && all(names(m) == c("x", "y"))))
}
## 'plot' and 'summary'-methods for marginals
plot.inla.marginal = function(x, y, ...)
{
m = inla.emarginal(function(xx) c(xx, xx^2), x)
xlab = paste("x (mean", format(m[1], digits=4),
", stdev", format(sqrt(max(0, m[2]-m[1]^2)), digits=4), ")")
plot(inla.smarginal(x, extrapolate=0), type = "l", xlab = xlab, ylab = "marginal density", ...)
return (invisible())
}
summary.inla.marginal = function(object, ...)
{
## option
silent=FALSE
inla.eval.dots(...)
m = inla.emarginal(function(xx) c(xx, xx^2), object)
q = inla.qmarginal(c(0.025, 0.25, 0.5, 0.75, 0.975), object)
s = sqrt(max(0, m[2]-m[1]^2))
if (!silent) {
cat("Properties of: ", "FIXME", "\n")
cat("+--------------+------------------\n")
cat("Mean ", format(m[1], digits=6), "\n")
cat("Stdev ", format(s, digits=6), "\n")
cat("Quantile 0.025", format(q[1], digits=6), "\n")
cat("Quantile 0.25 ", format(q[2], digits=6), "\n")
cat("Quantile 0.5 ", format(q[3], digits=6), "\n")
cat("Quantile 0.75 ", format(q[4], digits=6), "\n")
cat("Quantile 0.975", format(q[5], digits=6), "\n")
}
return (invisible(list(mean = m[1], sd = s, "quant0.025" = q[1],
"quant0.25" = q[2], "quant0.5" = q[3], "quant0.75"=q[4], "quant0.975" = q[5])))
}
`inla.plot.inla.marginals` = function(x, ...)
{
stop("NOT IN USE FOR THE MOMENT")
## input here is a list of marginals
n = length(x)
if (n == 1) {
## length is 1, so here we plot the marginal itself
plot(x[[1]], ...)
} else {
xx = sapply(x, summary, silent=TRUE)
stopifnot(dim(xx)[1] == 7)
stopifnot(dim(xx)[2] == n)
xx.min = min(as.numeric(xx))
xx.max = max(as.numeric(xx))
ylim = range(pretty(c(xx.min, xx.max)))
lab = c("mean", "quant0.5", "quant0.025", "quant0.975", "quant0.25", "quant0.75")
xval = 1:n
plot(xval, rep(ylim[2]*1000, n),
ylim = ylim, main = "FIXME",
xlab = "x", ylab = inla.paste(lab, sep=", "), ...)
polygon(c(xval, rev(xval)), c(xx[lab[3], ], rev(xx[lab[4], ])), col="lightgray", border=NA, ...)
polygon(c(xval, rev(xval)), c(xx[lab[5], ], rev(xx[lab[6], ])), col="gray", border=NA, ...)
## mean
j=1; lines(xval, xx[lab[j], ], lwd=1, lty=j, ...)
j=1; points(xval, xx[lab[j], ], pch=20, cex=0.4)
## median
j=2; lines(xval, xx[lab[j], ], lwd=1, lty=j, ...)
return (invisible())
}
}
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.