Correcting bias in samples selected by a prediction

Question

Here is the scenario:

I'm trying to find as many golden tickets as I can, so that I can sell them to kids that want to go on a tour of Wonka's chocolate factory.

Fortunately, I have a machine that can tell me the probability that a ticket is present within a chocolate bar without having to open it. I can use this machine on chocolate bars before I must buy them, so I can use it to select the chocolate bars most likely to have tickets.

Of course, the machine isn't accurate, so the probability it gives me will have some error, but this error will be normally distributed around the real probabilities, so it will even out across multiple chocolate bars. This means that if I tested every chocolate bar, I could determine the overall probability of finding a winning ticket by averaging the predictions produced by the machine.

So I test 10,000 chocolate bars (its Costco), and pick the 100 with the highest probabilities as predicted by the machine.

Now, to test the machine's accuracy, I average its predictions for these 100 bars, and I get 0.2 - so there is a 1 in 5 chance that any given chocolate bar will have a ticket - so I expect to find 20 in this group of 100.

The problem? The average of the predictions is above the actual rate at which I discover tickets in this group of 100.

In fact, if I repeat this experiment many times, the average of the predictions is always higher than reality.

I understand why this is, because when I select those chocolate bars with the highest predicted probabilities, I'm more likely to select those with positive error, and less likely to select those with negative error.

My question is: is there any way I can correct for this bias, so that I can accurately estimate the number of chocolate bars in the group of 100?

Answer 1

To try to make sense of this, we need a model. As Didier said, the probability that a given bar has a ticket is really either 1 or 0. But I'll suppose that in the process of determining whether bar number $j$ gets a ticket, a continuous random variable $X_j$ is used, with values in $[0,1]$. The $X_j$ are independent and identically distributed, with density $f$; the conditional probability that bar $j$ has a ticket, given all the $X_k$, is $X_j$. The machine tells you $Z_j = X_j + Y_j$, where $Y_j$ has density $g$ and is independent of the $X_k$ and the other $Y_k$'s. The conditional probability that bar $j$ has a ticket, given $Z_j = z$, is then $\frac{\int_0^1 dx \ x f(x) g(z - x)}{\int_0^1 dx \ f(x) g(z-x)}$. I doubt that we can say much more without knowing something about $f$.

Answer 2

There is an extra layer of probability in the question which I think obscures the issue without adding anything. Suppose $X_1,...,X_{10000}$ are IID random variables (whose values are numbers which are not necessarily probabilities of anything). Suppose $Y_1,...,Y_{10000}$ are IID standard Gaussians. Let $X_i + Y_i=Z_i$ be observable, and let $a_1,...,a_{100}$ be the indices with the highest values among the $Z_i$. The question is how much larger $(Z_{a_1} + ... + Z_{a_{100}})/100$ is over $(X_{a_1} + ... + X_{a_{100}})/100$ (perhaps conditioned on the values observed). This problem seems fundamental, and I hope a statistician can say how this is typically addressed.

Even if the $X_i$ are constant, this is still interesting, and it shows that you need to know the values $100$ and $10000$, with a good approximation coming from the ratio $r=1/100$. If you simply know that you have noisy versions of the top $100$, you don't know how strong the bias is. (In concrete terms, you need to know how many experiments were not published to estimate the bias toward type I statistical errors in what is reported.)

The bias in the case that the $X_i$ are constant is an upper bound for the bias in general. The proof is simply that $\sum Y_{a_i}$ is at most the sum of the greatest $100$ values among the $Y_i$.

If we take the limit as the size of the sample goes to infinity, fixing the ratio between the sample and the population, then the bias can be computes in terms of the distribution for $X_i$. Let $T$ be the boundary of the top $r$ of the distribution of each $Z_i$. The bias is $1/r \int_{-\infty}^\infty \int_{T-x}^\infty y ~d\mu_Y d\mu_X$ which can be computed explicitly in some cases.

For example, if $X_i$ is distributed as $N(0,\sigma^2)$, then each $Z_i$ is distributed as $N(0,\sigma^2+1)$ so $T = \sqrt{\sigma^2+1}~ \Phi^{-1}(1-r)$, and the bias is $\frac{\exp(-T^2/(2\sigma^2+2)}{r~\sqrt{2\pi (\sigma^2 +1)}} = \frac{\exp(\Phi^{-1}(1-r)^2/2)}{r \sqrt{2 \pi(\sigma^2+1)}}$.

For comparison, the average value of $Z_{a_i}$ is $\frac{\exp(\Phi^{-1}(1-r)^2/2)\sqrt{\sigma^2+1}}{r \sqrt{2 \pi}},$ or $\sigma^2+1$ times as great, so $1/(\sigma^2+1)$ of the average observed value is the bias from the noise $Y_{a_i}$ and $\sigma^2/(\sigma^2+1)$ of the average observed value comes from the $X_{a_i}$.