library(bootstrap)
str(cholost)
summary(cholost)
postscript("Boot-cholost-data.ps")
plot(y ~ z, data=cholost, xlab="Compliance", 
            ylab="Improvement",pch=19)
abline(lm(y ~ z, data=cholost), col="red")
dev.off()


n<-nrow(cholost)

q<-2
phi1<-2*pi/50
phi2<-2*pi/66


r<-order(x)
x<-cholost$z[r]
y<-cholost$y[r]


X1<-matrix(0,nrow=n,ncol=2*q+2)
X1[,1]<-1
X1[,2]<-x
X1[,3]<-cos(x*phi1)
X1[,4]<-sin(x*phi1)
X1[,5]<-cos(x*phi2)
X1[,6]<-sin(x*phi2)

beta.hat1<-solve(t(X1)%*%X1 ) %*% t(X1) %*% y

q<-3
phi1<-2*pi/50
phi2<-2*pi/66
phi3<-2*pi/75

X2<-matrix(0,nrow=n,ncol=2*q+2)
X2[,1]<-1
X2[,2]<-x
X2[,3]<-cos(x*phi1)
X2[,4]<-sin(x*phi1)
X2[,5]<-cos(x*phi2)
X2[,6]<-sin(x*phi2)
X2[,7]<-cos(x*phi3)
X2[,8]<-sin(x*phi3)


beta.hat2<-solve(t(X2)%*%X2 ) %*% t(X2) %*% y

xv<-seq(0,100,by=0.1)
Xv<-matrix(0,nrow=length(xv),ncol=2*q+2)
Xv[,1]<-1
Xv[,2]<-xv
Xv[,3]<-cos(xv*phi1)
Xv[,4]<-sin(xv*phi1)
Xv[,5]<-cos(xv*phi2)
Xv[,6]<-sin(xv*phi2)
Xv[,7]<-cos(xv*phi3)
Xv[,8]<-sin(xv*phi3)


X.fit1<-Xv[,1:6] %*% beta.hat1
X.fit2<-Xv %*% beta.hat2

postscript("Boot-cholost-fits123.ps")
plot(y ~ z, data=cholost, xlab="Compliance", 
            ylab="Improvement",pch=19)
abline(lm(y ~ z, data=cholost), col="red",lwd=2)

lines(xv,X.fit1,col="red",lty=2,lwd=2)
lines(xv,X.fit2,col="red",lty=3,lwd=2)
legend(0,100,c("q=1","q=2","q=3"),lty=c(1,2,3),col="red",lwd=2)
dev.off()


###################################################

q<-2
phi1<-2*pi/50
phi2<-2*pi/66

r<-order(x)
x<-cholost$z[r]
y<-cholost$y[r]
index<-c(1:n)

B<-200

xv<-seq(0,100,by=0.1)
Xv<-matrix(0,nrow=length(xv),ncol=2*q+2)
Xv[,1]<-1
Xv[,2]<-xv
Xv[,3]<-cos(xv*phi1)
Xv[,4]<-sin(xv*phi1)
Xv[,5]<-cos(xv*phi2)
Xv[,6]<-sin(xv*phi2)

postscript("Boot-cholost-boot1.ps")
plot(y ~ z, data=cholost, xlab="Compliance", 
            ylab="Improvement",pch=19)
dev.off()
postscript("Boot-cholost-boot2.ps")
plot(y ~ z, data=cholost, xlab="Compliance", 
            ylab="Improvement",pch=19)

for(b in 1:B){

	ind<-sort(sample(index,size=n,rep=T))
	xboot<-x[ind]
	yboot<-y[ind]

	Xboot<-matrix(0,nrow=n,ncol=2*q+2)
	Xboot[,1]<-1
	Xboot[,2]<-xboot
	Xboot[,3]<-cos(xboot*phi1)
	Xboot[,4]<-sin(xboot*phi1)
	Xboot[,5]<-cos(xboot*phi2)
	Xboot[,6]<-sin(xboot*phi2)

	beta.hat.boot<-solve(t(Xboot)%*%Xboot ) %*% t(Xboot) %*% yboot

	X.fitboot<-Xv%*% beta.hat.boot

	lines(xv,X.fitboot,col="red",lwd=1)


}
dev.off()

###################################################
postscript("Boot-cholost-boot3.ps")
plot(y ~ z, data=cholost, xlab="Compliance", 
            ylab="Improvement",pch=19)

for(b in 1:B){

	ind<-sort(sample(index,size=n,rep=T))
	xboot<-x[ind]
	yboot<-y[ind]

	X.fitboot<-predict(loess(yboot~xboot),xv)

	lines(xv,X.fitboot,col="red",lwd=1)


}
dev.off()


postscript("Boot-cholost-boot4.ps")
plot(y ~ z, data=cholost, xlab="Compliance", 
            ylab="Improvement",pch=19)

for(b in 1:B){

	ind<-sort(sample(index,size=n,rep=T))
	xboot<-x[ind]
	yboot<-y[ind]

	X.fitboot<-predict(loess(yboot~xboot,span=0.375),xv)

	lines(xv,X.fitboot,col="red",lwd=1)


}
dev.off()