Getting maximum likelihood to compute AIC, BIC and DIC

[jpmvferreira]

Hi there,

I'm using arviz to analyze the output of a fit performed by pystan3, and I wish to compute the value of both AIC, BIC and DIC.

I've seen a few discussions, namely this one over at the Stan forums, where it's argued that both AIC and BIC are not baeysian quantities, as they rely on the maximum likelihood instead of working with the full distribution. Instead, they recommend using WAIC and LOO, both of which are already implemented in arviz.

However, in my field, which is cosmology, AIC, BIC and DIC are often used and, although this might be due to the natural delay between different fields start communicating with each other, I would like to have them present in my research.

As such, I would like to know if it's possible using arviz, or even in Stan directly (although this is not the Stan forums), to obtain the maximum value for the likelihood such that I can compute AIC, BIC and DIC.

[jpmvferreira]

I am also curious as to why they're not shipped with arviz unlike WAIC and LOO. Is it because of the reasons raised previously?
Naturally, any discussion regarding using or not the 3 previously mentioned criteria is more that welcome, but having a background in physics most of the things I read are taken as a "fact", rather than something which I can argue or properly understand, unfortunately.

[ahartikainen]

Currently we have functionality only for mcmc, so we are missing metrics for e.g. optimize and variational inference.

AIC, DIC etc are defined on max likelihood, and there is no guarantee that mcmc ever samples near that.

You could use optimize method to find maximum likelihood value with pystan3 (use scipy optimizer + gradient and lp info) or use optimize in cmdstanpy

But to use optimize you need to be careful that the correct density is optimized and tbh I don't remember the correct procedure right now (cc @OriolAbril @avehtari )

[jpmvferreira]

Okay so you suggestion is to either define the likelihood in Python and use scipy, or work with the cmdstanpy package along side pystan in order to find the maximum likelihood, correct?
I also saw that you are working with pystan3, so for you to suggest this approach, I assume it is because this feature not available there yet, correct?

Also, based on the discussion over at the Stan forums entitled How to get best fit value of the parameters from pystan.StanModel.sampling result?, you opposed the usage of MLE.
Would MAP make any difference? Or your concerns are regarding the fact that we're using a single value whereas we should be working with the full distribution, which I assume is the case for WAIC and LOO-CV.
If so, considering that there's a bunch of cosmologists using AIC, BIC and DIC often (at least from what I can tell), would you advise against? Even considering that (again, at least from what I can tell) our a posteriori distributions look fairly well behaved?

[ahartikainen]

You still can use Stan model, just plug it in some solver

https://discourse.mc-stan.org/t/speed-of-evaluating-gradients-of-log-probabilities-in-pystan-2-x-vs-3/22303/10

[ahartikainen]

Stan outputs penalized MLE with optimize, calculating MAP is harder (user needs to make sure that optimize is going for MAP)

Point estimates against likelihood in higher dimensions is not always the best idea. If you can do mcmc, you can use posterior mean as "point estimate" but I would also recommend using the rest of the mcmc sample to do your decisions.

[jpmvferreira]

So, if I understood correctly, I can use httpstan, as in the post you have linked, to obtain the MLE using optimize, yet you recommend that I use the posterior mean instead?
The posterior mean, just to see if I have things straight, is the mean value of θ, correct? Does that happen to be the most likely value for the marginalized 1D constrain?

[avehtari]

1. DIC uses posterior mean and not maximum likelihood (but not recommended as LOO is more reliable)
2. If priors are used, you can't get maximum likelihood from MCMC
3. If priors are not used (which is not recommended), you could approximate the maximum from MCMC in low dimensions (but that would require more work than using LOO)
4. If you are doing something Bayesian, you should forget AIC
5. If you are using priors (which is recommended), you should forget BIC
6. If you are making predictions using only the posterior mean (which is not usually recommended) you could use DIC
7. As ArviZ is "Exploratory analysis of Bayesian models", I think we shouldn't add AIC and BIC to confuse the users to think that it's ok to use them. DIC is not worth adding as LOO and WAIC are available.
8. If you want to use maximum likelihood or posterior mean, even then using cross-validation (LOO or K-fold-CV) would be better
9. If needed I can provide references which you can cite when telling in the article why AIC, BIC and DIC would be invalid to use (even if bunch of people have been using them)
10. If you tell more about your model, I'm happy to help more to check different assumptions and find the best solution for you (it is possible that it really doesn't matter what you use). Since you are using pystan this discussion could happen in Stan discourse
11. Sorry that you get the same recommendations from ArviZ and Stan developers, but we do have some overlap 😂

[jpmvferreira]

Thanks for all of the tips provided!
They all seem to hint that I shouldn't be using AIC, BIC or DIC, and should use LOO (and WAIC as well?) instead.

9. I will accept your offer, and I would like to know a little bit more as to what I should be using, or not, and why.

10. Regarding the model, it's rather simple, and I'll give you a brief overview and share the .stan model file.
    In short, we have gravitational wave events in which we measure the redshift (z) and the luminosity distance (dL(z)) where the later comes with an error associated with it.
    We now take a simple model which has two free parameters, h (≥ 0) and Omega_m (0 ≤ Omega_m ≤ 1), that says that the luminosity distance function should be given by (if you have the dark theme click on the image below)
    
    where
    
    We then assume a Gaussian likelihood and wide Gaussian priors on each parameter.

11. It's great to have recommendation from people that work in the field, overlap means more experience (or bias, hehe)!

Just to illustrate my point, I will show you the following article: https://arxiv.org/abs/2104.15123
Here, the authors considered a given cosmological model, and then using different sets of data, each of which constrained a different set of parameters for that model, computed the AIC, BIC and DIC, and compared it to a base model (ΛCDM).
This was essentially my intention, to be able to develop a model, constrain it, and then comparing it with others found in the literature.
This analysis is somewhat standard in my field (take "standard" with a big grain of salt as I'm studying it for less than a year), so being aware of when to use and not to use these selection criteria is of great value, and also how to obtain them.

[OriolAbril]

It might be worth adding that AIC, DIC, WAIC and LOO are all trying to estimate the same quantity, the differences lie on how they approximate that quantity and what assumptions are needed for these approximations to hold. I'd recommend starting with https://link.springer.com/article/10.1007/s11222-013-9416-2/ and then https://link.springer.com/article/10.1007/s11222-016-9696-4 (both also on arxiv)

[avehtari]

Hi,
9. I recommend to ask in Stan discourse instead of github (even if this repo is called discussions). Discourse has also support for equations. When asking you can tag me so I notice it (although it can still take days before I have time to answer, but there are many others who can answer these questions, too). Did Oriol's answer help?
10. For such a simple model it might be sufficient to look at the posterior directly. LOO is useful if the parameters have such posterior correlations that interpretation of the marginals is not easy. Of course, LOO can be used also to estimate the predictive performance using a domain specific utility/loss that might help to communicate the relevance of the model.