Can i theoretically calculate the aic & bic from log_likelihood in the generated quantities? (apart from loo &waic)
# load libraries
library(rstan)
library(loo)
# Liverpool's Data Premier League scores 2018-2019
liverpool = c(4,2,1,2,2,3,1,0,1,4,1,2,3,1,3,4,3,2,4,5,1,1,4,1,1,3,0,5,0,4,2,2,3,2,2,5,3,2)
length(liverpool)
mean(liverpool)
var(liverpool)
# estimating Gamma parameters given the mean and the var of the data
estGammaParams <- function(mu, var) {
alpha <- mu^2 / var
beta <- mu / var
return(params = list(alpha = alpha, beta = beta))
}
# gamma parameters to pass in gamma prior
estGammaParams(2.25,1.125)
# 1st model of Poisson-Gamma Model
liv ="data {
int N;
int y[N];
}parameters {
real<lower=0> lambda;
}model {
lambda ~ gamma(4.5,2);
y ~ poisson(lambda);
}
generated quantities {
vector[N] log_lik;
for (n in 1:N) {
log_lik[n] = poisson_lpmf(y[N] | lambda);
}} "
# data as a list ALWAYS
stan.dat <- list(N = length(liverpool),y = liverpool)
# fit the sampling model to HMC
fit <- stan(model_code = liv,
data = stan.dat,
iter = 2000,
chains = 4,
seed = 12345,
control = list(adapt_delta=0.9))
# print posterior mean of lambda parameter
print(fit,pars = "lambda")
# plot parameter
plot(fit,pars = "lambda")
# extract log_lik from posterior samples
log_lik_1 <- extract_log_lik(fit,
parameter_name = "log_lik",
merge_chains = FALSE)
N = length(liverpool)
k = 1
deviance = -2*sum(log_lik_1);deviance
Aic = deviance + 2 * k ;Aic
Bic = deviance + log(N) * k;Bic
r_eff1 <- relative_eff(exp(log_lik_1))
loo_1 <- loo(log_lik_1, r_eff = r_eff1, cores = 4)
print(loo_1)
plot(loo_1)