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

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