SLR2_2.py

###########
# Created to run a spares linear regression analysis.
# Methods that perform elastic-net (lasso) regression,
# uses bootstrap sampling for penalty selection,
# uses bootstrap residual for determining statistic distributions
# uses permutation analysis to determine null distribution,
# estimates importance by a single step forward and backward
# selection. 1) Full feature matrix used for Hyperparameter selecion,
# 1b) once selected regression is performed on full dataset with full feature 
# matrix, 2) do bsr on selected varriables to id stats, 3) repeate from 
# the (1) for permutations.

# Created 20120330 RAT
# Editied 20120418 RAT
# Editied 20120731 RAT (added a simplified method)
#
###########
import elasticNetLinReg as enet
from glmnet import glmnet
import numpy as np
import math
import gpdPerm
import cvTools as st
from scipy.sparse import lil_matrix


def run(X,y,name,nPerms=1000,nSamp=100,alphaList=np.array([1]),reduceX=True,reselect=True):
	""" Run the full estimation / permutation routine 
	and then save the data to a file.
	"""
	if nPerms>0:
		solution, enm = permModel(X,y,nSamp,alphaList,nPerms,reselect=reselect,reduceX=reduceX)
	else:
		solution, enm = estModel(X,y,nSamp,alphaList,estImp=True,reduceX=reduceX)

	# we use the stDev latter so lets just get it now
	solution['sdY'] = np.sqrt(np.var(y))
	nObs,nRegs = X.shape
	sdX = np.zeros(nRegs)
	sdX[solution['indices']] = np.sqrt(np.var(X[:,solution['indices']],0))
	solution['sdX'] = sdX


	# we have it all, lets print it	
	
	f = open(name,'w')
	
	enm.lambdas[0].tofile(f,sep="\t")
	f.write("\n")
	
	enm.alpha.tofile(f,sep="\t")
	f.write("\n")

	enm.intercept[0].tofile(f,sep="\t")	
	f.write("\n")

	solution['aveErr'].tofile(f,sep="\t")
	f.write("\n")

	solution['sdErr'].tofile(f,sep="\t")
	f.write("\n")

	solution['aveNullErr'].tofile(f,sep="\t")
	f.write("\n")

	solution['sdNullErr'].tofile(f,sep="\t")
	f.write("\n")
	
	solution['sdY'].tofile(f,sep="\t")
	f.write("\n")	
	
	ind = solution['indices']

	if len(ind)>0:
		ind.tofile(f,sep="\t")
		f.write("\n")
	
		solution['sdX'][ind].tofile(f,sep="\t")
		f.write("\n")

		# need to get the coefs
		coef = np.zeros(nRegs)
		# note due to the bs residual, we may have 
		# non zero means for coefs with a zero global 
		# value, so lets set that up.
		coef[enm.indices] = enm.coef
		coef[ind].tofile(f,sep="\t")
		f.write("\n")
		
		solution['medCoef'][ind].tofile(f,sep="\t")
		f.write("\n")

		solution['sdCoef'][ind].tofile(f,sep="\t")
		f.write("\n")

		solution['pSup'][ind].tofile(f,sep="\t")
		f.write("\n")

		solution['errOut'][ind].tofile(f,sep="\t")
		f.write("\n")

		solution['errIn'][ind].tofile(f,sep="\t")
		f.write("\n")

		
		if nPerms > 0:
			solution['p'][ind].tofile(f,sep="\t")
			f.write("\n")



	f.close()


