blog - steve t. piantadosi

You are an octopus who loves sushi. Every friday night, you visit the local sushi restaurant and chow down. Your favorite sushi is the Spider Roll and least favorite is Tako. But, you are the kind of octopus who really likes variety—a totally random mix of sushi.

Unfortunately, you go to a sushi restaurant with a conveyor belt of sushi, where each roll moves past you in succession. As each roll passes, you can decide whether to pick it up (and put down one if your arm is currently holding one) or ignore it and wait for the next. Your problem is to find some way of choosing a sample of sushi rolls uniformly, such that each sushi roll that passes you has an equal probability of ending up being in one of your arms at the end. Sadly, you don’t know how long the sushi belt will continue to run. What to do?

One thing you might think to do is as each sushi roll passes you, you flip a coin. If the coin is heads, you pick up the current roll with an arm, and throw out that arm’s sushi. If the coin is tails, you ignore the current roll. Done dumbly, this, of course, biases you against keeping the earlier sushi, since they have a greater chance of being replaced if they’re even picked up. So you figure you can change the coin weight as the sushi roll past, and maybe, with some luck, you can end up with a uniform sample. But, if you don’t know how many sushi rolls will go by, you might have trouble deciding how quickly the coin weight should change…

But there’s a better way. It’s called reservoir sampling. Here’s one way to think about it, before you think about it: suppose that we assigned a random number chosen iid, uniformly from [0,1], to each sushi roll on the conveyor belt—all the ones you have seen and will see.

Now suppose that to take a uniform random sample of 8 of them (one for each arm), you take the rolls with the 8 lowest random numbers. Each roll has an equal probability of having its number in the lowest 8, so your sample of the lowest 8 is uniform.

What have you gotten? Well, something very nice. By thinking of it in those terms, rather than coin-flipping-accept-reject terms, you can make a pretty online algorithm. Each time you see a new sushi roll, you draw a random number. You hold on to the 8 sushi rolls with the lowest numbers, and if the new roll is lower than the largest of your 8, you add it to your 8, and drop the largest. At the end of this process, you are holding in your 8 arms the 8 lowest samples from all of them. And you didn’t need to ever hold on to more than 8 rolls. Cool.

(If you aren’t really an octopus, maybe you’re doing NLP or psycholinguistics and trying to draw N random sentences containing the word “octopus” from a corpus, without holding all occurrences of “octopus” in memory.)

Of course, that algorithm I gave isn’t great. In particular, you have to look at each sushi roll. But of course you could draw random numbers for the next several rolls and not even look at them if they are too large. Amazingly, you can sample 8 rolls uniformly, from N total in less than time. Thank this guy. In fact, drawing n samples uniformly from N is doable in time.

Pretty, clever algorithm.

Here’s a thought which recently came up: optimization problems with many possible solutions are more likely to have a single best solution than optimization problems with few solutions. A problem with 10 possible solutions is more likely to have multiple “best” solutions than one with 1000 possible solutions. Why might you think this?

Well, you might think about it in the following way: suppose that value or “goodness” or utility of each possible solution to an optimization problem is sampled from some underlying (discrete) distribution. So the optimal solution is the max of this distribution of samples. This assumption probably isn’t true of many real kinds of optimization problems, but let’s go with it for a minute.

If this assumption were true, then as your sample size increases, the best solution takes on a more and more extreme value in the underlying distribution. Analogously, the tallest person in a room with 5 people will be shorter than the tallest person in a room with 25 people. But if the max gets more and more extreme, its value is lower and lower probability (for any distribution, in the limit), and so you are less and less likely to have sampled multiple solutions with that maximum value. So the max should tend to be unique as the problem gets more and more complex. Analogously, if you grab 10 people and measure their height, you are more likely to find two people tied for tallest (rounded to, say, the nearest inch) than if you grabbed 1000 people.

Wait, that’s not obviously true. Actually, it’s not even right.

Here’s a simulation where you draw N samples from a poison distribution (lambda=0.1). The X axis is log N. The red line shows the probability of finding a tie for the maximum value. (Code here.)

What’s surprising, maybe, is that the probability of a tie for first place oscillates (it’s so surprising I thought that R had a bug in it’s random number generation until I wrote it in another language). In fact, for many distributions the limiting probability probability of a tie for the maximum value doesn’t exist. This was originally investigated for geometric distributions; for these distributions, you can picture the game is that N people flip a coin. If a person gets a tail, they are eliminated. The winner is the last person to be eliminated, but there is some probability that on the next flip, everyone left will be eliminated and thus tie as winners. What does the probability that there is a tie for winner, as a function of N? Well, it oscillates.

