What is a distribution?

This article uses object-oriented programming to explore of one of the most useful concepts in statistics, distributions.  The code is in a Jupyter notebook.

You can read a static version of the notebook on nbviewer.

OR

You can run the code in a browser by clicking this link and then selecting distribution.ipynb from the list.

The following is a summary of the material in the notebook, which you might want to read before you dive into the code.

Random processes and variables

One of the recurring themes of my books is the use of object-oriented programming to explore mathematical ideas. Many mathematical entities are hard to define because they are so abstract. Representing them in Python puts the focus on what operations each entity supports — that is, what the objects can do — rather than on what they are.

In this article, I explore the idea of a probability distribution, which is one of the most important ideas in statistics, but also one of the hardest to explain. To keep things concrete, I'll start with one of the usual examples: rolling dice.

• When you roll a standard six-sided die, there are six possible outcomes — numbers 1 through 6 — and all outcomes are equally likely.
• If you roll two dice and add up the total, there are 11 possible outcomes — numbers 2 through 12 — but they are not equally likely.  The least likely outcomes, 2 and 12, only happen once in 36 tries; the most likely outcome happens 1 times in 6.
• And if you roll three dice and add them up, you get a different set of possible outcomes with a different set of probabilities. 

What I've just described are three random number generators.   The output from each generator is a random variable. And each random variable has probability distribution, which is the set of possible outcomes and the corresponding set of probabilities.

Representing distributions

There are many ways to represent a probability distribution. The most obvious is a probability mass function, or PMF, which is a function that maps from each possible outcome to its probability. And in Python, the most obvious way to represent a PMF is a dictionary that maps from outcomes to probabilities.

So is a Pmf a distribution? No. At least in my framework, a Pmf is one of several representations of a distribution. Other representations include the cumulative distribution function (CDF) and the characteristic function (CF).

These representations are equivalent in the sense that they all contain the same information; if I give you any one of them, you can figure out the others.

So why would we want different representations of the same information? The fundamental reason is that there are many operations we would like to perform with distributions; that is, questions we would like to answer. Some representations are better for some operations, but none of them is the best for all operations.

Here are some of the questions we would like a distribution to answer:

• What is the probability of a given outcome?
• What is the mean of the outcomes, taking into account their probabilities?
• What is the variance of the outcome?  Other moments?
• What is the probability that the outcome exceeds (or falls below) a threshold?
• What is the median of the outcomes, that is, the 50th percentile?
• What are the other percentiles?
• How can get generate a random sample from this distribution, with the appropriate probabilities?
• If we run two random processes and choose the maximum of the outcomes (or minimum), what is the distribution of the result?
• If we run two random processes and add up the results, what is the distribution of the sum?

Each of these questions corresponds to a method we would like a distribution to provide. But there is no one representation that answers all of them easily and efficiently.

As I demonstrate in the notebook, the PMF representation makes it easy to look up an outcome and get its probability, and it can compute mean, variance, and other moments efficiently.

The CDF representation can look up an outcome and find its cumulative probability efficiently.  And it can do a reverse lookup equally efficiently; that is, given a probability, it can find the corresponding value, which is useful for computing medians and other percentiles.

The CDF also provides an easy way to generate random samples, and a remarkably simple way to compute the distribution of the maximum, or minimum, of a sample.

To answer the last question, the distribution of a sum, we can use the PMF representation, which is simple, but not efficient.  An alternative is to use the characteristic function (CF), which is the Fourier transform of the PMF.  That might sound crazy, but using the CF and the Convolution Theorem, we can compute the distribution of a sum in linearithmic time, or O(n log n).

If you are not familiar with the Convolution Theorem, you might want to read Chapter 8 of Think DSP.

So what's a distribution?

The Pmf, Cdf, and CharFunc are different ways to represent the same information. For the questions we want to answer, some representations are better than others. So how should we represent the distribution itself?

In my implementation, each representation is a mixin; that is, a class that provides a set of capabilities. A distribution inherits all of the capabilities from all of the representations. Here's a class definition that shows what I mean:

class Dist(Pmf, Cdf, CharFunc):
    
    def __init__(self, d):
        """Initializes the Dist.
        
        Calls all three __init__ methods.
        """
        Pmf.__init__(self, d)
        Cdf.__init__(self, *compute_cumprobs(d))
        CharFunc.__init__(self, compute_fft(d))

When you create a Dist, you provide a dictionary of values and probabilities.
Dist.__init__ calls the other three __init__ methods to create the Pmf, Cdf, and CharFunc representations. The result is an object that has all the attributes and methods of the three representations.

From a software engineering point of view, that might not be the best design, but it is meant to illustrate what it means to be a distribution.

In short, if you give me any representation of a distribution, you have told me everything I need to answer questions about the possible outcomes and their probabilities. Converting from one representation to another is mostly a matter of convenience and computational efficiency.

Conversely, if you are trying to find the distribution of a random variable, you can do it by computing whichever representation is easiest to figure out.

So that's the idea.  If you want more details, take a look at the notebook by following one of the links at the top of the page.

Bayesian Statistics for Undergrads

Yesterday Sanjoy Mahajan and I led a workshop on teaching Bayesian statistics for undergraduates.  The participants were college teachers from around New England, including Norwich University in Vermont and Wesleyan University in Connecticut, as well as our neighbors, Babson College and Wellelsey College.

The feedback we got was enthusiastic, and we hope the workshop will help the participants design new classes that make Bayesian methods accessible to their students.

Materials from the workshop are in this GitHub repository.  And here are the slides:



The goal of the workshop is to show that teaching Bayesian statistics to undergrads is possible and desirable.  To show that it's possible, we presented three approaches:

1. A computational approach, based on my class at Olin, Computational Bayesian Statistics, and the accompanying book, Think Bayes.  This material is appropriate for students with basic programming skills, although a lot of it could adapted for use with spreadsheets.
2. An analytic approach, based on Sanjoy's class, called Bayesian Inference.  This material is appropriate for students who are comfortable with mathematics including calculus.
3. We also presented core material that does not depend on programming or advanced math --really just arithmetic.

Why Bayes?

Reasons the participants gave for teaching Bayes included:

1. Some of them work and teach in areas like psychology and biology where the limitations of classical methods have become painfully apparent, and interest in alternatives is high.
2. Others are interested in applications like business intelligence and data analytics where Bayesian methods are a hot topic.
3. Some participants teach introductory classes that satisfy requirements in quantitative reasoning, and they are looking for material to develop students' ability to reason with and about uncertainty.

I think these are all good reasons.  At the introductory level, Bayesian methods are a great opportunity for students who might not be comfortable with math to gradually build confidence with mathematical methods as tools for better thinking.

Bayes's theorem provides a divide-and-conquer strategy for solving difficult problems by breaking them into smaller, simpler pieces.  And many of the classic applications of Bayes's theorem -- like interpreting medical tests and weighing courtroom evidence -- are real-world problems where careful thinking matters and mistakes have consequences!

For students who only take a few classes in mathematics, I think Bayesian statistics is a better choice than calculus, which the vast majority of students will never use again; and better than classical statistics, which (based on my observation) often leaves students more confused about quantitative reasoning than when they started.

At the more advanced level, Bayesian methods are appealing because they can be applied in a straightforward way to real-world decision making processes, unlike classical methods, which generally fail to answer the questions we actually want to answer.

For example, if we are considering several hypotheses about the world, it is useful to know the probability that each is true.  You can use that information to guide decision making under uncertainty.  But classical statistical inference refuses to answer that question, and under the frequentist interpretation of probability, you are not even allowed to ask it.

As another example, the result you get from Bayesian statistics is generally a posterior distribution for a parameter, or a joint distribution for several parameters.  From these results, it is straightforward to compute a distribution that predicts almost any quantity of interest, and this distribution encodes not only the most likely outcome or central tendency; it also represents the uncertainty of the prediction and the spread of the possible outcomes.