def permModelSimple2(X,y,nSamp=100,alphaList=np.array([1]),nPermsMin=1000,nPermsMax=100000,reselect=True,estImp=True):
	"""Fits the data to linear model using specified elastic net param 
	(defulat is 1, ie LASSO).  The penalty is specified by bootstrap permutations
	with nSamp (reselect determines if the penalty should be re-estimated over 
	the permutations or if the full model value should be used).  Permutations
	are done to find the permutation coef (key = 'medPermCoef') which is used 
	to estimate the p-value (key = 'p').  
	NOTE: in this version we do not calculate the standard error estimate over the 
	permutations, therfore we do not scale the coef, so the test statistic is simply
	the coefficent itself. 
	"""
	## ok this is a cut and paste job, 
	## some varriable names are not great (ie I still use the name tStat when its just
	## the abs of the coef and not really the tStat, but I think all is correct
	## in the technical sense that it does what I think it does.
	nObs,nRegs = X.shape
	solution, enm = estModel(X,y,nSamp,alphaList,estImp=estImp)
	medCoef = solution['medCoef']
	aveCoef = solution['aveCoef']
	sdCoef = solution['sdCoef']
	indices = solution['indices']
	solution['coef'] = np.zeros(nRegs)
	solution['coef'][enm.indices] = enm.coef
	lam = enm.lambdas[0]
	alpha = enm.alpha
	done = False
	p = np.ones(nRegs)
	medPermCoef = np.zeros(nRegs)
	nPerms = nPermsMin
	if len(indices)>0:
		while (nPerms<=nPermsMax and not done):
			tStat = np.zeros(nRegs)
			tStat[enm.indices] = np.abs(enm.coef)
			tStatPerm = lil_matrix((nRegs,nPerms))
			for i in range(nPerms):
				# permute the response 
				# *** probably should keep track to avoid repeats, future???
				yPerm = np.random.permutation(y)
				if reselect:
					enmPerm = select(X,yPerm,nSamp,alphaList) 
				else:
					enmPerm = enet.fit(X,yPerm,alpha,lambdas=[lam])[0]

				indPerm = enmPerm.indices
				
				
				if len(indPerm)>0:
					tmp = np.abs(enmPerm.coef)
					# more crzy shift cuz the dif from 1-d array in np and scipy
					tStatPerm[indPerm,i] = np.array(tmp,ndmin=2).T 
			#np.savetxt('tStat.dat',tStat)
			#np.savetxt('tStatPerm.dat',np.array(tStatPerm.todense()))		


				
			
			done = True 
			for i in range(nRegs):
				# even more confusion for scpy and np arrays
				# gdpPerm is expecting a vector which is diffrent 
				# from an nx1 matrix (apperently) 
				curTStatPerm = np.array(tStatPerm[i,:].todense())[0,:]
				medPermCoef[i] = np.median(curTStatPerm)
				p[i] = gpdPerm.est(tStat[i],curTStatPerm)
				# use standard permutation if this fails 
				if np.isnan(p[i]) or p[i] == 0:
					done = False
					tmp = np.sum(curTStatPerm>=tStat[i])+1
					p[i] = float(tmp)/(float(nPerms))
					if p[i]>1.0:p[i]=1.0
						
			nPerms = nPerms*10

	solution['p'] = p
	solution['medPermCoef'] = medPermCoef
	return solution, enm


def permModelSimple(X,y,nSamp=100,alphaList=np.array([1]),nPerms=1000,reselect=True):
	"""Fits the data to linear model using specified elastic net param 
	(defulat is 1, ie LASSO).  The penalty is specified by bootstrap permutations
	with nSamp (reselect determines if the penalty should be re-estimated over 
	the permutations or if the full model value should be used).  Permutations
	are done to find the permutation coef (key = 'medPermCoef') which is used 
	to estimate the p-value (key = 'p').  
	NOTE: in this version we do not calculate the standard error estimate over the 
	permutations, therfore we do not scale the coef, so the test statistic is simply
	the coefficent itself. 
	"""
	## ok this is a cut and paste job, 
	## some varriable names are not great (ie I still use the name tStat when its just
	## the abs of the coef and not really the tStat, but I think all is correct
	## in the technical sense that it does what I think it does.
	nObs,nRegs = X.shape
	solution, enm = estModel(X,y,nSamp,alphaList,estImp=True)
	medCoef = solution['medCoef']
	aveCoef = solution['aveCoef']
	sdCoef = solution['sdCoef']
	indices = solution['indices']
	solution['coef'] = np.zeros(nRegs)
	solution['coef'][enm.indices] = enm.coef
	lam = enm.lambdas[0]
	alpha = enm.alpha
	p = np.ones(nRegs)
	medPermCoef = np.zeros(nRegs)
	if len(indices)>0:
		tStat = np.zeros(nRegs)
		tStat[enm.indices] = np.abs(enm.coef)
		tStatPerm = lil_matrix((nRegs,nPerms))
		for i in range(nPerms):
			# permute the response 
			# *** probably should keep track to avoid repeats, future???
			yPerm = np.random.permutation(y)
			if reselect:
				enmPerm = select(X,yPerm,nSamp,alphaList) 
			else:
				enmPerm = enet.fit(X,yPerm,alpha,lambdas=[lam])[0]

			indPerm = enmPerm.indices
			
			
			if len(indPerm)>0:
				tmp = np.abs(enmPerm.coef)
				# more crzy shift cuz the dif from 1-d array in np and scipy
				tStatPerm[indPerm,i] = np.array(tmp,ndmin=2).T 
		#np.savetxt('tStat.dat',tStat)
		#np.savetxt('tStatPerm.dat',np.array(tStatPerm.todense()))		


			
		p = np.ones(nRegs)
		for i in range(nRegs):
			# even more confusion for scpy and np arrays
			# gdpPerm is expecting a vector which is diffrent 
			# from an nx1 matrix (apperently) 
			curTStatPerm = np.array(tStatPerm[i,:].todense())[0,:]
			medPermCoef[i] = np.median(curTStatPerm)
			p[i] = gpdPerm.est(tStat[i],curTStatPerm)
			# use standard permutation if this fails 
			if np.isnan(p[i]) or p[i] == 0:
				tmp = np.sum(curTStatPerm>=tStat[i])+1
				p[i] = float(tmp)/(float(nPerms))
				if p[i]>1.0:p[i]=1.0

	solution['p'] = p
	solution['medPermCoef'] = medPermCoef
	return solution, enm



