1-the-sock-drawer

最近在看《Fifty challenging problems in probability with solutions》^[1]，这是一本比较经典且有趣的概率习题集，本书不仅提供了50多道题目以及详细的解答，更是在书中阐述了许多概率思维，而且此书也是许多量化职位面试必备的参考书，因此非常值得一读。我想通过公众号记录下自己的阅读和思考的过程，以此加深自己的理解与记忆，与大家分享。

1. 问题

1. The Sock Drawer
   A drawer contains red socks and black socks. When two socks are drawn at random, the probability that both are red is .
   (a) How small can the number of socks in the drawer be?
   (b) How small if the number of black socks is even?

一个抽屉里有红色和黑色两种袜子，从中随中抽出两只袜子，都是红色的概率为 。那么这个抽屉里最少有多少只袜子？如果黑色袜子是偶数，抽屉里最少有多少只袜子？

2. 解答

这道题的解法比较暴力，直接硬算即可。
我们假设抽屉里分别有 只红袜子和 只黑袜子，根据题目一定有 。

所以，我们有

从而，直接求解上述关于 的方程即可，在求解之前可以简化一下这个问题。因为 ，一定有

求解上述两个不等式可得

这样我们得到了一个 和 之间的不等关系，这可以帮助我们关于缩小 的搜索范围。即先令 从 开始，求解上述不等式找到可能的 的取值，再逐个验证(1)式是否成立。

下面get_prob 函数是计算给定 和 时，抽出两只袜子均为红色的概率。find_rb函数用来循环查找给定 的范围时所有可能的解。

def get_prob(r=2, b=1):
    return r / (r + b) * (r - 1) / (r + b  - 1)

def find_rb(b_max=1000):
    for b in range(1, b_max):
        low = np.ceil((np.sqrt(2) + 1) * b).astype(int)
        high = np.ceil((np.sqrt(2) + 1) * b + 1).astype(int)
        for r in range(low, high):
            if get_prob(r, b) == 0.5:
                print('r={}, b={}, True'.format(r, b))
            else:
                continue

find_rb(1000)

我们令 ，找到如下四组可能的解：
r=3, b=1, True
r=15, b=6, True
r=85, b=35, True
r=493, b=204, True

至此，我们便可以回答原书中的两个问题:
(a) 抽屉里至少有 只袜子，其中 只为红色， 只为黑色。
(b) 若抽屉里黑色袜子为偶数，则抽屉里至少有 只袜子，其中 只为红色， 只为黑色。
反正总有袜子落单，哈哈。

3. 数值模拟

下面我们跑一下模拟实验，看一下我们上面找到的四组解是否确实概率为 。当然，直接带入(1)式可以直接验证啦。这里主要是为了画图直观感受一下，事件的频率确实是收敛到真实的概率 的。

def sample_socks(r=1, b=1, n=100000):
    drawer = np.zeros(r + b)  
    drawer[: r] = 1  # 1为红色，0为黑色
    res = []
    for i in range(n):
        samples = np.random.choice(drawer, 2, replace=False)  
        if samples.sum() == 2:
            res.append(1)
        else:
            res.append(0)
    return np.cumsum(res) / np.arange(1, n+1)

y1 = sample_socks(r=3, b=1)
y2 = sample_socks(r=15, b=6)
y3 = sample_socks(r=85, b=35)
y4 = sample_socks(r=493, b=204)

plt.plot(y1, label='r=3, b=1')
plt.plot(y2, label='r=15, b=6')
plt.plot(y3, label='r=85, b=35')
plt.plot(y4, label='r=493, b=204')
plt.ylim(0.4, 0.6)
plt.legend()
plt.show()

运行上述代码，我们可以得到如下图象：

参考

[1]

Mosteller F. Fifty challenging problems in probability with solutions[M]. Courier Corporation, 1987.: Book