## Source:   Chapter10.r (with modifications) in 
## http://www.stat.colostate.edu/computationalstatistics/CompStatR.zip

############################################################################
### EXAMPLE 10.6 WHALE MIGRATION (LOGSPLINE)
############################################################################
# whalemigration.dat = contains data
# y = observed data
############################################################################
## NOTES
# The following example uses the {polspline} package. See the ?logspline 
# help page for more information. 
############################################################################
## INITIAL VALUES
whalemigration.dat = read.table("whalemigration.dat")
y = whalemigration.dat$V1

library(polspline) #you may need to install this package

## MAIN, Logspline Density Estimation
fit1 = oldlogspline(y)                    # 1992 version
fit1 = oldlogspline.to.logspline(fit1,y)  # 1992 to 1997 version
fit2 = oldlogspline(y,nknots=6,delete=F)
fit2 = oldlogspline.to.logspline(fit2,y)
fit3 = oldlogspline(y,nknots=11,delete=F)
fit3 = oldlogspline.to.logspline(fit3,y)

## OUTPUT
fit1      # OPTIMAL FIT BASED ON BIC

## PLOTS (reproduce Figure 10.7)
hist(y,breaks=20,freq=FALSE, ylim=c(0,0.015), xlim=c(940,1420), main="Logspline density estimate", xlab="hours since midnight, April 5")
plot(fit1,add=T,lwd=2)
plot(fit2,add=T,lty=2, lwd=2)
plot(fit3,add=T,lty=3, lwd=2)

points(fit2$knots,rep(0,6),pch=21,cex=2,bg="white")
points(fit1$knots,rep(0,7),pch=21,cex=2,bg="black")
legend(940,0.015,legend=c("density with 7 knots", "density with 6 knots", "density with 11 knots"), lty=c(1,2,3), lwd=c(2,2,2), cex=0.6)