中国有嘻哈决赛冠军平票的概率问题(代码实现)

Posted on | In | |

上一篇博客讨论了“中国有嘻哈”决赛平票的概率问题，分简单情况和复杂情况给出了计算方法和结果，但是留了一个小尾巴，就是复杂情况中的$f(x)$到底怎么算的。由于我只是为了快速算出结果，实现$f(x)$的代码非常丑陋，因此在上一篇博客中没有好意思把代码贴出来。感谢“Tianheng”同学的回复，给了我分享代码的勇气，如果我的代码有问题，欢迎讨论。下面就简述下我的代码实现，代码基于Python。

回顾上一篇博客的计算公式：$$\frac{\sum_{r=0}^{100} C_{100}^{r} f(125-r)}{2^{100}\times 51^{3}}$$
我们需要组合函数$C_n^k$和整数规划函数$f(x)$。

组合函数

组合函数我直接从网上找了一段代码：

1
2
3
4
5
6
7
8
# 组合函数
import operator
import math
def c(n,k):
    if k == 0 or k == n:
        return 1
    else:
        return reduce(operator.mul, range(n - k + 1, n + 1)) / reduce(operator.mul, range(1, k +1))

整数规划函数

接下来就是求解 $$x+y+z = v, x_{min} \leq x \leq x_{max}, y_{min} \leq y \leq y_{max}, z_{min} \leq z \leq z_{max}$$ 这个整数规划问题，实际上我们只需要统计可行解的数目即可。原谅我一下子没想出来3个变量的这个问题怎么求解，(感觉直接穷举x,y,z有点太暴力了)，我就先从简单点的问题开始考虑，如何求解包含2个变量的这种问题，也就是$$x+y = v, x_{min} \leq x \leq x_{max}, y_{min} \leq y \leq y_{max}$$ 这个问题。

两变量问题：

两个变量的问题比较容易，在给定$v$和$y$的取值范围后，$x$的实际可取值范围是可以计算出来的，$x$不能小于$x_{min}$，也不能小于$v-y_{max}$；$x$不能大于$x_{max}$，也不能大于$v-y_{min}$。在确定$x$后，因为$y=v-x$，实际上同时也确定了$y$，因此，只需要得出$x$的实际可取值范围，就可以知道有多少个可行解。下面是代码实现（请忽略我丑陋的函数和变量命名）：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
def solution2(inputs, value):
    """
    find solutions of x + y = value
    inputs is a list of tuples,
    [(x_min, x_max), (y_min, y_max)]
    each tuple gives the domain of each variable
    return is a tuple, indicates the valid range of the first variable
    """
    x_min = max(inputs[0][0], value - inputs[1][1])
    x_max = min(inputs[0][1], value - inputs[1][0])
    if x_min > x_max:
        return False
    else:
        return (x_min, x_max)

这里程序的返回结果不是可行解的个数，而是$x$的可取值范围，可行解的个数就是$x_{max} - x_{min} + 1$，这么写是为了方便写三变量问题代码，

三变量问题：

有了两变量问题的解决办法，三变量问题就可以比较容易的解决了，我们可以把三变量问题先分解为两变量问题，让$t=y+z$，先求解$x$和$t$的两变量问题，求出$x$的真实取值范围，然后对于每个可能的$x$的取值，求解$y+z=v-x$这个两变量问题。下面是代码实现（同样请忽略我丑陋的函数和变量命名）：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
def solution3(inputs, value):
    """
    find the count of solutions of x + y + z = value
    """
    x_min = inputs[0][0]
    x_max = inputs[0][1]
    yz_min = inputs[1][0] + inputs[2][0]
    yz_max = inputs[1][1] + inputs[2][1]
    sx = solution2([(x_min,x_max), (yz_min,yz_max)], value)
    if not sx:
        return 0
    nx_min, nx_max = sx 
    counts = 0
    for x in range(nx_min, nx_max + 1):
        c = solution2([(inputs[1][0], inputs[1][1]), (inputs[2][0], inputs[2][1])], value-x)
        if c:
            counts += c[1] - c[0] + 1
    return counts

函数solution3就是我们需要的$f(x)$。

求解完整问题

在有了组合函数和整数规划函数，就可以根据公式：$$\frac{\sum_{r=0}^{100} C_{100}^{r} f(125-r)}{2^{100}\times 51^{3}}$$

计算概率了，代码也比较简单：

1
2
3
4
5
6
7
8
domain = (0, 50)
inputs = [domain,domain,domain]
counts = 0
for r in range(0, 101):
    counts += c(100,r) * solution3(inputs, 125 - r)
prob = counts *1.0 / (math.pow(2,100) * 51**3)
print prob

函数给出的结果是：0.0145193025331