121. 碟子游戏奖金(Disc game prize fund)

一开始包里装有一个红色碟子和一个蓝色碟子。在一场概率游戏中,每一轮玩家从包中取出一个碟子,记录其颜色,随后将碟子放回包中,并在包中加入一个红色碟子,再进行下一轮。

玩家需要付一英磅来玩这个游戏,如果他们在游戏结束时拿出的蓝色碟子比红色碟子更多,则获得胜利。

如果游戏进行四轮,那么玩家获胜的概率是\(11/120\),因此游戏所设定的最高奖金为十英磅,否则可能会造成亏损。注意奖金必须是整数,而且包含了玩家付出用于玩游戏的一英磅,也就是说玩家实际上赢得的数额是九英磅。

如果游戏进行十五轮,求游戏应该设定的最高奖金。

分析:首先我们先不管如何求玩家获胜的概率,假设我们已经求出了这个概率,那我们作为庄家应该如何来设定最高奖金避免出现亏损?看一下题目中的例子,在四轮游戏中玩家获胜的概率是\(11/120\)。假设我们设定的最高奖金为\(x\),那我们支付奖金的期望为\(\frac{11}{120}x\),为了避免亏损,庄家收到的一英磅至少要等于支付奖金的期望,也就是\(\frac{11}{120}x=1\Rightarrow x=\frac{120}{11}\approx10.91\),又由于奖金必须是整数,所以我们要把这个数向下取整得到十英磅。所以一般地,如果玩家获胜的概率为\(p\),则庄家设定的最高奖金应为\(\lfloor 1/p\rfloor\)

现在的问题是我们应该如何计算玩家获胜的概率,我采用了两种思路,分别是递归求解的思路和生成函数的思路,下面分别来看:

一、递归求解思路

还是来看题目中的例子。总共进行四轮游戏,每进行一轮加一个红色盘子,所马每轮游戏的盘子总数分别是\(2,3,4,5\),其中蓝色盘子都只有一个,剩下的就是红色盘子。假如进行四轮游戏,只要在玩家取出的蓝色盘子比红色盘子多的时候玩家才会赢,所以玩家必须取出三个或四个蓝色盘子,我们需要求出在四轮游戏中玩家获得三个和四个蓝色盘子的概率。假设我设玩家在\(n\)轮游戏中取出\(k\)个蓝色盘子的概率为\(p(n,k)\),这可分为两种情况:

因此,我们可以得到\(p(n,k)\)的递归表达式为:

$$ p(n,k)=p(n-1,k)*\frac{n}{n+1}+p(n-1,k-1)*\frac{1}{n+1} $$

再来看一些边界条件,主要有两个:

因此,完整的状态转移方程为:

$$ p(n,k)=\left\{\begin{matrix} \frac{1}{n+1} & k=0\\ \frac{1}{(n+1)!} & n=k\\ p(n-1,k)*\frac{n}{n+1}+p(n-1,k-1)*\frac{1}{n+1} & else \end{matrix}\right. $$

在总共有\(n\)轮的游戏中,玩家至少有取出\(\lceil{n/2}\rceil\)个蓝色盘子才能赢得游戏,也就是在总共有十五轮的游戏中,至少要取出八个蓝色盘子,所以只需要求出\(\sum_{i=8}^{15}p(15,i)\)就可以了。

二、生成函数思路

从上面的推理我们可以明显看出,获胜概率的分母为\((n+1)!\),对于十五轮的游戏,获胜概率的分母为\(16!\)。对于分子,我采用了生成函数的思路。假设我用\(x\)表示蓝色盘子,用\(y\)表示红色盘子,在第一轮游戏中红色与蓝色盘子各有一个,我把它表示成为\((x+y)\),第二轮有一个蓝色盘子两个红色盘子,可以表示成为\((x+2y)\)。那么进行四轮游戏,我们可以把每轮的表示相乘,有:

$$ (x+y)(x+2y)(x+3y)(x+4y)=\displaystyle x^{4} + 10 x^{3} y + 35 x^{2} y^{2} + 50 x y^{3} + 24 y^{4} $$

取出四个和三个蓝色盘子的情况分别对应着\(x\)指数为4和3的情况,它们的系数则是总的可能性数,两个系数加起来的结果正好是11。对于总共十五轮的游戏,也是一样的道理,有:

$$ \prod_{i=1}^{15}(x+iy)=\displaystyle x^{15} + 120 x^{14} y + 6580 x^{13} y^{2} + 218400 x^{12} y^{3} + 4899622 x^{11} y^{4} + 78558480 x^{10} y^{5} + \cdots $$

我们只需要关心其中\(x\)项指数大于等于八的情况下系数就可以了,把这些系数加总起来即为玩家获胜的总的可能性数。

三、代码实现

对于第一种思路,为了避免浮点误差,我使用了标准库\(fractions\)对概率用分数形式进行计算。对于第二种思路,我使用了\(python\)的符号计算库\(sympy\),求出生成函数多项式的相应系数。具体代码如下:

# approach 1, time cost = 490 µs ± 8.63 µs

import sympy as sp
from math import prod,factorial

def main(n=15):
    x = sp.var('x')
    expr = prod([(x+i) for i in range(1,n+1)])
    p = sp.Poly(expr)
    res = factorial(n+1)//sum(p.coeffs()[:n//2+1])
    return res

# approach 2, time cost = 1.01ms

from fractions import Fraction as f
from functools import lru_cache
from math import factorial

@lru_cache(maxsize=None)
def prob(n,k):
    if k == 0:
        return f(1,n+1)
    elif n == k:
        return f(1,factorial(n+1))
    else:
        return prob(n-1,k)*f(n,n+1)+prob(n-1,k-1)*f(1,n+1)

def main(n=15):
    p = sum([prob(n,i) for i in range(n//2+1,n)])
    return 1//p