################################################################################################
#####R CODE FOR VIRGINIA TECH DEPARTMENT OF STATISTICS TECHNICAL REPORT NO. 05-7           #####
#####"ROBUST PARAMETER DESIGN: A SEMI-PARAMETRIC APPROACH"                                 #####
#####STEPHANIE M. PICKLE, TIMOTHY J. ROBINSON, JEFFREY B. BIRCH, CHRISTINE M. ANDERSON-COOK#####
################################################################################################

#The following code will perform the parametric, nonparametric, and semi-parametric approach to the printing ink data#

#Read Ink Data in from a file 
#Note that data is coded such that each regressor is between 0 and 1
ink.data.coded<-read.table('FileName',sep=",",header=TRUE)

#Create separate vectors for each column of data matrix
x1<-ink.data.coded[,1]			#x1=speed
x2<-ink.data.coded[,2]			#x2=pressure
x3<-ink.data.coded[,3]			#x3=distance
X.design<-cbind(x1,x2,x3)		#design matrix
ybar<-ink.data.coded[,4]		#mean response
s<-ink.data.coded[,5]			#standard deviation
t<-ink.data.coded[,6]			#transformation of standard deviation = log(s^2+1)

#############################################################################################################################
#############################################################################################################################
#############################################################################################################################

#############################
#####PARAMETRIC APPROACH#####
#############################

#The following code will estimate responses parametrically using OLS for the transformed variance and EWLS for the mean#

#######################
#Preliminary Functions#
#######################
#Leave one out Cross-validation 
yhat.minusi.p.fn<-function(y,res,H){
	n<-length(y)
	h<-diag(H)	#grabs diagonals from H
	yhat.minusi<-rep(0,n)
	for(i in 1:n){
		yhat.minusi[i]<-y[i]-(res[i]/(1-h[i]))
	}	
	return(yhat.minusi)
}
obj.p.fn<-function(beta.mean,beta.t,theta,x0){
	#This function will find the estimated MSE=(yhat-theta)^2+varhat
	#at the location x0 based on the estimated ols functions for the
	#transformed variance and wls functions for mean
	#here the mean model is second order 
	#and the t model is first order
	x1.x2<-x0[1]*x0[2]
	x1.x3<-x0[1]*x0[3]
	x2.x3<-x0[2]*x0[3]
	x1.x1<-x0[1]*x0[1]
	x2.x2<-x0[2]*x0[2]
	x3.x3<-x0[3]*x0[3]
	modelx0.mean<-c(1,x0,x1.x2,x1.x3,x2.x3,x1.x1,x2.x2,x3.x3)
	modelx0.t<-c(1,x0)
	yhat<-beta.mean%*%modelx0.mean
	bias<-yhat-theta
	that<-beta.t%*%modelx0.t
	varhat<-exp(that)-1
	msehat<-bias^2+varhat	
	msehat
}
ga.p.fn<-function(beta.mean,beta.t,theta){
	#This function will perform a Simple Genetic Algorithm	
	#For the objective function estimated MSE=(yhat-theta)^2+var
	#Where yhat is the estimate of the mean and theta is the target value
	#And var is the estimate of the variance
	#var will be found with ols and ybar will be found with wls

	#Stopping Criteria
	maxiter<-10000			#set maximum number of iterations
	satiter<-1000			#set number of iterations during which the best results keep saturating
	epsln<-1e-8			#smallest gain worth recognizing

	#GA Parameters
	m<-50				#population size at each generation
	mutrate<-.2			#mutation rate
	crossrate<-.9			#crossover rate
	zerorate<-.2			#zero rate
	extremerate<-.2			#extreme rate
	keep<-2				#keep top two from each generation to remain unchanged
	num.parents<-m-keep		#number of population that will be used for mating
	mate<-ceiling((num.parents)/2)	#number of matings
	
	#Create initial population and initialize 'best result' holders
	xrange<-matrix(c(0,1,0,1,0,1),2,3)
	k<-ncol(xrange)
	elite<-matrix(0,keep,k)		#matrix for elite chromosomes to remain unchanged
	cross<-matrix(0,num.parents,k)	#matrix for crossover	
	iter<-0				#generation (iteration) counter
	stopcode<-0
	inarow<-0			#number of iterations with function's value consecutively less than epsln
	bestfun<-1e30			#essential positive infinity
	bestx<-matrix(0,1,k)		
	f<-rep(0,m)			#initialize function value for each of m vectors
	G<-matrix(runif(m*k,0,1),m,k)	#initial generation with population size m

	#####Iterate through generations#####
	while(stopcode==0){
		iter<-iter+1		#increments counter
		if(iter>maxiter) stopcode<-2	#loop will exit on stopcode=2 for exceeding the maximum number of iterations

		#Evaluate current generation
		for(i in 1:m){
			f[i]<-obj.p.fn(beta.mean,beta.t,theta,G[i,])
			}

		mat<-cbind(G,f)			#create matrix of x locations and corresponding function values
		newmat<-mat[order(mat[,4]),]	#sort matrix of by function values increasing
		minf<-newmat[1,4]		#optimal function (minimum)
		bf<-min(minf, bestfun)
		fgain<-bestfun-bf	#fgain is always non-negative it measure the change
		if(fgain>epsln)	inarow<-0 else inarow<-(inarow+1)

		if(fgain>0){
			bestfun<-bf
			bestx<-newmat[1,1:3]		#optimal x location
			}
 	
		if(inarow>satiter) stopcode<-1	#loop will exit on stopcode=1 for best result having been achieved

		#Select elite to remain unchanged
		elite<-newmat[1:keep,1:3]
	
		#Select parents for mating
		cross<-newmat[(keep+1):m,1:3]
		
		#Do Crossover
		for(i in 1:mate){
			z<-rbinom(1,1,crossrate)
			if(z==1){
				x<-ceiling(runif(1,0,(k-1)))
				temp<-cross[i,(x+1):k]
				cross[i,(x+1):k]<-cross[i+mate,(x+1):k]
				cross[i+mate,(x+1):k]<-temp
				}
			}
	
		#Do Mutation
		M<-matrix(runif((num.parents*k),0,1),num.parents,k)
		for(i in 1:num.parents){
			for(j in 1:k){
				zz<-rbinom(1,1,mutrate)
				if(zz==1) cross[i,j]<-M[i,j]
				}
			}
		#Do Zero Gene Operator
		for(i in 1:num.parents){
			for(j in 1:k){
				zzz<-rbinom(1,1,zerorate)
				if(zzz==1) cross[i,j]<-0
				}
			}
		#Do Extreme Values Operator
		Extreme<-matrix(1,num.parents,k)
		for(i in 1:num.parents){
			for(j in 1:k){
				zzzz<-rbinom(1,1,extremerate)
				if(zzzz==1) cross[i,j]<-Extreme[i,j]
				}
			}
		G<-rbind(elite,cross)
	}
	
	if(iter<maxiter) status<-0
	else status<-1

	result<-c(bestx,bestfun,iter)
	return(result)
}


