Table of Contents for

Data Analysis with Open Source Tools

The central formula for Bernoulli trials gives the probability of observing k Successes in N trials with Success probability p, and it is also known as the Binomial distribution (see Figure 9-1):

This should make good sense: we need to obtain k Successes, each occurring with probability p, and N – k Failures, each occurring with probability 1 – p. The term:

consisting of a binomial coefficient is combinatorial in nature: it gives the number of distinct arrangements for k successes and N – k failures. (This is easy to see. There are N! ways to arrange N distinguishable items: you have N choices for the first item, N – 1 choices for the second, and so on. However, the k Successes are indistinguishable from each other, and the same is true for the N – k Failures. Hence the total number of arrangements is reduced by the number of ways in which the Successes can be rearranged, since all these rearrangements are identical to each other. With k Successes, this means that k! rearrangements are indistinguishable, and similarly for the N – k failures.) Notice that the combinatorial factor does not depend on p.

Figure 9-1. The Binomial distribution: the probability of obtaining k Successes in N trials with Success probability p.

This formula gives the probability of obtaining a specific number k of Successes. To find the expected number of Successes μ in N Bernoulli trials, we need to average over all possible outcomes:

This result should come as no surprise. We use it intuitively whenever we say that we expect “about five Heads in ten tosses of fair coin” (N = 10, p = 1/2) or that we expect to obtain “about ten aces in sixty tosses of a fair die” (N = 60, p = 1/6).

Another result that can be worked out exactly is the standard deviation:

The standard deviation gives us the range over which we expect the outcomes to vary. (For example, assume that we perform m experiments, each consisting of N tosses of a fair coin. The expected number of Successes in each experiment is Np, but of course we won’t obtain exactly this number in each experiment. However, over the course of the m experiments, we expect to find the number of Successes in the majority of them to lie between and .

Notice that σ grows more slowly with the number of trials than does μ (σ ~ versus μ ~ N). The relative width of the outcome distribution therefore shrinks as we conduct more trials.

The primary reason why I place so much emphasis on the concept of Bernoulli trials is that it lends itself naturally to the development of mean-field models (see Chapter 8). Suppose we try to develop a model to predict the staffing level required for a call center to deal with customer complaints. We know from experience that about one in every thousand orders will lead to a complaint (hence p = 1/1000). If we shipped a million orders a day, we could use the Binomial distribution to work out the probability to receive 1, 2, 3,..., 999,999, 1,000,000 complaints a day and then work out the required staffing levels accordingly—a daunting task! But in the spirit of mean-field theories, we can cut through the complexity by realizing that we will receive “about Np = 1,000” complaints a day. So rather than working with each possible outcome (and its associated probability), we limit our attention to a single expected outcome. (And we can now proceed to determine how many calls a single person can handle per day to find the required number of customer service people.) We can even go a step further and incorporate the uncertainty in the number of complaints by considering the standard deviation, which in this example comes out to . (Here I made use of the fact that 1 – p is very close to 1 for the current value of p.) The spread is small compared to the expected number of calls, lending credibility to our initial approximation of replacing the full distribution with only its expected outcome. (This is a demonstration for the observation we made earlier that the width of the resulting distribution grows much more slowly with N than does the expected value itself. As N gets larger, this effect becomes more drastic, which means that mean-field theory gets better and more reliable the more urgently we need it! The tough cases can be situations where N is of moderate size—say, in the range of 10,..., 100. This size is too large to work out all outcomes exactly but not large enough to be safe working only with the expected values.)

Having seen this, we can apply similar reasoning to more general situations. For example, notice that the number of orders shipped each day will probably not equal exactly one million—instead, it will be a random quantity itself. So, by using N = 1,000,000 we have employed the mean-field idea already. It should be easy to generalize to other situations from here.

Figure 9-2. The Gaussian probability density.

Imagine you have a source of data points that are distributed according to some common distribution. The data could be numbers drawn from a uniform random-number generator, prices of items in a store, or the body heights of a large group of people.

Now assume that you repeatedly take a sample of n elements from the source (n random numbers, n items from the store, or measurements for n people) and form the total sum of the values. You can also divide by n to get the average. Notice that these sums (or averages) are random quantities themselves: since the points are drawn from a random distribution, their sums will also be random numbers.

Note that we don’t necessarily know the distributions from which the original points come, so it may seem it would be impossible to say anything about the distribution of their sums. Surprisingly, the opposite is true: we can make very precise statements about the form of the distribution according to which the sums are distributed. This is the content of the Central Limit Theorem.

The Central Limit Theorem states that the sums of a bunch of random quantities will be distributed according to a Gaussian distribution. This statement is not strictly true; it is only an approximation, with the quality of the approximation improving as more points are included in each sample (as n gets larger, the approximation gets better). In practice, though, the approximation is excellent even for quite moderate values of n.

This is an amazing statement, given that we made no assumptions whatsoever about the original distributions (I will qualify this in a moment): it seems as if we got something for nothing! After a moment’s thought, however, this result should not be so surprising: if we take a single point from the original distribution, it may be large or it may be small—we don’t know. But if we take many such points, then the highs and the lows will balance each other out “on average.” Hence we should not be too surprised that the distribution of the sums is a smooth distribution with a central peak. It is, however, not obvious that this distribution should turn out to be the Gaussian specifically.

We can now state the Central Limit Theorem formally. Let {x_i} be a sample of size n, having the following properties:

Then the sample average is distributed according to a Gaussian with mean μ and standard deviation . The approximation improves as the sample size n increases. In other words, the probability of finding the value x for the sample mean becomes Gaussian as n gets large:

Notice that, as for the binomial distribution, the width of the resulting distribution of the average is smaller than the width of the original distribution of the individual data points. This aspect of the Central Limit Theorem is the formal justification for the common practice to “average out the noise”: no matter how widely the individual data points scatter, their averages will scatter less.

On the other hand, the reduction in width is not as fast as one might want: it is not reduced linearly with the number n of points in the sample but only by . This means that if we take 10 times as many points, the scatter is reduced to only percent of its original value. To reduce it to 10 percent, we would need to increase the sample size by a factor of 100. That’s a lot!

Finally, let’s take a look at the Central Limit Theorem in action. Suppose we draw samples from a uniform distribution that takes on the values 1, 2,..., 6 with equal probability—in other words, throws of a fair die. This distribution has mean μ = 3.5 (that’s pretty obvious) and standard deviation (not as obvious but not terribly hard to work it out, or you can look it up).

We now throw the die a certain number of times and evaluate the average of the values that we observe. According to the Central Limit Theorem, these averages should be distributed according to a Gaussian distribution that becomes narrower as we increase the number of throws used to obtain an average. To see the distribution of values, we generate a histogram (see Chapter 2). I use 1,000 “repeats” to have enough data for a histogram. (Make sure you understand what is going on here: we throw the die a certain number of times and calculate an average based on those throws; and this entire process is repeated 1,000 times.)

The results are shown in Figure 9-3. In the upper-left corner we have thrown the die only once and thus form the “average” over only a single throw. You can see that all of the possible values are about equally likely: the distribution is uniform. In the upper-right corner, we throw the dice twice every time and form the average over both throws. Already a central tendency in the distribution of the average of values can be observed! We then continue to make longer and longer averaging runs. (Also shown is the Gaussian distribution with the appropriately adjusted width: , where n is the number of throws over which we form the average.)

I’d like to emphasize two observations in particular. First, note how quickly the central tendency becomes apparent—it only takes averaging over two or three throws for a central peak to becomes established. Second, note how well the properly scaled Gaussian distribution fits the observed histograms. This is the Central Limit Theorem in action.

The most predominant feature of the Gaussian density function is the speed with which it falls to zero as |x| (the absolute value of x—see Appendix B) becomes large. It is worth looking at some numbers to understand just how quickly it does decay. For x = 2, the standard Gaussian with zero mean and unit variance is approximately p(2, 0, 1) = 0.05 .... For x = 5, it is already on the order of 10^–6; for x = 10 it’s about 10^–22; and not much further out, at x = 15, we find p(15, 0, 1) ≈ 10^–50. One needs to keep this in perspective: the age of the universe is currently estimated to be about 15 billion years, which is about 4 · 10¹⁷ seconds. So, even if we had made a thousand trials per second since the beginning of time, we would still not have found a value as large or larger than x = 10!

Figure 9-3. The Central Limit Theorem in action. Distribution of the average number of points when throwing a fair die several times. The boxes show the histogram of the value obtained; the line shows the distribution according to the Central Limit Theorem.

Although the Gaussian is defined for all x, its weight is so strongly concentrated within a finite, and actually quite small, interval (about [–5, 5]) that values outside this range will not occur. It is not just that only one in a million events will deviate from the mean by more than 5 standard deviations: the decline continues, so that fewer than one in 10²² events will deviate by more than 10 standard deviations. Large outliers are not just rare—they don’t happen!

This is both the strength and the limitation of the Gaussian model: if the Gaussian model applies, then we know that all variation in the data will be relatively small and therefore “benign.” At the same time, we know that for some systems, large outliers do occur in practice. This means that, for such systems, the Gaussian model and theories based on it will not apply, resulting in bad guidance or outright wrong results. (We will return to this problem shortly.)

It is the combination of two properties that makes the Gaussian probability distribution so common and useful: because of the Central Limit Theorem, the Gaussian distribution will occur whenever we we dealing with averages; and because so much of the Gaussian’s weight is concentrated in the central region, almost any expression can be approximated by concentrating only on the central region, while largely disregarding the tails.

As we will discuss in Chapter 10 in more detail, the first of these two arguments has been put to good use by the creators of classical statistics: although we may not know anything about the distribution of the actual data points, the Central Limit Theorem enables us to make statements about their averages. Hence, if we concentrate on estimating the sample average of any quantity, then we are on much firmer ground, theoretically. And it is impressive to see how classical statistics is able to make rigorous statements about the extent of confidence intervals for parameter estimates while using almost no information beyond the data points themselves! I’d like to emphasize these two points again: through clever application of the Central Limit Theorem, classical statistics is able to give rigorous (not just intuitive) bounds on estimates—and it can do so without requiring detailed knowledge of (or making additional assumptions about) the system under investigation. This is a remarkable achievement!

The price we pay for this rigor is that we lose much of the richness of the original data set: the distribution of points has been boiled down to a single number—the average.

The second argument is not so relevant from a conceptual point, but it is, of course, of primary practical importance: we can actually do many integrals involving Gaussians, either exactly or in very good approximation. In fact, the Gaussian is so convenient in this regard that it is often the first choice when an integration kernel is needed (we have already seen examples of this in Chapter 2, in the context of kernel density estimates, and in Chapter 4, when we discussed the smoothing of a time series).

The basic idea goes like this: we want to evaluate an integral of the form:

We know that the Gaussian is peaked around x = 0, so that only nearby points will contribute significantly to the value of the integral. We can therefore expand f(x) in a power series for small x. Even if this expansion is no good for large x, the result will not be affected significantly because those points are suppressed by the Gaussian. We end up with a series of integrals of the form

which can be performed exactly. (Here, a_n is the expansion coefficient from the expansion of f(x).)

We can push this idea even further. Assume that the kernel is not exactly Gaussian but is still strongly peaked:

where the function g(x) has a minimum at some location (otherwise, the kernel would not have a peak at all). We can now expand g(x) into a Taylor series around its minimum (let’s assume it is at x = 0), retaining only the first two terms: g(x) ≈ g(0) + g″ (0)x²/2 + ···. The linear term vanishes because the first derivative g′ must be zero at a minimum. Keeping in mind that the first term in this expansion is a constant not depending on x, we have transformed the original integral to one of Gaussian type:

which we already know how to solve.

This technique goes by the name of Laplace’s method (not to be confused with “Gaussian integration,” which is something else entirely).

Given that the Central Limit Theorem is a rigorously proven theorem, what could possibly go wrong? After all, the Gaussian distribution guarantees the absence of outliers, doesn’t it? Yet we all know that unexpected events do occur.

There are two things that can go wrong with the discussion so far:

So what should you do when you encounter a situation described by a power-law distribution? The most important thing is to stop using classical methods. In particular, the mean-field approach (replacing the distribution by its mean) is no longer applicable and will give misleading or incorrect results.

From a practical point of view, you can try segmenting the data (and, by implication, the system) into different groups: the majority of data points at small values (on the lefthand side in Figure 9-5), the set of data points in the tail of the distribution (for relatively large values), and possibly even a group of data points making up the intermediate regime. Each such group is now more homogeneous, so that standard methods may apply. You will need insight into the business domain of the data, and you should exercise discretion when determining where to make those cuts, because the data itself will not yield a natural “scale” or other quantity that could be used for this purpose.

There is one more practical point that you should be aware of when working with power-law distributions: the form ~ 1/x^β is only valid “asymptotically” for large values of x. For small x, this rule must be supplemented, since it obviously cannot hold for x → 0 (we can’t divide by zero). There are several ways to augment the original form near x = 0. We can either impose a minimum value x_min of x and consider the distribution only for values larger than this. That is often a reasonable approach because such a minimum value may exist naturally. For example there is an obvious “minimum” number of pages (i.e., one page) that a website visitor can view and still be considered a “visitor.” Similar considerations hold for the population of a city and the copies of books sold—all are limited on the left by x_min = 1. Alternatively, the behavior of the observed distribution may be different for small values. Look again at Figure 9-5: for values less than about 5, the curve deviates from the power-law behavior that we find elsewhere.

Depending on the shape that we require near zero, we can modify the original rule in different ways. Two examples stand out: if we want a flat peak for x = 0, then we can try a form like ~ 1/(a + x^β) for some a > 0, and if we require a peak at a nonzero location, we can use a distribution like ~ exp(–C/x)/x^β (see Figure 9-6). For specific values of β, two distributions of this kind have special names:

The expectation value E(f) of a function f(x), which in turn depends on some random quantity x, is nothing but the weighted average of that function in which we use the probability density p(x) of x as the weight function:

Of particular importance are the expectation values for simple powers of the variable x, the so called moments of the distribution:

Figure 9-6. The Lévy distribution for several values of the parameter c.

The first expression must always equal 1, because we expect p(x) to be properly normalized. The second is the familiar mean, as the weighted average of x. The last expression is used in the definition of the standard deviation:

For power-law distributions, which behave as ~ 1/x^β with β > 1 for large x, some of these integrals may not converge—in this case, the corresponding moment “does not exist.” Consider the kth moment (C is the normalization constant C = E(1) = ∫ p(x) dx):

Unless β – k > 1, this integral does not converge at the upper limit of integration. (I assume that the integral is proper at the lower limit of integration, through a lower cutoff x_min or another one of the methods discussed previously.) In particular, if β < 2, then the mean and all higher moments do not exist; if β < 3, then the standard deviation does not exist.

We need to understand that this is an analytical result—it tells us that the distribution is ill behaved and that, for instance, the Central Limit Theorem does not apply in this case. Of course, for any finite sample of n data points drawn from such a distribution, the mean (or other moment) will be perfectly finite. But these analytical results warn us that, if we continue to draw additional data points from the distribution, then their average (or other moment) will not settle down: it will grow as the number of data points in the sample grows. Any summary statistic calculated from a finite sample of points will therefore not be a good estimator for the true (in this case: infinite) value of that statistic. This poses an obvious problem because, of course, all practical samples contain only a finite number of points.

Power-law distributions have no parameters that could (or need) be estimated—except for the exponent, which we know how to obtain from a double logarithmic plot. There is also a maximum likelihood estimator for the exponent:

where x₀ is the smallest value of x for which the asymptotic power-law behavior holds.

If you want to dig deeper into the theory of heavy-tail phenomena, you will find that it is a mess. There are two reasons for that: on the one hand, the material is technically hard (since one must make do without two standard tools: expectation values and the Central Limit Theorem), so few simple, substantial, powerful results have been obtained—a fact that is often covered up by excessive formalism. On the other hand, the “colorful” and multi disciplinary context in which power-law distributions are found has led to much confusion. Similar results are being discovered and re-discovered in various fields, with each field imposing its own terminology and methodology, thereby obscuring the mathematical commonalities.

The unexpected and often almost paradoxical consequences of power-law behavior also seem to demand an explanation for why such distributions occur in practice and whether they might all be expressions of some common mechanisms. Quite a few theories have been proposed toward this end, but none has found widespread acceptance or proved particularly useful in predicting new phenomena—occasionally grandiose claims to the contrary notwithstanding.

At this point, I think it is fair to say that we don’t understand heavy-tail phenomena: not when and why they occur, nor how to handle them if they do.

The geometric distribution (see Figure 9-7):

p(k, p) = p(1 – p)^k–1 with k = 1, 2, 3,...

is a special case of the binomial distribution. It can be viewed as the probability of obtaining the first Success at the kth trial (i.e., after observing k – 1 failures). Note that there is only a single arrangement of events for this outcome, hence the combinatorial factor is equal to one. The geometric distribution has mean μ = 1/p and standard deviation .

The binomial distribution gives us the probability of observing exactly k events in n distinct trials. In contrast, the Poisson distribution describes the probability of finding k events during some continuous observation interval of known length. Rather than being characterized by a probability parameter and a number of trials (as for the binomial distribution), the Poisson distribution is characterized by a rate λ and an interval length t.

The Poisson distribution p(k, t, λ) gives the probability of observing exactly k events during an interval of length t when the rate at which events occur is λ (see Figure 9-8):

Figure 9-8. The Poisson distribution: .

Because t and λ only occur together, this expression is often written in a two-parameter form as p(k, υ) = e^–υ υ^k/k!. Also note that the term e^–λt does not depend on k at all—it is merely there as a normalization factor. All the action is in the fractional part of the equation.

Let’s look at an example. Assume that phone calls arrive at a call center at a rate of 15 calls per hour (so that λ = 0.25 calls/minute). Then the Poisson distribution p(k, 1, 0.25) will give us the probability that k = 0, 1, 2,... calls will arrive in any given minute. But we can also use it to calculate the probability that k calls will arrive during any 5-minute time period: p(k, 5, 0.25). Note that in this context, it makes no sense to speak of independent trials: time passes continuously, and the expected number of events depends on the length of the observation interval.

We can collect a few results. Mean μ and standard deviation σ for the Poisson distribution are given by:

Notice that only a single parameter (λt) controls both the location and the width of the distribution. For large λ, the Poisson distribution approaches a Gaussian distribution with μ = λ and . Only for small values of λ (say, λ < 20) are the differences notable. Conversely, to estimate the parameter λ from observations, we divide the number k of events observed by the length t of the observation period: λ = k/t. Keep in mind that when evaluating the formula for the Poisson distribution, the rate λ and the length t of the interval of interest must be of compatible units. To find the probability of k calls over 6 minutes in our call center example above, we can either use t = 6 minutes and λ = 0.25 calls per minute or t = 0.1 hours and λ = 15 calls per hour, but we cannot mix them. (Also note that 6 · 0.25 = 0.1 · 15 = 1.5, as it should.)

The Poisson distribution is appropriate for processes in which discrete events occur independently and at a constant rate: calls to a call center, misprints in a manuscript, traffic accidents, and so on. However, you have to be careful: it applies only if you can identify a rate at which events occur and if you are interested specifically in the number of events that occur during intervals of varying length. (You cannot expect every histogram to follow a Poisson distribution just because “we are counting events.”)

Some quantities are inherently asymmetrical. Consider, for example, the time it takes people to complete a certain task: because everyone is different, we expect a distribution of values. However, all values are necessarily positive (since times cannot be negative). Moreover, we can expect a particular shape of the distribution: there will be some minimum time that nobody can beat, then a small group of very fast champions, a peak at the most typical completion time, and finally a long tail of stragglers. Clearly, such a distribution will not be well described by a Gaussian, which is defined for both positive and negative values of x, is symmetric, and has short tails!

The log-normal distribution is an example of an asymmetric distribution that is suitable for such cases. It is related to the Gaussian: a quantity follows the log-normal distribution if its logarithm is distributed according to a Gaussian.

The probability density for the log-normal distribution looks like this:

(The additional factor of x in the denominator stems from the Jacobian in the change of variables from x to log x.) You may often find the log-normal distribution written slightly differently:

This is the same once you realize that log(x/μ) = log(x) – log(μ) and make the identification . The first form is much better because it expresses clearly that μ is the typical scale of the problem. It also ensures that the argument of the logarithm is dimensionless (as it must be).

Figure 9-9. The log-normal distribution.

Figure 9-9 shows the log-normal distribution for a few different values of σ. The parameter σ controls the overall “shape” of the curve, whereas the parameter μ controls its “scale.” In general, it can be difficult to predict what the curve will look like for different values of the parameters, but here are some results (the mode is the position of the peak).

Values for the parameters can be estimated from a data set as follows:

The log-normal distribution is important as an example of a standard statistical distribution that provides an alternative to the Gaussian model for situations that require an asymmetrical distribution. That being said, the log-normal distribution can be fickle to use in practice. Not all asymmetric point distributions are described well by a log-normal distribution, and you may not be able to obtain a good fit for your data using a log-normal distribution. For truly heavy-tail phenomena in particular, you will need a power-law distribution after all. Also keep in mind that the log-normal distribution approaches the Gaussian as σ becomes small compared to μ (i.e., σ/μ ≪ 1), at which point it becomes easier to work with the familiar Gaussian directly.

Many additional distributions have been defined and studied. Some, such as the gamma distribution, are mostly of theoretical importance, whereas others—such as the chi-square, t, and F distributions—are are at the core of classical, frequentist statistics (we will encounter them again in Chapter 10). Still others have been developed to model specific scenarios occurring in practical applications—especially in reliability engineering, where the objective is to make predictions about likely failure rates and survival times.

I just want to mention in passing a few terms that you may encounter. The Weibull distribution is used to express the probability that a device will fail after a certain time. Like the log-normal distribution, it depends on both a shape and a scale parameter. Depending on the value of the shape parameter, the Weibull distribution can be used to model different failure modes. These include “infant mortality” scenarios, where devices are more likely to fail early but the failure rate declines over time as defective items disappear from the population, and “fatigue death” scenarios, where the failure rate rises over time as items age.

Yet another set of distributions goes by the name of extreme-value or Gumbel distributions. They can be used to obtain the probability that the smallest (or largest) value of some random quantity will be of a certain size. In other words, they answer the question: what is the probability that the largest element in a set of random numbers is precisely x?

Quite intentionally, I don’t give formulas for these distributions here. They are rather advanced and specialized tools, and if you want to use them, you will need to consult the appropriate references. However, the important point to take away here is that, for many typical scenarios involving random quantities, people have developed explicit models and studied their properties; hence a little research may well turn up a solution to whatever your current problem is.