# LAPS_dens: Bayesian density estimation In JOPS: Practical Smoothing with P-Splines

 LAPS_dens R Documentation

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

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.