#######################
#Create model matrices#
#######################
ones<-rep(1,length(x1))			#intercept column
X.model.1st<-cbind(ones,X.design)	#first-order model matrix

#TwoWay Interactions
x1.x2<-x1*x2
x1.x3<-x1*x3
x2.x3<-x2*x3

#Quadratic Terms
x1.x1<-x1^2
x2.x2<-x2^2
x3.x3<-x3^2

X.model.2nd<-cbind(X.model.1st,x1.x2,x1.x3,x2.x3,x1.x1,x2.x2,x3.x3)	#second-order model matrix

##########################################
#FIRST-ORDER OLS FOR TRANSFORMED VARIANCE#
##########################################
t.ols<-lm(t~x1+x2+x3)
betas.t.ols<-t.ols$coefficients
betas.t.ols
fit.t.ols<-t.ols$fitted.values
fit.t.ols
res.t.ols<-t.ols$residuals
res.t.ols
sum.t.ols<-summary(t.ols)
sum.t.ols
anova(t.ols)
par(mfrow=c(2,2))
plot(t.ols,pch=16)
mtext("Residual Plots for First-Order Model For Variance Transformation Using OLS",outer=T,line=-2)

#Hat matrix for OLS
H.t.ols<-X.model.1st%*%(solve(t(X.model.1st)%*%X.model.1st))%*%t(X.model.1st)
H.t.ols

#Get variance weights for means model
varhat.ols<-exp(fit.t.ols)-1
varhat.ols
wts.t.ols<-1/varhat.ols
wts.t.ols
W.t.ols<-diag(wts.t.ols)

#Leave one out Cross-validation 
fit.minusi.t.ols<-yhat.minusi.p.fn(t,res.t.ols,H.t.ols)
fit.minusi.t.ols
res.minusi.t.ols<-t-fit.minusi.t.ols
res.minusi.t.ols

############################
#SECOND-ORDER EWLS FOR MEAN#
############################
mean.ewls<-lm(ybar~x1+x2+x3+x1.x2+x1.x3+x2.x3+x1.x1+x2.x2+x3.x3,weights=wts.t.ols)
betas.mean.ewls<-mean.ewls$coefficients
betas.mean.ewls
fit.mean.ewls<-mean.ewls$fitted.values
fit.mean.ewls
res.mean.ewls<-mean.ewls$residuals
res.mean.ewls
sum.mean.ewls<-summary(mean.ewls)
sum.mean.ewls
anova(mean.ewls)
par(mfrow=c(2,2))
plot(mean.ewls,pch=16)
mtext("Residual Plots for Second-Order Model For Mean Using EWLS",outer=T,line=-2)

#Hat matrix for EWLS
H.mean.ewls<-X.model.2nd%*%(solve(t(X.model.2nd)%*%W.t.ols%*%X.model.2nd))%*%t(X.model.2nd)%*%W.t.ols
H.mean.ewls

##################
#OPTIMIZE WITH GA#
##################
#Mean target is 500
ga.p.fn(betas.mean.ewls,betas.t.ols,500)

#############################################################################################################################
#############################################################################################################################
#############################################################################################################################

################################
#####NONPARAMETRIC APPROACH#####
################################

#The following code will estimate responses nonparametrically using LLR for the transformed variance and EWLLR for the mean#