def permModel(X,y,nSamp=100,alphaList=np.array([1]),nPerms=1000,reselect=True,reduceX=False):
	"""This does it all folks, runs the full model 
	with bs samp for model selection, at the selected model 
	runs the bs res on the full feature set X
	to estimate the sdDev and means for various values 
	(ie coef), while doing that it uses 10 fold CV at each 
	sample to estimate prediction errors.  It then does all
	of that nPerms time to do a full permutation analysis,
	and uses these values to estimate tStats and finaly 
	p-values for the coef (this is imporved, when possible,
	by using the GPD approx.  Also we estimate importance 
	scores for the various coef, but this is not done over
	the random permutations.  If reselect then the bs 
	selection process is reimplimented at each permutation
	otherwise model selected on non permuted values is 
	used.
	"""
	nObs,nRegs = X.shape
	solution, enm = estModel(X,y,nSamp,alphaList,estImp=True)
	medCoef = solution['medCoef']
	aveCoef = solution['aveCoef']
	sdCoef = solution['sdCoef']
	indices = solution['indices']
	solution['coef'] = np.zeros(nRegs)
	solution['coef'][enm.indices] = enm.coef
	lam = enm.lambdas[0]
	alpha = enm.alpha
	p = np.ones(nRegs)
	if len(indices)>0:
		sdCoef[sdCoef<1E-21] = 1E-21
		tStat = np.abs(medCoef/sdCoef)
		tStatPerm = lil_matrix((nRegs,nPerms))
		for i in range(nPerms):
			# permute the response 
			# *** probably should keep track to avoid repeats, future???
			yPerm = np.random.permutation(y)
			if reselect:
				solutionPerm, _= estModel(X,yPerm,nSamp,alphaList,estErr=False,reduceX=reduceX)
			else:
				solutionPerm, _= estModel(X,yPerm,nSamp,alphaList,estErr=False,reduceX=reduceX,params=(lam,alpha))
			medCoefPerm = solutionPerm['medCoef']
			sdCoefPerm = solutionPerm['sdCoef']
			indPerm = solutionPerm['indices']
			
			if len(indPerm)>0:
				sdCoefPerm[sdCoefPerm<1E-21] = 1E-21
				tmp = np.abs(medCoefPerm/sdCoefPerm)
				# more crzy shift cuz the dif from 1-d array in np and scipy
				tStatPerm[indPerm,i] = np.array(tmp[indPerm],ndmin=2).T 
		#np.savetxt('tStat.dat',tStat)
		#np.savetxt('tStatPerm.dat',np.array(tStatPerm.todense()))		


		# no values should have 2 in the end,
		# this will let us know if something goes wrong
		p = np.ones(nRegs)
		for i in range(nRegs):
			# even more confusion for scpy and np arrays
			# gdpPerm is expecting a vector which is diffrent 
			# from an nx1 matrix (apperently) 
			curTStatPerm = np.array(tStatPerm[i,:].todense())[0,:]
			p[i] = gpdPerm.est(tStat[i],curTStatPerm)
			# use standard permutation if this fails 
			if np.isnan(p[i]) or p[i] < 1E-21:
				tmp = np.sum(curTStatPerm>=tStat[i])+1
				p[i] = float(tmp)/(float(nPerms))
				if p[i]>1.0:p[i]=1.0

	solution['p'] = p
	return solution, enm