In 1995 Baryshnikov, Eisenberg, and Stengle proved the general theorem you really want: the limiting probability of a tie for the max exists if and only if

where P is the underlying discrete distribution on . If the limit doesn’t exist, you get something that looks like the plot above. So, for many distributions, the probability of a tie for the maximum oscillates as the sample sizes increases. Bizarre, I think.

But here is one way to understand it. Imagine having dawn a sample, finding the max, and slowing adding more samples. As you add more and more samples, you will get the max again before sampling something larger, assuming that the max is sufficiently higher probability than the numbers that follow it (which is one way to read what happens when the above limit is not 0). Then with the new max, you start all over, getting the new max again before finding a new new max. So your probability of getting the max more than once oscillates with the sample size.

This might all be irrelevant to actual optimization problems since I’m not sure which ones have spaces of possible solutions which can be thought of as draws from some underlying distribution. But I thought it was odd anyway that the probability of a tie for the max of N samples oscillates as N increases. So yeah.

Here’s the most amazing thing I read recently (in Mar-Apr edition of American Scientist–probably the best science writing around–which attributed this puzzle to David Blackwell).

Suppose that you write down two numbers, each on a slip of paper. You shuffle them up and put them, face-down, on a table. I am free to flip over the first one. Now I say whether or not the first number is larger than the second number. If I am right, you give me $100. If I am wrong, I give you $100.

First question: is this a fair bet for you? You should try for a minute to prove that this is a fair bet (and may “succeed”). Informally, you might think the first number doesn’t give me any information about the second since you could write down anything, so it seems like there is no strategy I can use to, on average, do better than breaking even.

Ahh our dumb intuition. Here’s how I can make money on this bet. This is the remarkable part. I flip over the first slip of paper. I then choose a number at random (say, from a normal distribution). If the random number is larger than the one I flipped over, I tell you that the hidden number is also larger than the one I flipped over. If the random number is smaller than the one I flipped over, I tell you that the hidden number is also smaller than the one I flipped over.

How could this possibly make me money on an even bet? Well suppose that you wrote down numbers X and Y. When my random number is less than both X and Y, then I will win 50% of the time since I will turn over the smaller number first half the time, and the larger number first half the time. In both of these cases, I will tell you the hidden number is the smaller of the two numbers, and so will be right half the time. Similarly, when the random number is larger than both X and Y, I will be right half the time.

But, when the random number is between X and Y, I will will all the time. To see this, suppose that X < Y. When I turn over X, my random number is larger than it, so I say that the unseen number, Y, is larger than the seen number X, and I am right. When I turn over Y first, my random number is smaller than it, so I say that the unseen number, X, is smaller than the seen number Y, and I am right.

So when I randomly choose a number between X and Y, I win all the time. Thus, if I sample from a distribution which assigns nonzero probability to every range, I will on average win more than half the time. Crazy, huh?

I played around with a fun problem this weekend with Celeste. We started to talk about descendants and ancestors. What are the odds that someone descended from me will be alive T years in the future? What are the odds I will have no descendants in 10000 years? What is the relationship between the number of kids I have and the odds my descendants will all die out?

There are some funny phenomena relevant to this. Mitochondrial Eve is one. The number of people descended from Genghis Kahn is another. In Knuth’s Volume 1, he talks about a “proof” by H. W. Watson that under certain assumptions, infinitely many people will be born in the future, but each family line will die with probability 1 (pg 383). This is of course logically inconsistent, and Knuth proves so: that in an infinitely tall tree, you can find an infinite path of descent. If people live forever, someone’s family will live forever. Which should be obvious.

Here’s a simple model for thinking about infinite family trees. Suppose that there is some bounded population P. To create people for the next generation, the following is repeated: two people are chosen at random, bred, and their (one) child inhabits the next generation. That child is a descendant of their two parents. This is repeated until the next generation is filled up with people. And so on. This neglects a lot of important things–like gender, natural selection, sexual selection, and the grossness of incest. But maybe its not such a bad start…

