Data <- read.delim("../../data/bmd.txt")
MData <- Data[Data$gender=="male",]
FData <- Data[Data$gender=="female",]

## Slide 4
xx <- seq(min(Data$age),max(Data$age),len=101)
f1 <- f2 <- numeric(101)
lam <- 0.5
for (i in 1:101)
  {
    Ind <- abs(FData$age-xx[i]) < lam
    f1[i] <- sum(FData$spnbmd*Ind)/sum(Ind)
    Ind <- abs(MData$age-xx[i]) < lam
    f2[i] <- sum(MData$spnbmd*Ind)/sum(Ind)
  }
par(mfrow=c(1,2))
plot(FData$age,FData$spnbmd,xlab="Age",ylab="Relative change in spinal BMD",pch=19,cex=0.7,col="gray",main="Females")
lines(xx,f1,type="s",lwd=2,col="red")
plot(MData$age,MData$spnbmd,xlab="Age",ylab="Relative change in spinal BMD",pch=19,cex=0.7,col="gray",main="Males")
lines(xx,f2,type="s",lwd=2,col="red")

## Slide 8
x <- c(  1,  3,  5,  6,6.2,6.6,  7,  8, 10)
y <- c(1.1,0.9,1.2,2.8,2.9,2.7,1.5,1.7,1.2)
n <- length(x)
xx <- seq(min(x),max(x),len=101)
plot(x,y,col="gray",pch=19,cex=0.7,bty="n",ylim=c(0,max(y)))
for (i in 1:n)
  {
    lines(xx,10*dnorm(xx,x[i],0.75)/10,col="blue")
  }
f <- numeric(length(xx))
for (i in 1:length(xx))
  {
    K <- dnorm(x,mean=xx[i],sd=.75)
    f[i] <- sum(y*K)/sum(K)
  }
lines(xx,f,lwd=2)

## Slide 10
xx <- seq(-3,3,len=101)
y1 <- dnorm(xx,sd=1)
y2 <- 3/4*(1-xx^2)*(abs(xx) <= 1)
y3 <- (1 - abs(xx)^3)^3 * (abs(xx) <= 1)##/1.157
df <- data.frame(x=rep(xx,3),Density=c(y1,y2,y3),kernel=c(rep("Gaussian",101),rep("Epanechnikov",101),rep("Tri-cube",101)))
tp1 <- xyplot(Density~x,df,group=kernel,auto.key=list(points=FALSE,lines=TRUE,columns=3),type="l",ylab="K(x)")
y1 <- dnorm(xx,sd=0.5)
y3 <- (1 - abs(xx)^3)^3 * (abs(xx) <= 1)/1.157
df <- data.frame(x=rep(xx,3),Density=c(y1,y2,y3),kernel=c(rep("Gaussian",101),rep("Epanechnikov",101),rep("Tri-cube",101)))
tp2 <- xyplot(Density~x,df,group=kernel,auto.key=list(points=FALSE,lines=TRUE,columns=3),type="l",ylab="K(x)")
trellis.par.set(superpose.line=list(col=c("red","green","blue"),lwd=3))
print(tp1,split=c(1,1,2,1),more=TRUE)
print(tp2,split=c(2,1,2,1))

## Slide 11
ep <- function(x){3/4*(abs(x)<1)*(1-x^2)}
xx <- seq(min(Data$age),max(Data$age),len=101)
f1 <- f2 <- numeric(101)
lam <- c(0.25,2,10)
par(mfcol=c(2,3))
for (j in 1:3)
  {
    for (i in 1:101)
      {
        K <- ep((FData$age-xx[i])/lam[j])
        f1[i] <- sum(FData$spnbmd*K)/sum(K)
        K <- ep((MData$age-xx[i])/lam[j])
        f2[i] <- sum(MData$spnbmd*K)/sum(K)
      }
    plot(FData$age,FData$spnbmd,xlab="Age",ylab="Relative change in spinal BMD",pch=19,cex=0.7,col="gray",main="Females")
    lines(xx,f1,type="l",lwd=2,col="red")
    plot(MData$age,MData$spnbmd,xlab="Age",ylab="Relative change in spinal BMD",pch=19,cex=0.7,col="gray",main="Males")
    lines(xx,f2,type="l",lwd=2,col="red")
  }

## Slide 12
xx <- seq(min(Data$age),max(Data$age),len=101)
f1 <- f2 <- numeric(101)
lam <- c(0.25,2,10)/2
par(mfcol=c(2,3))
for (j in 1:3)
  {
    for (i in 1:101)
      {
        K <- dnorm((FData$age-xx[i])/lam[j])
        f1[i] <- sum(FData$spnbmd*K)/sum(K)
        K <- dnorm((MData$age-xx[i])/lam[j])
        f2[i] <- sum(MData$spnbmd*K)/sum(K)
      }
    plot(FData$age,FData$spnbmd,xlab="Age",ylab="Relative change in spinal BMD",pch=19,cex=0.7,col="gray",main="Females")
    lines(xx,f1,type="l",lwd=2,col="red")
    plot(MData$age,MData$spnbmd,xlab="Age",ylab="Relative change in spinal BMD",pch=19,cex=0.7,col="gray",main="Males")
    lines(xx,f2,type="l",lwd=2,col="red")
  }