Hello, Stan friends,
I am trying to fit a logistic regression model using L1 and L2 regularization (elastic-net). I have followed the Stan reference manual, and a basic model that does not include QR reparametrization is working. However, when the QR reparametrization is included, the model does not compile. Below is a full example for R showing the two models. Any suggestions and feedback will be greatly appreciated.
Thanks,
Leandro
library(rstan)
X <- matrix(rnorm(10000), 1000, 10)
factor_variable <- sample(c(1,2), 1000, replace = TRUE)
y <- sample(c(0,1), 1000, replace = TRUE)
model <- "data {
int<lower=0> N;
int<lower=0> M;
matrix[N, M] X;
int factor_variable[N];
int<lower=0,upper=1> y[N];
}
parameters {
vector[M] beta;
vector[2] beta_factor_variable;
real<lower=0> lambda1;
real<lower=0> lambda2;
}
model {
beta_factor_variable ~ normal(0, 5);
beta ~ normal(0, 1);
y ~ bernoulli_logit(beta_factor_variable[factor_variable] + X * beta);
target += -lambda1 * fabs(beta);
target += -lambda2 * dot_self(beta);
}
generated quantities {
vector[M] beta_elastic_net;
beta_elastic_net = (1 + lambda2) * beta;
}"
model_qr <- "
data {
int<lower=0> N;
int<lower=0> M;
matrix[N, M] X;
int factor_variable[N];
int<lower=0,upper=1> y[N];
}
transformed data {
// Compute, thin, and then scale QR decomposition
matrix[N, M] Q;
matrix[M, M] R;
matrix[M, M] R_inv;
Q = qr_Q(X)[, 1:M] * sqrt(N - 1);
R = qr_R(X)[1:M, ] / sqrt(N - 1);
R_inv = inverse(R);
}
parameters {
vector[M] beta_tilde;
vector[2] beta_factor_variable;
real<lower=0> lambda1;
real<lower=0> lambda2;
}
transformed parameters {
vector[M] beta = R_inv * beta_tilde;
}
model {
beta_factor_variable ~ normal(0, 5);
// I use target notation here for beta (instead of ~normal)
// to avoid a warning message:
// Left-hand side of sampling statement (~)
// may contain a non-linear transform of a parameter or local variable.
target += normal_lpdf(beta | 0, 1);
y ~ bernoulli_logit(Q * beta_tilde + beta_factor_variable[factor_variable]);
target += -lambda1 * fabs(beta);
target += -lambda2 * dot_self(beta);
}
generated quantities {
vector[M] beta_elastic_net;
beta_elastic_net = (1 + lambda2) * beta;
}
"
data <- list(N = nrow(X), M = ncol(X), X= X,
factor_variable= factor_variable,
y = y)
# Basic model
fit <- stan(model_code = model, data = data, cores = 3,
chains = 3, iter = 1000, warmup = 500, thin = 10)
# QR model
fit_qr <- stan(model_code = model_qr, data = data, cores = 3,
chains = 3, iter = 1000, warmup = 500, thin = 10)
This is my sessioninfo():
R version 3.5.0 (2018-04-23)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Arch Linux
Matrix products: default
BLAS: /usr/lib/libblas.so.3.8.0
LAPACK: /usr/lib/liblapack.so.3.8.0
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] rstan_2.17.3 StanHeaders_2.17.2 ggplot2_3.0.0