Probability and Statistics

Why fitting a logistic is nearly impossible from early data

Posted on by John

Nothing grows exponentially forever. What appears to be an exponential curve often turns out to be some sort of S curve, such as a logistic curve.logistic curve with extrapolations

Suppose you’re collecting data on the left side of the curve. If there’s even a small amount of error in your data, you won’t be able to predict the asymptotic value with any accuracy. But if you have data on both sides of the inflection point, you can make a good prediction of the limiting value.

I’ve written about this before, explaining that the problem is hard, but I didn’t say why it’s hard. Here I’d like to give an idea why it’s hard.

Suppose you want to fit a logistic equation

y(t) = \frac{L}{1 + \exp(-k(t - t_0))}

to three distinct values of t and the corresponding values of y. There is a unique solution, but in general you cannot find a solution in closed form. However, if the values of t are evenly spaced

y_2 - y_1 = y_1 - y_0 = h

there is a method [1] to solve for the parameters L, k, and t0. For this post we’re only interested in the limiting value L, and it can be found by

L = \frac{y_1^2(y_0 + y_2) - 2y_0 y_1 y_2}{y_1^2 - y_0 y_2}

independent of h.

To find out how small changes in the y‘s change the estimate of L, we take the partial derivatives of L with respect to the y‘s and find

\frac{\partial L}{\partial y_0} = \frac{\partial L}{\partial y_2} = \frac{y_1^2\, h^2}{\left(y_1^2 - y_0 y_2\right)^2}

and

\frac{\partial L}{\partial y_1} = \frac{-2\, y_0 y_2\, h^2}{\left(y_1^2 - y_0 y_2\right)^2}
All three derivatives have the same expression in the denominator: y1² − y0 y2.

If the function y(t) were an exponential, this expression would be exactly zero [2]. The function y(t) is not exactly exponential, but it is approximately exponential when the t‘s are in the left or right tail of the logistic curve. The further out in either tail the t‘s are, the closer the expression is to zero.

So when all the t‘s come from the same side of the inflection point, y(t) is nearly exponential the partial derivatives are huge and so the fitted value of L is extremely sensitive to changes in the y‘s.

As a concrete example, set L = k = 1 and t0 = 0. Evaluate y(t) at −2, −1.5, and −1. Then the values of y are

y0 = 0.11920292
y1 = 0.18242552
y2 = 0.26894142

If you forecast L using exactly these three values you’ll get L = 1.

But if you change y0 to 0.12374097, the forecasted value of L is infinite. Values of y0 in the interval [0.11920292, 0.12374097] predict values of K in [1, ∞].

[1] Raymond Pearl and Lowell J. Reed. On the Rate of Growth of the Population of the United States Since 1790 and its Mathematical Representation. Proceedings of the National Academy of Sciences of the United States of America, Vol. 6, No. 6 (Jun. 15, 1920), pp. 275-288

[2] exp(x + h)² = exp(x)² exp(h)² = exp(x) exp(x + 2h)

 

Coming soon

Posted on by John

There’s a pizza shop near my home with a sign out front that says “Coming Soon.” When I drove by it this morning I thought about how you would model the time until an event happens that is “coming soon.”

Suppose I look at the sign one day and guess how many days until the pizza shop will open. When I drive by a week later and guess again, should my guess be smaller? You might argue that the shop will open some day, fixed in time but unknown to me, and so every day I’m one day closer to the eventual opening.

You might model the pizza shop opening like radioactive decay and say that the estimated number of days until it opens is always the same until the day it actually opens.

Now I think this shop has been “coming soon” for over a year. So instead of decreasing, every day I increase my estimate of the time until the shop opens. Something has gone wrong that the owners didn’t expect when they put up the sign.

Maybe the reasonable thing would be for estimated days until opening to decrease over time, but only up to a point. After some point, the longer a business has been “coming soon” the less like that it is coming soon, or coming at all.

This brings up an interesting point about modeling. There are two probability distributions at work: the probability that the shop will eventually open, and the time until opening assuming it eventually opens.

