Bayesian density estimation with P-splines and Laplace approximation.
LAPS_dens(B, P, y, loglambdas, tol = 1e-05, mon = FALSE)A list with elements:
P-spline coefficients of length n.
weights from the Laplace approximation, which sum to 1 and are
the same length as loglambdas.
a vector of length m of expected values.
covariance matrix (m by m) of log(mu).
the penalty parameter.
the effective model dimension.
matrix (m by n) with B-spline basis, see bbase().
penalty matrix (n by n).
vector (length m) of counts, usually a histogram.
a vector of values of logarithms of lambda to explore.
convergence tolerance (relative change in coefficients), default 1e-5.
TRUE or FALSE to monitor the iteration history (default FALSE).
Paul Eilers
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.
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.
# 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)
Run the code above in your browser using DataLab