rm(list=ls(all=T))

library(mgcv)
library(gss)


################################################################################
# Function for preliminary detection of outliers
################################################################################
fn.outliers <- function(X,Y,a1=2.5)
{
 x1<- X[,1];x2<- X[,2];x3<- X[,3];x4<- X[,4] # Must be corresponding to the dimension of X
 temp.model<- ssanova(Y~x1+x2+x3+x4,alpha=1)
 temp.X<- data.frame(x1=x1,x2=x2,x3=x3,x4=x4)
 temp.fitted<- predict(temp.model,newdata=temp.X)
 temp.resid<- Y-temp.fitted
 temp.resid<- (temp.resid-mean(temp.resid))/sd(temp.resid)
 n.out <- length(which(abs(temp.resid)>a1))- 2*(1-pnorm(a1))*length(Y)
 if (n.out< 0)
 n.out<- 0
 return(round(n.out))
}
################################################################################
#  Outlier Shifting Smoothing Splines ANOVA
################################################################################

out.shift.ssanova<-function(X,Y,lambda,a1=2.5,weights=NULL){
 n<-length(Y)
 n.out <-fn.outliers(X,Y,a1)
 x1<- X[,1];x2<- X[,2];x3<- X[,3];x4<- X[,4] # Must be corresponding to the dimension of X

 temp.ss<- ssanova(Y ~ x1+x2+x3+x4,alpha=lambda,method="m",weights=weights)
 NEW.data<- data.frame(x1=x1,x2=x2,x3=x3,x4=x4)
 fit.est<- predict(temp.ss,NEW.data)
 sigma.est<- mad(Y-fit.est)
 gamma.est<-rep(0,n)
 Y.old<- Y

 if (n.out>0) { lambda.gamma<- sigma.est*qnorm((2*n-n.out)/(2*n))}
 if (n.out==0) {lambda.gamma<- sigma.est*qnorm((2*n-1e-06)/(2*n))}

 tol<- 100000
 n.iter <- 0
 while(tol>1e-5 & n.iter<100 )
{
fit.pre<- fit.est
nonzero.index<-which(abs(Y.old-fit.est)>=lambda.gamma)
gamma.est[nonzero.index]<- (Y.old-fit.est)[nonzero.index]
Y.new<- Y.old - gamma.est
model.update <- ssanova(Y.new~ x1+x2+x3+x4,alpha=lambda,weights = weights,method="m")
fit.est<- predict(model.update,NEW.data)
tol<- mean((fit.pre-fit.est)^2)
n.iter<- n.iter+1
}

Y.fit<- predict(model.update,NEW.data)

object<- list(fit=Y.fit, iter = n.iter, sigma.est = sigma.est, gamma.est = gamma.est,model=model.update,n.out=n.out)
}

################################################################################
#   Outlier Shifting Generalized Additive Model with  \tilde{lambda}_gamma    ##
################################################################################

out.shift.gam<-function(X,Y,a1=2.5,weights=NULL)
{
 n<-length(Y)
 n.out <-fn.outliers(X,Y,a1)
 x1<- X[,1]; x2<- X[,2]; x3<- X[,3]; x4<- X[,4] # Should be corresponding to the dimension of X

 temp.gam<- gam(Y ~ s(x1,bs="ts")+ s(x2,bs="ts")+s(x3,bs="ts")+s(x4,bs="ts"),weights=weights,method="REML")
 fit.est<- temp.gam$fitted
 sigma.est<- mad(Y-fit.est)
 gamma.est<-rep(0,n)
 Y.old<- Y
 lambda.gamma<- sigma.est*qnorm((2*n-n.out)/(2*n))

 tol<- 100000
 n.iter <- 0
 while(tol>1e-5 & n.iter<100 )
{
fit.pre<- fit.est
nonzero.index<-which(abs(Y.old-fit.est)>=lambda.gamma)
gamma.est[nonzero.index]<- (Y.old-fit.est)[nonzero.index]
Y.new<- Y.old - gamma.est
model.update <- gam(Y.new ~ s(x1,bs="ts")+ s(x2,bs="ts")+s(x3,bs="ts")+s(x4,bs="ts"),weights=weights,method="REML")
fit.est<- model.update$fitted
tol<- mean((fit.pre-fit.est)^2)
n.iter<- n.iter+1
}

Y.fit<- model.update$fitted

object<- list(fit=Y.fit, iter = n.iter, sigma.est = sigma.est, gamma.est = gamma.est,model=model.update)
}

##################################
# Simulation Begins from here.   #
##################################

f1<-function(x){ temp<- 4*sin(pi*x)
return(temp) }
f2<-function(x){ temp<- exp(2*x)-2.76
return(temp) }
f3<-function(x){ temp<- x^11*(10*(1-x))^6+10^4*x^3*(1-x)^10-1.4
return(temp) }
f4<- function(x) {
 return(rep(-1,length(x)))
 }

n.sam = 200
n.rep=500
con.prop=0.05
n.con = n.sam*con.prop