When the sign first goes up saying the business is coming soon, there’s some change that it is in fact not coming. Maybe you’re optimistic and think this probability is small, but it would seem unreasonable to think the probability is zero. That means the expected number of days until opening is always infinite. If there’s a probability ε that the shop never opens, the expected time to opening is

ε × ∞ + (1 − ε) × something = ∞.

Probability of correcting errors

Posted on by John

Error correcting codes are most simply described in terms of the errors they can certainly correct. For example, the Hadamard code used for the Mariner 9 probe to Mars encoded each 6-bit pixel to a 32-bit codeword in such a way that the original pixel could be recovered if no more than 7 bits were corrupted in transit.

What is the probability that a pixel could be repaired if corrupted? That depends on your probability model. We will assume that the probability of each bit being flipped is p and that errors are independent.

(Are errors independent, i.e. if a bit flips, is the next bit more or less likely to flip? That would depend on context.)

It’s straight-forward to calculate the probability that 7 or fewer bits out of 32 flip; this is the cumulative distribution of a binomial random variable. The following Python code will return the probability of k or fewer successes out of n trials, each with probability of success p:

    from scipy.stats import binom
    print(binom.cdf(k, n, p))

For example, if there is a 10% chance that each bit will flip, there’s a 98.8% chance that 7 or fewer bits out of 32 will flip.

However this only gives a lower bound on the probability of correcting an error. If eight bits flip in transit, we cannot tell with certainty which codeword was sent, but there will be a couple possibilities that stand out. We’ll have to guess, but we’ve narrowed down the possibilities. With even more flipped bits, there’s always a chance of recovering the original data. Still, the lower bound captures most of the probability of recovery.

Now suppose you’re given a desired error recovery rate and have to determine what value of p it can sustain. For example, someone might say they want a 98.8% chance of recovering a pixel correctly, and you could come back and say p must be less than or equal to 0.1. This would be a conservative answer because as discussed above, p = 0.1 gives a pixel recovery probability of something more than 98.8, though it’s messy to calculate how much more.

You could solve for p by trial and error, or you could use some more sophisticated math to compute p directly. Given a probability F, you can solve for p such that the probability of up to k successes out of n trials using the inverse of the regularized incomplete beta function.

from scipy.special import betaincinv
p = 1 - betaincinv(n - k, k + 1, F)

Calculating F given n, k, and p could be a homework exercise in an introductory probability course. Solving for p given F, n, and k either requires some numerical programming or special functions and so would be a more challenging problem.

Related posts

Posterior variance

Posted on by John

A few days ago I wrote a post entitled Does additional data always reduce posterior variance?. In a nutshell, the answer is no, not always.

That led the previous post which looked at posterior means for three Bayesian models, showing how the posterior mean is a weighted average of the prior mean and the mean of the new data. The weights are precisions, which means something different for each model.

For the beta-binomial model, variance may increase when seeing unexpected data (details here), but precision always increases.

For the normal-normal model precision is the reciprocal of variance. Every new data point makes precision go up and posterior variance go down.

The Poisson-gamma model may be the most interesting. As stated in the previous post, if data has a Poisson distribution with parameter λ, and λ has a gamma(α0, β0) prior distribution, then the posterior distribution on λ after observing k events over time t has a gamma(α0 + k, β0 + t) posterior distribution. Therefore the posterior variance is

(α0 + k) / (β0 + t)².

Note the posterior variance is an increasing function of k and a decreasing function of t. This means that the posterior variance increases every time an event is observed, and it decreases quadratically between observations.

Here’s an illustration. I simulated data from a Poisson process with λ and used a gamma(1, 1) prior on λ. Here’s a plot of the posterior variance.

Online (one-pass) algorithms

Posted on by John

Canonical example

The sample variance of a set of numbers is defined in terms of the sum of the squared distances from each point to the mean.

s^2 = \frac{1}{n-1}\sum_{i=1}^n (x_i -\bar{x})^2