#######################
#Preliminary Functions#
#######################
#Get Kernel Weights at x0
kern.fn<-function(X,y,b,x0){
	#This function will find the normal kernel weights at the location x0
	#Given the design matrix X, the observed values y, and the bandwidth b
	n<-length(y)			#number of obervations
	p<-ncol(X)			#number of regressors
	u<-matrix(0,n,p)		#matrix of arguement for kernel function with n rows and p columns
	kern.norm<-matrix(0,n,p)	#matrix for kernel function with n rows and p columns
	K.norm<-rep(0,n)		#vector of mult. norm. kernel
	wt<-rep(0,n)			#vector of kernel weights with n rows
	for(j in 1:p){
		for(i in 1:n){
			u[i,j]<-(x0[j]-X[i,j])/b	
			kern.norm[i,j]<-exp(-(u[i,j]^2))	#normal kernel
		}
	}
	K.norm<-kern.norm[,1]*kern.norm[,2]*kern.norm[,3]	#multiplicative normal kernel for p=3 regressors	
	for(i in 1:n){
			wt[i]<-K.norm[i]/sum(K.norm)		#kernel weights
	}
	return(wt)
}
#Get LPR Estimate at x0
lpr.fn<-function(X,y,b,d,x0){
	#This function will find the Local Polynomial Regression estimate at the location x0
	#Using the normal kernel for the weights
	#Given the design matrix X, the observed values y, the bandwidth b, and the degree d
	#If d = 0 a constant will be fit; i.e., Kernel Regression
	#If d = 1 a simple linear regression line will be fit; i.e., Local Linear Regression
	wt<-kern.fn(X,y,b,x0)	#kernel weight vector
	yhat<-0			
	if(d==0){
		model<-lm(y~1,weight=wt)
		beta<-model$coefficients
		beta0<-beta[1]
		beta1<-0
		beta2<-0
		beta3<-0
	}else{model<-lm(y~X,weight=wt)
		beta<-model$coefficients
		beta0<-beta[1]
		beta1<-beta[2]
		beta2<-beta[3]
		beta3<-beta[4]
	}

	yhat<-beta0+beta1*x0[1]+beta2*x0[2]+beta3*x0[3]
	return(yhat)
}
#Get Weighted LPR Estimate at x0
lpr.w.fn<-function(X,y,b,d,x0,W){
	#This function will find the Weigthed Local Polynomial Regression estimate at the location x0
	#Using the normal kernel for the weights
	#Given the design matrix X, the observed values y, the bandwidth b, and the degree d
	#and the estimated weight matrix from variance model W
	#If d = 0 a constant will be fit; i.e., Kernel Regression
	#If d = 1 a simple linear regression line will be fit; i.e., Local Linear Regression
	k<-kern.fn(X,y,b,x0)	#kernel weight vector
	k.half<-sqrt(k)
	K.half<-diag(k.half)
	W.star<-K.half%*%W%*%K.half
	wt<-diag(W.star)
	yhat<-0			
	if(d==0){
		model<-lm(y~1,weight=wt)
		beta<-model$coefficients
		beta0<-beta[1]
		beta1<-0
		beta2<-0
		beta3<-0
	}else{model<-lm(y~X,weight=wt)
		beta<-model$coefficients
		beta0<-beta[1]
		beta1<-beta[2]
		beta2<-beta[3]
		beta3<-beta[4]
	}
	yhat<-beta0+beta1*x0[1]+beta2*x0[2]+beta3*x0[3]
	return(yhat)
}
#Fit original data points with LPR
fit.lpr.fn<-function(X,y,b,d){
	#This function will fit the original data points by local polynomial regression with a degree (d) polynomial
	#Given the design matrix X, the observed y values, and the bandwidth b
	n<-length(y)	#number of observations
	X0<-X		#fit at original data points
	fit<-rep(0,n)	
	for(i in 1:n){
		fit[i]<-lpr.fn(X,y,b,d,X0[i,])
	}
	return(fit)
}
#Fit original data points with weighted LPR
fit.lpr.w.fn<-function(X,y,b,d,W){
	#This function will fit the original data points by local polynomial regression with a degree (d) polynomial
	#Given the design matrix X, the observed y values, and the bandwidth b
	#and the estimated weight matrix from variance model W
	n<-length(y)	#number of observations
	X0<-X		#fit at original data points
	fit<-rep(0,n)	
	for(i in 1:n){
		fit[i]<-lpr.w.fn(X,y,b,d,X0[i,],W)
	}
	return(fit)
}
#Leave one out Cross-validation 
fit.lpr.minusi.fn<-function(X,y,b,d){
	#This function will fit the original data points using minus i technique for LPR
	#Given the design matrix X, the observed values y, the bandwidth b, the degree of polynomial d
	n<-length(y)
	X0<-X
	fit.minusi<-rep(0,n)
	for(i in 1:n){
		fit.minusi[i]<-lpr.fn(X[-i,],y[-i],b,d,X0[i,])
	}
	return(fit.minusi)
}
#Weighted Leave one out Cross-validation 
fit.lpr.w.minusi.fn<-function(X,y,b,d,wt){
	#This function will fit the original data points using minus i technique for Weigthed LPR
	#Given the design matrix X, the observed values y, the bandwidth b, the degree of polynomial d
	#and the estimated weight matrix from variance model W
	num<-length(y)
	X0<-X
	fit.w.minusi<-rep(0,num)
	W.minusi<-array(0,dim=c(num-1,num-1,num))
	for(i in 1:num){
		W.minusi[,,i]<-diag(wt[-i])
	}
	for(i in 1:num){	
		fit.w.minusi[i]<-lpr.w.fn(X[-i,],y[-i],b,d,X0[i,],W.minusi[,,i])
	}
	return(fit.w.minusi)
}
#Get row of Hat Matrix for x0
Hrow.fn<-function(X,y,b,d,x0){
	#This function will find the row of the Hat matrix corresponding to x0 such that Hrow*y=yhat at x0
	#Given the design matrix X, the observed values y, the bandwidth b, the degree d of the LPR estimate
	#If d = 0 kernel regression
	#If d = 1 LLR	
	n<-length(y)			#number of observations
	k<-ncol(X)			#number of regressors
	p<-k+1				#number of parameters
	ones<-rep(1,n)
	modelX<-cbind(ones,X)
	modelx0<-c(1,x0)
	Hrow<-rep(0,n)
	wt.kern<-kern.fn(X,y,b,x0)	#vector of kernel weights with n rows
	kern.mat<-diag(wt.kern)		#diagonal matrix of kernel weights
	XpWX<-t(modelX)%*%kern.mat%*%modelX
	inv<-solve(XpWX)
	if(d==0){
		Hrow<-wt.kern
		}else{
			Hrow<-modelx0%*%inv%*%t(modelX)%*%kern.mat
		}
	return(Hrow)
}
#Get row of weighted Hat Matrix for x0
Hrow.w.fn<-function(X,y,b,d,x0,W){
	#This function will find the row of the Hat matrix corresponding to x0 such that Hrow*y=yhat at x0
	#Given the design matrix X, the observed values y, the bandwidth b, the degree d of the LPR estimate
	#and the estimated weight matrix from variance model W
	#If d = 0 kernel regression
	#If d = 1 LLR	
	n<-length(y)			#number of observations
	k<-ncol(X)			#number of regressors
	p<-k+1				#number of parameters
	ones<-rep(1,n)
	modelX<-cbind(ones,X)
	modelx0<-c(1,x0)
	Hrow.w<-rep(0,n)
	wt.kern<-kern.fn(X,y,b,x0)	#vector of kernel weights with n rows
	wt.kern.half<-sqrt(wt.kern)
	kern.half.mat<-diag(wt.kern.half)	#diagonal matrix of half kernel weights
	W.mat<-kern.half.mat%*%W%*%kern.half.mat 	#weight matrix
	XpWX<-t(modelX)%*%W.mat%*%modelX
	inv<-solve(XpWX)
	if(d==0){
		Hrow.w<-diag(W.mat)
		}else{
			Hrow.w<-modelx0%*%inv%*%t(modelX)%*%W.mat
		}
	return(Hrow.w)
}
#Hat matrix for lpr
H.lpr.fn<-function(X,y,b,d){
	n<-length(y)		#number of observations
	X0<-X			#fit at original data points
	H<-matrix(0,n,n)	#create H matrix
	for(i in 1:n){
		H[i,]<-Hrow.fn(X,y,b,d,X0[i,])
	}
	return(H)
}
#Hat matrix for weighted lpr
H.lpr.w.fn<-function(X,y,b,d,W){
	n<-length(y)		#number of observations
	X0<-X			#fit at original data points
	H.w<-matrix(0,n,n)	#create H matrix
	for(i in 1:n){
		H.w[i,]<-Hrow.w.fn(X,y,b,d,X0[i,],W)
	}
	return(H.w)
}
#Get Trace of square matrix
trace.fn<-function(X){
	#This function will find the trace of a square matrix
	n<-nrow(X)
	diag<-rep(0,n)
	for(i in 1:n){
		diag[i]<-X[i,i]
	}
	trace<-sum(diag)
	return(trace)
}
#Get SSE
sse.fn<-function(y,yhat){
	#This function will find SSE given observed values y and fitted values yhat
	res<-y-yhat
	sse<-t(res)%*%res
	return(sse)
}
#Get weighted SSE
sse.w.fn<-function(y,yhat,W){
	#This function will find weighted SSE given observed values y, fitted values yhat
	#and the estimated weight matrix from variance model W
	res<-y-yhat
	sse.w<-t(res)%*%W%*%res
	return(sse.w)
}
#Get PRESS** for one bandwidth value
press_star_star.b.fn<-function(X,y,b,d){
	#This function will find the PRESS** value corresponding to a given bandwidth value(b)
	#Given the design matrix X, the observed values y, and the degree (d) of the polynomial for LPR
	n<-length(y)			#number of observations
	k<-ncol(X)			#number of regressors
	p<-k+1				#number of parameters
	X0<-X				#fit at original data points
	H<-matrix(0,n,n)		#create H matrix
	Htrace<-0			#create variable for trace of H values
	yhat<-rep(0,n)			#create vector for yhat values
	press<-rep(0,n)			#create vector for PRESS values
	press_star_star<-0		#create variable for PRESS** value
	for(i in 1:n){
		H[i,]<-Hrow.fn(X,y,b,d,X0[i,])
	}
	Htrace<-trace.fn(H)
	yhat.minusi<-fit.lpr.minusi.fn(X,y,b,d)
	yhat<-fit.lpr.fn(X,y,b,d)
	res<-y-yhat
	SSEb<-t(res)%*%res
	if(d==0){
		MSE<-(summary(lm(y~1))$sigma)^2
		dfE<-summary(lm(y~1))$df[2]
		SSEmax<-MSE*dfE
	}else{MSE<-(summary(lm(y~X))$sigma)^2
		dfE<-summary(lm(y~X))$df[2]
		SSEmax<-MSE*dfE
	}
	for(i in 1:n){
		press[i]<-(y[i]-yhat.minusi[i])^2
	}
	press_star_star<-sum(press)/(n-Htrace+(n-p)*((SSEmax-SSEb)/SSEmax))
	return(press_star_star)
}
#Get weighted PRESS** for one bandwidth value
press_star_star.b.w.fn<-function(X,y,b,d,wt){
	#This function will find the PRESS** value corresponding to a given bandwidth value(b)
	#Given the design matrix X, the observed values y, and the degree (d) of the polynomial for LPR
	#and the estimated weights wt
	n<-length(y)			#number of observations
	k<-ncol(X)			#number of regressors
	p<-k+1				#number of parameters
	X0<-X				#fit at original data points
	H.w<-matrix(0,n,n)		#create H matrix
	Htrace<-0			#create variable for trace of H values
	yhat<-rep(0,n)			#create vector for yhat values
	press.w<-rep(0,n)		#create vector for PRESS values
	press_star_star.w<-0		#create variable for PRESS** value
	W<-diag(wt)			#Create weight matrix
	for(i in 1:n){
		H.w[i,]<-Hrow.w.fn(X,y,b,d,X0[i,],W)
	}
	Htrace<-trace.fn(H.w)
	yhat.w<-fit.lpr.w.fn(X,y,b,d,W)
	res.w<-y-yhat.w	
	SSEb.w<-t(res.w)%*%W%*%res.w
	yhat.minusi.w<-fit.lpr.w.minusi.fn(X,y,b,d,wt)	
	if(d==0){
		MSE<-(summary(lm(y~1,weights=wt))$sigma)^2
		dfE<-summary(lm(y~1,weights=wt))$df[2]
		SSEmax<-MSE*dfE
	}else{MSE<-(summary(lm(y~X,weights=wt))$sigma)^2
		dfE<-summary(lm(y~X,weights=wt))$df[2]
		SSEmax<-MSE*dfE
	}
	for(i in 1:n){
		press.w[i]<-((y[i]-yhat.minusi.w[i])^2)*wt[i]
	}
	press_star_star.w<-sum(press.w)/(n-Htrace+(n-p)*((SSEmax-SSEb.w)/SSEmax))
	return(press_star_star.w)
}
#Get Optimal Bandwidth from Candidate List use with coded x's (0,1)
band.press_star_star.fn<-function(X,y,band,d){
	#This function will find the optimal bandwidth that minimizes PRESS**
	#given a candidate list of bandwidth values (band)
	m<-length(band)
	press_star_star.value<-rep(1e10,m)
	tol<-rep(0,m)
	press_star_star.value[1]<-press_star_star.b.fn(X,y,band[1],d)
	limit<-0.01
	for(i in 2:m){
		press_star_star.value[i]<-press_star_star.b.fn(X,y,band[i],d)
		tol[i]<-abs(press_star_star.value[i]-press_star_star.value[i-1])/press_star_star.value[i-1]
		if(tol[i]<=limit){break}
	}
	mat<-cbind(band,press_star_star.value,tol)	#create matrix of all b and corresponding PRESS** values
	newmat<-mat[order(mat[,2]),]		#sort matrix of b and PRESS** values by PRESS** increasing
	bopt<-newmat[1,1]			#optimal bandwidth
	press_star_star_opt<-newmat[1,2]	#corresponding optimal PRESS** (minimum)
	return(bopt,press_star_star_opt)	
}
#Get weighted PRESS** Optimal Bandwidth from Candidate List use with coded x's (0,1)
band.press_star_star.w.fn<-function(X,y,band,d,wt){
	#This function will find the optimal bandwidth that minimizes PRESS**
	#given a candidate list of bandwidth values (band)
	#and the estimated weights wt
	m<-length(band)
	press_star_star.w.value<-rep(1e10,m)
	tol<-rep(0,m)
	press_star_star.w.value[1]<-press_star_star.b.w.fn(X,y,band[1],d,wt)
	limit<-0.01
	for(i in 2:m){
		press_star_star.w.value[i]<-press_star_star.b.w.fn(X,y,band[i],d,wt)
		tol[i]<-abs(press_star_star.w.value[i]-press_star_star.w.value[i-1])/press_star_star.w.value[i-1]
		if(tol[i]<=limit){break}
	}
	mat<-cbind(band,press_star_star.w.value,tol)	#create matrix of all b and corresponding PRESS** values
	newmat<-mat[order(mat[,2]),]		#sort matrix of b and PRESS** values by PRESS** increasing
	bopt<-newmat[1,1]			#optimal bandwidth
	press_star_star_opt<-newmat[1,2]	#corresponding optimal PRESS** (minimum)
	return(bopt,press_star_star_opt)	
}
#Get Plot of PRESS** vs Bandwidth from Candidate List
band.press_star_star.plot.fn<-function(X,y,band,d){
	#This function plot PRESS** values vs bandwidth values
	#given a candidate list of bandwidth values (band)
	m<-length(band)
	press_star_star.value<-rep(0,m)
	for(i in 1:m){
		press_star_star.value[i]<-press_star_star.b.fn(X,y,band[i],d)
	}
	mat<-cbind(band,press_star_star.value)	#create matrix of all b and corresponding PRESS** values
	newmat<-mat[order(mat[,2]),]		#sort matrix of b and PRESS** values by PRESS** increasing
	bopt<-newmat[1,1]			#optimal bandwidth
	press_star_star_opt<-newmat[1,2]	#corresponding optimal PRESS** (minimum)
	plot(mat[,1],mat[,2],xlab="Bandwidth",ylab="PRESS**",type="l")	
}
#Get Plot of Weighted PRESS** vs Bandwidth from Candidate List
band.press_star_star.plot.w.fn<-function(X,y,band,d,wt){
	#This function plot PRESS** values vs bandwidth values
	#given a candidate list of bandwidth values (band)
	#and the estimated weights wt
	m<-length(band)
	press_star_star.w.value<-rep(0,m)
	for(i in 1:m){
		press_star_star.w.value[i]<-press_star_star.b.w.fn(X,y,band[i],d,wt)
	}
	mat<-cbind(band,press_star_star.w.value)	#create matrix of all b and corresponding PRESS** values
	newmat<-mat[order(mat[,2]),]		#sort matrix of b and PRESS** values by PRESS** increasing
	bopt<-newmat[1,1]			#optimal bandwidth
	press_star_star_opt<-newmat[1,2]	#corresponding optimal PRESS** (minimum)
	plot(mat[,1],mat[,2],xlab="Bandwidth",ylab="PRESS**",type="l")	
}
obj.llr.fn<-function(X,y,theta,b1,d1,t,b2,d2,x0,W){
	#This function will find the estimated MSE=(yhat-theta)^2+varhat
	#Given the data (X), the mean response (y), the target value for
	#the mean (theta), the bandwidth for the estimation of the mean (b1), the 
	#degree of the lpr for the mean (d1), the transformed variance (t), 
	#the bandwidth for the estimation of the transformed variance (b2), the
	#degree of the lpr for the transformed variance, and the location of
	#estimation (x0)

	yhat<-lpr.w.fn(X,y,b1,d1,x0,W)
	bias<-yhat-theta
	that<-lpr.fn(X,t,b2,d2,x0)
	varhat<-exp(that)-1
	msehat<-bias^2+varhat
	msehat
}

