# NOT RUN {
set.seed(1)
par(mfrow = c(1, 1))
x = rnorm(1000)
xx = seq(-4, 4, 0.01)
y = dnorm(xx)
# Plenty of histogram bins (100)
breaks = seq(-4, 4, length.out=101)
# P-spline fitting with cubic B-splines, 2nd order penalty and 10 internal segments
# CV search for penalty coefficient.
fit = fpsdengpd(x, useq = seq(0, 3, 0.1), fixedu = TRUE,
lambdaseq = 10^seq(-5, 5, 0.25), breaks = breaks,
xrange = c(-4, 4), nseg = 10, degree = 3, ord = 2)
hist(x, freq = FALSE, breaks = breaks, xlim = c(-6, 6))
lines(xx, y, col = "black") # true density
# P-splines+GPD
with(fit, lines(xx, dpsdengpd(xx, beta, nbinwidth,
u = u, sigmau = sigmau, xi = xi, design = design.knots),
lwd = 2, col = "red"))
abline(v = fit$u, col = "red", lwd = 2, lty = 3)
# P-splines density estimate
with(fit, lines(xx, dpsden(xx, beta, nbinwidth, design = design.knots),
lwd = 2, col = "blue", lty = 2))
# vertical lines for all knots
with(fit, abline(v = design.knots, col = "red"))
# internal knots
with(fit, abline(v = design.knots[(degree + 2):(length(design.knots) - degree - 1)], col = "blue"))
# boundary knots (support of B-splines)
with(fit, abline(v = design.knots[c(degree + 1, length(design.knots) - degree)], col = "green"))
legend("topright", c("True Density","P-spline density","P-spline+GPD"),
col=c("black", "blue", "red"), lty = c(1, 2, 1))
legend("topleft", c("Internal Knots", "Boundaries", "Extra Knots", "Threshold"),
col=c("blue", "green", "red", "red"), lty = c(1, 1, 1, 2))
# }
# NOT RUN {
# }
Run the code above in your browser using DataLab