# Comparison Study Between Gamma Kernel Regression and Chaubey et al. (2010)'s Method # Chaubey's Method: Exponeential Design set.seed(9999); S=matrix(rep(0,12),nrow=3) for(k in seq(200)) { Smse=matrix(rep(0,12),nrow=3) ri=1; for(n in c(100,150,200)) { ci=1; for(model in c(1,2,3,4)) { # Model and Data x=rexp(n,1); y1=x+2*exp(-16*x^2)+rnorm(n,0,0.5); y2=x+2*exp(-16*x^2)+rexp(n,sqrt(2)/0.5)-rexp(n,sqrt(2)/0.5); y3=sin(2*x)+2*exp(-16*x^2)+rnorm(n,0,0.5); y4=sin(2*x)+2*exp(-16*x^2)+rexp(n,sqrt(2)/0.5)-rexp(n,sqrt(2)/0.5); y=y1*(model==1)+y2*(model==2)+y3*(model==3)+y4*(model==4) # Estimated MSE Using Leave-one-Out CV Criterion mse=function(aa) { vn=aa; en=aa^2; xij=matrix(kronecker(x,vn^2*(x+en),"/"),nrow=n); Xiy=exp(-xij)%*%diag(y*x^{1/vn^2-1}) Xi=exp(-xij)%*%diag(x^{1/vn^2-1}) Asum=apply(Xiy,1,sum)-x^{1/vn^2-1}*exp(-x/(vn^2*(x+en)))*y; Bsum=apply(Xi,1,sum)-x^{1/vn^2-1}*exp(-x/(vn^2*(x+en))); mean((y-Asum/Bsum)^2) } # Searching for optimal smoothing parameter a=seq(0.001,1,length=200); fmse=rep(0,200); for(i in seq(200)) { fmse[i]=mse(a[i]) } Smse[ri,ci]=min(fmse[fmse!="NaN"]) ci=ci+1; } ri=ri+1; } S=S+Smse cat(k,"MSE=",S/k,"\n") } S=S/200 dimnames(S)=list(c(100,150,200),c("model 1","model 2","model 3","model 4")) S #===================================================================================================== # Comparison Study Between Gamma Kernel Regression and Chaubey et al. (2010)'s Method # Song's Method: Exponential Design set.seed(9999); S=matrix(rep(0,12),nrow=3) for(k in seq(200)) { Smse=matrix(rep(0,12),nrow=3) ri=1; for(n in c(100,150,200)) { ci=1; for(model in c(1,2,3,4)) { x=rexp(n,1); y1=x+2*exp(-16*x^2)+rnorm(n,0,0.5); y2=x+2*exp(-16*x^2)+rexp(n,sqrt(2)/0.5)-rexp(n,sqrt(2)/0.5); y3=sin(2*x)+2*exp(-16*x^2)+rnorm(n,0,0.5); y4=sin(2*x)+2*exp(-16*x^2)+rexp(n,sqrt(2)/0.5)-rexp(n,sqrt(2)/0.5); y=y1*(model==1)+y2*(model==2)+y3*(model==3)+y4*(model==4) smse=function(h) { xij=matrix(kronecker(x,x/h,"^"),nrow=n); Xiy=xij%*%diag(y*exp(-x/h)); Xi=xij%*%diag(exp(-x/h)); Asum=apply(Xiy,1,sum)-x^{x/h}*exp(-x/h)*y; Bsum=apply(Xi,1,sum)-x^{x/h}*exp(-x/h); mean((y-Asum/Bsum)^2); } hseq=seq(0.001,1,length=200); fsmse=rep(0,length(hseq)) i=1; for(h in hseq) { fsmse[i]=smse(h) i=i+1 } Smse[ri,ci]=min(fsmse[fsmse!="NaN"]) ci=ci+1; } ri=ri+1; } S=S+Smse } S=S/200 dimnames(S)=list(c(100,150,200),c("model 1","model 2","model 3","model 4")) S