Project Euler 744

发表于 分类于 阅读次数： 本文字数： 4.1k 阅读时长 ≈ 4 分钟

Project Euler 744

题目

What? Where? When?

“What? Where? When?” is a TV game show in which a team of experts attempt to answer questions. The following is a simplified version of the game.

It begins with \(2n+1\) envelopes. \(2n\) of them contain a question and one contains a RED card.

In each round one of the remaining envelopes is randomly chosen. If the envelope contains the RED card the game ends. If the envelope contains a question the expert gives their answer. If their answer is correct they earn one point, otherwise the viewers earn one point. The game ends normally when either the expert obtains n points or the viewers obtain n points.

Assuming that the expert provides the correct answer with a fixed probability \(p\), let \(f(n,p)\) be the probability that the game ends normally (i.e. RED card never turns up).

You are given (rounded to 10 decimal places) that

\(f(6,\dfrac{1}{2})=0.2851562500,f(10,\dfrac{3}{7})=0.2330040743,f(10^4,0.3)=0.2857499982\).

Find \(f(10^{11},0.4999)\). Give your answer rounded to \(10\) places behind the decimal point.

解决方案

假设整个过程中 RED 从未出现，那么游戏就是一串独立伯努利试验：成功概率为 \(p\)，“失败”（观众得分）概率为 \(q\)。

令随机变量 \(X\) 表示在正常结束时，总共抽到的问题数（回合数）。那么,至少要抽 \(n\) 次问题，才能有一方到 \(n\) 分，因此 \(X \ge n\) ；且最晚不会超过 \(2n-1\) 次问题，因为到第 \(2n-1\) 次之前，两边最多都是 \(n-1\) 分，第 \(2n-1\) 次一定决出胜负。因此\(X \in \{n,n+1,\dots,2n-1\}.\)

如果游戏一共进行了 \(m\) 轮问题（即抽到 \(m\) 个问题信封），意味着前 \(m\) 次抽取都没有抽到 RED。

由于 \(2n+1\) 个信封中只有 1 个 RED，不放回抽取，前 \(m\) 次都避开 RED 的概率为\(\dfrac{2n}{2n+1}\cdot\dfrac{2n-1}{2n}\cdot\dfrac{2n-2}{2n-1}\cdots\dfrac{2n+1-m}{2n+2-m}= \dfrac{2n+1-m}{2n+1}= 1-\dfrac{m}{2n+1}.\)

若专家在第 \(m\) 回合赢得第 \(n\) 分，则第 \(m\) 回合必须答对（概率 \(p\)），并且前 \(m-1\) 回合中专家得了 \(n-1\) 分（观众得 \(m-n\) 分）

前 \(m-1\) 回合中安排 \(n-1\) 个“专家得分”的位置数为\(\dbinom{m-1}{n-1}.\)于是专家恰好在第 \(m\) 回合获胜的概率为\(\dbinom{m-1}{n-1} p^n q^{m-n}.\)

同理，观众在第 \(m\) 回合赢得第 \(n\) 分需要：第 \(m\) 回合专家答错（概率 \(q\)），前 \(m-1\) 回合观众已有 \(n-1\) 分，即专家已有 \(m-n\) 分。因此概率为\(\dbinom{m-1}{n-1} q^n p^{m-n}.\)

最终，因此在 RED 不出现 的前提下，正常结束回合数为 \(m\) 的概率为\(\dbinom{m-1}{n-1}\Big(p^n q^{m-n}+q^n p^{m-n}\Big).\)

再乘上“前 \(m\) 次没抽到 RED”的概率，得到：

\[ \Pr\{X=m\}= \left(1-\frac{m}{2n+1}\right) \binom{m-1}{n-1}\Big(p^n q^{m-n}+q^n p^{m-n}\Big). \]

对 \(m=n,\dots,2n-1\) 求和即得

\[ f(n,p)=\sum_{m=n}^{2n-1}\Pr\{X=m\}= \sum_{m=n}^{2n-1} \left(1-\frac{m}{2n+1}\right) \binom{m-1}{n-1}\Big(p^n q^{m-n}+q^n p^{m-n}\Big). \]

将 \(m=n+i\)（其中 \(i=0,\dots,n-1\)）代入，可写成更对称的形式：

\[ f(n,p)= \sum_{i=0}^{n-1} \frac{n+1-i}{2n+1} \binom{n+i-1}{n-1}\Big(p^n q^{i}+q^n p^{i}\Big). \]

这就是一个关于 \(p\) 的多项式。继续把上式稍作变形。利用恒等式\(m\dbinom{m-1}{n-1}=n\dbinom{m}{n},\)可得

\[ \begin{aligned} f(n,p) &=1-\sum_{m=n}^{2n-1}\frac{m}{2n+1}\binom{m-1}{n-1}\Big(p^n q^{m-n}+q^n p^{m-n}\Big)\\ &=1-\frac{n}{2n+1}\sum_{m=n}^{2n-1}\binom{m}{n}\Big(p^n q^{m-n}+q^n p^{m-n}\Big)\\ &=1-\frac{n}{2n+1}\sum_{i=0}^{n-1}\binom{n+i}{n}\Big(p^n q^{i}+q^n p^{i}\Big). \end{aligned} \]