So it would seem that you first need to calculate the mean, then go back and compute the squared differences from the mean. And yet sample variance can be computed in one pass through the data.

You’ll find two equivalent equations in statistics books: the one described above and another based on the sum of the data points and the sum of the data points squared.

s^2 = \frac{1}{n(n-1)}\left(n\sum_{i=1}^n x_i^2 -\left(\sum_{i=1}^n x_i\right)^2\right)

While this equation is theoretically correct, it is numerically unstable. Code that directly implements this equation can return a negative value for a quantity that is theoretically positive. I’ve seen this happen with real data, causing a program to crash when taking the square root of the variance to get the standard deviation.

However, there is an algorithm that computes mean and variance in one pass that is accurate and numerically stable. This algorithm was developed by B. P. Welford in 1962. I discuss Welford’s algorithm and give code for implementing it here.

Online algorithms

Welford’s algorithm is known in computer science as an “online” algorithm. This term was coined well before the Internet. For example, see the paper [1] from 1965.

But of course now “online” means something else, and so the technical and colloquial uses of “online algorithm” have split. Technical literature uses the phrase to describe the kinds of algorithms in this post. Most people would take “online algorithm” to mean code that runs on a remote server. You may see “streaming algorithm” as a contemporary technical term, but I’d still search on “online algorithm” to find papers.

Computing higher moments online

Welford’s algorithm computes the first two moments, mean and variance, of a data set online. It is also possible to compute skewness and kurtosis online, as well as higher moments.

Online regression

Simple linear regression is closely related to calculating mean and variance, and it is possible to compute simple regression coefficients online. I have some old notes on this here.

This post was motivated by an email asking me about multiple regression. It is also possible to compute multiple regression coefficients online, but I haven’t done this. I found a couple references, [2] and [3], but I have not read them. There is a simple procedure for two predictor variables but I believe things get a little more complicated with three or more predictors, requiring a recursive least squares algorithm.

Related posts

The notion of online algorithms is closely related to the notion of a fold in functional programming. Here are several posts on computing things with folds.

[1] One-Tape, Off-Line Turing Machine Computations by F. C. Hennie. Information and Control. 8, 553-578 (1965). Available here. In this paper Hennie writes “In an on-line computation the input data are supplied to the machine, one symbol at a time, at a special input terminal. … In an off-line computation all of the input symbols are written on one of the machine’s tapes prior to the start of the computation.

[2] Arthur Albert and Robert W. Sittler, “A Method for Computing Least Squares Estimators that Keep Up with the Data,” Journal of the Society for Industrial and Applied Mathematics, Series A: Control, 3(3), 384–417, 1965. DOI: 10.1137/0303026.

[3] Petre Stoica and Per Ashgren. Exact initialization of the recursive least-squares algorithm. Int. J. Adapt. Control Signal Process. 2002; 16:219&ndashh;230.

Mills ratio and tail thickness

Posted on by John

The Mills ratio [1] is the ratio of the CCDF to the PDF. That is, for a random variable X, the Mills ratio at x is the complementary cumulative distribution function divided by the density function. If the density function of X is f, then

m(x) = \frac{\int_x^\infty f(x)\, dx}{f(x)}

The Mills ratio highlights an important difference between the Student t distribution and the normal distribution.

Introductory statistics classes will say things like “you can approximate a t distribution with a normal if it has more than 30 degrees of freedom.” That may be true, depending on the application. A t(30) distribution and a normal distribution behave similarly in the middle but not in the tails.

Mills ratio plot for t(30) vs normal

The Mills ratio for a t distribution with ν degrees of freedom is asymptotically x/ν,  while the Mills ratio for a standard normal distribution is asymptotically 1/x. Note that increasing ν does make the Mills function smaller, but it still eventually grows linearly whereas the Mills function of a normal distribution decays linearly.

In general, the Mills ratio is a decreasing function for thin-tailed distributions and an increasing function for fat-tailed distributions. The exponential distribution is in the middle, with constant Mills function.

Related posts

[1] Named after John P. Mills, so there’s no apostrophe before the s.