ga.llr.fn<-function(X,y,theta,b1,d1,t,b2,d2,W){
	#This function will perform a Simple Genetic Algorithm	
	#For the objective function estimated MSE=(yhat-theta)^2+var
	#Given the data (X), the mean response (y), the target value for
	#the mean (theta), the bandwidth for the estimation of the mean (b1), the 
	#degree of the lpr for the mean (d1), the transformed variance (t), 
	#the bandwidth for the estimation of the transformed variance (b2), the
	#degree of the lpr for the transformed variance

	#Stopping Criteria
	maxiter<-10000			#set maximum number of iterations
	satiter<-1000			#set number of iterations during which the best results keep saturating
	epsln<-1e-8			#smallest gain worth recognizing
	
	#GA Parameters
	m<-50				#population size at each generation
	mutrate<-.2			#mutation rate
	crossrate<-.9			#crossover rate
	zerorate<-.2			#zero rate
	extremerate<-.2			#extreme rate
	keep<-2				#keep top two from each generation to remain unchanged
	num.parents<-m-keep		#number of population that will be used for mating
	mate<-ceiling((num.parents)/2)	#number of matings
	
	#Create initial population and initialize 'best result' holders
	xrange<-matrix(c(0,1,0,1,0,1),2,3)
	k<-ncol(xrange)
	elite<-matrix(0,keep,k)		#matrix for elite chromosomes to remain unchanged
	cross<-matrix(0,num.parents,k)	#matrix for crossover	
	iter<-0				#generation (iteration) counter
	stopcode<-0
	inarow<-0			#number of iterations with function's value consecutively less than epsln
	bestfun<-1e30			#essential positive infinity
	bestx<-matrix(0,1,k)		
	f<-rep(0,m)			#initialize function value for each of m vectors
	G<-matrix(runif(m*k,0,1),m,k)	#initial generation with population size m

	#####Iterate through generations#####
	while(stopcode==0){
		iter<-iter+1		#increments counter
		if(iter>maxiter) stopcode<-2	#loop will exit on stopcode=2 for exceeding the maximum number of iterations

		#Evaluate current generation
		for(i in 1:m){
			f[i]<-obj.llr.fn(X,y,theta,b1,d1,t,b2,d2,G[i,],W)
			}

		mat<-cbind(G,f)			#create matrix of x locations and corresponding function values
		newmat<-mat[order(mat[,4]),]	#sort matrix of by function values increasing
		minf<-newmat[1,4]		#optimal function (minimum)
		bf<-min(minf, bestfun)
		fgain<-bestfun-bf	#fgain is always non-negative it measure the change
		if(fgain>epsln)	inarow<-0 else inarow<-(inarow+1)

		if(fgain>0){
			bestfun<-bf
			bestx<-newmat[1,1:3]		#optimal x location
			}
 	
		if(inarow>satiter) stopcode<-1	#loop will exit on stopcode=1 for best result having been achieved

		#Select elite to remain unchanged
		elite<-newmat[1:keep,1:3]
	
		#Select parents for mating
		cross<-newmat[(keep+1):m,1:3]
		
		#Do Crossover
		for(i in 1:mate){
			z<-rbinom(1,1,crossrate)
			if(z==1){
				x<-ceiling(runif(1,0,(k-1)))
				temp<-cross[i,(x+1):k]
				cross[i,(x+1):k]<-cross[i+mate,(x+1):k]
				cross[i+mate,(x+1):k]<-temp
				}
			}
	
		#Do Mutation
		M<-matrix(runif((num.parents*k),0,1),num.parents,k)
		for(i in 1:num.parents){
			for(j in 1:k){
				zz<-rbinom(1,1,mutrate)
				if(zz==1) cross[i,j]<-M[i,j]
				}
			}
		#Do Zero Gene Operator
		for(i in 1:num.parents){
			for(j in 1:k){
				zzz<-rbinom(1,1,zerorate)
				if(zzz==1) cross[i,j]<-0
				}
			}
		#Do Extreme Values Operator
		Extreme<-matrix(1,num.parents,k)
		for(i in 1:num.parents){
			for(j in 1:k){
				zzzz<-rbinom(1,1,extremerate)
				if(zzzz==1) cross[i,j]<-Extreme[i,j]
				}
			}
		G<-rbind(elite,cross)
	}
	
	if(iter<maxiter) status<-0
	else status<-1

	result<-c(bestx,bestfun,iter)
	return(result)
}

