require(fields) #provides sreg = Smoothing spline regression

# define data : either simulated or real data
experiment <- "simulated";
# experiment <- "realData";
if(experiment == "simulated"){ # simulated data
  f <- function(t) {t*(1-t)*sin(2*pi*t)}
  s<-0.08; n<-50;
  t<-runif(n); t<-sort(t);
  theta<-f(t);
  y<-theta+s*rnorm(n);
} else { #real data: New-Jersey lottery payoff as a function of the number
  njlot <- read.table("njlottery.txt", header=T, sep="\t");
  tu <- njlot$number; 
  yu <- njlot$payoff; 
  t<- sort(tu); y <- yu[order(tu)];
  n<- length(njlot$number);
}

# plot a spline regression 
tt <- seq(0,max(tu), 0.01);
smo <- sreg(t, y, df=15); # Smoothing spline regression
yfit0 <- predict(smo, tt); 
plot(t, y, type='p', lwd=3)
lines(tt, yfit0, type='l', col='blue', lwd=8)

# build a "bootstrap confidence tube" for the spline regression
N<- 100;
for (k in 1:N){
  bK <- sample(1:n, replace=T); bK <- sort(bK);
  bt <- t[bK];
  by <- y[bK]; #bKo <- order(bK); by <- y[bKo]; bt <- bt[bKo];
  smo <- sreg(bt, by, df=15);
  yfit <- predict(smo, tt); 
  #lines(bt, by, type='b', lwd=3, col='cyan')
  lines(tt, yfit, type='l', col='magenta', lwd=1)
}

lines(tt, yfit0, type='l', col='blue', lwd=8)

if(experiment == "simulated"){ # plot true function 
  lines(tt, f(tt), type='l', col='red', lwd=2)
}