Hi, I am using Hierarchical structure for linear model with two data samples and 5 (excluding the intrinsic scatter) parameters to be determined out which 4 are common between the two samples. I want to calculate the BIC or the WAIC for this (any goodness of fit criteria). I am using the default NUT sampler.
But when I use pm.waic(trace), I get the following error:
TypeError Traceback (most recent call last)
in
18 logp = trace[‘model_logp’]
19 print(np.amax(logp))
—> 20 pm.waic(trace)
~/anaconda3/lib/python3.7/site-packages/pymc3/stats/init.py in wrapped(*args, **kwargs)
22 )
23 kwargs[new] = kwargs.pop(old)
—> 24 return func(*args, **kwargs)
25
26 return wrapped
~/anaconda3/lib/python3.7/site-packages/arviz/stats/stats.py in waic(data, pointwise, scale)
1162 )
1163 if “log_likelihood” not in inference_data.sample_stats:
-> 1164 raise TypeError(“Data must include log_likelihood in sample_stats”)
1165 log_likelihood = inference_data.sample_stats.log_likelihood
1166
TypeError: Data must include log_likelihood in sample_stats
Since this is the trace objects, It contains ‘model_logp’ but I am still getting this error. Any help or any way to calculate BIC.
1 Like
What versions of PyMC3 and ArviZ do you have installed?
To begin with, I’d recommend updating both to the latest version and then explicitely converting the trace to InferenceData (more details on trace to inference data conversion here):
import arviz as az
...
idata = az.from_pymc3(trace, model=<...>) # model arg is not necessary but highly recommended
print(idata)
# make sure that the log_likelihood group is present
az.waic(idata)
Stats and plots are delegated to ArviZ. Everything should still work when passing traces directly, but it is recommended to explicitly convert the results for better performance and debugging.
2 Likes
In addition to updating to latest ArviZ, I’d also recommend reading Calculating WAIC for models with multiple likelihood functions. There are many cases where multiple ways to calculate waic are correct and therefore ArviZ cannot automatically select one of the options for you as all of them are correct, choosing which one to use depends on what question should be answered.
2 Likes
Hi @OriolAbril,
Thank much for your answers.
I have updated pymc3 and arviz but my model contains two likelihood functions in hierarchical model so I am reading your post for the solution.
However, on a side note, how can one get log-likelihood values in case of one likelihood? I think I can combine my two likelihoods into one as their sum and use that hardcoded customized likelihood function with DensityDist. Since sampler stats contain only `model_logp, what is the way to get likelihood values.
Thanks again.
Pointwise log likelihood values can be obtained by evaluating logp_elemwise for all observed variables, here is how it is done in ArviZ. I would not recommend this approach though, adding a couple lines to sum or concat the multiple likelihoods will most surely be easier to implement and to understand when coming back to the code.
The same is probably true when trying to rewrite the model so it has only one observed variable, not only you are bound to end up with harder to understand code but you could also end up having a multimodal distribution and making sampling harder and slower.
We have improved waic and loo since I wrote that post, so now things should be a little bit simpler. var_name argument has been added. Therefore, the last example can be simplified to:
az.waic(idata, var_name="away_points")
And the first example could now be:
idata.log_likelihood["log_lik_sum"] = log_lik.home_points + log_lik.away_points
az.waic(idata, var_name="log_lik_sum")
This makes clearer what is waic being calculated for and makes it easier to calculate waic for different approaches without having to repeat any calculation, only changing the var_name is enough.
If you were interested in the log likelihood value (which is not useful to calculate waic but can have other uses) in addition to the pointwise ones you should read this other answer
2 Likes
HI @OriolAbril,
Thanks a lot for your detailed reply. I think this should solve my problem. I can calculate the WAIC like the first example summing the two observables.
1 Like