##############################
#LLR FOR TRANSFORMED VARIANCE#
##############################
#Find optimal bandwidth from candidate list based on coded x's
par(mfrow=c(1,1))
band<-seq(.1,1,.01)
band.press_star_star.plot.fn(X.design,t,band,1)
band<-seq(.3,1,.01)
bopt_t<-band.press_star_star.fn(X.design,t,band,1)$bopt
bopt_t

#Find LLR Estimates for transformed variance
fit.t.llr<-fit.lpr.fn(X.design,t,bopt_t,1)
fit.t.llr
res.t.llr<-t-fit.t.llr
res.t.llr

#Leave one out Cross-validation 
fit.minusi.t.llr<-fit.lpr.minusi.fn(X.design,t,bopt_t,1)
fit.minusi.t.llr
res.minusi.t.llr<-t-fit.minusi.t.llr
res.minusi.t.llr

H.t.llr<-H.lpr.fn(X.design,t,bopt_t,1)
H.t.llr

#Get variance weights for means model 
varhat.llr<-exp(fit.t.llr)-1
varhat.llr
wts.t.llr<-1/varhat.llr
wts.t.llr
W.t.llr<-diag(wts.t.llr)

################################
#####Weigthted LLR FOR MEAN#####
################################
#Find optimal bandwidth from candidate list based on coded x's
par(mfrow=c(1,1))
band<-seq(.1,1,.01)
band.press_star_star.plot.w.fn(X.design,ybar,band,1,wts.t.llr)
band<-seq(.3,1,.01)
bopt_mean<-band.press_star_star.w.fn(X.design,ybar,band,1,wts.t.llr)$bopt
bopt_mean

