Lasso算法与AIC、bic

sgn1<-matrix(rep(0,6),1,6) for (k1 in 1:6) {

if (abs(beta0[k1,1]<1e-4)) beta0[k1,1]<-1e-4 #如果beta0的绝对值太小，赋予定值 sgn1[1,k1]<-1/abs(beta0[k1,1]) }

U4<-t(X)%*%X+200*lam*diag(c(sgn1))

P<-X%*%solve(U4)%*%t(X) #U3是倒数第二次迭代的beta生成的绝对值倒数矩阵， tr<-sum(diag(P)) #求P得迹

t1<-t(Y-X%*?ta0)%*%(Y-X%*?ta0) G[1,k]<-t1/((1-tr/200)^2) }

f<-which(G==min(G)) #求使判断方程最小的下角标

beta<-Bet[,f] #经过有限次迭代和GCV方程比较后最后得到的beta

for(j in 1:3) #计算beta前三个和后三个变量分别等于0的个数 {

if(beta[j]==0) incorr<-incorr+1 if(beta[j+3]==0) corr<-corr+1 }

#二进制表示全子集 Index <- function(index) {

index <- index + c(1,0,0,0,0,0)

for(i in 1:6) #逢二进一 {

if(index[i] == 2) {

index[i] <- 0

index[i+1] <- index[i+1] + 1 } }

Index <- index }

index <- matrix(0,64,6) #计算AIC、BIC aic <- matrix(0,63,1) bic <- matrix(0,63,1) for (i in 1:63) {

index[i+1,] <- Index(index[i,]) # 变量的全子集

p <- sum(index[i+1,]) # 选入变量的个数

X1 <- X[,which(index[i+1,]==1)] # 选入的子集重新构造矩阵X1

aic[i] <- 200*log(sum((Y-X1%*%(solve(t(X1)%*%X1)%*%t(X1)%*%Y))^2)) + 2*p

bic[i] <- log((sum((Y-X1%*%(solve(t(X1)%*%X1)%*%t(X1)%*%Y))^2))/200) + (p/200)*log(200) }

aic1<-rbind(index[which.min(aic)+1,]) #输出使AIC最小的子集 bic1<-rbind(index[which.min(bic)+1,]) #输出使BIC最小的子集

for(j in 1:3) #在AIC和BIC算法中计算beta前三个和后三个变量分别等于0的个数 {

if(aic1[j]==0) incorra<-incorra+1 if(aic1[j+3]==0) corra<-corra+1 if(bic1[j]==0) incorrb<-incorrb+1 if(bic1[j+3]==0) corrb<-corrb+1 }

# Stepwise算法 x1 <- X[,1] x2 <- X[,2] x3 <- X[,3] x4 <- X[,4] x5 <- X[,5] x6 <- X[,6]

index1 <- matrix(0,1,6) # Regression

sol <- data.frame(x1,x2,x3,x4,x5,x6,Y)

lm <- lm(Y ~ x1+x2+x3+x4+x5+x6, data=sol) lmstep <- step(lm, trace=0) form <- formula(lmstep) vari <- as.character(form) vari <- vari[3] if(grepl(\ index1[1] <- 1 if(grepl(\ index1[2] <- 1 if(grepl(\ index1[3] <- 1 if(grepl(\ index1[4] <- 1 if(grepl(\ index1[5] <- 1 if(grepl(\ index1[6] <- 1

for(j in 1:3) #计算stepwise算法中beta前三个和后三个变量分别等于0的个数 {

if(index1[j]==0) incorrs<-incorrs+1 if(index1[j+3]==0) corrs<-corrs+1 } }

correct<-corr/1000 #输出lasso的结果 correct

incorrect<-incorr/1000 incorrect

correcta<-corra/1000 #输出AIC的结果 correcta

incorrecta<-incorra/1000 incorrecta

correctb<-corrb/1000 #输出BIC的结果 correctb

incorrectb<-incorrb/1000 incorrectb

corrects<-corrs/1000 #输出stepwise的结果 corrects

incorrects<-incorrs/1000 incorrects