Drawing without replacement with a twist

Question

Imagine the following scenario: You have an urn with $r$ red balls and $b$ blue balls. The red balls have two hooks, the blue balls have one hook. Consequently, if $r=1$ and $b=1$, you would be twice as likely to fish a red ball rather than a blue ball.

What is the expected fraction of blue balls that you fish without replacement ($f_{bf}$) as function of the sample size ($n=r+b$) and the fraction of blue balls in the urn ($f_{bu} = b/n$)?

If I were to draw all balls, then $f_{bf}=f_{bu}$. If I were to draw one ball only, then $f_{bf}=\frac{f_{bu}}{2-f_{bu}}$. What if I draw somewhere in between 1 and $n$?

I can simulate the expected fraction (see below), but I am looking for how to calculate the formula. How would I go about that?

Answer 1

Although you do not mention the words, I think you are looking for the expected fraction of blue balls drawn without replacement.

In general there is not a closed form, though there are approximations. You should read about Wallenius' noncentral hypergeometric distribution

Answer 2

This past question is related.

If $r$ and $b$ are not too large, you can compute the distribution recursively in $O(rb)$ steps.

Let $p(x,y)$ be the probability that at some point there are $x$ red balls and $y$ blue balls left. If $x \lt r, y \lt b$ then

$$p(x,y) = \frac{2x+2}{2x+y+2}p(x+1,y) + \frac{y+1}{2x+y+1}p(x,y+1)$$

Replace the term with $0$ if $x=r$ xor $y=b$, and $p(r,b)=1$.

To compute the average for $1000$ draws from $1000$ red and $1000$ blue balls using this formula, $618.05$ red, took about $0.05$ seconds. Here is some C# code:

        double[,] p = new double[r+1, b+1];

        for (Int32 x = r; x>=0; x--)
            for (Int32 y = b; y >= 0; y--)
            {
                if (x == r && y == b)
                {
                    p[x, y] = 1.0;
                }
                else
                {
                    p[x,y] += (x==r)? 0 : (2*(double)x+2.0)/(2*x+2.0+y) * p[x+1,y];
                    p[x,y] += (y==r)? 0 : ((double)y+1.0)/(2*x+y+1.0) * p[x,y+1];
                }
            }