#Find weighted LPR Estimates for mean
fit.mean.llr.w<-fit.lpr.w.fn(X.design,ybar,bopt_mean,1,W.t.llr)
fit.mean.llr.w
res.mean.llr.w<-ybar-fit.mean.llr.w
res.mean.llr.w

#Weighted Leave one out Cross-validation 
fit.minusi.w.mean.llr<-fit.lpr.w.minusi.fn(X.design,ybar,bopt_mean,1,wts.t.llr)
fit.minusi.w.mean.llr
res.minusi.w.mean.llr<-ybar-fit.minusi.w.mean.llr
res.minusi.w.mean.llr

H.mean.llr.w<-H.lpr.w.fn(X.design,ybar,bopt_mean,1,W.t.llr)
H.mean.llr.w

##################
#OPTIMIZE WITH GA#
##################
#Mean target is 500
ga.llr.fn(X.design,ybar,500,bopt_mean,1,t,bopt_t,1,W.t.llr)

#############################################################################################################################
#############################################################################################################################
#############################################################################################################################

##################################
#####SEMI-PARAMETRIC APPROACH#####
##################################

#The following code will estimate responses semi-parametrically using MRR1 for the transformed variance and MRR2 for the mean#

#######################
#Preliminary Functions#
#######################
lambda.mrr1.fn<-function(y,yhat.p,yhat.np,yhat.minusi.p,yhat.minusi.np){
	#This will find asymptotically optimal data driven lambda for MRR1
	diff.minusi<-yhat.minusi.np-yhat.minusi.p
	res.p<-y-yhat.p
	diff<-yhat.np-yhat.p
	numerator<-t(diff.minusi)%*%res.p
	denominator<-t(diff)%*%diff
	lambda.mrr1<--(numerator/denominator)
	return(lambda.mrr1)
}
lambda.mrr2.fn<-function(y,yhat.p,rhat.llr){
	#This will find asymptotically optimal data driven lambda for MRR2
	res.p<-y-yhat.p	
	numerator<-t(rhat.llr)%*%res.p
	denominator<-t(rhat.llr)%*%rhat.llr
	lambda.mrr2<-numerator/denominator
	opt<-min(lambda.mrr2,1)
	return(opt)
}
obj.dmrr.fn<-function(X,y,beta.mean,theta,b1,d1,l1,t,beta.t,b2,d2,l2,x0){
	#This function will find the estimated MSE=(yhat-theta)^2+varhat
	#Given the data (X), the mean response (y), the coeffiecents for means wls model (beta.mean)
	#the target value for the mean (theta), the bandwidth for the mean residuals (b1), the 
	#degree of the lpr for the mean (d1), the lambda for mean (l1) the transformed variance (t),
	#the coefficients for variance ols model 
	#the bandwidth for the estimation of the transformed variance (b2), the
	#degree of the lpr for the transformed variance (d2), the lambda for variance (l2),
	#and the location of estimation (x0)
	x1.x2<-x0[1]*x0[2]
	x1.x3<-x0[1]*x0[3]
	x2.x3<-x0[2]*x0[3]
	x1.x1<-x0[1]*x0[1]
	x2.x2<-x0[2]*x0[2]
	x3.x3<-x0[3]*x0[3]
	modelx0.mean<-c(1,x0,x1.x2,x1.x3,x2.x3,x1.x1,x2.x2,x3.x3)
	modelx0.t<-c(1,x0)
	yhat.p<-beta.mean%*%modelx0.mean
	that.p<-beta.t%*%modelx0.t
	
	res.yhat<-y-yhat.p
	reshat.llr<-lpr.fn(X,res.yhat,b1,d1,x0)
	that.llr<-lpr.fn(X,t,b2,d2,x0)
	
	yhat.mmrr<-yhat.p+l1*reshat.llr
	that.vmrr<-l2*that.llr+(1-l2)*that.p

	bias<-theta-yhat.mmrr
	varhat<-exp(that.vmrr)-1
	msehat<-bias^2+varhat
	msehat

}

