R Functions :Forward Selection Based on Robust Multivariate Correlation

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

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

#R_ functions :Forward Selection Based on Robust Multivariate Correlation

#Author: Hassan S. Uraibi and Habshah Midi

#E-mail:

#Date: Aug , 10, 2013

#Version: 0.1

#Copyright (C) 2013 Hassan S. Uraibi and Habshah Midi

#

######################################################################## Robust Forward Selection Based on RFCH correlations Matrix

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

RFS.RFCH.M<-function(Xtr,ytr,Xts,yts)

{

Aic<-c()

H<-c()

Ch<-c()

TT<-c()

pe<-pvalue<-c()

active <- c()

library(MASS)

library(robustbase)

library(robust)

temp.rob <- standardize.rob(Xtr,ytr)

X <- temp.rob$X

y <- temp.rob$y

temp.rob.ts <- standardize.rob(Xts,yts)

X.ts <- temp.rob.ts$X

y.ts <- temp.rob.ts$y

data<-data.frame(y,X)

fit.rob<-rlm(y~.,data,method="MM",psi="psi.bisquare")

weight<- fit.rob$ "w"

prob<-weight/sum(weight)

XR<-X[prob!=0,]

YR<-y[prob!=0]

inactive <- c(seq(1, (k <- ncol(XR))))

names(inactive)<-1:length(inactive)

Z1<-as.matrix(data.frame(XR,YR))

b<-dim(Z1)[2]

R<-matrix(0,b,b) # construction a robust corrlation matrix

Cov<-as.matrix(Rfch(x=Z1)$rcovf)

for(i in 1:b){

for(j in i:b){

c<-Cov[i,j]/(sqrt(Cov[i,i])*sqrt(Cov[j,j]))

R[i,j]<-c

R[j,i]<-c}}

R<-round(R,4)

N<-ncol(R)

RR<-R

diag(RR)<-0

RY<-RR[N,]

pos<-which.max(abs(RY))

active <- c(active,as.numeric(pos))

inactive <- inactive[inactive != pos]

MG<-RY[pos]

x<-as.matrix(Z1[,pos])

n<-dim(x)[1]

day<-data.frame(x,YR)

aic<-AIC(rlm(YR~.,day,method="MM"))

Aic<-c(Aic,aic )

n<-dim(x)[1]

s<-MG^2

S<-1-s

Fcal<-((n-2)*s)/S

df1<-1

df2<-(n-2)

Ftab<-qf(.95, df1, df2)

Ch<-c(Ch,Fcal)

TT<-c(TT,Ftab)

pval<-pf(Fcal,1,df2,lower.tail=F)

Pval<- 2 * min(pval, 1 - pval)

pvalue<-c(pvalue,Pval)

data.ts<-data.frame(y.ts,X.ts)

fit.rob.ts<-rlm(y.ts~.,data.ts,method="MM",psi="psi.bisquare")

weight.ts<- fit.rob.ts$ "w"

prob.ts<-weight.ts/sum(weight.ts)

XR.ts<-X.ts[prob.ts!=0,]

YR.ts<-y.ts[prob.ts!=0]

p.clas<-prediction(XR.ts,YR.ts,XR,YR,active)

pe<-c(pe,p.clas)

g<-c(inactive,b)

m<-R

for( j in 1:length(g)){

par<- P.cor.(m,g, active)

diag(par)<-0

b2<-dim(par)[2]

py<-par[,b2]

#MG<- max(abs(py))

PO<-which.max(abs(py))

MG<- py[PO]

names(py)<-g

bo<-as.numeric(names( py[PO]))

g<- g[-PO]

active <- c(active,bo)

p.clas2<-prediction(XR.ts,YR.ts,XR,YR,active)

pe<-c(pe,p.clas2)

x<-X[,active]

n<-length(y)

aic2<-AIC.MM.has(x,y,n)$aic

Aic<-c(Aic,aic2 )

H<-active

fp<-F.p(H,n,MG,S)

F.p2<-fp$F.p1

Ftab2<-fp$Ftab1

pval<-pf(F.p2,length(H),fp$df2,lower.tail=F)

Pval<- 2 * min(pval, 1 - pval)

pvalue<-c(pvalue,Pval)

Ch<-c(Ch,F.p2)

TT<-c(TT,Ftab2)

if(F.p2<Ftab2){

active<-active[-length(H)]

Aic<-Aic[-length(Aic)]

pe<-pe[-length(pe)]

break("\'F.p2\' is not significant!\n")

}

if(g[1]==dim(X)[2]+1)break

S<-fp$S

}

dat<-data.frame(y,X[,active])

Pval.cof<- summary(lmRob(y~.,dat))$coef[,4][-1]

names(Pval.cof)<-active

sel<-which(Pval.cof<(0.05))

Active<-as.numeric(names(sel))

select<-colnames(X)[Active]

list( Active=Active,active=active,pe=pe,fc=Ch,ft=TT,pval=round(pvalue,4))}

