Lab: Density Estimation, Splines and other Smoothers


Getting Started

To Start R: Click on the R icon


Example 1: Fitting a Univariate Density

Let's take a look at the Old Faithful geyser dataset:

> help(faithful)
> attach(faithful)

Let's make a histogram of the eruption times:

> hist(eruptions)

Let's computes kernel density estimates:

> help(density)

What are the choices of kernels and bandwidths?

Let's tryout the normal or Gaussian kernel with the default plug-in bandwidth:

> temp <- density(eruptions)

What's in temp? Can you plot it?

Let's try some other bandwiths with the same kernel. Look at the code in geyser.r

You can copy paste the lines of code or at the R prompt:
> source("http://www.rohan.sdsu.edu/~babailey/stat672/geyser.r")

What do you conclude?


Example 2: Other Kernels

> help(density)

Try some other kernels for example, the rectangular kernel:

> plot(density(eruptions, kernel="rectangular"))


Function Estimation:

Example 3: Smoothing Splines

> help(smooth.spline)

How is the the smoothing parameter lambda determined? (the default).

Let's look at the dataset mcycle:

The dataset is in the MASS package:

> help(mcycle)
> attach(mcycle)

If you can not get the data set:

> mcycle <- read.table("http://www.rohan.sdsu.edu/~babailey/stat672/mcycle.dat", header=T)
> attach(mcycle)

There is a help file available: mcycle.help

Let's fit a cubic smoothing spline:

> fit1 <- smooth.spline(times, accel)

What's in the fit object? Let's plot the data and overlay the fitted values:

> plot(times, accel, main = "data(mcycle) & smoothing spline")
> lines(fit1, col = "blue")

Can you predict the accel from times?


Example 4: Other Smoothers

Let's check out the kernel regression estimates:

> help(ksmooth)

Fit a normal kernel with with the bandwith=2:

> fit2 <- ksmooth(times, accel, kernel="normal", bandwidth=2)
> lines(fit2, col = "red")

Can you make a plot with different bandwidths=2, 5, and 10?

You could do polynomical regression,

> lines(times,fitted(lm(accel~poly(times,3))))
> lines(times,fitted(lm(accel~poly(times,6))))

There are other smoothers, you check them out!

> help(lowess)
> help(supsmu)
> help(ppr)


Source a file and creating graphics output:

Let's make geyser.r into program that creates a pdf file for output

help(pdf)