lambda.ssanova = lambda.ssanova.w02 = lambda.ssanova.w04 = lambda.ssanova.w05 = lambda.ssanova.w06 = lambda.ssanova.w08 = rep(0,n.rep)
lambda.ssanova.S = lambda.ssanova.w02.S = lambda.ssanova.w04.S = lambda.ssanova.w05.S = lambda.ssanova.w06.S = lambda.ssanova.w08.S = rep(0,n.rep)
gamma.ssanova = gamma.ssanova.w02= gamma.ssanova.w04= gamma.ssanova.w05= gamma.ssanova.w06= gamma.ssanova.w08 = matrix(0,nrow=n.rep,ncol=n.sam)
outlier.mat = matrix(0,nrow=n.rep,ncol=n.sam)
fitted.ssanova = fitted.ssanova.w02 = fitted.ssanova.w04 = fitted.ssanova.w05 = fitted.ssanova.w06 = fitted.ssanova.w08 =matrix(0,nrow=n.rep,ncol=n.sam)
fitted.ssanova.S = fitted.ssanova.w02.S = fitted.ssanova.w04.S = fitted.ssanova.w05.S = fitted.ssanova.w06.S = fitted.ssanova.w08.S =matrix(0,nrow=n.rep,ncol=n.sam)

abs.err.ssanova = abs.err.ssanova.w02= abs.err.ssanova.w04 =abs.err.ssanova.w05 =abs.err.ssanova.w06 =abs.err.ssanova.w08 = rep(0,n.rep) # absolute error. Cannot find for each term
abs.err.ssanova.S=abs.err.ssanova.w02.S=abs.err.ssanova.w04.S=abs.err.ssanova.w05.S=abs.err.ssanova.w06.S=abs.err.ssanova.w08.S=  rep(0,n.rep)

sq.err.ssanova = sq.err.ssanova.w02= sq.err.ssanova.w04=sq.err.ssanova.w05=sq.err.ssanova.w06=sq.err.ssanova.w08 =  rep(0,n.rep) # square error by term
sq.err.ssanova.S = sq.err.ssanova.w02.S = sq.err.ssanova.w04.S=sq.err.ssanova.w05.S=sq.err.ssanova.w06.S=sq.err.ssanova.w08.S =  rep(0,n.rep)


sigma.ssanova = sigma.ssanova.w02 = sigma.ssanova.w04 = sigma.ssanova.w05 = sigma.ssanova.w06 = sigma.ssanova.w08 = rep(0,n.rep)
sigma.ssanova.S = sigma.ssanova.w02.S = sigma.ssanova.w04.S = sigma.ssanova.w05.S = sigma.ssanova.w06.S = sigma.ssanova.w08.S = rep(0,n.rep)
################################################################################
lambda.gam = lambda.gam.w02 = lambda.gam.w04 = lambda.gam.w05 = lambda.gam.w06 = lambda.gam.w08 = rep(0,n.rep)
lambda.gam.S = lambda.gam.w02.S = lambda.gam.w04.S = lambda.gam.w05.S = lambda.gam.w06.S = lambda.gam.w08.S = rep(0,n.rep)
gamma.gam = gamma.gam.w02= gamma.gam.w04= gamma.gam.w05= gamma.gam.w06= gamma.gam.w08 = matrix(0,nrow=n.rep,ncol=n.sam)
fitted.gam = fitted.gam.w02 = fitted.gam.w04 = fitted.gam.w05 = fitted.gam.w06 = fitted.gam.w08 =matrix(0,nrow=n.rep,ncol=n.sam)
fitted.gam.S = fitted.gam.w02.S = fitted.gam.w04.S = fitted.gam.w05.S = fitted.gam.w06.S = fitted.gam.w08.S =matrix(0,nrow=n.rep,ncol=n.sam)

abs.err.gam = abs.err.gam.w02= abs.err.gam.w04 =abs.err.gam.w05 =abs.err.gam.w06 =abs.err.gam.w08 = matrix(0,nrow=n.rep,ncol=4) # absolute error by term
abs.err.gam.S=abs.err.gam.w02.S=abs.err.gam.w04.S=abs.err.gam.w05.S=abs.err.gam.w06.S=abs.err.gam.w08.S= matrix(0,nrow=n.rep,ncol=4)

sq.err.gam = sq.err.gam.w02= sq.err.gam.w04=sq.err.gam.w05=sq.err.gam.w06=sq.err.gam.w08 = matrix(0,nrow=n.rep,ncol=4) # square error by term
sq.err.gam.S = sq.err.gam.w02.S = sq.err.gam.w04.S=sq.err.gam.w05.S=sq.err.gam.w06.S=sq.err.gam.w08.S = matrix(0,nrow=n.rep,ncol=4)

abs.err = abs.err.w02 = abs.err.w04 = abs.err.w05=abs.err.w06=abs.err.w08=rep(0,n.rep)
abs.err.S = abs.err.w02.S = abs.err.w04.S = abs.err.w05.S=abs.err.w06.S=abs.err.w08.S=rep(0,n.rep)
sq.err = sq.err.w02 = sq.err.w04 = sq.err.w05=sq.err.w06=sq.err.w08=rep(0,n.rep)
sq.err.S = sq.err.w02.S = sq.err.w04.S = sq.err.w05.S=sq.err.w06.S=sq.err.w08.S=rep(0,n.rep)