Mahanalobis disteance function

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

mahalanobis2<-function (x, center, cov, inverted = FALSE, ...)

{

x <- if (is.vector(x))

matrix(x, ncol = length(x))

else as.matrix(x)

x <- sweep(x, 2, center)

if (!inverted)

cov <- ginv(cov, ...)

retval <- rowSums((x %*% cov) * x)

names(retval) <- rownames(x)

retval

}

# Partial correlation function #

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

partial.r<-function (m, x, y)

{

cl <- match.call()

if (dim(m)[1] != dim(m)[2]) {

n.obs <- dim(m)[1]

m <- cor(m, use = "pairwise")

}

if (!is.matrix(m))

m <- as.matrix(m)

nm <- dim(m)[1]

t.mat <- matrix(0, ncol = nm, nrow = nm)

xy <- c(x, y)

numx <- length(x)

numy <- length(y)

nxy <- numx + numy

for (i in 1:nxy) {

t.mat[i, xy[i]] <- 1

}

reorder <- t.mat %*% m %*% t(t.mat)

reorder[abs(reorder) > 1] <- NA

X <- reorder[1:numx, 1:numx]

Y <- reorder[1:numx, (numx + 1):nxy]

phi <- reorder[(numx + 1):nxy, (numx + 1):nxy]

phi.inv <- solve(phi)

X.resid <- X - Y %*% phi.inv %*% t(Y)

#X.resid <- cov2cor(X.resid)

#colnames(X.resid) <- rownames(X.resid) <- colnames(m)[x]

class(X.resid) <- c("psych", "partial.r")

return(X.resid)

}

# Function of square root of classical prediction error#

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

prediction<-function(Xts,yts,Xtr,ytr,active){

xtr<-as.matrix(Xtr[,active])

colnames(xtr)<-active

data<-data.frame(xtr,ytr)

cof.tr<-as.matrix(lm(ytr~.,data)$coef[-1])

xts<-as.matrix(Xts[,active])

colnames(xts)<-active

yhad<-xts%*% cof.tr

resid<-yts-yhad

rmse<-sqrt(sum(resid^2)/(length(yts)-length(active)-1))

return(rmse)}

# Function of square of robust prediction error #

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

R.prediction<-function(Xts,yts,Xtr,ytr,active){

xtr<-as.matrix(Xtr[,active])

colnames(xtr)<-active

data<-data.frame(xtr,ytr)

cof.tr<-as.matrix(rlm(ytr~.,data,method="MM")$coef[-1])

xts<-as.matrix(Xts[,active])

colnames(xts)<-active

data.2<-data.frame(xts,yts)

rfit<-rlm(yts~.,data.2,method="MM")

rfit$w

w<- rfit $"weights"

ywts<-yts*w

yhad<-xts%*% cof.tr

resid<-ywts-yhad

rmse<-sqrt(sum(resid^2)/(length(yts)-length(active)-1))

return(rmse)}

