记一道有趣的应用题

我亲爱的高中同学发来一个应用题目，好像困扰了她和室友一阵子。
主要是用Python写了一个小程序，模拟了一下概率。写的时候挺兴奋的，姑且把解题记录记在这里。

题目为：
在一个抽奖活动中，顾客有$4$次抽球机会，在有$10$种颜色、每种颜色各有$7$个，共计$70$个彩球中，每次不放回地随机取出$4$个，如果这四个球中存在两个颜色相同的球，则可以多抽取一次.

但一个人总共抽取的次数不超过$7$次。一旦抽到特定的含有“中奖”字样的小球时活动结束，顾客得到奖品。
每人次的票价是$29.9$元，奖品的成本约$35$元。

求问，为不亏本，含有“中奖”字样的小球占全部球的最大比例。

题目分解

实际上要求算每人次的中奖概率，所以可以把问题分解成先求每人所摸球次数的平均期望，因为每人不放回地每次摸$4$个球，所以可以通过前者算出每人能摸到的球数的平均期望，进而调整中奖小球的比例。

数学方法求算第一次摸球概率

不妨先计算第一次摸球的概率。
一开始直接联想到了超几何分布，不过超几何分布那都是高中学的了，该忘的已经忘得差不多。这种问题并不复杂，也肯定有人研究，我们不需要研究它，只是要找到前人的研究，然后把他的成果拿来用就好了。终于在网上搜索半天，边学边推导，发现了一个有意思的关键词叫做“多元超几何分布”。顺藤摸瓜找到了相关文献。
因为是多种颜色，不放回摸球，所以满足多元超几何分布的前提条件，可以直接代公式。
这张图片截自论文《多维超几何分布高阶混合矩的算法》:

我们计算摸到四个球颜色各不相同的情况。
公式中“$r$等品”对应本题中的颜色，另有：$N=70$, $n=4$.
那么概率为
$$C_{10}^{4} \times P_{(X_1=1,X_2=1,X_3=1,X_4=1)}=C_{10}^{4} \times\frac{C_{7}^{1}\times C_{7}^{1}\times C_{7}^{1}\times C_{7}^{1}}{C_{70}^{4}} =\frac{10\times 9\times 8\times 7\times 7^4}{70\times 69\times 68\times 67}\approx 54.99 \% $$
也就是说，顾客在第一次取球后没能获得奖励的取球次数的概率约为$55\%.$
但这只是计算了第一次取球，并未考虑不放回的多次取球对于样本总体的影响。

其他小伙伴的解法

大家想到了有趣的方法:

计算抽到四个不同颜色球的概率。
给定条件：
总共有$7$种颜色；每种颜色有$10$个球；总共$70$个球。

计算步骤：

第一个球：可以是任何颜色，概率为 $1$

第二个球：必须是与第一个球不同的颜色 $概率 = \frac{60}{69}$ (因为还剩$60$个不同颜色的球，总共剩$69$个球)

第三个球：必须是与前两个球不同的颜色 $概率 = \frac{50}{68}$

第四个球：必须是与前三个球不同的颜色 $概率 = \frac{40}{67}$

将这些概率相乘：
$$P_{(四个球都不同颜色)} = 1 \times \frac{60}{69} \times \frac{50}{68} \times \frac{40}{67} ≈ 0.2581$$
因此，抽到四个不同颜色球的概率约为$25.81\%$，或者说大约是$\frac{1}{4}$的几率。

她面对的题干好像和我的不大一样，大概要归咎于某个传话筒，但如果题干一致，这种方法得出的算式和我先前引用的文献里的这个一样，也是得到正确结果的一种好方法。

程序模拟方法模拟第一次摸球概率

接下来考虑采用Python模拟法计算。代码如下：

import random
import numpy as np  

i=1
j=1
win = 0
colour = 0
num = 0
arrmap = np.zeros((10, 7))

while (j <= 100 ):
    while (i <= 4 ):
        colour = random.randint(1, 10)
        num = random.randint(1, 7)
        if  arrmap[colour-1][num-1] != 1 :
            arrmap[colour-1][num-1] = 1
        else:
            i=i-1
            print ('本次取球重复')
        i= i+1
        # print (arrmap)

    row_sums = np.sum(arrmap, axis=1)
    i=1
    j=j+1
    arrmap = np.zeros((10, 7))
    print (row_sums)
    for x in row_sums:
        if x > 1 :
            win = win + 1
            print("此时成功")
 
