MVSDIST in SciPy

My pathway to probability theory was a little tortured. Like most people, I sat through my first college-level "Statistics" class fairly befuddled. I was good at math and understood calculus pretty well. As a result, I did well in the course, but didn't feel that I really understood what was going on. I took a course that used as its text this book by Papoulis. Now the text is a great reference for me, but at the time I didn't really understand the point of most of the more theoretical ideas. It wasn't until later after I had studied measure theory, and understood more of the implications of the set-theory studies of George Cantor that I began to see the significance of a Borel algebra and why some of the complexity was necessary from a foundational perspective.

I still believe, however, that diving into the details of measure theory is over-kill for introducing probability theory. I've been convinced by E.T. Jaynes that probability theory is an essential component of any education and as such should be presented in multiple ways at multiple times and certainly not thrown at you as "just an application of measure theory" the way it sometimes is in math courses. I think this is improving, but there is still work to do.

What typically still happens is that people get their "taste" of probability theory (or worse, their taste of "statistics") and then move on not ever really assimilating the lessons in their life. The trouble is everyone must deal with uncertainty. Our brains are hard-wired to deal with it --- often in ways that can be counter-productive. At its core, probability theory is just a consistent and logical way to deal with uncertainty using real numbers. In fact, it can be argued that it is the only way to deal with uncertainty.

I've done a bit of experimentation over the years and dealt with a lot of data (MRI, ultrasound, electrical impedance data). In probability theory, I found a framework for understanding what the data really tells me which led me to spend several years studying inverse problems. There are a lot of problems that can be framed as inverse problems. Basically, inverse problem theory can be applied to any problem where you have data and you want to understand what the data tells you. To apply probability theory to solve an inverse problem you have to have some model that determines how what you want to know leads to the data you've got. Then, you basically invert the model. Bayes' theorem provides a beautiful framework for this inversion.

The result of this Bayesian approach to inverse problems, though, is not just a number. It is explicitly a probability density function (or probability mass function). In other words, the result of a proper solution to an inverse problem is a random variable, or probability distribution. Seeing the result of any inverse problem as a random variable changes the way you think about drawing conclusions from data.

Think about the standard problem of fitting a line to data. You plug-and-chug using a calculator or a spreadsheet (or a function call in NumPy), and you get two numbers (the slope and intercept). If you properly understand inverse problems as requiring the production of a random variable, then you will not be satisfied with just these numbers. You will want to know, how certain am I about these numbers. How much should I trust them? What if I am going to make a decision on the basis of these numbers? (Speaking of making decisions, someday, I would like to write about how probability theory is also under-utilized in standard business financial projections and business decisions).

Some statisticians when faced with this regression problem will report the "goodness" of fit and feel satisfied, but as one who sees the power and logical simplicity of Bayesian inverse theory, I'm not satisfied by such an answer. What I want is the joint probability distribution for slope and intercept based on the data. A lot of common regression techniques do not provide this. I'm not going to go into details regarding the historical reasons for why this is. You can use google to explore some of the details if you are interested. A lot of it comes down to the myth of objectivity and the desire to eliminate the need for a prior which Bayesian inverse theory exposes.

As an once very active contributor to SciPy (now an occasional contributor who is still very interested in its progress), I put in the scipy.stats package a few years ago a little utility for estimating the mean, standard deviation, and variance from data that expresses my worldview a little bit. I recently updated this utility and created a function called mvsdist. This function finally returns random variable objects (as any good inverse problem solution should!) for the Mean, Standard deviation, and Variance derived from a vector of data. The assumptions are 1) the data were all sampled from a random variable with the same mean and variance, 2) the standard deviation and variance are "scale" parameters, and 3) non-informative (improper) priors.

The details of the derivation are recorded in this paper. Any critiques of this paper are welcome as I never took the time to try and get formal review for it (I'm not sure where I would have submitted it for one --- and I'm pretty sure there is a paper out there that already expresses all of this, anyway).

It is pretty simple to get started playing with mvsdist (assuming you have SciPy 0.9 installed). This function is meant to be called any time you have a bunch of data and you want to "compute the mean" or "find the standard deviation." You collect the data into a list or NumPy array of numbers and pass this into the mvsdist function:

>>> from scipy.stats import mvsdist
>>> data = [9, 12, 10, 8, 6, 11, 7]
>>> mean, var, std = mvsdist(data)

This returns three distribution objects which I have intentionally named mean, var, and std because they represent the estimates of mean, variance, and standard-deviation of the data. Because they are estimates, they are not just numbers, but instead are (frozen) probability distribution objects. These objects have methods that let you evaluate the probability density function: .pdf(x), compute the cumulative density function: .cdf(x), generate random samples drawn from the distribution: .rvs(size=N), determine an interval that contains some percentage of the random draws from this distribution: .interval(alpha), and calculate simple statistics: .stats(), .mean(), .std(), .var().

In this case, consider the following example:

>>> mean.interval(0.90)
(7.4133999449331132, 10.586600055066887)
>>> mean.mean()
9.0
>>> mean.std()
0.99999999999999989


>>> std.interval(0.90)
(1.4912098929401241, 4.137798046658852)
>>> std.mean()
2.4869681414837035
>>> std.std()
0.90276766847572409

Notice that once we have the probability distribution we can report many things about the estimate that provide for not only the estimate itself, but also any question we might have regarding the uncertainty in the estimate. Often, we may want to visualize the probability density function as is shown below for the standard deviation estimate and the mean estimate