# Robust scale function #

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

standardize.rob<-function(X, y)

{

X <- apply(X,2,robstand)

y <- robstand(y)

list(X=X, y=y)

}

robstand<-function(x) # robust standardize

{

if ((m<-mad(x))==0){ (x-mean(x))/sd(x) }

else {(x-median(x))/m }

}

# Special function for simulation #

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

LARS_RFCH<-function(X,y, csteps = 10)

{

x<-data.frame(y,X)

# works for p = 1

zx <- x

x <- as.matrix(x)

p <- dim(x)[2]##get the DGK estimator

Covs <- var(x)

mns <- apply(x, 2, mean)## concentrate

for(i in 1:csteps) {

md2 <- mahalanobis(x, mns, Covs)

medd2 <- median(md2)

if(p > 1) {

mns <- apply(x[md2 <= medd2, ], 2, mean)

}

if(p == 1) {

mns <- mean(zx[md2 <= medd2])

}

covs <- var(x[md2 <= medd2, ])

if(sum(abs(diag(covs)))-sum(abs(diag(Covs)))==0)break

Covs<-covs

}

XDGK<-x[md2 <= medd2, ]

ptm <-proc.time()

LARS_DGK<-lars_cor(X=XDGK[,-1], y=XDGK[,1], m=ncol(XDGK[,-1]))

covd <- covs

mnd <- mns##get the MB estimator

covv <- diag(p)

med <- apply(x, 2, median)

md2 <- mahalanobis(x, center = med, covv)

medd2 <- median(md2) ##get the location criterion

lcut <- medd2

## get the start

if(p > 1) {

mns <- apply(x[md2 <= medd2, ], 2, mean)

}

if(p == 1) {

mns <- mean(zx[md2 <= medd2])

}

Covs <- var(x[md2 <= medd2, ])

## concentrate

for(i in 1:csteps) {

md2 <- mahalanobis(x, mns, covs)

medd2 <- median(md2)

if(p > 1) {

mns <- apply(x[md2 <= medd2, ], 2, mean)

}

if(p == 1) {

mns <- mean(zx[md2 <= medd2])

}

covs <- var(x[md2 <= medd2, ])

if(sum(abs(diag(covs)))-sum(abs(diag(Covs)))==0)break

Covs<-covs

}

covm <- covs

mnm <- mns

XMB<-x[md2 <= medd2, ]

LARS_MB<-lars_cor(X=XMB[,-1], y=XMB[,1],m=ncol(XMB[,-1]))

indx<-c(0,1)

##get FCH attractor

covf <- covm

mnf <- mnm

val <- mahalanobis(t(mnd), med, covv)

if(val < lcut) {##crit = square root of det(cov)

critd <- prod(diag(chol(covd)))

critm <- prod(diag(chol(covm)))

if(critd < critm) {

covf <- covd

mnf <- mnd

IND<-indx[1]

}

IND<-indx[2]

}

## get FCH estimator

chisqm <- qchisq(0.5, p)

rd2 <- mahalanobis(x, mnf, covf)

const <- median(rd2)/chisqm

covf <- const * covf

##reweight the above FCH estimator (mnf,covf) to get the RFCH estimator

## (rmnf,rcovf)

rd2 <- mahalanobis(x, mnf, covf)

up <- qchisq(0.975, p)

if(p > 1) {

rmnf <- apply(x[rd2 <= up, ], 2, mean)

}

if(p == 1){

rmnf = mean(zx[rd2 <= up])

}

rcovf <- var(x[rd2 <= up, ])

rd2 <- mahalanobis(x, rmnf, rcovf)

const <- median(rd2)/chisqm

rcovf <- const * rcovf## reweight again

rd2 <- mahalanobis(x, rmnf, rcovf)

if(p > 1){

rmnf <- apply(x[rd2 <= up, ], 2, mean)

}

if(p == 1){

rmnf = mean(zx[rd2 <= up])

}

rcovf <- var(x[rd2 <= up, ])

XRFCH<-x[rd2 <= up, ]

LARS_RFCH<-lars_cor(X=XRFCH[,-1], y=XRFCH[,1],m=ncol(XRFCH[,-1]))

list( loc.rf= rmnf,rcovf=rcovf,LARS_RFCH= LARS_RFCH, LARS_MB= LARS_MB,LARS_DGK=LARS_DGK,XDGK=XDGK,XRFCH=XRFCH,XMB=XMB)

}

