File: soft.R

package info (click to toggle)
r-cran-regsem 1.6.2+dfsg-1
  • links: PTS, VCS
  • area: main
  • in suites: bullseye, sid
  • size: 496 kB
  • sloc: cpp: 263; ansic: 15; makefile: 2
file content (66 lines) | stat: -rw-r--r-- 1,967 bytes parent folder | download
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66


soft <- function(par,lambda,type,step,e_alpha,gamma){
  if(type=="lasso"){
    lambda <- lambda*step
   # print(par)
      ret.val <- sign(par)*max(abs(par)-lambda,0)
     # print(ret.val)
  }else if(type=="enet"){
    #http://www.stat.washington.edu/courses/stat527/s13/readings/zouhastie05.pdf ****p. 305****
      lambda <- lambda*step
      lambda2 <- e_alpha*(lambda)
      lambda1 <- (1-e_alpha)*lambda
      ret.val <- sign(par)*(max(abs(par)-(lambda1),0)/(1+lambda2))

    #  ret.val <- (sign(par)*max(abs(par)-step*lambda,0))/(1+lambda2)

    #  if(abs(par) < e_alpha*lambda){
     #   ret.val <- 0
     # }else{
        # might be missing max(0,lambda)
    #    ret.val <- (sign(par)*(abs(par)-e_alpha*lambda))/(1+(1-e_alpha)*lambda)
    #    ret.val <- sign(par)*(max(abs(par)-(step*lambda)/2,0)/(1+step*lambda))
    #  }

  }else if(type=="alasso"){
    # ftp://ftp.stat.math.ethz.ch/Teaching/buhlmann/advanced-comput-statist/notes1.pdf
    ret.val <- sign(par)*max(abs(par)-(step*lambda)/(2*abs(par)),0)


  }else if(type=="scad"){
    lambda <- lambda*step
    gamma <- gamma*step
    #stop("currently not supported")



    if(abs(par) <= 2*lambda){
      ret.val <- sign(par)*max(abs(par)-lambda,0)

    #  ret.val <- sign(par)*max(abs(par)-lambda*gamma,0)
    }else if(abs(par) > 2*lambda & abs(par) <= lambda*gamma){
      ret.val <- ((gamma - 1)/(gamma - 2)) * sign(par)*max(abs(par)-((lambda*gamma)/(gamma-1)),0)

     # ret.val <- ((gamma-1)*par - sign(par)*gamma*lambda)/(gamma-2)
    }else if(abs(par) > lambda*gamma){
      ret.val <- par

    }

  }else if(type=="mcp"){

    lambda <- lambda*step
    gamma <- gamma*step
    #print(lambda*gamma)
    #stop("currently not supported")

    if(abs(par) <= lambda * gamma){
      ret.val <- (gamma/(gamma-1)) * sign(par)*max(abs(par)-lambda,0)
    }else if(abs(par) > lambda*gamma){
      ret.val <- par
    }

  }
  ret.val
}