{{ message }}

Instantly share code, notes, and snippets.

# dirkschumacher/armacmp_logreg.md

Last active Aug 4, 2019
```library(armacmp)
# Arnold, T., Kane, M., & Lewis, B. W. (2019). A Computational Approach to Statistical Learning. CRC Press.
# logistic regression using the Newton-Raphson
log_reg <- armacmp(function(X, y) {
beta <- rep.int(0, ncol(X))
for (i in seq_len(25)) {
b_old <- beta
alpha <- X %*% beta
p <- 1 / (1 + exp(-alpha))
W <- p * (1 - p)
XtX <- crossprod(X, diag(W) %*% X)
score <- t(X) %*% (y - p)
delta <- solve(XtX, score)
beta <- beta + delta
}
return(beta)
})

log_reg_r <- function(X, y) {
beta <- rep.int(0, ncol(X))
for (i in seq_len(25)) {
b_old <- beta
alpha <- X %*% beta
p <- 1 / (1 + exp(-alpha))
W <- as.numeric(p * (1 - p))
XtX <- crossprod(X, diag(W) %*% X)
score <- t(X) %*% (y - p)
delta <- solve(XtX, score)
beta <- beta + delta
}
return(beta)
}

n <- 1000 ; p <- 50
true_beta <- rnorm(p)
X <- cbind(1, matrix(rnorm(n * (p - 1)), ncol = p - 1))
y <- runif(n) < plogis(X %*% true_beta)

bench::mark(
as.numeric(log_reg(X, y)),
as.numeric(log_reg_r(X, y))
)
#> Warning: Some expressions had a GC in every iteration; so filtering is
#> disabled.
#> # A tibble: 2 x 6
#>   expression                       min  median `itr/sec` mem_alloc `gc/sec`
#>   <bch:expr>                  <bch:tm> <bch:t>     <dbl> <bch:byt>    <dbl>
#> 1 as.numeric(log_reg(X, y))    54.19ms  65.1ms    14.6      10.8KB     0
#> 2 as.numeric(log_reg_r(X, y))    1.35s   1.35s     0.741   212.1MB     2.96```

Created on 2019-08-03 by the reprex package (v0.3.0)