ga.dmrr.fn<-function(X,y,beta.mean,theta,b1,d1,l1,t,beta.t,b2,d2,l2){
	#This function will perform a Simple Genetic Algorithm	
	#For the objective function estimated MSE=(yhat-theta)^2+var
	#Given the data (X), the mean response (y), the coeffiecents for means wls model (beta.mean)
	#the target value for the mean (theta), the bandwidth for the mean residuals (b1), the 
	#degree of the lpr for the mean (d1), the lambda for mean (l1) the transformed variance (t),
	#the coefficients for variance ols model 
	#the bandwidth for the estimation of the transformed variance (b2), the
	#degree of the lpr for the transformed variance (d2), the lambda for variance (l2)

	#Stopping Criteria
	maxiter<-10000			#set maximum number of iterations
	satiter<-1000			#set number of iterations during which the best results keep saturating
	epsln<-1e-8			#smallest gain worth recognizing
	
	#GA Parameters
	m<-50				#population size at each generation
	mutrate<-.2			#mutation rate
	crossrate<-.9			#crossover rate
	zerorate<-.2			#zero rate
	extremerate<-.2			#extreme rate
	keep<-2				#keep top two from each generation to remain unchanged
	num.parents<-m-keep		#number of population that will be used for mating
	mate<-ceiling((num.parents)/2)	#number of matings
	
	#Create initial population and initialize 'best result' holders
	xrange<-matrix(c(0,1,0,1,0,1),2,3)
	k<-ncol(xrange)
	elite<-matrix(0,keep,k)		#matrix for elite chromosomes to remain unchanged
	cross<-matrix(0,num.parents,k)	#matrix for crossover	
	iter<-0				#generation (iteration) counter
	stopcode<-0
	inarow<-0			#number of iterations with function's value consecutively less than epsln
	bestfun<-1e30			#essential positive infinity
	bestx<-matrix(0,1,k)		
	f<-rep(0,m)			#initialize function value for each of m vectors
	G<-matrix(runif(m*k,0,1),m,k)	#initial generation with population size m

	#####Iterate through generations#####
	while(stopcode==0){
		iter<-iter+1		#increments counter
		if(iter>maxiter) stopcode<-2	#loop will exit on stopcode=2 for exceeding the maximum number of iterations

		#Evaluate current generation
		for(i in 1:m){
			f[i]<-obj.dmrr.fn(X,y,beta.mean,theta,b1,d1,l1,t,beta.t,b2,d2,l2,G[i,])
			}

		mat<-cbind(G,f)			#create matrix of x locations and corresponding function values
		newmat<-mat[order(mat[,4]),]	#sort matrix of by function values increasing
		minf<-newmat[1,4]		#optimal function (minimum)
		bf<-min(minf, bestfun)
		fgain<-bestfun-bf	#fgain is always non-negative it measure the change
		if(fgain>epsln)	inarow<-0 else inarow<-(inarow+1)

		if(fgain>0){
			bestfun<-bf
			bestx<-newmat[1,1:3]		#optimal x location
			}
 	
		if(inarow>satiter) stopcode<-1	#loop will exit on stopcode=1 for best result having been achieved

		#Select elite to remain unchanged
		elite<-newmat[1:keep,1:3]
	
		#Select parents for mating
		cross<-newmat[(keep+1):m,1:3]
		
		#Do Crossover
		for(i in 1:mate){
			z<-rbinom(1,1,crossrate)
			if(z==1){
				x<-ceiling(runif(1,0,(k-1)))
				temp<-cross[i,(x+1):k]
				cross[i,(x+1):k]<-cross[i+mate,(x+1):k]
				cross[i+mate,(x+1):k]<-temp
				}
			}
	
		#Do Mutation
		M<-matrix(runif((num.parents*k),0,1),num.parents,k)
		for(i in 1:num.parents){
			for(j in 1:k){
				zz<-rbinom(1,1,mutrate)
				if(zz==1) cross[i,j]<-M[i,j]
				}
			}
		#Do Zero Gene Operator
		for(i in 1:num.parents){
			for(j in 1:k){
				zzz<-rbinom(1,1,zerorate)
				if(zzz==1) cross[i,j]<-0
				}
			}
		#Do Extreme Values Operator
		Extreme<-matrix(1,num.parents,k)
		for(i in 1:num.parents){
			for(j in 1:k){
				zzzz<-rbinom(1,1,extremerate)
				if(zzzz==1) cross[i,j]<-Extreme[i,j]
				}
			}
		G<-rbind(elite,cross)
	}
	
	if(iter<maxiter) status<-0
	else status<-1

	result<-c(bestx,bestfun,iter)
	return(result)
}

