LAPS_dens: Bayesian density estimation

In JOPS: Practical Smoothing with P-Splines

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)

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).

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.

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.

Author(s)

Paul Eilers

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)

JOPS documentation built on Sept. 8, 2023, 5:42 p.m.