histosmooth {EpiLPS}

Histogram smoothing with Laplacian-P-splines

Description

This function provides a smooth density estimate to a histogram using Laplacian-P-splines. The B-spline basis is computed on the midpoints of the histogram bins. The default number of (cubic) B-splines is 30 and a third-order penalty is specified. The negative binomial distribution is used to model the number of observations falling in each bin.

Usage

histosmooth(x, xl = min(x), xr = max(x), K = 30)

Arguments

x

A vector of real numbers from which the histogram will be constructed.

xl

The left bound for the domain of x which also coincides with the left bound of the domain of the B-spline basis. The default is taken to be the minimum of x.

xr

The right bound for the domain of x which also coincides with the right bound of the domain of the B-spline basis. The default is taken to be the maximum of x.

K

Number of B-splines in the basis.

Value

A list containing the left (xl) and right (xr) bounds of the domain of the estimated density, the binwidth and a function to be evaluated between xl and xr.

Author(s)

Oswaldo Gressani [e-mail]

References

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

Gressani, O., Wallinga, J., Althaus, C. L., Hens, N. and Faes, C. (2022). EpiLPS: A fast and flexible Bayesian tool for estimation of the time-varying reproduction number. Plos Computational Biology, 18(10): e1010618.

Examples

# Old Faithful geyser application
data(eruptions)
x <- eruptions
ffit <- histosmooth(x, xl = 1, xr = 6)
tt <- seq(ffit$xl, ffit$xr, length = 500)
dtt <- tt[2] - tt[1]
graphics::hist(x, breaks = seq(ffit$xl, ffit$xr, by = ffit$binwidth),
     freq = FALSE, ylim = c(0, 0.8), main = "Old Faithful Geyser",
    xlab = "Eruption time (minutes)")
densfit <- sapply(tt, ffit$fdens)
densfit <- densfit / (sum(densfit * dtt))
graphics::lines(tt, densfit, col = "red", lwd = 2)