Sigmas and Student

Posted on by John

I saw something yesterday saying that the Japanese bond market had experienced a six standard deviation move. This brought to mind a post I’d written eight years ago.

All probability statements depend on a model. And if you’re probability model says an event had a probability six standard deviations from the mean, it’s more likely that your model is wrong than that you’ve actually seen something that rare. I expand on this idea here.

How likely is it that a sample from a random variable will be six standard deviations from its mean? If you have in mind a normal (Gaussian) distribution, as most people do, then the probability is on the order of 1 chance in 10,000,000. Six sigma events are not common for any distribution, but they’re not unheard of for distributions with heavy tails.

Let X be a random variable with a Student t distribution and ν degrees of freedom. When ν is small, i.e. no more than 2, the tails of X are so fat that the standard deviation doesn’t exist. As ν → ∞ the Student t distribution approaches the normal distribution. So in some sense this distribution interpolates between fat tails and thin tails.

What is the probability that X takes on a value more than six standard deviations from its mean at 0, i.e. what does the function

f(ν) = Prob(X > 6σ)

look like as a function of ν where σ² = ν/(ν − 2) is the variance of X?

As you’d expect, the limit of f(ν) as ν → ∞ is the probability of a six-sigma event for a normal distribution, around 10−7 as mentioned above. Here’s a plot of f(ν) for ν > 3. Notice that the vertical axis is on a log scale, i.e. the probability decreases exponentially.

What you might not expect is that f(ν) isn’t monotone. It rises to a maximum value before it decays exponentially. In hindsight this makes sense. As ν → 2+ the variance becomes infinite, and the probability of being infinitely far from the mean is 0. Here’s a plot of f(ν) between 2 and 3.

So six sigma probabilities for a Student t distribution rise from 0 up to a maximum of around 10−3 then decrease exponentially, then asymptotically approach a value around 10−7.

Related posts

The middle binomial coefficient

Posted on by John

The previous post contained an interesting observation:

\sqrt{52!} \approx 26! \, 2^{26}

Is it true more generally that

\sqrt{(2n)!} \approx n! \, 2^n

for large n? Sorta, but the approximation gets better if we add a correction factor.

If we square both sides of the approximation and move the factorials to one side, the question becomes whether

\frac{(2n)!}{(n!)^2} = \binom{2n}{n} \approx 4^n

Now the task becomes to estimate the middle coefficient in when we apply the binomial theorem to (x + y)2n.

A better approximation for the middle binomial coefficient is

\binom{2n}{n} \approx \frac{4^n}{\sqrt{\pi n}}

Now the right hand side is the first term of an asymptotic series for the left. The ratio of the two sides goes to 1 as n → ∞.

We could prove the asymptotic result using Stirling’s approximation, but it’s more fun to use a probability argument.

Let X be a binomial random variable with distribution B(2n, 1/2). As n grows, X converges in distribution to a normal random variable with the same mean and variance, i.e. with μ = n and σ² = n/2. This says for large n,

\text{Prob}(X = 1/2) = \binom{2n}{n} 4^{-n} \approx \frac{1}{\sqrt{2\pi\sigma^2}} = \frac{1}{\sqrt{\pi n}}

The argument above only gives the first term in the asymptotic series for the middle coefficient. If you want more terms in the series, you’ll need to use more terms in Stirling’s series. If we add a couple more terms we get

\binom{2n}{n} = \frac{4^n}{\sqrt{\pi n}} \left(1 - \frac{1}{8n} + \frac{1}{128n^2} + {\cal O}\left(\frac{1}{n^3}\right) \right)

Let’s see how much accuracy we get in estimating 52 choose 26.

from scipy.special import binom
from numpy import pi, sqrt

n = 26
exact = binom(2*n, n)
approx1 = 4**n/sqrt(pi*n)
approx2 = approx1*(1 - 1/(8*n))
approx3 = approx1*(1 - 1/(8*n) + 1/(128*n**2))

for a in [approx1, approx2, approx3]:
    print(exact/a)