print (win , "为取到同色球的次数")
# print ('循环正常退出')

其主要思路是，用一个二维数组模拟球的选中状态，行表示颜色，列表示球的序号，通过random产生两个随机数，定位每一次取到的球，如果此球已经取到，则再取一次。然后，通过numpy自带的行列求和，求出每行（也就是每种颜色）的和，如果大于1，则计数器win加一。每次循环重置一次数组。通过调节j的循环次数，增大实验次数。
经我测试检验，在循环$1000$次时，计数器值为$453$，符合数学方法计算得到的结果，也就是存在两球颜色相同的概率约为$45\%$，四个球颜色各不相同的概率约为$55\%.$

但是，由于二维数组求和，再判断其是否大于1的方法只适用于一次取球，要想计算多次取球，考虑第二次取球时，第一次取得的球还未放回，还需要更改计算逻辑。

程序模拟方法模拟多次摸球概率

因为并不知道用数学方法怎么计算，而程序模拟法能够看到一点儿希望，那不如优化一下程序。代码如下：

import random
import numpy as np  

people = 1
fetch_time_total = 0
bonus_total = 0
win_total = 0

while(people <= 1000):
    print("----------以下是第",people,"人取球")
    people += 1
    health_point = 4
    fetch_time = 1
    arrmap = np.zeros((10, 7))
    bonus = 0
    win = False
    
    i = 1
    while (i <= 3 ):
            colour = random.randint(1, 10)
            num = random.randint(1, 7)
            if  arrmap[colour-1][num-1] == 0 :
                arrmap[colour-1][num-1] = 777
            else:
                i -= 1
                print ('本次设置中奖球重复')
            i += 1

    while (fetch_time <= 7 and health_point >= 1 ):
        print('-----以下是第',fetch_time,'次取球，生命值为',health_point)
        get_bonus = False
        i = 1
        while (i <= 4 ):
                colour = random.randint(1, 10)
                num = random.randint(1, 7)
                if  arrmap[colour-1][num-1] == 0 :
                    arrmap[colour-1][num-1] = fetch_time
                else:
                    if arrmap[colour-1][num-1] == 777 :
                        win = True
                        win_total += 1
                        break
                    else:
                        i -= 1
                        print ('本次取球重复')
                i += 1
        # print(arrmap)
        if win :
            print ("中奖！")
            print("中奖次数",win_total)
            break

        for x in arrmap :
            count = 0
            for y in x:
                if y == fetch_time :
                    count += 1
                if count > 1 :
                    get_bonus = True
                    break
            if get_bonus :
                bonus += 1
                if bonus > 3 :
                    bonus -= 1
                    health_point -= 1
                    print("bonus达到上限")
                health_point += 1
                bonus_total += 1
                print("获得bonus,现在bonus值为",bonus)
                break
            # else:
            #     print("未获得bonus")

        health_point -= 1
        fetch_time += 1
        fetch_time_total = fetch_time_total + fetch_time - 1

print("共1000人取球","总重复取球数为",fetch_time_total,"总bonus为",bonus_total,"总中奖次数为",win_total)

最终用$1000$次模拟，分别计算了中奖小球在$1, 2, 3$ 时的中奖概率，分别为$33\%, 57\%, 72\%.$
这里的代码思路是，在每人抽球之前随机几个位置赋个高值当成中奖标记，判断颜色是否相同的时候顺便判断一下。然后正式摸球，摸球时如果摸到中奖球，则改变win的布尔值，进入下一个人的循环.

如果以票价为变量计算，则：

当中奖小球数为1，$\frac{0.33\times 35}{1-0.33}=17.24$
当中奖小球数为2，$\frac{0.57\times 35}{1-0.57}=46.39$
所以最好设定中奖小球数为$1.$

感叹

程序写得其实并不怎么优雅。不过实际上这算是我第一个Python小项目了，以前光说要学，但只是写个hello world罢了。
现在的经历使我找回了高中写C++程序的时光。
偶尔这样玩玩还是很有意思的。