lars_cor<-function(X, y, m=ncol(X))

{

ptm <- proc.time()

X <- scale(X)

y <- scale(y)

inactive <- c(seq(1, (k <- ncol(X))))

active <- c()

signvec<-c()

#up.x <- c()

actmat<-diag(rep(1,k))

up.cor <- rep(0, k)

for(j in 1:k) {

d<-dim(as.matrix(X[,j]))[2]

up.cor[j] <- cor(y,x=X[,j])

}

ycorvec<-up.cor

up.maxcor <- max(abs(up.cor))

newsub <- which(abs(up.cor) == max(abs(up.cor)))

newsign <- sign(up.cor[newsub])

active <- c(active, newsub)

inactive <- inactive[inactive != newsub]

signvec<-c(signvec,newsign)

#up.x <- cbind(up.x, newsign * X[, newsub])

for(j in inactive)

{actmat[newsub,j]<-actmat[j,newsub]<- cor(y=X[,newsub],x=X[,j]) }

ones<-rep(1,length(active))

G <- diag(ones)

inv.G <- solve(G)

up.A <- as.double((ones %*% inv.G %*% ones)^(-0.5))

up.w <- up.A * (inv.G %*% ones)

#up.u <- up.x %*% up.w

up.a <- rep(0, k)

for(j in inactive){

corvec <- actmat[j,active]*signvec

up.a[j] <- sum(corvec * up.w)}

up.gam <- 9999999

for(i in inactive) {

#cat("i=",i,"\n")

if((aa <- (up.maxcor - up.cor[i])/(up.A - up.a[i])) > 0 & aa <

up.gam) {

up.gam <- aa

#cat("what ", up.gam, "\n")

newsub <- i

newsign <- 1

}

if((bb <- (up.maxcor + up.cor[i])/(up.A + up.a[i])) > 0 & bb <

up.gam) {

up.gam <- bb

#cat("what ", up.gam, "\n")

newsub <- i

newsign <- -1

}

}

up.gam <- as.double(up.gam)

#cat("current gamma =", up.gam, "\n")

#up.mean <- up.mean + up.gam*up.u

#cat("current mean =", up.mean, "\n")

#up.resid <- y - up.mean

up.maxcor <- up.maxcor - up.gam * up.A

up.cor <- up.cor - up.gam * up.a

#cat("step", "\t", j, "\n")

for(l in 1:(m - 2)) {

# cat(l,"\n")

active <- c(active, newsub)

signvec<-c(signvec,newsign)

inactive <- inactive[inactive != newsub]

#up.x <- cbind(up.x, newsign * X[, newsub])

for(j in inactive){

actmat[newsub,j]<-actmat[j,newsub]<- cor(y=X[,newsub],x=X[,j])}

ones <- rep(1, (d = length(active)))

G <- (signvec%*%t(signvec))*actmat[active,active]

d1=dim(G)[1]

check <- eigen(G,symmetric=T)

lamb <- check$values

if (min(lamb)<0)

{

evec <- check$vectors

lamb <- (apply(X[,active]%*%evec,2,mad))^2

G <- matrix(0,ncol=d1,nrow=d1)

for (o in 1:d1)

{G <- G + lamb[o]*(evec[,o]%*%t(evec[,o]))}

}

inv.G <- solve(G)

# inv.G=solve(G, tol = 1e-25)

up.A <- as.double((ones %*% inv.G %*% ones)^(-0.5))

up.w <- up.A * (inv.G %*% ones)

#up.u <- up.x %*% up.w

up.a <- rep(0, k)

for(j in inactive){

corvec <- actmat[j,active]*signvec

up.a[j] <- sum(corvec * up.w)}

up.gam <- 9999999

for(i in inactive) {

#cat("i=",i,"\n")

if((aa <- (up.maxcor - up.cor[i])/(up.A - up.a[i])) >

0 & aa < up.gam) {

up.gam <- aa

#cat("what ", up.gam, "\n")

newsub <- i

newsign <- 1

}

if((bb <- (up.maxcor + up.cor[i])/(up.A + up.a[i])) >

0 & bb < up.gam) {

up.gam <- bb

#cat("what ", up.gam, "\n")

newsub <- i

newsign <- -1

}

}

up.gam <- as.double(up.gam)

#cat("current gamma =", up.gam, "\n")

#up.mean <- up.mean + up.gam*up.u

#cat("current mean =", up.mean, "\n")

#up.resid <- y - up.mean

up.maxcor <- up.maxcor - up.gam * up.A

up.cor <- up.cor - up.gam * up.a

}

active <- c(active, newsub)

s.t <- (proc.time()-ptm)[1]

list(active=active,time= s.t)

}