#######################################
#####MRR1 FOR TRANSFORMED VARIANCE#####
#######################################
#Find optimal mixing parameter
lambda.opt_t<-c(lambda.mrr1.fn(t,fit.t.ols,fit.t.llr,fit.minusi.t.ols,fit.minusi.t.llr))
lambda.opt_t

#Find VMRR Estimates
fit.t.vmrr<-lambda.opt_t*fit.t.llr+(1-lambda.opt_t)*fit.t.ols
fit.t.vmrr
res.t.vmrr<-t-fit.t.vmrr
res.t.vmrr

#Find H_vmrr
H.t.vmrr<-lambda.opt_t*H.t.llr+(1-lambda.opt_t)*H.t.ols
H.t.vmrr

#Get variance weights for MMRR
varhat.t.vmrr<-exp(fit.t.vmrr) - 1
varhat.t.vmrr
wts.t.vmrr<-1/varhat.t.vmrr
wts.t.vmrr
W.t.vmrr<-diag(wts.t.vmrr)

#######################
#####MRR2 FOR MEAN#####
#######################
#LLR FOR Residuals#
#Find optimal bandwidth from candidate list based on coded x's
par(mfrow=c(1,1))
band<-seq(.1,1,.01)
band.press_star_star.plot.fn(X.design,res.mean.ewls,band,1)
band<-seq(.3,1,.01)
bopt_res<-band.press_star_star.fn(X.design,res.mean.ewls,band,1)$bopt
bopt_res

#Find LLR Estimates for residuals
fit.res.llr<-fit.lpr.fn(X.design,res.mean.ewls,bopt_res,1)
fit.res.llr
res.res.llr<-res.mean.ewls-fit.res.llr
res.res.llr

#Find H_res
H.res.llr<-H.lpr.fn(X.design,res.mean.ewls,bopt_res,1)
H.res.llr

#Semi-parametric Fits Using Means Model Robust Regression 2#
#Find optimal mixing parameter
lambda.opt_mean<-c(lambda.mrr2.fn(ybar,fit.mean.ewls,fit.res.llr))
lambda.opt_mean

#Find MMRR Estimates 
fit.mean.mmrr<-fit.mean.ewls+lambda.opt_mean*fit.res.llr
fit.mean.mmrr
res.mean.mmrr<-ybar-fit.mean.mmrr
res.mean.mmrr

#Find H mmrr
I<-diag(1,27,27)
H.mean.mmrr<-H.mean.ewls+lambda.opt_mean*H.res.llr%*%(I-H.mean.ewls)
H.mean.mmrr

##################
#OPTIMIZE WITH GA#
##################
#Mean target is 500
ga.dmrr.fn(X.design,ybar,betas.mean.ewls,500,bopt_res,1,lambda.opt_mean,t,betas.t.ols,bopt_t,1,lambda.opt_t)


#############################################################################################################################
#############################################################################################################################
#############################################################################################################################