Automatic Parameterization of Scipy's Random Number Generators

I am working on a new version of Model-Builder, which will incorporate Uncertainty analysis features based on my BIP package.

One key aspect of doing uncertainty analysis, is the ability to attribute prior distributions to the model's parameters, and programatically, this is easily accomplished through the use of scipy.stats random number generators which provide a large variety of both continuous and discrete distributions. However, since Model-Builder is a GUI-based software, I need to offer a unified interface for the parameterization of such distribution.

This should not be as hard as it seems because in the stats package, all distributions are either a subclass of stats.rv_continuous or stats.rv_discrete, and thus [almost] all of them can be parameterized by calling the class constructor with a tuple of shape parameters, the length of which will vary from distribution to distribution but is available as the attribute numargs. This is exemplified in the doc string of all the distributions as the following pseudo-doctest:

>>> import matplotlib.pyplot as plt
>>> numargs = norm.numargs
>>> [ ] = [0.9,]*numargs
>>> rv = norm()

Display frozen pdf

>>> x = np.linspace(0,np.minimum(rv.dist.b,3))
>>> h=plt.plot(x,rv.pdf(x))

So I set out to test if this would actually work for all generators. If it worked for at least the most important ones, it would mean I could just have the user provide me with a -tuple of numbers, and pass them to whichever distribution the user chose. To run this test I came up with the following script:

# -*- coding: utf-8 -*-
import pylab as P
import numpy as np
f1 = P.figure()
f2 = P.figure()
a1 = f1.add_subplot(111)
a1.set_title("Continuous Distributions")
a2 = f2.add_subplot(111)
a2.set_title("Discrete Distributions")
for d in dir(st):
   obj = getattr(st, d)
   if isinstance(obj, st.rv_continuous):
       try:
           numargs = obj.numargs
           shape_pars = [0.9]*numargs
           rv = obj(shape_pars)
           x = np.linspace(0,np.minimum(rv.dist.b,3))
           h=a1.plot(x,rv.pdf(x),  label=obj.name)
           a1.legend()
           print "+++>", obj.name
       except:
           print "--->", obj.name
   if isinstance(obj, st.rv_discrete):
       try:
           numargs = obj.numargs
           shape_pars = [0.9]*numargs
           rv = obj(shape_pars)
           x = np.linspace(np.maximum(-3, rv.dist.a),np.minimum(rv.dist.b,3))
           h=a2.plot(x,rv.pmf(x),  label=obj.name)
           a2.legend()
           print "+++>", obj.name
       except:
           print "--->", obj.name
P.show()

When I ran this simple script, I not got the plots of each probability distribution that conformed with the parameterization via the shape tuple, but also a report which did (+++>) or did not(--->) conform. Here are the results:

+++> alpha
---> anglit
---> arcsine
+++> bernoulli
---> beta
---> betaprime
---> binom
---> boltzmann
+++> bradford
---> burr
---> cauchy
+++> chi
+++> chi2
---> cosine
+++> dgamma
+++> dlaplace
+++> dweibull
+++> erlang
---> expon
+++> exponpow
---> exponweib
---> f
+++> fatiguelife
+++> fink
+++> foldcauchy
+++> foldnorm
+++> frechet_l
+++> frechet_r
+++> gamma
---> gausshyper
---> genexpon
+++> genextreme
---> gengamma
+++> genhalflogistic
+++> genlogistic
+++> genpareto
+++> geom
---> gilbrat
+++> gompertz
---> gumbel_l
---> gumbel_r
---> halfcauchy
---> halflogistic
---> halfnorm
---> hypergeom
---> hypsecant
+++> invgamma
+++> invnorm
+++> invweibull
---> johnsonb
---> johnsonsu
+++> ksone
---> kstwobign
---> laplace
---> levy
---> levy_l
---> levy_stable
+++> loggamma
---> logistic
+++> loglaplace
+++> lognorm
+++> logser
+++> lomax
---> maxwell
---> mielke
+++> nakagami
---> nbinom
---> ncf
---> nct
---> ncx2
---> norm
+++> pareto
+++> planck
+++> poisson
+++> powerlaw
---> powerlognorm
+++> powernorm
---> randint
---> rayleigh
+++> rdist
+++> recipinvgauss
---> reciprocal
+++> rice
---> semicircular
+++> t
+++> triang
+++> truncexpon
---> truncnorm
+++> tukeylambda
---> uniform
+++> vonmises
---> wald
+++> weibull_max
+++> weibull_min
+++> wrapcauchy
+++> zipf

From the list above, it is clear that although the majority passed the test, many important ones did not (n, beta, uniform, etc.). This don't mean those cannot be used or won't work, but only that they cannot be parameterized in the same way as the rest. The norm distribution, for instance, has numargs = 0, (WTF?)

Although initially I was pretty excited about the possibility of handling every probability distribution in the same standard way, in the end I was disappointed that the Scipy people, after having come up with a such nice design for the stats package, did not think about running their own doc-tests (or making them runnable in the first place).

I am still looking for a solution to this. If you have an idea of how to elegantly work around this issue, please comment below.

Object Oriented Bayesian Inference

A couple of days ago, while working on a chapter of my book, I came across a manuscript from Kerman and Gelman, from Columbia University talking about a fully Bayesian computing environment. They argue that the basic building block for such a system, would be a "random variable" data type.

I liked their paper. Especially because I was already thinking along those lines for a while. So I decided to give it a shot. Well it was easy, since I already had a nice library for calculating likelihood functions, I didn't have to implement that. I came down to about a hundred lines of Python code (not counting the likelihood library)

I was pretty impressed with it, I had a object that I could specify from a prior distribution, that I could throw data at it and it would return a posterior distribution, a likelihood function, samples from the posterior, etc.

Pretty nice! It also allows the recursive inference typical of Bayesian statistics by which, if available, the posterior becomes the prior when new data is available.

It can still be extended by inhering from basic types such a numpy array. This would allow for this object to be used as a full-fledged type and participate in all sorts of mathematical operations...

Since the code is a bit long and blogger sucks at posting code, I am simply adding a UML representation of the class. If you are interested in discussing about this, getting access to the code or make suggestions, please email me or just leave a comment.