Rlars<-function(X, y, m=ncol(X))

{

temp <- standardize.rob(X, y)

X <- temp$X

y <- temp$y

inactive <- c(seq(1, (k <- ncol(X))))

active <- c()

signvec<-c()

#up.x <- c()

actmat<-diag(rep(1,k))

up.cor <- rep(0, k)

for(j in 1:k) {

up.cor[j] <- bi.hubcor(X[, j], y)

}

ycorvec<-up.cor

up.maxcor <- max(abs(up.cor))

newsub <- which(abs(up.cor) == max(abs(up.cor)))

newsign <- sign(up.cor[newsub])

active <- c(active, newsub)

inactive <- inactive[inactive != newsub]

signvec<-c(signvec,newsign)

#up.x <- cbind(up.x, newsign * X[, newsub])

for(j in inactive)

{actmat[newsub,j]<-actmat[j,newsub]<-bi.hubcor(X[,newsub],X[, j])}

ones<-rep(1,length(active))

G <- diag(ones)

inv.G <- solve(G)

up.A <- as.double((ones %*% inv.G %*% ones)^(-0.5))

up.w <- up.A * (inv.G %*% ones)

#up.u <- up.x %*% up.w

up.a <- rep(0, k)

for(j in inactive){

corvec <- actmat[j,active]*signvec

up.a[j] <- sum(corvec * up.w)}

up.gam <- 9999999

for(i in inactive) {

#cat("i=",i,"\n")

if((aa <- (up.maxcor - up.cor[i])/(up.A - up.a[i])) > 0 & aa <

up.gam) {

up.gam <- aa

#cat("what ", up.gam, "\n")

newsub <- i

newsign <- 1

}

if((bb <- (up.maxcor + up.cor[i])/(up.A + up.a[i])) > 0 & bb <

up.gam) {

up.gam <- bb

#cat("what ", up.gam, "\n")

newsub <- i

newsign <- -1

}

}

up.gam <- as.double(up.gam)

#cat("current gamma =", up.gam, "\n")

#up.mean <- up.mean + up.gam*up.u

#cat("current mean =", up.mean, "\n")

#up.resid <- y - up.mean

up.maxcor <- up.maxcor - up.gam * up.A

up.cor <- up.cor - up.gam * up.a

#cat("step", "\t", j, "\n")

for(l in 1:(m - 2)) {

# cat(l,"\n")

active <- c(active, newsub)

signvec<-c(signvec,newsign)

inactive <- inactive[inactive != newsub]

#up.x <- cbind(up.x, newsign * X[, newsub])

for(j in inactive){

actmat[newsub,j]<-actmat[j,newsub]<-bi.hubcor(X[,newsub],X[,j])}

ones <- rep(1, (d = length(active)))

G <- (signvec%*%t(signvec))*actmat[active,active]

d1=dim(G)[1]

check <- eigen(G,symmetric=T)

lamb <- check$values

if (min(lamb)<0)

{

evec <- check$vectors

lamb <- (apply(X[,active]%*%evec,2,mad))^2

G <- matrix(0,ncol=d1,nrow=d1)

for (o in 1:d1)

{G <- G + lamb[o]*(evec[,o]%*%t(evec[,o]))}

}

inv.G <- solve(G)

# inv.G=solve(G, tol = 1e-25)

up.A <- as.double((ones %*% inv.G %*% ones)^(-0.5))

up.w <- up.A * (inv.G %*% ones)

#up.u <- up.x %*% up.w

up.a <- rep(0, k)

for(j in inactive){

corvec <- actmat[j,active]*signvec

up.a[j] <- sum(corvec * up.w)}

up.gam <- 9999999

for(i in inactive) {

#cat("i=",i,"\n")

if((aa <- (up.maxcor - up.cor[i])/(up.A - up.a[i])) >

0 & aa < up.gam) {

up.gam <- aa

#cat("what ", up.gam, "\n")

newsub <- i

newsign <- 1

}

if((bb <- (up.maxcor + up.cor[i])/(up.A + up.a[i])) >

0 & bb < up.gam) {

up.gam <- bb

#cat("what ", up.gam, "\n")

newsub <- i

newsign <- -1

}

}

up.gam <- as.double(up.gam)

#cat("current gamma =", up.gam, "\n")

#up.mean <- up.mean + up.gam*up.u

#cat("current mean =", up.mean, "\n")

#up.resid <- y - up.mean

up.maxcor <- up.maxcor - up.gam * up.A

up.cor <- up.cor - up.gam * up.a

}

active <- c(active, newsub)

list(active=active)

}

