# NOT RUN {
## Old Faithful eruptions data histogram and clsd density
library(MASS)
data(faithful)
attach(faithful)
model <- clsd(eruptions)
ylim <- c(0,max(model$density,hist(eruptions,breaks=20,plot=FALSE)$density))
plot(model,ylim=ylim)
hist(eruptions,breaks=20,freq=FALSE,add=TRUE,lty=2)
rug(eruptions)
summary(model)
coef(model)
## Simulated data
set.seed(42)
require(logspline)
## Example - simulated data
n <- 250
x <- sort(rnorm(n))
f.dgp <- dnorm(x)
model <- clsd(x)
## Standard (cubic) estimate taken from the logspline package
## Compute MSEs
mse.clsd <- mean((fitted(model)-f.dgp)^2)
model.logspline <- logspline(x)
mse.logspline <- mean((dlogspline(x,model.logspline)-f.dgp)^2)
ylim <- c(0,max(fitted(model),dlogspline(x,model.logspline),f.dgp))
plot(model,
ylim=ylim,
sub=paste("MSE: logspline = ",format(mse.logspline),", clsd = ",
format(mse.clsd)),
lty=3,
col=3)
xer <- model$xer
lines(xer,dlogspline(xer,model.logspline),col=2,lty=2)
lines(xer,dnorm(xer),col=1,lty=1)
rug(x)
legend("topright",c("DGP",
paste("Cubic Logspline Density (package `logspline', knots = ",
model.logspline$nknots,")",sep=""),
paste("clsd Density (degree = ", model$degree, ", segments = ",
model$segments,", penalty = ",round(model$penalty,2),")",sep="")),
lty=1:3,
col=1:3,
bty="n",
cex=0.75)
summary(model)
coef(model)
## Simulate data with known bounds
set.seed(42)
n <- 10000
x <- runif(n,0,1)
model <- clsd(x,lbound=0,ubound=1)
plot(model)
# }
# NOT RUN {
# }
# NOT RUN {
<!-- %% End dontrun -->
# }
Run the code above in your browser using DataCamp Workspace