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.

Converting C code to Python

This may sound like a weird thing to do, but actually, I have craved for something like that every time I have to read C code. Converting C code to Python, can not only help us understand code more easily, but also turn non perfomance-critical code easier to maintain. C can be easy enough to read if well-written and formatted (indented) adequately, however having to move back and forth between .h and .c files drives me mad sometimes.

You can imagine my joy when I came across "ctopy", create by none other than Eric Raymond himself!

The first problem I ran into was the fact that ctopy was not available on my Ubuntu. No problem, a quick search through Eric's web site lead me to it:
http://www.catb.org/~esr/ctopy/

Naturally, ctopy is written in python and my immediate impulse was to download and test it.

For test code I decided to download a small C program from the "computer language shoutout" site:

http://shootout.alioth.debian.org/u64q/program.php?test=mandelbrot&lang=gcc&id=6

Unfortunately even for this simple program, ctopy didn't do a good job. I am not including the results here. but you can easily replicate them with the command:

cat mandelbrot.c | indent |./ctopy

The man page of ctopy recommends to pipe the C code through indent, because it relies on the indentation of the C code to do the translation.

Well as Raymond himself says in the beggining of ctopy's man page, its a quick and dirty translator, which requires a human to finish the job.

Anyway, that's how far my dream of a "Google Translate" for code went. I decided to post the links in case someone decides to continue where Raymond stopped.

Utopia

Interesting post, But hoping to invert those proportions, i.e., making the smarter guys the 80% is harder that having a society where everyone earns more money than they can spend... Pretty next to impossible...;-)

in reference to: Coding Horror: The Two Types of Programmers (view on Google Sidewiki)