standardize.rob<-function(X, y)

{

X <- apply(X,2,robstand)

y <- robstand(y)

list(X=X, y=y)

}

# Following function adjusts for 0-1 variable

robstand<-function(x)

{

if ((m<-mad(x))==0){ (x-mean(x))/sd(x) }

else {(x-median(x))/m }

}

bi.hubcor<-function(x, y, a = 2, p = 0.95)

{

robcov <- diag(c(1,1))

robcov[1,2] <- simp.hubercor(x,y, a)

robcov[2,1] <- robcov[1,2]

xy <- cbind(x, y)

n <- nrow(xy)

d <- qchisq(p, 2)

invcov<-ginv(robcov)

D <- rowSums((xy %*% invcov) * xy)

#D <- diag(xy %*% ginv(robcov) %*% t(xy))

ind <- D>d

xy[ind,] <- xy[ind,]*sqrt(d/D[ind])

x <- xy[, 1]

y <- xy[, 2]

cor(x,y)

}

simp.hubercor<-function(x, y, a=2)

{

indx<-abs(x)>a

x[indx]<-a*sign(x[indx])

indy<-abs(y)>a

y[indy]<-a*sign(y[indy])

n <- length(x)

aa <- x*y

ind <- aa < 0

n1 <- sum(ind)

n2 <- n - n1

if(n1 <= n2)

{

a1 <- a*sqrt(n1/n2)

indx<-abs(x[ind])>a1

x[ind][indx] <- a1*sign(x[ind][indx])

indy<-abs(y[ind])>a1

y[ind][indy] <- a1*sign(y[ind][indy])

}

else

{

ind <- aa > 0

a1 <- a*sqrt(n2/n1)

indx<-abs(x[ind])>a1

x[ind][indx] <- a1*sign(x[ind][indx])

indy<-abs(y[ind])>a1

y[ind][indy] <- a1*sign(y[ind][indy])

}

cor(x, y)

}