It is not always easy to solve an inverse problem by providing the full probability distribution object (especially in multiple dimensions). But, when it's possible, it really does provide for a more thorough understanding of the problem. I'm very interested in SciPy growing more of these kinds of estimator approaches where possible.

Zen of NumPy

While I was on-site working for a client, one of the developers I worked with would begin each day with a brief discussion of one of the tenets from the "Zen of Python." For those who are not familiar with this little pearl of Python goodness. You can find the "Zen of Python" as an Easter egg inside a Python distribution:

>>> import this
The Zen of Python, by Tim Peters

Beautiful is better than ugly.
Explicit is better than implicit.
Simple is better than complex.
Complex is better than complicated.
Flat is better than nested.
Sparse is better than dense.
Readability counts.
Special cases aren't special enough to break the rules.
Although practicality beats purity.
Errors should never pass silently.
Unless explicitly silenced.
In the face of ambiguity, refuse the temptation to guess.
There should be one-- and preferably only one --obvious way to do it.
Although that way may not be obvious at first unless you're Dutch.
Now is better than never.
Although never is often better than *right* now.
If the implementation is hard to explain, it's a bad idea.
If the implementation is easy to explain, it may be a good idea.
Namespaces are one honking great idea -- let's do more of those!

The Zen of Python is often quoted from one Python user to another in trying to communicate something of the essence of what makes programming in Python different. While we were discussing one of the points, one of my co-workers suggested that there should be a "Zen of NumPy". This isn't the first time I've heard that suggestion. Actually David Morrill (author of Traits) was the first person who suggested there should be a book about the "Zen of NumPy." I totally agree with him. The only problem is that everybody involved with NumPy has apparently been too busy to write one :-)

With this idea in my mind, when it came time to give a talk on NumPy at the New York Python Meetup group in Manhattan, I decided to create a first-draft of the Zen of NumPy. The phrases are included on one slide in the deck shared here.

I'm interested in feedback on these before proposing them for placement as

numpy.this

Here is my attempt at a "Zen of NumPy"

Strided is better than scattered
Contiguous is better than strided
Descriptive is better than imperative (use data-types)
Array-oriented is often better than object-oriented
Broadcasting is a great idea -- use where possible
Vectorized is better than an explicit loop
Unless it’s complicated --- then use numexpr, weave, or Cython
Think in higher dimensions

I think there are useful edits as well as more statements that could be added to this list. Your feedback is welcome.

A New Blog

Lately I have been finding a need to have a voice --- an authentic voice. A voice, which I occasionally expressed in the days when I had the time to be more active on open source mailing lists (SciPy, NumPy, and even Python itself). When I was younger, I didn't have as many endearing entanglements to the future that depend on my present. As a result, I could spend much time pursuing efforts that gave me a tremendous sense of accomplishment.

For as long as I can remember, I have been driven by discovery. Much to their annoyance, I would constantly ask my parents and 9 siblings "Why?" I used to be quite proud of myself as they would relate these stories of my inquisitive childhood at family gatherings. My particular combination of infused biochemistry that led to my knowledge addiction certainly drove most pursuits during my formative years, and this has had a strong impact in my life.

During my nearly 40 years, however, I have encountered an impressive cadre of awe-inspiring people each uniquely different. This has led me to conclude that it is not the particular current physical emergence that I find myself in. Rather, it is the particular use I am making of it. Do I pursue an agenda that barely extends beyond my internal neurobiology, or do I use my combination of skills and knowledge to seek a wider consistency that can harmonize with a beautifully complex society.

Earlier tonight, I listened to technology leaders and entrepreneurs tell their view of what society would be like if their respective companies were wildly successful. I listened to this message in a stunning lecture hall in Peterhouse at Cambridge University. While they each brought a distinct perspective, their unifying message was the power of technology to change the world.

Search for "Silicon valley comes to Cambridge" in a few days to get a summary and possibly even video of the talks. Megan Smith from Google (www.google.org) spoke of the power of big data to solve social injustices such as the sexual exploitation of children. Reid Hoffman, co-founder of LinkedIn, spoke of the power of inter-connectedness to solve big problems by bringing the right people together quickly.

Other people spoke and gave interesting perspectives including Mike Lynch, founder of Autonomy, who gave a wonderful talk on the importance of meaningful interaction with data so that our lives are enhanced and not enslaved by the explosion of data and technology. He also gave a tribute to Thomas Bayes. By looking on his site, I noticed that he gives similar props to Claude Shannon. I'm already impressed. These are two thinkers who were able to present important concepts that remain under-appreciated.

I do think it's important what people think. The ideas we carry in our heads are critical. It is these ideas which drive our necessarily individual pursuits and can lead to disharmony. I like to pass along useful information, colored of course by my own experiences and perspective in the simple and perhaps naive hopes that sustainable, lasting solutions can be discovered.

Most of my posts will be technical, as I am hoping to use this forum as a way to write about the thoughts I am having in my own attempt to hone and pare them. In particular, most of these posts will be about technology that I am involved with or have some exposure to. Upcoming posts include "The Zen of NumPy", "7 Heresies of Technical Computing", and "What I've learned from SciPy and Open Source"

If you happen to come across these musings, your feedback is welcome. I would love to hear about your experiences with any thoughts that are covered in my posts.

-Travis