LASSO算法

徐国盛 浙江财经大学

  ci = 100
  beta0 = c(rep(0,15))
  AA = matrix(0,100,length(beta0))
  for(j in 79:100)  
  {  
    betam = c(7,-5,rep(0,13))
    betam <- as.matrix(betam)
    n = 150
    v = 2
    xm1 = matrix(0,n,v)
    for(i in 1:n)
    {
      xm1[i,] =0.2*rbinom(2,1,0.5)
    }
    xm2 = matrix(0,n,length(betam)-v)
    for(i in 1:n)
    {
      xm2[i,] = rnorm(length(betam)-v,0,10)
    }
    xm =  cbind(xm1,xm2)
    
    u = 0
    for(i in 1:n)
    {
      u[i] =as.numeric(xm[i,]%*%betam)
    }
    #ym = rgamma(n,1)
    #ym = (ym - mean(ym))
    #
    ym = 0
    for(i in 1:n)
    {
      ym[i] = rgamma(1,exp(u[i]))
    }
    y <- as.matrix(ym-mean(ym))
    x <- xm 
    ##########################################################
    likely<- function(beta){
      ll = 0;
      for(i in 1:nrow(x)){
        t = as.numeric(x[i,]%*%as.matrix(beta))
        ll[i] = as.numeric(y[i,])*t-exp(t)
      }
      l = -sum(ll)
      return(l)
    }
    beta0 <- optim(betam,likely)$par
    ##############################################
    #求lammda值，满足BIC准则
    
    n = nrow(x) 
    lammda = 0
    for(i in 1:ncol(x)){
      lammda[i] = as.numeric(log(n)/(n*abs(beta0[i])))
    }
    #求惩罚项之和
    
    pena <- function(beta){
      spena = 0
      for(i in 1:ncol(x)){
        spena[i] =lammda[i]*beta[i]^2/abs(beta0[i])
      }
      spena = sum(spena)
      return(spena)
    }
    #构造目标函数
    
    likely2 <- function(beta){
      
      l2 = likely(beta) + nrow(x)*pena(beta)
      return(l2)
    }
    #进行迭代计算，求得收敛的beta。
    
    iterate <- function(beta){
      std = 1
      while(std > 0.000001){
        root = optim(beta,likely2)$par
        temp = as.matrix(root - beta)
        std = crossprod(temp,temp)
        beta0 = root
        beta = root
        # write.table(beta,file = "betam2.txt",append = T,eol = "\r\n")
        # write.table(std,file = "stdm2.txt",append = T,eol = "\r\n")
      }
      answer = beta
      return(answer)
    }
    AA[j,] = iterate(beta0)
    write.table(AA[j,],"AA.txt",append=T,eol="\r\n")#linux下为eol="\n"
  }