JOPS (version 0.1.19)

LAPS_dens: Bayesian density estimation

Description

Bayesian density estimation with P-splines and Laplace approximation.

Usage

LAPS_dens(B, P, y, loglambdas, tol = 1e-05, mon = FALSE)

Value

A list with elements:

alpha

P-spline coefficients of length n.

weights

weights from the Laplace approximation, which sum to 1 and are the same length as loglambdas.

mu

a vector of length m of expected values.

Cov

covariance matrix (m by m) of log(mu).

lambda

the penalty parameter.

ed

the effective model dimension.

Arguments

B

matrix (m by n) with B-spline basis, see bbase().

P

penalty matrix (n by n).

y

vector (length m) of counts, usually a histogram.

loglambdas

a vector of values of logarithms of lambda to explore.

tol

convergence tolerance (relative change in coefficients), default 1e-5.

mon

TRUE or FALSE to monitor the iteration history (default FALSE).

Author

Paul Eilers

Details

The B-spline basis should be based on the midpoints of the histogram bins. See the example below. This function is based on the paper of Gressani and Lambert (2018) and code input by Oswaldo Gressani.

References

Eilers, P.H.C. and Marx, B.D. (2021). Practical Smoothing, The Joys of P-splines. Cambridge University Press.

Gressani, O. and Lambert, P. (2018). Fast Bayesian inference using Laplace approximations in a flexible promotion time cure model based on P-splines. Computational Statistics and Data Analysis 124, 151-167.

Examples

# Smoothing a histogram of Old Faithful eruption durations
data(faithful)
durations = faithful[, 1]  # Eruption length

# Histogram with narrow bin widths
bw = 0.05
hst = hist(durations, breaks = seq(1, 6, by = bw), plot = TRUE)
x = hst$mids
y = hst$counts

# B-spline basis matrices, for fitting and plotting
nseg = 30
B = bbase(x, nseg = nseg)
xg = seq(min(x), max(x), by = 0.01)
Bg = bbase(xg, nseg = nseg)
n = ncol(B)

# Penalty matrix
D2 = diff(diag(n), diff = 2)
P2 = t(D2) %*% D2

# Fit the model
loglambs = seq(-1, 2, by = 0.05)
laps2 = LAPS_dens(B, P2, y, loglambs, mon = FALSE)
fhat2 = exp(Bg %*% laps2$alpha)
lines(xg, fhat2, col = "blue", lwd = 2)