As noted yesterday, I figured out how to sample PyMC chains in parallel recently. I’ve also been working on some plotting functions to assess chain convergence and autocorrelation, goodness of model fit, and to illustrate a fitted (multilevel GRT) model.

In assessing chain convergence, you can calculate the Gelman-Rubin statistic , which is, in essence, an F-like statistic providing an ANOVA-like test of any overall differences between the values in your chains. PyMC has a number of useful diagnostic tools built in, including a gelman_rubin method. It also has methods for calculating and plotting autocorrelation functions and histograms of chains, but for some reason, with the model I’m dealing with, the out-of-the-box plotting functions aren’t working well for me. So, I wrote my own, and for each parameter, it produces a figure like this (which is illustrating , autocorrelation (using the statsmodels package time series analysis acf function), and each of three chain’s distribution of values for a parameter called mu_mx):

I’m pretty happy with it, even if there are still some kinks to be worked out (e.g., the x-axis label for the top left panel is clipped by the bottom left panel; for some parameters, the tick labels under the histogram overlap).

I also wrote a function to plot observed and predicted (identification-confusion) response probabilities for the model I’m working with:

The vertical error bars indicate 95% highest density intervals, with the circles indicating the mean posterior value (on the y-axis) and the observed value (on the x-axis). The closer the symbols to the diagonal line, the better the fit.

Finally, I wrote a function to illustrate the fitted model for each individual subject (the small panels) and the group-level model (the big panel). These two figures are for models fit to the non-speech auditory perception data presented in my 2009 paper with Lentz and Townsend (see my CV for more details, if you’re so inclined). This one illustrates the multilevel GRT model fit to the frequency-by-duration, broadband noise stimuli:

And this one illustrates the model fit to the pitch-by-timbre (F0-by-spectral-prominence-location) data:

Again, there are probably some ways these could be improved, but overall I’m quite happy with how they look. I’m using my recently-acquired ellipse-plotting knowledge for its intended purpose (for the purpose *I* intended for it, anyway), which is nice, and it’s very easy to make this kind of multi-panel plot with matplotlib, which is also nice.

It feels good to have this model functioning in a way that others could, in principle, use. The two papers I’ve published using this model presented analyses done with WinBUGS, and the only way I could figure out how to get WinBUGS to fit this model was by feeding it a giant array of pre-calculated bivariate normal CDF values and using trilinear interpolation. The model works, but it’s ungainly. I don’t think I would want anyone else to mess with it, and I don’t imagine others would particularly want to do so. This isn’t a recipe for reproducible research.

I never could figure out how to get this model to work in JAGS, which is a superior, cross-platform version of BUGS. JAGS would just hang when I tried to feed it the giant CDF array (BUGS would always *seem* like it was hanging, but if I went away for long enough, it would have a model fit for me – no such luck with JAGS).

I haven’t tried very hard (at all, really) to get this model working with Stan, though I did request multivariate normal CDFs as a feature (it’s on the to-do list, so maybe I’ll come back to work on this model in Stan later on). Given that it’s functioning in PyMC, I don’t feel much need to get it working in Stan right now.