First, two observations. If the population size is bounded, then it is eventually going to be the case that everyone is a descendant of me, or nobody is. This is because everyone being my descendant, and nobody being my descendant are fixed points—once either is true, it is true for all of time. And each generation there is some probability of either of these happening, so given enough time, one must come true. The second thing to notice is that the expected number of descendants of me at the next time step is 1-(1-p)^2 = p(2-p), where p is the proportion of people who are my descendants at the current time step. This is because the probability of getting someone in the next generation who is not a descendant of me is the probability that neither parent was a descendant of me, or (1-p)(1-p). (This also shows that p=0 and p=1 are the only fixed points). Thus, at each generation the number of my descendants should grow by a factor of (2-p). Things are looking good.

But of course there is a lot of variability, especially when you only have a few kids initially. So how many kids should someone have in order to be pretty sure your genes will take over? Since I know more about perl than stochastic processes, I wrote a little script to figure this out. If anyone wants to solve this analytically, I’d like to know, because I can’t run the perl script very quickly on a population of size 6 billion.

Here are some results, showing the probability that everyone is eventually a descendant of me for various population sizes and number of kids I have, averaged over 100 runs:

So, you don’t need many kids to eventually take have a good odds of taking over. At least when the population is relatively small. But in the larger population, probably most people’s kids reproduce themselves more readily than in this model, so you don’t need many kids to eventually take over. I wonder if anyone has attempted to model the distribution of family names along these lines…

Here’s a problem I recently had: suppose you are given a sample of some joint distribution on X and Y. That is you are given a bunch of (X,Y) pairs. On thing you may want to know is whether the distribution P( Y | X ) depends on X. You can’t always test this by doing a regression–sometimes the expected of Y may not depend on X, but the distribution of Ys might. Sometimes P(Y | X) may be entirely different types of distributions for various X.

But here’s one fun technique: it is easy to compute the marginal densities P(Y) and P(X) using, say, density estimation. Then we can do a 2D density estimation where we weight each data point (x,y) by 1 / ( P(x) P(y) ).

Weighting the density estimation like this doesn’t recover a function for the actual density. Instead, it gives you a value P(Y | X) / P(Y). To see this, note that the density we get out should approximate

P(X,Y) / ( P(x) P(y) ) = P(X) P(Y | X) / ( P(x) P(y) ) = P(Y | X) / P(Y)

Thus, if P(Y | X) is independent of X–and therefore equal to P(Y)–you will recover a constant function. So if the function looks non-constant, you know that there is some kind of interaction. This thing P(Y | X) / P(Y) can be interpreted as how much the probability changes by conditioning on X. That is, how much the probability of that specific Y value depends on whether or not you know X. If you take logs, you get how much the suprisal of Y changes by conditioning on X.

An example is shown below:

This was generated by taking X to be distributed uniformly in [0,1]. Y is then chosen to be a normal distribution with SD 1, centered at 10*X. I ran this density estimation technique, and plotted the log densities. As things get more red, they get more negative. A script to generate it (and can be modified to play with other data) can be found here. If you squint a bit, you can see how this represents that distribution’s nonindependence: for example, if you don’t know X, a high Y is very possible and maybe as likely as a small Y. But if you know X is small, a high Y is very unlikely. So, when X is small and Y is large, you get red, meaning that by conditioning on X, the probability in that region goes down a lot. Hopefully this can help visualize how they are nonindependent.

Something I see sometimes is people running experiments and repeatedly checking p-values to see when they are significant. I don’t know what this is called, so I’ll call it p-value sniffing–you sniff the p-values every so often, and stop testing more subjects when things smell good for your hypothesis. Everyone knows that p-value sniffing is bad, but I wondered how much it really hurt, so I made a quick R script (here) to play around with it.

In the script, you can set N to be the number of subjects you run. There is a vector called sniff_at which sees which values you sniff p-values at. The script works by pretending to run a t-test experiment in which the null hypothesis is true, and sniffing p-values at each of the sniff_at values. If it finds anything significant at any of these sniffed values, it stops the experiment and rejects the nully hypothesis. The script simulates a bunch of runs, and prints out the average number of times the null hypothesis generates data sets which are significantly different.

Therefore, you can interpret the output of this program as something like adjusted p-values: how often do you get a large enough difference of means under the null hypothesis and sniffing experimental procedure. If you didn’t sniff at all (until your one sniff at the end of the experiment) you should get the correct significance value.

Here is some output for N=40 and N=20, and two different significance levels:

You can see that when we just test at 40, we get about the right significance levels of p=0.0518 and p=0.0098. But with more sniffing, your stats go to crap, but not as fast as I would have thought. It’s certainly not as bad as doing multiple comparisons, but still not good!