inv_logit <- function(input, param){
  return(1 / (1 + exp(-input%*%param)))
}

library(rstan)
rstan_options(auto_write=T)
options(mc.cores = parallel::detectCores())

n <- 100
p <- 10

xrange <- c(-5, 5)
x <- matrix(runif(n*p, min = xrange[1], max=xrange[2]), nrow=n, ncol=p)

true_beta <- rnorm(p, mean=0, sd=1)

bernoulli_realization <- runif(n)

bernoulli_prob <- inv_logit(x, true_beta)

y <- as.numeric(bernoulli_prob > bernoulli_realization)

hy_center <- 0
hy_sd <- 10

logistic_data <- list(n = n, p=p, x=x, y=y, hy_center=hy_center, hy_sd=hy_sd)
logistic_fit <- stan(file="logistic_regression.stan", data=logistic_data, chains=1, iter=2000)

logistic_param <- extract(logistic_fit, permuted=F)
parameter_num <- dim(logistic_param)[1]
predicted_bernoulli_prob <- numeric(n)
for(l in 1:parameter_num){
  predicted_bernoulli_prob <- predicted_bernoulli_prob + inv_logit(x, logistic_param[l,1,1:p])/parameter_num
}
print(predicted_bernoulli_prob)

logisitic_model <- stan_model("logistic_regression.stan")
gradient_logistic <- optimizing(logisitic_model, data=logistic_data)