Given a predictive distribution, you can answer whatever questions are relevant to the domain, like the probability of exceeding some bound, or the range of values most likely to contain the true value (another question classical inference refuses to answer).  And it is straightforward to feed the entire distribution into other analyses, like risk-benefit analysis and other kinds of optimization, that directly guide decision making.

I mention these advantages in part to address one of the questions that came up in the workshop.  Several of the participants are currently teaching traditional introductory statistics classes, and they would like to introduce Bayesian methods, but are also required to cover certain topics in classical statistics, notably null-hypothesis significance testing (NHST).

So they want to know how to design a class that covers these topics and also introduces Bayesian statistics.  This is an important challenge, and I was frustrated that I didn't have a better answer to offer at the workshop.  But with some time to organize my thoughts, I have a two suggestions:

Avoid direct competition

I don't recommend teaching a class that explicitly compares classical and Bayesian statistics.  Pedagogically, it is likely to be confusing.  Strategically, it is asking for intra-departmental warfare.  And importantly, I think it misrepresents Bayesian methods, and undersells them, if you present them as a tool-for-tool replacement for classical methods.

The real problem with classical inference is not that it gets the wrong answer; the problem is that is asks the wrong questions.  For example, a fundamental problem with NHST is that it requires a binary decision: either we reject the null hypothesis or we fail to reject it (whatever that means).  An advantage of the Bayesian approach is that it helps us represent and work with uncertainty; expressing results in terms of probability is more realistic, and more useful, than trying to cram the world into one of two holes.

If you use Bayesian methods to compute the probability of a hypothesis, and then apply a threshold to decide whether the theory is true, you are missing the point.  Similarly, if you compute a posterior distribution, and then collapse it to a single point estimate (or even an interval), you are throwing away exactly the information that makes Bayesian results more useful.

Bayesian methods don't do the same things better; they do different things, which are better.  If you want to demonstrate the advantages of Bayesian methods, do it by solving practical problems and answering the questions that matter.

As an example, this morning my colleague Jon Adler sent me a link to this paper, Bayesian Benefits for the Pragmatic Researcher, which is a model of what I am talking about. 

Identify the goals

As always, it is important to be explicit about the learning goals of the class you are designing. Curriculum problems that seems impossible can sometimes be simplified by unpacking assumptions about what needs to be taught and why.  For example, if we think about why NHST is a required topic, we get some insight into how to present it: if you want to make sure students can read papers that report p-values, you might take one approach; if you imagine they will need to use classical methods, that might require a different approach.

For classical statistical inference, I recommend "The New Statistics", an approach advocated by Geoff Cumming (I am not sure to what degree it is original to him).  The fundamental idea of is that statistical analysis should focus on estimating effect sizes, and should express results in terms that emphasize practical consequences, as contrasted with statistical significance.

If "The New Statistics" is what we should teach, computational simulation is how.  Many of the ideas that take the most time, and seem the hardest, in a traditional stats class, can be taught much more effectively using simulation.  I wrote more about this just last week, in this post, There is Still Only One Test, and there are links there to additional resources.

But if the goal is to teach classical statistical inference better, I would leave Bayes out of it.  Even if it's tempting to use a Bayesian framework to explain the problems with classical inference, it would be more likely to confuse students than help them.

If you only have space in the curriculum to teach one paradigm, and you are not required to teach classical methods, I recommend a purely Bayesian course.  But if you have to teach classical methods in the same course, I suggest keeping them separated.

I experienced a version of this at PyCon this year, where I taught two tutorials back to back: Bayesian statistics in the morning and computational statistical inference in the afternoon.  I joked that I spent the morning explaining why the afternoon was wrong.  But the reality is that they two topics hardly overlap at all.  In the morning I used Bayesian methods to formulate real-world problems and answer practical questions.  In the afternoon, I helped people understand classical inference, including its limitations, and taught them how to do it well, if they have to.

I think a similar balance (or compromise?) could work in the undergraduate statistic curriculum at many colleges and universities.

There is still only one test