def estModel(XFull,y,nSamp=100,alphaList=np.array([1]),estErr=True,estImp=False,reduceX=False,params=[]):
	"""Estimate a mean and standard deviation
	for an elastic net model using bootstrap 
	residual.
	Note: Bootstrap resampling is used to select
	model parameters, then the bs res at these 
	params is used on the full feature set X
	to calculate means and standard errors.
	Note: if estErr then 10 fold CV is used to estimate 
	the prediction error at each iteration of the bs.
	This is ten extra iterations at each bs res 
	sample, but reduces the bias in prediction error.
	The mean and sdDev of the CV error is then reported.
	Note: If params are passed then we assume its a tuple
	with the (lambda,alpha) model parameters.  In this case 
	model selection is bipassed. and these params are used.
	"""

	nObs,nRegsFull = XFull.shape
	# select full model values
	if len(params)==2:
		lam,alpha = params
		enm = enet.fit(XFull,y,alpha,lambdas=[lam])[0]
	else:
		enm = select(XFull,y,nSamp,alphaList)
	lam = enm.lambdas[0]
	yHat = enm.predict(XFull)
	intercept = enm.intercept[0]
	globalCoef =enm.coef[np.abs(enm.coef)>1E-21]
	coefIndex = enm.indices[np.abs(enm.coef)>1E-21]
	alpha = enm.alpha

	# now is when we reduce the x if we need too!
	if reduceX:
		nRegs = len(coefIndex)
		if nRegs > 0:
			X = XFull[:,coefIndex]
			nObs, _ = X.shape
	else:
		X = XFull
		nRegs = nRegsFull

	# get the bootstrap residual response samples
	res = y - yHat
	resCent = res-np.mean(res)
	ySample = np.zeros((nObs,nSamp))
	for i in range(nSamp):
		resSample = st.sampleWR(resCent)
		ySample[:,i] = yHat+resSample

	if nRegs > 0:
	
		# residual bs time
		if estErr:
			sumErr = 0
			sumSqErr = 0
			sumNullErr = 0
			sumSqNullErr = 0

		sc = np.zeros(nRegs)
		sSqc = np.zeros(nRegs)
		ac = lil_matrix((nRegs,nSamp))
		sumSup = np.zeros(nRegs)
		

		for i in range(nSamp):
			# cv to get the errors
			if estErr:
				err,tmpEnm,tmpallVals = fitSampling(X,ySample[:,i],alpha,10,method='cv',lambdas=[lam])
				sumErr = err.mErr[0] + sumErr
				sumSqErr = err.mErr[0]**2 + sumSqErr
				# cv over this thing to get the null model errors
				nullErr,a = fitSamplingNull(ySample[:,i],10, method='cv')
				sumNullErr = sumNullErr + nullErr
				sumSqNullErr = sumSqNullErr + nullErr**2

			# need the coef
			# they change so we need to map the back to the original
			tmpEnm = enet.fit(X,ySample[:,i], alpha,lambdas=[lam])
			sc[tmpEnm.indices] = sc[tmpEnm.indices] + tmpEnm.coef[:,0]
			sSqc[tmpEnm.indices] = sSqc[tmpEnm.indices] + tmpEnm.coef[:,0]**2
			if len(tmpEnm.indices)>0:
				ac[tmpEnm.indices,i] = tmpEnm.coef
			# find supports 
			occur = np.zeros(len(tmpEnm.coef[:,0]))
			occur[abs(tmpEnm.coef[:,0])>1E-25] = 1.0
			sumSup[tmpEnm.indices] = sumSup[tmpEnm.indices] + occur
				

		# get averages and variances
		if estErr:
			aveErr = sumErr/nSamp
			sdErr = np.sqrt(sumSqErr/nSamp - aveErr**2)
			aveNullErr = sumNullErr/nSamp
			sdNullErr = np.sqrt(sumSqNullErr/nSamp - aveNullErr**2)

		aveCoef = sc/nSamp
		sdCoef = np.sqrt(sSqc/nSamp - aveCoef**2)
		#some crazy stuff here becase of the way scipy mat is shaped
		medCoef = np.array(np.median(ac.todense(),1))[:,0]
		pSup = sumSup/nSamp
		indices = np.arange(nRegs)[np.abs(medCoef)>1E-21]
		# put it in a dict for simplicity 
		solution = {}
		if estErr:
			solution['aveErr'] = aveErr
			solution['sdErr'] = sdErr
			solution['aveNullErr'] = aveNullErr
			solution['sdNullErr'] = sdNullErr
		if reduceX:
			# need to go back to the original indicies 
			solution['aveCoef'] = np.zeros(nRegsFull)
			solution['sdCoef'] = np.zeros(nRegsFull)
			solution['medCoef'] = np.zeros(nRegsFull)
			solution['pSup'] = np.zeros(nRegsFull)

			solution['aveCoef'][coefIndex] = aveCoef
			solution['sdCoef'][coefIndex] = sdCoef
			solution['medCoef'][coefIndex] = medCoef
			solution['pSup'][coefIndex] = pSup
			solution['indices'] = coefIndex[indices]
		else:
			solution['aveCoef'] = aveCoef
			solution['sdCoef'] = sdCoef
			solution['medCoef'] = medCoef
			solution['pSup'] = pSup
			solution['indices'] = indices
		
		nRegsHat = len(indices)
		if nRegsHat>0 and estImp:
			Xhat = X[:,indices]
			# lets do the leave one out importance deal
			errOutHat = np.zeros(nRegsHat) 
			if nRegsHat>1:
				for j in range(nRegsHat):
					Xprime = np.delete(Xhat,j,axis=1)

					# residual bs time
					sumErr = 0
					sumSqErr = 0
					
					for i in range(nSamp):
						# cv to get the errors
						err,tmpenm,tmpallVals = fitSampling(Xprime,ySample[:,i],alpha,10,method='cv',lambdas=[lam])
						sumErr = err.mErr[0] + sumErr
						sumSqErr = err.mErr[0]**2 + sumSqErr

					errOutHat[j] = sumErr/nSamp

			elif nRegsHat==1:
				errOutHat[0] = aveNullErr

			# lets do leave only one
			errInHat = np.zeros(nRegsHat) 
			for j in range(nRegsHat):
				Xprime = np.zeros((nObs,1))
				Xprime[:,0] = Xhat[:,j]

				# residual bs time
				sumErr = 0
				sumSqErr = 0
				
				for i in range(nSamp):
					# cv to get the errors
					err,tmpenm,tmpallVals = fitSampling(Xprime,ySample[:,i],alpha,10,method='cv',lambdas=[lam])
					sumErr = err.mErr[0] + sumErr
					sumSqErr = err.mErr[0]**2 + sumSqErr

				errInHat[j] = sumErr/nSamp

			errOut = np.zeros(nRegs)
			errOut[indices] = errOutHat
			solution['errOut'] = errOut
			errIn = np.zeros(nRegs)
			errIn[indices] = errInHat
			solution['errIn'] = errIn



	else:
			solution = {}
			if estErr:
				sumNullErr = 0
				sumSqNullErr = 0
			
				for i in range(nSamp):
					# cv over this thing to get the null model errors
					nullErr,a = fitSamplingNull(ySample[:,i],10, method='cv')
					sumNullErr = sumNullErr + nullErr
					sumSqNullErr = sumSqNullErr + nullErr**2
			
				# get averages and variances
				aveNullErr = sumNullErr/nSamp
				sdNullErr = np.sqrt(sumSqNullErr/nSamp - aveNullErr**2)
				aveErr = aveNullErr
				sdErr = sdNullErr
				solution['aveErr'] = aveErr
				solution['sdErr'] = sdErr
				solution['aveNullErr'] = aveNullErr
				solution['sdNullErr'] = sdNullErr

			solution['aveCoef'] = np.zeros(nRegsFull)
			solution['sdCoef'] = np.zeros(nRegsFull)
			solution['medCoef'] = np.zeros(nRegsFull)
			solution['pSup'] = np.zeros(nRegsFull)
			solution['indices'] = np.array([])

		
	
			

	return solution, enm 







