"pdensity"<-
function(data, xmin = NULL, xmax = NULL, m = 10 * log(length(data)), ndx = 40, 
	q = 3, d = 3, lambda = 10)
{
#Function pdensity: density estimation with P-splines
#
#Input: data=vector of observations.
#Input: xmin=left boundry of domain.
#Input: xmax=right boundry of domain.
#Input: m=number of histogram bins.
#Input: ndx=number of intervals for B-splines.
#Input: q=degree of B-splines.
#Input: d=order of differences in the penalty.
#Input: lambda=smoothness parameter >=0.
#
#Result: a plot of the raw and the smoothed histogram.
#Output: AIC=deviance+2*trace(Hat).
#
#Reference: Eilers, P.H.C. and Marx, B.D. (1996). Flexible smoothing with B-splines and 
#           penalties (with comments and rejoinder). Statistical Science, 11(2): 89-121.
#
# (c) 1994 Paul Eilers 
#
	if(missing(xmin)) {
		xmin <- min(data)
	}
	if(missing(xmax)) {
		xmax <- max(data)
	}
	delta <- 0.005 * (xmax - xmin)
	x <- seq(xmin + delta, xmax - delta, length = m)
	d1 <- data[xmin < data]
	d2 <- d1[d1 <= xmax]
	y <- tabulate(cut(d2, x), m)
	dx <- (xmax - xmin)/ndx
	knots <- seq(xmin - q * dx, xmax + q * dx, by = dx)
	n <- ndx + q
	p <- diag(n)
	for(j in 1:d) {
		p <- diff(p)
	}
	p <- sqrt(lambda) * p
	b <- spline.des(knots, x, q + 1, 0 * x)$design
	nix <- rep(0, n - d)
	znew <- log(y + mean(y)/10)
	z <- 0 * znew
	it <- 0
	repeat {
		it <- it + 1
		if(it > 15)
			break
		dz <- max(abs(z - znew))
		if(dz < 1e-005)
			break
		print(c(it, dz))
		z <- znew
		mu <- exp(z)
		u <- (y - mu)/mu + z
		f <- lsfit(rbind(b, p), c(u, nix), wt = c(mu, (nix + 1)), 
			intercept = F)
		znew <- b %*% as.vector(f$coef)
	}
	plot(x, y, type = "s")
	lines(x + (0.5 * (xmax - xmin))/m, mu)
	dev <- 2 * sum(y[y > 0] * log(y[y > 0]/mu[y > 0]))
	h <- hat(f$qr)[1:m]
	aic <- dev + 2 * sum(h)
	aic
}