Thanks for the tip! I looked into it and turns out Turing interfaces with AdvancedVI.jl, which has a built in function for extracting the ELBO, detailed in its doc.
The doc example only passes in a custom log joint probability function to get the maximized ELBO, but it seems like passing in a Turing model works too, fortunately.
using Turing, AdvancedVI, Turing: Variational
d = 2; n = 100;
observations = randn((d, n));
advi = ADVI(10, 10000)
## turing model
@model function model(y, d)
μ ~ MvNormal(ones(d))
y ~ MvNormal(μ, ones(d))
end;
q = vi(model(observations, d), advi)
#I'm not sure yet what the last number does in elbo()
AdvancedVI.elbo(advi, q, model(observations, d), 10000)
The contributors seem to be interested in redesigning it at some point, so the solution above might not be applicable to future readers.