# Simulation Study #

#######################################################################d=30

n=500

p=9

roh=0

sigma=5

epsilon=0.1

B=200

library(MASS)

library(mvtnorm)

library(robust)

n.zero.fs<-n.zero.fws<-n.zero.rfc<-c()

zero.fs<-zero.fws<-zero.rfc<-c()

len.fs<-len.fws<-len.rfc<-c()

mse.fs<-mse.fws<-mse.rfc<-c()

for( j in 1:B){

beta <- c(c(p:1),c(rep(0,d-p))) # coefficients

std<-1:p

Sigma <- roh^t(sapply(1:d, function(i, j) abs(i-j), 1:d))

X <- rmvnorm(n, sigma=Sigma) # predictor matrix

colnames(X)<-1:d

e <- rnorm(n) # error terms

i <- 1:ceiling(epsilon*n) # observations to be contaminated

e[i] <- e[i] + 100 # vertical outliers

y <- c(X%*% beta + sigma * e) # response

#X[i,(1:p)] <- X[i,(1:p)]+ 100 # bad leverage points

samp<-sample(1:length(y),.6*length(y),rep=FALSE)

Xtr<-X[samp,]

Xts<-X[-samp,]

ytr<-y[samp]

yts<-y[-samp]

samp0<-sample(1:length(samp),.1*length(samp),rep=FALSE)

Xtr[samp0,c(1:10)]<-100

fs<-FS.M(Xtr,ytr,Xts,yts)

Active.fs<-fs$active

Nzero.fs<-as.numeric(which(std%in%Active.fs[1:p]))

n.zero.fs<-c(n.zero.fs,length(Nzero.fs))

nois.fs<-as.numeric(which(!Active.fs%in%std))

zero.fs<-c(zero.fs,length(nois.fs))

L.fs<-length(Active.fs)

len.fs<-c(len.fs,L.fs)

mse.cl<-fs$pe[L.fs]*L.fs

mse.fs<-c(mse.fs,mse.cl)

fws<-RFS.Winso.M(Xtr,ytr,Xts,yts)

Active.fws<-fws$active

Nzero.fws<-as.numeric(which(std%in%Active.fws[1:p]))

n.zero.fws<-c(n.zero.fws,length(Nzero.fws))

nois.fws<-as.numeric(which(!Active.fws%in%std))

zero.fws<-c(zero.fws,length(nois.fws))

L.fws<-length(Active.fws)

len.fws<-c(len.fws,L.fws)

mse.fw<-fws$pe[L.fws]*L.fws

mse.fws<-c(mse.fws,mse.fw)

rfc<-RFS.RFCH.M(Xtr,ytr,Xts,yts)

Active.rfc<-rfc$Active

Nzero.rfc<-as.numeric(which(std%in%Active.rfc[1:p]))

n.zero.rfc<-c(n.zero.rfc,length(Nzero.rfc))

nois.rfc<-as.numeric(which(!Active.rfc%in%std))

zero.rfc<-c(zero.rfc,length(nois.rfc))

L.rfc<-length(Active.rfc)

len.rfc<-c(len.rfc,L.rfc)

mse.rf<-rfc$pe[L.rfc]*L.rfc

mse.rfc<-c(mse.rfc,mse.rf)

}

Non_zero<-rbind(mean(n.zero.fs),mean(n.zero.fws),mean(n.zero.rfc))

Zero<-rbind(mean(zero.fs),mean(zero.fws),mean(zero.rfc))

length<-rbind(mean(len.fs),mean(len.fws),mean(len.rfc))

MSE<-rbind(mean(mse.fs)/mean(n.zero.fs),mean(mse.fws)/mean(n.zero.fws),mean(mse.rfc)/mean(n.zero.rfc))

data.frame(Non_zero,Zero,length,MSE)