library(splines) # loads the library of functions to compute B-splines

x <- c(1:1000)/1000
X1 <- 1
for (k in 1:5) X1 <- cbind(X1,cos(2*k*pi*x),sin(2*k*pi*x))
X2 <- cbind(1,bs(x,df=10))

plot(x, X2[,2])
for(i in 3:11) points(x, X2[,i])

leverage <- function(x,y,w,r=10,m=100) {
  qrx <- qr(x)
  qry <- qr(y)
  n <- nrow(x)
  lev <- NULL
  for (i in 1:m) {
    vx <- ifelse(runif(n)>0.5,1,-1)
    vx[-w] <- 0
    v0x <- qr.fitted(qrx,vx)
    v0y <- qr.fitted(qry,vx)
    fx <- v0x
    for (j in 2:r) {
      v0x[-w] <- 0
      v0x <- qr.fitted(qrx,v0x)
      fx <- fx + v0x/j
    }
    levx <- c(lev,sum(vx*fx))
    fy <- v0y
    for (j in 2:r) {
      v0y[-w] <- 0
      v0y <- qr.fitted(qry,v0y)
      fx <- fy + v0y/j
    }
    levy <- c(lev,sum(vx*fy))
  }
  std.errx <- exp(-mean(levx))*sd(levx)/sqrt(m)
  levx <- 1 - exp(-mean(levx))
  std.erry <- exp(-mean(levy))*sd(levy)/sqrt(m)
  levy <- 1 - exp(-mean(levy))
  r <- list(levx=levx,std.errx=std.errx,levy=levy,std.erry=std.erry)
  r
}

lev <- leverage(X2, X1, 50)
lev

for (i in 1:1000) {
  lev <- leverage(X2, i)
  print(lev)
}