This prints

0.9952041409266293
1.0000118903997048
1.0000002776131290

and so we see substantial improvement from each additional term. This isn’t always the case with asymptotic series. We’re guaranteed that for a fixed number of terms, the relative error goes to zero as n increases. For a fixed n, we do not necessarily get more accuracy by including more terms.

Related posts

Perfect and imperfect shuffles

Posted on by John

Take a deck of cards and cut it in half, placing the top half of the deck in one hand and the bottom half in the other. Now bend the stack of cards in each hand and let cards alternately fall from each hand. This is called a rifle shuffle.

Random shuffles

Persi Diaconis proved that it takes seven shuffles to fully randomize a desk of 52 cards. He studied videos of people shuffling cards in order to construct a realistic model of the shuffling process.

Shuffling randomizes a deck of cards due to imperfections in the process. You may not cut the deck exactly in half, and you don’t exactly interleave the two halves of the deck. Maybe one card falls from your left hand, then two from your right, etc.

Diaconis modeled the process with a probability distribution on how many cards are likely to fall each time. And because his model was realistic, after seven shuffles a deck really is well randomized.

Perfect shuffles

Now suppose we take the imperfection out of shuffling. We do cut the deck of cards exactly in half each time, and we let exactly one card fall from each half each time. And to be specific, let’s say the first card will always fall from the top half of the deck. That is, we do an in-shuffle. (See the next post for a discussion of in-shuffles and out-shuffles.) A perfect shuffle does not randomize a deck because it’s a deterministic permutation.

To illustrate a perfect in-shuffle, suppose you start with a deck of these six cards.

A23456

Then you divide the deck into two halves.

A23 456

Then after the shuffle you have the following.

4A5263

Incidentally, I created the images above using a font that included glyphs for the Unicode characters for playing cards. More on that here. The font produced black-and-white images, so I edited the output in GIMP to turn things red that should be red.

Coming full circle

If you do enough perfect shuffles, the deck returns to its original order. This could be the basis for a magic trick, if the magician has the skill to repeatedly perform a perfect shuffle.

Performing k perfect in-shuffles will restore the order of a deck of n cards if

2k = 1 (mod n + 1).

So, for example, after 52 in-shuffles, a deck of 52 cards returns to its original order. We can see this from a quick calculation at the Python REPL:

>>> 2**52 % 53
1

With slightly more work we can show that less than 52 shuffles won’t do.

>>> for k in range(1, 53):
    ... if 2**k % 53 == 1: print(k)
52

The minimum number of shuffles is not always the same as the size of the deck. For example, it takes 4 shuffles to restore the order of a desk of 14 cards.

>>> 2**4 % 15
1

Shuffle code

Here’s a function to perform a perfect in-shuffle.

def shuffle(deck):
    n = len(deck)
    return [item for pair in zip(deck[n//2 :], deck[:n//2]) for item in pair]

With this you can confirm the results above. For example,

n = 14
k = 4
deck = list(range(n))
for _ in range(k):
    deck = shuffle(deck)
print(deck)

This prints 0, 1, 2, …, 13 as expected.

Related posts

Trying to fit exponential data

Posted on by John

The first difficulty in trying to fit an exponential distribution to data is that the data may not follow an exponential distribution. Nothing grows exponentially forever. Eventually growth slows down. The simplest way growth can slow down is to follow a logistic curve, but fitting a logistic curve has its own problems, as detailed in the previous post.

Suppose you are convinced that whatever you’re wanting to model follows an exponential curve, at least over the time scale that you’re interested in. This is easier to fit than a logistic curve. If you take the logarithm of the data, you now have a linear regression problem. Linear regression is numerically well-behaved and has been thoroughly explored.

There is a catch, however. When you extrapolate a linear regression, your uncertainly region flares out as you go from observed data to linear predictions based on the observed data. Your uncertainty grows linearly. But remember that we’re not working with the data per se; we’re working with the logarithm of the data. So on the original scale, the uncertainty flares out exponentially.