def select(X,y,nSamp=100,alphaList=np.array([1])):
	"""Select an elastic net model based 
	on a resampling method and return that
	model.
	"""
	nObs,nRegs = X.shape
	sdY = np.sqrt(np.var(y))
	# selection via bootstrap
	bestMin = 1E10
	for a in alphaList:
		tmpErr,tmpEnm,allVals = fitSampling(X,y,a,nSamp,method='bs')
		tmpErrV = tmpErr.mErr
		tmpMin = np.min(tmpErrV)
		
		if tmpMin < bestMin:
			bestMin = tmpMin
			modelIndex = np.argmin(tmpErrV)
			enm = tmpEnm
			err = tmpErr
			alpha = a
	
	# important values
		
	return enm[modelIndex]





def fitSampling(regressors, response, alpha, nSamp, method='cv', 
		memlimit=None, largest=None, **kwargs):
	"""Performs an elastic net constrained linear regression,
	see fit, with selected sampleing method to estimate errors
	using nSamp number of sampleings.
	methods:
	'cv'	cross validation with nSamp number of folds
	'bs'	bootstrap 
	'bs632'	boostrap 632 (weighted average of bs and training error)
	Returns a TrainingError object (cvTools) and an 
	ENetModel object for the full fit (err,enm).
	Function requires cvTools
	"""
	
	nObs,nRegs = regressors.shape
	# get the full model fit 
	fullEnm = enet.fit(regressors, response, alpha, memlimit,
                largest, **kwargs)
	# get the lambda values determined in the full fit (going to force these lambdas for all cv's)
	lam = fullEnm.lambdas
	# the lambdas may have been user defined, don't want it defined twice 
	if kwargs.has_key('lambdas'):
		del kwargs['lambdas']

	# lets partition the data via our sampling method
	if method=='cv':
		t,v = st.kFoldCV(range(nObs),nSamp,randomise=True)
	elif (method=='bs') or (method=='bs632'):
		t,v = st.kRoundBS(range(nObs),nSamp)
	else:
		raise ValueError('Sampling method not correct')

	# lets consider many versions of errors
	# with our error being mean squared error
	# we want the epected mean squared error
	# and the corisponding variance over the diffrent versions
	nModels = len(lam)
	smse = np.zeros(nModels)
	sSqmse = np.zeros(nModels)
	allVals = np.zeros((nModels,nSamp))

		# loop through the folds
	for i in range(nSamp):
		# get the training values
		X = regressors[t[i]]
		y = response[t[i]]
		enm =  enet.fit(X, y, alpha, memlimit,
                	largest, lambdas=lam, **kwargs)
			# get the validation values
		Xval = regressors[v[i]]
		Yval = response[v[i]]
		nVal = float(len(Yval))
		# get the predicted responses from validation regressors
		Yhat = enm.predict(Xval)
				# what is the mean squared error?
		# notice the T was necassary to do the subtraction
		# the rows are the models and the cols are the observations
		mse = np.sum((Yhat.T-Yval)**2,1)/nVal
		# sum the rows (errors for given model)
		smse = smse + mse
		sSqmse = sSqmse + mse**2
		allVals[:,i] = mse
		
	# now it is time to average and send back
	# I am putting the errors in a container 
	nSampFlt = float(nSamp)
	meanmse = smse/nSampFlt
	varmse = sSqmse/nSampFlt - meanmse**2
	if method=='bs632':
		yhat = fullEnm.predict(regressors)
		resubmse = np.sum((yhat.T-response)**2,1)/float(nObs)
		meanmse = 0.632*meanmse+(1-0.632)*resubmse
		
	err = enet.ENetTrainError(lam,nSamp,meanmse,varmse,[0],[0],alpha)
	err.setParamName('lambda')

	fullEnm.setErrors(err.mErr)
	
	return err, fullEnm, allVals 

