### 5.Simulation.R

##########################
########################## Define model y = f + e
##########################

################# 1. Choose mean function f
m <- 3; meanf <- c("heav", "dopp", "ppoly", "fg12")

################# 2. Define error term e
# if phi <- 0.0 then generate Gaussian noise
# if phi <- 0.5 then generate AR(0.5) noise
# if phi <- "t" then generate t(3) noise
phi <- 0.5 
if(phi != 0.0 && phi != "t") vscale <- "level" else vscale <- "independent"

################# 3. Define the model y = f + e
# with sample size=n and signal-to-noise-ratio=snr
snr <- 7        
nlevel <- 10            
n <- 2^nlevel
testfun <- get(paste(meanf[m]))
testdata <- testfun(n)
sigma <- testdata$sdf / snr
x <- testdata$x

################# 4. Define cross-validation scheme
bsize <- 4; kfold <- 8

################# 5. Define cross-validation levels
optlevel <- 3

################# 6. Comparing the MSE of EBayes, cv.Nason and cv.level
# Define the number of replication, noofcase
noofcase <- 10 # 100
mse <- NULL
for (i in 1:noofcase) {
    if(phi != 0.0 && phi != "t")
        error <- as.numeric(arima.sim(n=n, list(ar=phi, ma=0.0))) * 
                               sqrt(1 - phi^2) * sigma
    if(phi == 0.0)
        error <- rnorm(n, 0, sigma)
    if(phi == "t")
        error <- rt(n, 3) / sqrt(3) * sigma
    y <- testdata$meanf + error
    ywd <- wd(y) 
    #If the wavethresh4 BETA is installed, run following
    #assign(paste("yn.", i, sep=""), wr(threshold(ywd, policy="cv")))
    assign(paste("ye.", i, sep=""), 
           wr(ebayesthresh.wavelet(ywd, vscale=vscale, smooth.levels=nlevel-3)))
    assign(paste("out.", i, sep=""), 
        cvwavelet(y=y, ywd=ywd, cv.optlevel=optlevel, cv.bsize=bsize, 
                      cv.kfold=kfold, impute.vscale=vscale))
                                     
    mse <- rbind(mse, 
                   c(mean((get(paste("ye.", i, sep=""))-testdata$meanf)^2),
                     #If the wavethresh4 BETA is installed, run following
                     #mean((get(paste("yn.", i, sep=""))-testdata$meanf)^2),
                     mean((get(paste("out.", i, sep=""))$yc-testdata$meanf)^2)))          
}

par(mfrow=c(1,1))    
boxplot(data.frame(mse), main=meanf[m], names=c("EBayes", "cv.level")) 
#If the wavethresh4 BETA is installed, run following
#boxplot(data.frame(mse), main=meanf[m], names=c("EBayes", "cv.Nason", "cv.level"))

################# 7. Plotting the reconstructions of EBayes, cv.Nason and cv.level
checking.plot <- T
par(mfrow=c(2,2), oma=c(0,0,0,0), mar=0.1+c(4,4,2.0,1)) 
for (i in 1:noofcase) {
 tmp <- get(paste("out.", i, sep=""))
 plot(x, tmp$y, main="Observation", xlab="", ylab="", cex=0.3)
 plot(x, get(paste("ye.", i, sep="")), main="EBayes", xlab="", ylab="", type="l")
 #If the wavethresh4 BETA is installed, run following
 #plot(x, get(paste("yn.", i, sep="")), main="cv.Nason", xlab="", ylab="", type="l")
 plot(x, tmp$yc, main="cv.level", xlab="", ylab="", type="l")
 if(checking.plot) locator(1)
}