Optimizing a stochastic "flip and prune" procedure for selecting a subset of coins

I place some number of coins, $(c_1, ..., c_N) \in C$ on a table, where each coin is originally tails up. Let's call the "tails" state $0$ and the "heads" state $1$. I then perform the following procedure some number of times $R$ with the goal of having only a specific subset of $M$ coins, $T \subset C$, left on the table when I'm finished:

(a) I sequentially select each coin and flip it to the "heads" or $1$ state with probability $p$, and do nothing with probability $(1-p)$.

(b) If all of the coins $c_i \in T$ were flipped to the "heads" state in (a), I remove (UPDATE) some fraction $q$ of the coins on the table in the "tails" state (perhaps with some Gaussian error). Otherwise I do nothing.

(c) I flip all coins on the table back to the "tails" or $0$ state.

As a function of $N$ and $M$, what value of $p$ minimizes the mean number of necessary rounds for the above procedure until I'm left with only my desired subset of coins $T$? With this optimized value of $p$, what probability distribution should I expect for the number of necessary rounds of the above procedure?

If there isn't a good general procedure for performing the minimization, what if we apply the restriction that $M = 1$, i.e. that we only want to select for a single coin $c_i$ to be left on the table?

Also, is it much easier to find a solution if one allows the value of $p$ to be adjusted each round with feedback for the number of leftover coins on the table?

• @DouglasZare Yes, you are correct. However, if there is some particularly "nice" way to perform the minimization for some metric other than the mean, I would certainly be interested to hear it. Apr 29 '13 at 6:03
• The distribution can be described as (the maximum of some geometric distributions) compounded with a negative binomial distribution. I believe the maximum of IID geometric distributions is surprisingly complicated, with a mean of about $-\log_{1-p}(M-N) + O(1)$, but there is a $\Theta(1)$ oscillating term which depends on the fractional part of $\log_{1-p}(M-N)$. Apr 29 '13 at 7:25

First, suppose $|T| = 0$. The number of times you toss every coin for each coin to turn up heads is the maximum of $N$ IID geometric random variables with success probability $1-p$. (There are two conventions for geometric random variables differing by $1$. Some people count failures before the first success, but here we need tosses up to and including the first success.) Let the average number of tosses needed be $a(N,p)$.

If we increase $N$ and $|T|$ by the same amount, then some of the tosses don't count because at least one of the specified subset doesn't come up heads. It takes an average of $p^{-|T|}$ tosses to get one which counts. The average number of tosses before you have selected $T$ is $a(n,p) p^{-|T|}$, where $n = N-|T|$. Indeed, the whole distribution can be viewed as a negative binomial distribution compounded with the maximum of IID geometric distributions: First take the maximum of IID geometric random variables $X$, and then sample from a negative binomial distribution with success rate $p^{|T|}$ and count $X$.

To optimize, we would like to know $a(n,p)$. The exact value is complicated in the way I commented, but it does have some useful estimates, and some exact formulas. See Mike Spivey's answers to "What is the expected maximum out of a sample (size N) from a geometric distribution?" A lot of useful information can be found in Bennett Eisenberg, 2008. "On the expectation of the maximum of IID geometric random variables." Statistics & Probability Letters 78 (2008) 135–143. This includes short proofs of many earlier results. Here is an exact formula:

$$a(n,p) = \sum_{k=1}^n (-1)^{k+1} \frac{n \choose k}{1-p^k}$$

Unfortunately, this alternating sum seems hard to analyze directly. Here is an easy estimate by viewing geometric random variables as exponential random variables with corresponding rates rounded up. The average maximum of IID exponentials is easy to calculate. Let $\lambda = -\log p, p = \exp(-\lambda)$.

$$\frac{1}{\lambda} H_n \lt a(n,p) \lt 1+\frac{1}{\lambda} H_n$$

where $H_n = \sum_{k=1}^n \frac{1}{k}$ is the expected maximum of IID standard exponential random variables. $H_n = \log n + \gamma +o(1)$ where $\gamma$ is Euler's gamma constant. In a complicated sense, $a(n,p)$ may be viewed as approximately the midpoint $\frac{1}{2} + \frac{1}{\lambda}H_n$, although the error term does not go to $0$ as $n \to \infty$. The absolute magnitude of the error is small, and has a Cesàro mean of $0$.

If we pretend that $a(n,p) = \frac{1}{2} + \frac{1}{\lambda}(\log n + \gamma)$, then we can optimize $a(n,p) p^{-|T|}$. Using Mathematica I find one maximum on $[0,1]$:

$$p = n \exp (\gamma - \sqrt{(\gamma + \log n)(\gamma +\log n + 2/|T|)})$$

For example, if $n=100, |T|=1$, this estimates the optimal value of $p$ as $0.399$. If $n=100, |T|=2$, then $p = 0.620$. It appears that the optimal value of $p$ is much more sensitive to $|T|$ than to $n$, it is roughly a value which makes the specified subset all heads with probability $1/e$, $p \sim 1-\frac{1}{|T|+1/2}$, but it is slightly higher for small values of $n$. Since the optimal $p$ is relatively insensitive to $n$, allowing yourself to change $p$ based on $n$ in future rounds should slightly decrease the optimal value of $p$ for any given $n$.

For small values of $n$ and $|T|$, we can optimize the exact value, and compare this with the approximation. For $n=|T|=10$, the optimum value of $p$ is $0.906333$. The approximation gives $p=0.906357$. These result in averages of $80.965470$ and $80.965473$ rounds of tosses, respectively, a difference of about $3\times 10^{-6}$.