def fitSamplingNull(response,nSamp, method='cv', 
		memlimit=None, largest=None, **kwargs):
	nObs = len(response)
	# lets partition the data via our sampling method
	if method=='cv':
		t,v = st.kFoldCV(range(nObs),nSamp,randomise=True)
	elif (method=='bs') or (method=='bs632'):
		t,v = st.kRoundBS(range(nObs),nSamp)
	else:
		raise ValueError('Sampling method not correct')

	smse = 0
	sSqmse = 0

	for i in range(nSamp):
		# get the training values
		
		y = response[t[i]]
		Yval = response[v[i]]
		nVal = float(len(Yval))
		mse = np.sum((Yval-np.mean(y))**2)/nVal
		# sum the rows (errors for given model)
		smse = smse + mse
		sSqmse = sSqmse + mse**2
		
	# now it is time to average and send back
	# I am putting the errors in a container 
	nSampFlt = float(nSamp)
	meanmse = smse/nSampFlt
	varmse = sSqmse/nSampFlt - meanmse**2
	if method=='bs632':
		yhat = fullEnm.predict(regressors)
		resubmse = np.sum((yhat.T-response)**2,1)/float(nObs)
		meanmse = 0.632*meanmse+(1-0.632)*resubmse
		
	return meanmse, varmse