knitr::opts_chunk$set( message = FALSE, warning = FALSE, include = TRUE ) library(Rcpp) suppressMessages(library(tidyverse)) suppressMessages(library(R.matlab)) suppressMessages(library(glmnet)) library(rbenchmark)
The purpose of this document is to verify my ADMM implementation matches that of Stephen Boyd https://stanford.edu/~boyd/papers/admm/logreg-l1/logreg.html
n <- 200 p <- 5 set.seed(815) beta <- c(rnorm(floor(p/2)), rep(0,ceiling(p/2))) b0 <- rnorm(1) # random intercept X <- matrix(rnorm(n*p), nrow = n, ncol = p) pr <- 1/(1+exp(-X %*% beta + b0)) y = rbinom(n,1,pr) y <- ifelse(y==1,1,-1) qplot(as.numeric(y==1), geom = "histogram")
Now save the data in format matlab can read:
writeMat("input.mat", X=X, y=y)
Note: Boyd's implementation uses $m \times \mu$ in place of $\lambda$ for some reason, where m is the number of observations.
#sourceCpp("src/lasso_logisticC.cpp") library(aeffp) admmlasso_logC(X,y, 3, 0.2)
n <- 200 p <- 50 set.seed(815) beta <- c(rnorm(10), rep(0,40)) b0 <- rnorm(1) # random intercept X <- matrix(rnorm(n*p), nrow = n, ncol = p) pr <- 1/(1+exp(-X %*% beta + b0)) y = rbinom(n,1,pr) y <- ifelse(y==1,1,-1) qplot(as.numeric(y==1), geom = "histogram")
benchmark("aeffp" = { mod1 <- admmlasso_logC(X,y,1) }, "glmnet" = { glmnet(X, y, family = "binomial", alpha = 1) }, replications = 10 )
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.