sigma.gam = sigma.gam.w02 = sigma.gam.w04 = sigma.gam.w05 = sigma.gam.w06 = sigma.gam.w08 = rep(0,n.rep)
sigma.gam.S = sigma.gam.w02.S = sigma.gam.w04.S = sigma.gam.w05.S = sigma.gam.w06.S = sigma.gam.w08.S = rep(0,n.rep)


for (i in 1:n.rep){

set.seed(i+100)
x1<-runif(n.sam);x2<-runif(n.sam);x3<-runif(n.sam);x4<-runif(n.sam)
error = epsilon = rnorm(n.sam,0,2)
index.out<- which(rank(-abs(epsilon))<n.con+1)
error[index.out]<- 2*epsilon[index.out]
outlier.mat[i,index.out]<-1

true.fct<- f1(x1)+f2(x2)+f3(x3)+f4(x4)
Y<- true.fct +error
X<- cbind(x1,x2,x3,x4)

################################################################################
#  Fit SS-ANOVA and outlier shifting SS-ANOVA
################################################################################

 model.ssanova<-  ssanova(Y ~ x1+x2+x3+x4,method="m")   # method="m" for REML.
  # Find weights
 temp.weight<- abs(summary(model.ssanova)$residuals)
 weight02 <- 1/(0.2*temp.weight+0.8*mean(temp.weight))
 weight04 <- 1/(0.4*temp.weight+0.6*mean(temp.weight))
 weight05 <- 1/(0.5*temp.weight+0.5*mean(temp.weight))
 weight06 <- 1/(0.6*temp.weight+0.4*mean(temp.weight))
 weight08 <- 1/(0.8*temp.weight+0.2*mean(temp.weight))
 # Fit Weighted SS-ANOVA
 model.ssanova.w02 <- ssanova(Y ~ x1+x2+x3+x4,alpha=1,weights=weight02,method="m")
 model.ssanova.w04 <- ssanova(Y ~ x1+x2+x3+x4,alpha=1,weights=weight04,method="m")
 model.ssanova.w05 <- ssanova(Y ~ x1+x2+x3+x4,alpha=1,weights=weight05,method="m")
 model.ssanova.w06 <- ssanova(Y ~ x1+x2+x3+x4,alpha=1,weights=weight06,method="m")
 model.ssanova.w08 <- ssanova(Y ~ x1+x2+x3+x4,alpha=1,weights=weight08,method="m")
 # record lambda selected by GCV
 lambda.ssanova[i] <- summary(model.ssanova)$penalty
 lambda.ssanova.w02[i]<- summary(model.ssanova.w02)$penalty
 lambda.ssanova.w04[i]<- summary(model.ssanova.w02)$penalty
 lambda.ssanova.w05[i]<- summary(model.ssanova.w02)$penalty
 lambda.ssanova.w06[i]<- summary(model.ssanova.w02)$penalty
 lambda.ssanova.w08[i]<- summary(model.ssanova.w02)$penalty
 # fitted value
 fitted.ssanova[i,] <- summary(model.ssanova)$fitted
 fitted.ssanova.w02[i,] <- summary(model.ssanova.w02)$fitted
 fitted.ssanova.w04[i,] <- summary(model.ssanova.w04)$fitted
 fitted.ssanova.w05[i,] <- summary(model.ssanova.w05)$fitted
 fitted.ssanova.w06[i,] <- summary(model.ssanova.w06)$fitted
 fitted.ssanova.w08[i,] <- summary(model.ssanova.w08)$fitted
 # estimate of scale parameter
 sigma.ssanova[i] = summary(model.ssanova)$sigma
 sigma.ssanova.w02[i] = summary(model.ssanova.w02)$sigma
 sigma.ssanova.w04[i] = summary(model.ssanova.w04)$sigma
 sigma.ssanova.w05[i] = summary(model.ssanova.w05)$sigma
 sigma.ssanova.w06[i] = summary(model.ssanova.w06)$sigma
 sigma.ssanova.w08[i] = summary(model.ssanova.w08)$sigma

# average absolute errors
 abs.err.ssanova[i]<- mean(abs(fitted.ssanova[i,]-true.fct))
 abs.err.ssanova.w02[i]<- mean(abs(fitted.ssanova.w02[i,]-true.fct))
 abs.err.ssanova.w04[i]<- mean(abs(fitted.ssanova.w04[i,]-true.fct))
 abs.err.ssanova.w05[i]<- mean(abs(fitted.ssanova.w05[i,]-true.fct))
 abs.err.ssanova.w06[i]<- mean(abs(fitted.ssanova.w06[i,]-true.fct))
 abs.err.ssanova.w08[i]<- mean(abs(fitted.ssanova.w08[i,]-true.fct))

# average squared errors
 sq.err.ssanova[i]<- mean((fitted.ssanova[i,]-true.fct)^2)
 sq.err.ssanova.w02[i]<- mean((fitted.ssanova.w02[i,]-true.fct)^2)
 sq.err.ssanova.w04[i]<- mean((fitted.ssanova.w04[i,]-true.fct)^2)
 sq.err.ssanova.w05[i]<- mean((fitted.ssanova.w05[i,]-true.fct)^2)
 sq.err.ssanova.w06[i]<- mean((fitted.ssanova.w06[i,]-true.fct)^2)
 sq.err.ssanova.w08[i]<- mean((fitted.ssanova.w08[i,]-true.fct)^2)
 

################################################################################  
 # Fit weighted Outlier shifting SS-ANOVA
 model.ssanova.shift<- out.shift.ssanova(X,Y,lambda=1)
 model.ssanova.shift.w02<- out.shift.ssanova(X,Y,lambda=1,weights=weight02)
 model.ssanova.shift.w04<- out.shift.ssanova(X,Y,lambda=1,weights=weight04)
 model.ssanova.shift.w05<- out.shift.ssanova(X,Y,lambda=1,weights=weight05)
 model.ssanova.shift.w06<- out.shift.ssanova(X,Y,lambda=1,weights=weight06)
 model.ssanova.shift.w08<- out.shift.ssanova(X,Y,lambda=1,weights=weight08)
 # record lambda selected by GCV
 lambda.ssanova.S[i] <- summary(model.ssanova.shift$model)$penalty
 lambda.ssanova.w02.S[i]<- summary(model.ssanova.shift.w02$model)$penalty
 lambda.ssanova.w04.S[i]<- summary(model.ssanova.shift.w04$model)$penalty
 lambda.ssanova.w05.S[i]<- summary(model.ssanova.shift.w05$model)$penalty
 lambda.ssanova.w06.S[i]<- summary(model.ssanova.shift.w06$model)$penalty
 lambda.ssanova.w08.S[i]<- summary(model.ssanova.shift.w08$model)$penalty
# record the estimated gamma
 gamma.ssanova[i,] <- model.ssanova.shift$gamma.est
 gamma.ssanova.w02[i,] <- model.ssanova.shift.w02$gamma.est
 gamma.ssanova.w04[i,] <- model.ssanova.shift.w04$gamma.est
 gamma.ssanova.w05[i,] <- model.ssanova.shift.w05$gamma.est
 gamma.ssanova.w06[i,] <- model.ssanova.shift.w06$gamma.est
 gamma.ssanova.w08[i,] <- model.ssanova.shift.w08$gamma.est
  # fitted value
 fitted.ssanova.S[i,] <- model.ssanova.shift$fit
 fitted.ssanova.w02.S[i,] <- model.ssanova.shift.w02$fit 
 fitted.ssanova.w04.S[i,] <- model.ssanova.shift.w04$fit 
 fitted.ssanova.w05.S[i,] <- model.ssanova.shift.w05$fit 
 fitted.ssanova.w06.S[i,] <- model.ssanova.shift.w06$fit 
 fitted.ssanova.w08.S[i,] <- model.ssanova.shift.w08$fit 
  # estimate of scale parameter
 sigma.ssanova.S[i] = summary(model.ssanova.shift$model)$sigma
 sigma.ssanova.w02[i] = summary(model.ssanova.shift.w02$model)$sigma
 sigma.ssanova.w04[i] = summary(model.ssanova.shift.w04$model)$sigma
 sigma.ssanova.w05[i] = summary(model.ssanova.shift.w05$model)$sigma
 sigma.ssanova.w06[i] = summary(model.ssanova.shift.w06$model)$sigma
 sigma.ssanova.w08[i] = summary(model.ssanova.shift.w08$model)$sigma
 # average absolute errors
 abs.err.ssanova.S[i]<- mean(abs(fitted.ssanova.S[i,]-true.fct))
 abs.err.ssanova.w02.S[i]<- mean(abs(fitted.ssanova.w02.S[i,]-true.fct))
 abs.err.ssanova.w04.S[i]<- mean(abs(fitted.ssanova.w04.S[i,]-true.fct))
 abs.err.ssanova.w05.S[i]<- mean(abs(fitted.ssanova.w05.S[i,]-true.fct))
 abs.err.ssanova.w06.S[i]<- mean(abs(fitted.ssanova.w06.S[i,]-true.fct))
 abs.err.ssanova.w08.S[i]<- mean(abs(fitted.ssanova.w08.S[i,]-true.fct))

 # average squared errors
 sq.err.ssanova.S[i]<- mean((fitted.ssanova.S[i,]-true.fct)^2)
 sq.err.ssanova.w02.S[i]<- mean((fitted.ssanova.w02.S[i,]-true.fct)^2)
 sq.err.ssanova.w04.S[i]<- mean((fitted.ssanova.w04.S[i,]-true.fct)^2)
 sq.err.ssanova.w05.S[i]<- mean((fitted.ssanova.w05.S[i,]-true.fct)^2)
 sq.err.ssanova.w06.S[i]<- mean((fitted.ssanova.w06.S[i,]-true.fct)^2)
 sq.err.ssanova.w08.S[i]<- mean((fitted.ssanova.w08.S[i,]-true.fct)^2)


################################################################################
#  Fit GAM and outlier shifting GAM
################################################################################

 model.gam<- gam(Y ~ s(x1,bs="ts")+ s(x2,bs="ts")+s(x3,bs="ts")+s(x4,bs="ts"),method="REML")

   # Find weights
 temp.weight<- abs(model.gam$residuals)
 weight02 <- 1/(0.2*temp.weight+0.8*mean(temp.weight))
 weight04 <- 1/(0.4*temp.weight+0.6*mean(temp.weight))
 weight05 <- 1/(0.5*temp.weight+0.5*mean(temp.weight))
 weight06 <- 1/(0.6*temp.weight+0.4*mean(temp.weight))
 weight08 <- 1/(0.8*temp.weight+0.2*mean(temp.weight))
 # Fit Weighted GAM
 model.gam.w02 <- gam(Y ~ s(x1,bs="ts")+ s(x2,bs="ts")+s(x3,bs="ts")+s(x4,bs="ts"),weights=weight02,method="REML")
 model.gam.w04 <- gam(Y ~ s(x1,bs="ts")+ s(x2,bs="ts")+s(x3,bs="ts")+s(x4,bs="ts"),weights=weight04,method="REML")
 model.gam.w05 <- gam(Y ~ s(x1,bs="ts")+ s(x2,bs="ts")+s(x3,bs="ts")+s(x4,bs="ts"),weights=weight05,method="REML")
 model.gam.w06 <- gam(Y ~ s(x1,bs="ts")+ s(x2,bs="ts")+s(x3,bs="ts")+s(x4,bs="ts"),weights=weight06,method="REML")
 model.gam.w08 <- gam(Y ~ s(x1,bs="ts")+ s(x2,bs="ts")+s(x3,bs="ts")+s(x4,bs="ts"),weights=weight08,method="REML")
 # record GCV values
 lambda.gam[i] <- summary(model.gam)$sp.criterion
 lambda.gam.w02[i]<- summary(model.gam.w02)$sp.criterion
 lambda.gam.w04[i]<- summary(model.gam.w02)$sp.criterion
 lambda.gam.w05[i]<- summary(model.gam.w02)$sp.criterion
 lambda.gam.w06[i]<- summary(model.gam.w02)$sp.criterion
 lambda.gam.w08[i]<- summary(model.gam.w02)$sp.criterion
 # fitted value
 fitted.gam[i,] <- model.gam$fitted
 fitted.gam.w02[i,] <- model.gam.w02$fitted
 fitted.gam.w04[i,] <- model.gam.w04$fitted
 fitted.gam.w05[i,] <- model.gam.w05$fitted
 fitted.gam.w06[i,] <- model.gam.w06$fitted
 fitted.gam.w08[i,] <- model.gam.w08$fitted
 # estimate of scale parameter
 sigma.gam[i] = summary(model.gam)$dispersion
 sigma.gam.w02[i] = summary(model.gam.w02)$dispersion
 sigma.gam.w04[i] = summary(model.gam.w04)$dispersion
 sigma.gam.w05[i] = summary(model.gam.w05)$dispersion
 sigma.gam.w06[i] = summary(model.gam.w06)$dispersion
 sigma.gam.w08[i] = summary(model.gam.w08)$dispersion

 # average absolute errors by the term
 abs.err.gam[i,] <- c(mean(abs(predict(model.gam,type="terms")[,1]+mean(f1(x1))-f1(x1))),
mean(abs(predict(model.gam,type="terms")[,2]+mean(f2(x2))-f2(x2))),mean(abs(predict(model.gam,type="terms")[,3]+mean(f3(x3))-f3(x3))),
mean(abs(predict(model.gam,type="terms")[,4]+mean(f4(x4))-f4(x4))))
 abs.err.gam.w02[i,] <- c(mean(abs(predict(model.gam.w02,type="terms")[,1]+mean(f1(x1))-f1(x1))),
mean(abs(predict(model.gam.w02,type="terms")[,2]+mean(f2(x2))-f2(x2))),mean(abs(predict(model.gam.w02,type="terms")[,3]+mean(f3(x3))-f3(x3))),
mean(abs(predict(model.gam.w02,type="terms")[,4]+mean(f4(x4))-f4(x4))))
 abs.err.gam.w04[i,] <- c(mean(abs(predict(model.gam.w04,type="terms")[,1]+mean(f1(x1))-f1(x1))),
mean(abs(predict(model.gam.w04,type="terms")[,2]+mean(f2(x2))-f2(x2))),mean(abs(predict(model.gam.w04,type="terms")[,3]+mean(f3(x3))-f3(x3))),
mean(abs(predict(model.gam.w04,type="terms")[,4]+mean(f4(x4))-f4(x4))))
 abs.err.gam.w05[i,] <- c(mean(abs(predict(model.gam.w05,type="terms")[,1]+mean(f1(x1))-f1(x1))),
mean(abs(predict(model.gam.w05,type="terms")[,2]+mean(f2(x2))-f2(x2))),mean(abs(predict(model.gam.w05,type="terms")[,3]+mean(f3(x3))-f3(x3))),
mean(abs(predict(model.gam.w05,type="terms")[,4]+mean(f4(x4))-f4(x4))))
 abs.err.gam.w06[i,] <- c(mean(abs(predict(model.gam.w06,type="terms")[,1]+mean(f1(x1))-f1(x1))),
mean(abs(predict(model.gam.w06,type="terms")[,2]+mean(f2(x2))-f2(x2))),mean(abs(predict(model.gam.w06,type="terms")[,3]+mean(f3(x3))-f3(x3))),
mean(abs(predict(model.gam.w06,type="terms")[,4]+mean(f4(x4))-f4(x4))))
 abs.err.gam.w08[i,] <- c(mean(abs(predict(model.gam.w08,type="terms")[,1]+mean(f1(x1))-f1(x1))),
mean(abs(predict(model.gam.w08,type="terms")[,2]+mean(f2(x2))-f2(x2))),mean(abs(predict(model.gam.w08,type="terms")[,3]+mean(f3(x3))-f3(x3))),
mean(abs(predict(model.gam.w08,type="terms")[,4]+mean(f4(x4))-f4(x4))))

 # average squared errors by the term
 sq.err.gam[i,] <- c(mean((predict(model.gam,type="terms")[,1]+mean(f1(x1))-f1(x1))^2),
mean((predict(model.gam,type="terms")[,2]+mean(f2(x2))-f2(x2))^2),mean((predict(model.gam,type="terms")[,3]+mean(f3(x3))-f3(x3))^2),
mean((predict(model.gam,type="terms")[,4]+mean(f4(x4))-f4(x4))^2))
 sq.err.gam.w02[i,] <- c(mean((predict(model.gam.w02,type="terms")[,1]+mean(f1(x1))-f1(x1))^2),
mean((predict(model.gam.w02,type="terms")[,2]+mean(f2(x2))-f2(x2))^2),mean((predict(model.gam.w02,type="terms")[,3]+mean(f3(x3))-f3(x3))^2),
mean((predict(model.gam.w02,type="terms")[,4]+mean(f4(x4))-f4(x4))^2))
 sq.err.gam.w04[i,] <- c(mean((predict(model.gam.w04,type="terms")[,1]+mean(f1(x1))-f1(x1))^2),
mean((predict(model.gam.w04,type="terms")[,2]+mean(f2(x2))-f2(x2))^2),mean((predict(model.gam.w04,type="terms")[,3]+mean(f3(x3))-f3(x3))^2),
mean((predict(model.gam.w04,type="terms")[,4]+mean(f4(x4))-f4(x4))^2))
 sq.err.gam.w05[i,] <- c(mean((predict(model.gam.w05,type="terms")[,1]+mean(f1(x1))-f1(x1))^2),
mean((predict(model.gam.w05,type="terms")[,2]+mean(f2(x2))-f2(x2))^2),mean((predict(model.gam.w05,type="terms")[,3]+mean(f3(x3))-f3(x3))^2),
mean((predict(model.gam.w05,type="terms")[,4]+mean(f4(x4))-f4(x4))^2))
 sq.err.gam.w06[i,] <- c(mean((predict(model.gam.w06,type="terms")[,1]+mean(f1(x1))-f1(x1))^2),
mean((predict(model.gam.w06,type="terms")[,2]+mean(f2(x2))-f2(x2))^2),mean((predict(model.gam.w06,type="terms")[,3]+mean(f3(x3))-f3(x3))^2),
mean((predict(model.gam.w06,type="terms")[,4]+mean(f4(x4))-f4(x4))^2))
 sq.err.gam.w08[i,] <- c(mean((predict(model.gam.w08,type="terms")[,1]+mean(f1(x1))-f1(x1))^2),
mean((predict(model.gam.w08,type="terms")[,2]+mean(f2(x2))-f2(x2))^2),mean((predict(model.gam.w08,type="terms")[,3]+mean(f3(x3))-f3(x3))^2),
mean((predict(model.gam.w08,type="terms")[,4]+mean(f4(x4))-f4(x4))^2))
  
 # average absolute errors of the fitted value
 abs.err[i] = mean(abs(fitted.gam[i,]-true.fct)) 
 abs.err.w02[i] = mean(abs(fitted.gam.w02[i,]-true.fct))
 abs.err.w04[i] = mean(abs(fitted.gam.w04[i,]-true.fct))
 abs.err.w05[i] = mean(abs(fitted.gam.w05[i,]-true.fct))
 abs.err.w06[i] = mean(abs(fitted.gam.w06[i,]-true.fct))
 abs.err.w08[i] = mean(abs(fitted.gam.w08[i,]-true.fct))
 
 # average squared errors of the fitted value
 sq.err[i] = mean((fitted.gam[i,]-true.fct)^2) 
 sq.err.w02[i] = mean((fitted.gam.w02[i,]-true.fct)^2)
 sq.err.w04[i] = mean((fitted.gam.w04[i,]-true.fct)^2)
 sq.err.w05[i] = mean((fitted.gam.w05[i,]-true.fct)^2)
 sq.err.w06[i] = mean((fitted.gam.w06[i,]-true.fct)^2)
 sq.err.w08[i] = mean((fitted.gam.w08[i,]-true.fct)^2)
 
################################################################################  
#  Fit weighted Outlier shifting GAM
 
 model.gam.shift<- out.shift.gam(X,Y)
 model.gam.shift.w02<- out.shift.gam(X,Y,weights=weight02)
 model.gam.shift.w04<- out.shift.gam(X,Y,weights=weight04)
 model.gam.shift.w05<- out.shift.gam(X,Y,weights=weight05)
 model.gam.shift.w06<- out.shift.gam(X,Y,weights=weight06)
 model.gam.shift.w08<- out.shift.gam(X,Y,weights=weight08)
 # record  selected GCV
 lambda.gam.S[i] <- summary(model.gam.shift$model)$sp.criterion
 lambda.gam.w02.S[i]<- summary(model.gam.shift.w02$model)$sp.criterion
 lambda.gam.w04.S[i]<- summary(model.gam.shift.w04$model)$sp.criterion
 lambda.gam.w05.S[i]<- summary(model.gam.shift.w05$model)$sp.criterion
 lambda.gam.w06.S[i]<- summary(model.gam.shift.w06$model)$sp.criterion
 lambda.gam.w08.S[i]<- summary(model.gam.shift.w08$model)$sp.criterion

# record the estimated gamma
 gamma.gam[i,] <- model.gam.shift$gamma.est
 gamma.gam.w02[i,] <- model.gam.shift.w02$gamma.est
 gamma.gam.w04[i,] <- model.gam.shift.w04$gamma.est
 gamma.gam.w05[i,] <- model.gam.shift.w05$gamma.est
 gamma.gam.w06[i,] <- model.gam.shift.w06$gamma.est
 gamma.gam.w08[i,] <- model.gam.shift.w08$gamma.est
  # fitted value
 fitted.gam.S[i,] <- model.gam.shift$fit
 fitted.gam.w02.S[i,] <- model.gam.shift.w02$fit 
 fitted.gam.w04.S[i,] <- model.gam.shift.w04$fit 
 fitted.gam.w05.S[i,] <- model.gam.shift.w05$fit 
 fitted.gam.w06.S[i,] <- model.gam.shift.w06$fit 
 fitted.gam.w08.S[i,] <- model.gam.shift.w08$fit 
  # estimate of scale parameter
 sigma.gam.S[i] = summary(model.gam.shift$model)$dispersion
 sigma.gam.w02.S[i] = summary(model.gam.shift.w02$model)$dispersion
 sigma.gam.w04.S[i] = summary(model.gam.shift.w04$model)$dispersion
 sigma.gam.w05.S[i] = summary(model.gam.shift.w05$model)$dispersion
 sigma.gam.w06.S[i] = summary(model.gam.shift.w06$model)$dispersion
 sigma.gam.w08.S[i] = summary(model.gam.shift.w08$model)$dispersion


 # average absolute errors by the term
 abs.err.gam.S[i,] <- c(mean(abs(predict(model.gam.shift$model,type="terms")[,1]+mean(f1(x1))-f1(x1))),
mean(abs(predict(model.gam.shift$model,type="terms")[,2]+mean(f2(x2))-f2(x2))),mean(abs(predict(model.gam.shift$model,type="terms")[,3]+mean(f3(x3))-f3(x3))),
mean(abs(predict(model.gam.shift$model,type="terms")[,4]+mean(f4(x4))-f4(x4))))
 abs.err.gam.w02.S[i,] <- c(mean(abs(predict(model.gam.shift.w02$model,type="terms")[,1]+mean(f1(x1))-f1(x1))),
mean(abs(predict(model.gam.shift.w02$model,type="terms")[,2]+mean(f2(x2))-f2(x2))),mean(abs(predict(model.gam.shift.w02$model,type="terms")[,3]+mean(f3(x3))-f3(x3))),
mean(abs(predict(model.gam.shift.w02$model,type="terms")[,4]+mean(f4(x4))-f4(x4))))
 abs.err.gam.w04.S[i,] <- c(mean(abs(predict(model.gam.shift.w04$model,type="terms")[,1]+mean(f1(x1))-f1(x1))),
mean(abs(predict(model.gam.shift.w04$model,type="terms")[,2]+mean(f2(x2))-f2(x2))),mean(abs(predict(model.gam.shift.w04$model,type="terms")[,3]+mean(f3(x3))-f3(x3))),
mean(abs(predict(model.gam.shift.w04$model,type="terms")[,4]+mean(f4(x4))-f4(x4))))
 abs.err.gam.w05.S[i,] <- c(mean(abs(predict(model.gam.shift.w05$model,type="terms")[,1]+mean(f1(x1))-f1(x1))),
mean(abs(predict(model.gam.shift.w05$model,type="terms")[,2]+mean(f2(x2))-f2(x2))),mean(abs(predict(model.gam.shift.w05$model,type="terms")[,3]+mean(f3(x3))-f3(x3))),
mean(abs(predict(model.gam.shift.w05$model,type="terms")[,4]+mean(f4(x4))-f4(x4))))
 abs.err.gam.w06.S[i,] <- c(mean(abs(predict(model.gam.shift.w06$model,type="terms")[,1]+mean(f1(x1))-f1(x1))),
mean(abs(predict(model.gam.shift.w06$model,type="terms")[,2]+mean(f2(x2))-f2(x2))),mean(abs(predict(model.gam.shift.w06$model,type="terms")[,3]+mean(f3(x3))-f3(x3))),
mean(abs(predict(model.gam.shift.w06$model,type="terms")[,4]+mean(f4(x4))-f4(x4))))
 abs.err.gam.w08.S[i,] <- c(mean(abs(predict(model.gam.shift.w08$model,type="terms")[,1]+mean(f1(x1))-f1(x1))),
mean(abs(predict(model.gam.shift.w08$model,type="terms")[,2]+mean(f2(x2))-f2(x2))),mean(abs(predict(model.gam.shift.w08$model,type="terms")[,3]+mean(f3(x3))-f3(x3))),
mean(abs(predict(model.gam.shift.w08$model,type="terms")[,4]+mean(f4(x4))-f4(x4))))

 # average squared errors by the term
 sq.err.gam.S[i,] <- c(mean((predict(model.gam.shift$model,type="terms")[,1]+mean(f1(x1))-f1(x1))^2),
mean((predict(model.gam.shift$model,type="terms")[,2]+mean(f2(x2))-f2(x2))^2),mean((predict(model.gam.shift$model,type="terms")[,3]+mean(f3(x3))-f3(x3))^2),
mean((predict(model.gam.shift$model,type="terms")[,4]+mean(f4(x4))-f4(x4))^2))
 sq.err.gam.w02.S[i,] <- c(mean((predict(model.gam.shift.w02$model,type="terms")[,1]+mean(f1(x1))-f1(x1))^2),
mean((predict(model.gam.shift.w02$model,type="terms")[,2]+mean(f2(x2))-f2(x2))^2),mean((predict(model.gam.shift.w02$model,type="terms")[,3]+mean(f3(x3))-f3(x3))^2),
mean((predict(model.gam.shift.w02$model,type="terms")[,4]+mean(f4(x4))-f4(x4))^2))
 sq.err.gam.w04.S[i,] <- c(mean((predict(model.gam.shift.w04$model,type="terms")[,1]+mean(f1(x1))-f1(x1))^2),
mean((predict(model.gam.shift.w04$model,type="terms")[,2]+mean(f2(x2))-f2(x2))^2),mean((predict(model.gam.shift.w04$model,type="terms")[,3]+mean(f3(x3))-f3(x3))^2),
mean((predict(model.gam.shift.w04$model,type="terms")[,4]+mean(f4(x4))-f4(x4))^2))
 sq.err.gam.w05.S[i,] <- c(mean((predict(model.gam.shift.w05$model,type="terms")[,1]+mean(f1(x1))-f1(x1))^2),
mean((predict(model.gam.shift.w05$model,type="terms")[,2]+mean(f2(x2))-f2(x2))^2),mean((predict(model.gam.shift.w05$model,type="terms")[,3]+mean(f3(x3))-f3(x3))^2),
mean((predict(model.gam.shift.w05$model,type="terms")[,4]+mean(f4(x4))-f4(x4))^2))
 sq.err.gam.w06.S[i,] <- c(mean((predict(model.gam.shift.w06$model,type="terms")[,1]+mean(f1(x1))-f1(x1))^2),
mean((predict(model.gam.shift.w06$model,type="terms")[,2]+mean(f2(x2))-f2(x2))^2),mean((predict(model.gam.shift.w06$model,type="terms")[,3]+mean(f3(x3))-f3(x3))^2),
mean((predict(model.gam.shift.w06$model,type="terms")[,4]+mean(f4(x4))-f4(x4))^2))
 sq.err.gam.w08.S[i,] <- c(mean((predict(model.gam.shift.w08$model,type="terms")[,1]+mean(f1(x1))-f1(x1))^2),
mean((predict(model.gam.shift.w08$model,type="terms")[,2]+mean(f2(x2))-f2(x2))^2),mean((predict(model.gam.shift.w08$model,type="terms")[,3]+mean(f3(x3))-f3(x3))^2),
mean((predict(model.gam.shift.w08$model,type="terms")[,4]+mean(f4(x4))-f4(x4))^2))

 
 # average absolute errors of the fitted value
 abs.err.S[i] = mean(abs(fitted.gam.S[i,]-true.fct)) 
 abs.err.w02.S[i] = mean(abs(fitted.gam.w02.S[i,]-true.fct))
 abs.err.w04.S[i] = mean(abs(fitted.gam.w04.S[i,]-true.fct))
 abs.err.w05.S[i] = mean(abs(fitted.gam.w05.S[i,]-true.fct))
 abs.err.w06.S[i] = mean(abs(fitted.gam.w06.S[i,]-true.fct))
 abs.err.w08.S[i] = mean(abs(fitted.gam.w08.S[i,]-true.fct))
 
 # average squared errors of the fitted value
 sq.err.S[i] = mean((fitted.gam.S[i,]-true.fct)^2) 
 sq.err.w02.S[i] = mean((fitted.gam.w02.S[i,]-true.fct)^2)
 sq.err.w04.S[i] = mean((fitted.gam.w04.S[i,]-true.fct)^2)
 sq.err.w05.S[i] = mean((fitted.gam.w05.S[i,]-true.fct)^2)
 sq.err.w06.S[i] = mean((fitted.gam.w06.S[i,]-true.fct)^2)
 sq.err.w08.S[i] = mean((fitted.gam.w08.S[i,]-true.fct)^2)
 


if (i/10==round(i/10)) print(i)

save(list=ls(all=T),file="REMLprop05.RData")

}