Calculating Bayes Factor using Stan

Hello Stan Community!

As we are all aware that, Stan gives an estimate of the log marginal likelihood (lp__) in the output along with the estimates of the other parameters.

My question is, can I use this estimate of lp__ to calculate Bayes Factors (exp(lp__Null - lp__Alternative) ? Any help is much appreciated!

Thank you!

Use the bridgesampler package.

1 Like

Thank you very much for your reply! Actually, I have been using the bridgesampling package and facing weird issues with it . Like, after fitting the model if I exit R and reload the .Rdata with the stanfit object having the posterior draws, the bridge_sampler won’t work and I have to refit and run the function in that same session. That is why I was wondering why I could not use the lp__ from Stan.

Anyway, thank you very much!

Anything using lp__ directly is not going to be a good estimator of a Bayes Factor.

The reloading thing is actually a rstan issue rather than a bridgesampling one. When you close your R session, the dynamic shared object that is loaded in R’s memory gets lost, and thus it is unable to transform the draws to the unrestricted parameter space where the bridge sampling happens.

1 Like

I see! Thank you very much for the clarification!

lp__ is the logarithm of the numerator of Bayes rule (possibly with constants dropped depending on how the model is coded). It’s true that the log of the integral of the numerator wrt the parameters is the log marginal likelihood but I hesitate to refer to lp__ as an estimate of the log marginal likelihood as it’s rarely a good idea to use it that way.


Thank you for the explanation of the lp__. I now understand why one should not use that as an estimate of the log marginal likelihood.

To add to this answer, the bridge_sampler method for stanfit objects has an argument stanfit_model which accepts a second empty (i.e., no samples) stanfit object created in the current R session. This allows to use \samples from a previous session. To create this empty object just call stan with chains = 0. This will simply compile the model which can then be used.


Thanks! This is really helpful and I always give a name to my stan-model. Thank you very much!