In 2011 I wrote an article called "There is Only One Test", where I explained that all hypothesis tests are based on the same framework, which looks like this:

Here are the elements of this framework:

1) Given a dataset, you compute a test statistic that measures the size of the apparent effect.  For example, if you are describing a difference between two groups, the test statistic might be the absolute difference in means.  I'll call the test statistic from the observed data 𝛿*.

2) Next, you define a null hypothesis, which is a model of the world under the assumption that the effect is not real; for example, if you think there might be a difference between two groups, the null hypothesis would assume that there is no difference.

3) Your model of the null hypothesis should be stochastic; that is, capable of generating random datasets similar to the original dataset.

4) Now, the goal of classical hypothesis testing is to compute a p-value, which is the probability of seeing an effect as big as 𝛿* under the null hypothesis.  You can estimate the p-value by using your model of the null hypothesis to generate many simulated datasets.  For each simulated dataset, compute the same test statistic you used on the actual data.

5) Finally, count the fraction of times the test statistic from simulated data exceeds 𝛿*.  This fraction approximates the p-value.  If it's sufficiently small, you can conclude that the apparent effect is unlikely to be due to chance (if you don't believe that sentence, please read this).

That's it.  All hypothesis tests fit into this framework.  The reason there are so many names for so many supposedly different tests is that each name corresponds to

1) A test statistic,

2) A model of a null hypothesis, and usually,

3) An analytic method that computes or approximates the p-value.

These analytic methods were necessary when computation was slow and expensive, but as computation gets cheaper and faster, they are less appealing because:

1) They are inflexible: If you use a standard test you are committed to using a particular test statistic and a particular model of the null hypothesis.  You might have to use a test statistic that is not appropriate for your problem domain, only because it lends itself to analysis.  And if the problem you are trying to solve doesn't fit an off-the-shelf model, you are out of luck.

2) They are opaque: The null hypothesis is a model, which means it is a simplification of the world.  For any real-world scenario, there are many possible models, based on different assumptions.  In most standard tests, these assumptions are implicit, and it is not easy to know whether a model is appropriate for a particular scenario.

One of the most important advantages of simulation methods is that they make the model explicit.  When you create a simulation, you are forced to think about your modeling decisions, and the simulations themselves document those decisions.

And simulations are almost arbitrarily flexible.  It is easy to try out several test statistics and several models, so you can choose the ones most appropriate for the scenario.  And if different models yield very different results, that's a useful warning that the results are open to interpretation.  (Here's an example I wrote about in 2011.)

More resources

A few days ago, I saw this discussion on Reddit.  In response to the question "Looking back on what you know so far, what statistical concept took you a surprising amount of effort to understand?", one redditor wrote

The general logic behind statistical tests and null hypothesis testing took quite some time for me. I was doing t-tests and the like in both work and classes at that time, but the overall picture evaded me for some reason. 

I remember the exact time where everything started clicking - that was after I found a blog post (cannot find it now) called something like "There is only one statistical test". And it explained the general logic of testing something and tied it down to permutations. All of that seemed very natural.

I am pretty sure they were talking about my article.  How nice!  In response, I provided links to some additional resources; and I'll post them here, too.

First, I wrote a followup to my original article, called "More hypotheses, less trivia", where I provided more concrete examples using the simulation framework.

Later in 2011 I did a webcast with O'Reilly Media where I explained the whole idea:

In 2015 I developed a workshop called "Computational Statistics", where I present this framework along with a similar computational framework for computing confidence intervals.  The slides and other materials from the workshop are here.

And I am not alone!  In 2014, John Rauser presented a keynote address at Strata+Hadoop, with the excellent title "Statistics Without the Agonizing Pain":

And for several years, Jake VanderPlas has been banging a similar drum, most recently in an excellent talk at PyCon 2016:

UPDATE: John Rauser pointed me to this excellent article, "The Introductory Statistics Course: A Ptolemaic Curriculum" by George W. Cobb.

UPDATE: Andrew Bray has developed an R package called "infer" to do computational statistical inference.  Here's an excellent talk where he explains it.