接下来观察第一个和：\(\displaystyle{S_p=\sum_{i=0}^{n-1}\binom{n+i}{n}p^n q^i.}\)

考虑做 \(2n\) 次伯努利试验（成功概率为 \(p\)），令 \(Y \sim \mathrm{B}(2n,p)\)。事件 \(Y \ge n+1\) 等价于第 \(n+1\) 次成功发生在第 \(2n\) 次或更早。

设第 \(n+1\) 次成功发生在第 \(n+i+1\) 次（即前面恰有 \(i\) 次失败），那么对 \(i=0,\dots,n-1\)，并且：前 \(n+i\) 次中恰好有 \(n\) 次成功：\(\dbinom{n+i}{n}\)，此时概率为\(p^n q^i\)，此外在第 \(n+i+1\) 次必须成功：需要再乘 \(p\)。

于是\(\displaystyle{\Pr\{Y\ge n+1\}=\sum_{i=0}^{n-1}\binom{n+i}{n}p^{n+1}q^i= pS_p.}\)所以有\(S_p=\dfrac{1}{p}\Pr\{\mathrm{B}(2n,p)\ge n+1\}.\)同理有\(\displaystyle{S_q=\sum_{i=0}^{n-1}\dbinom{n+i}{n}q^n p^i=\frac{1}{q}\Pr\{\mathrm{B}(2n,q)\ge n+1\}.}\)

因此

\[ f(n,p)= 1-\frac{n}{2n+1}\left( \frac{\Pr\{\mathrm{B}(2n,p)\ge n+1\}}{p}+ \frac{\Pr\{\mathrm{B}(2n,q)\ge n+1\}}{q} \right). \]

再用正则化不完全 Beta 函数表示二项分布尾概率：\(\Pr\{\mathrm{B}(2n,p)\ge n+1\}=I_p(n+1,n),\)

最终闭式为

\[ f(n,p)= 1-\frac{n}{2n+1}\left( \frac{I_p(n+1,n)}{p}+ \frac{I_q(n+1,n)}{q} \right), q=1-p. \]

本题 \(p=0.4999 < 1/2\)，且 \(n=10^{11}\) 极大。在上面的闭式\(f(n,p)\)中，关键在于估计 \(I_p(n+1,n)\) 与 \(I_q(n+1,n)\) 的量级。

令随机变量 \(Y\sim\mathrm{B}(2n,p)\)，则\(\mathbb{E}[Y]=2np.\)当 \(p<\dfrac12\) 时，有 \(\mathbb{E}[Y] < n\)，而我们关心的事件是 \(Y\ge n+1\)，即\(\dfrac{Y}{2n}\ge \dfrac{n+1}{2n}=\dfrac12+\dfrac{1}{2n}.\) 这要求 \(\dfrac{Y}{2n}\) 从其极限均值 \(p<\dfrac12\) 向右偏移到接近 \(\dfrac12\) 的位置，是一个典型的大偏差事件。

由大数定律，\(\dfrac{Y}{2n}\to p\)，因此该事件的概率趋于 \(0\)，从而得到极限结论

\[ \lim_{n\to\infty} I_p(n+1,n)=\lim_{n\to\infty}\Pr\big(\mathrm{B}(2n,p)\ge n+1\big)=0,p<\dfrac12. \]

同理，因为 \(q=1-p>\dfrac12\)，则\(I_q(n+1,n)=\Pr\big(\mathrm{B}(2n,q)\ge n+1\big)\to 1.\)

由于\(n\)已经趋近于非常大的数，因此经过化简，可以得到

\[\lim_{n\to\infty} f(n,p)=1-\frac{1}{2\max(p,1-p)}, 0<p<1.\]

代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
from math import comb

N = 10 ** 11
P = 0.4999
def f_bruteforce(n: int, p: float) -> float:
    q = 1.0 - p
    s = 0.0
    for m in range(n, 2 * n):
        s += (1.0 - m / (2 * n + 1.0)) * comb(m - 1, n - 1) * (
            (p ** n) * (q ** (m - n)) + (q ** n) * (p ** (m - n))
        )
    return s


def f(n: int, p: float) -> float:
    if n >= 10 ** 7:
        return 1 - 1 / (2 * max(p, 1 - p))
    else:
        return f_bruteforce(n, p)

print(f"{f(N, P):.10f}")
1
2
3
4
5
6
7
8
9
10
11
from scipy.special import betainc

N = 10 ** 11
P = 0.4999
def f(n: int, p: float):
    q = 1.0 - p
    Ip = betainc(n + 1, n, p)
    Iq = betainc(n + 1, n, q)
    return 1.0 - (n / (2 * n + 1.0)) * (Ip / p + Iq / q)

print(f"{f(N,P):.10f}")