突发奇想

看到一个视频（源链接：https://www.bilibili.com/video/BV1MHzHY2EQ6）：甲和乙猜拳，甲连赢 3 次算赢，累计输 10 次算输，问甲获胜的概率是多少？

每一次猜拳，都有三种可能性，每一种可能性是 $\frac{1}{3}$。假设 $P(L, W)$ 代表在已经 L 局并且连赢了 W 局的情况下，甲获胜的概率，那么有以下递推式：

$$ \left\{ \begin{aligned} & P(l, 2) = \frac{1}{3} + \frac{1}{3} P(l, 0) + \frac{1}{3} P(l + 1, 0)\\ & P(l, 1) = \frac{1}{3} P(l, 2) + \frac{1}{3} P(l, 0) + \frac{1}{3} P(l + 1, 0)\\ & P(l, 0) = \frac{1}{3} P(l, 1) + \frac{1}{3} P(l, 0) + \frac{1}{3} P(l + 1, 0) \end{aligned} \right. $$

可以解得

$P(l, 0) = \frac{13}{14} P(l + 1, 0) + \frac{1}{14}$

代码如下：

from fractions import Fraction
from random import Random

def solve_number(times: int, L: int, W: int) -> Fraction:
    """这是计算机模拟代码"""

    wins = 0
    random = Random()
    for i in range(times):
        win_times = 0
        loss_times = 0
        while True:
            result = random.randint(1, 3)
            if result == 1:
                win_times += 1
            elif result == 2:
                win_times = 0
            else:
                win_times = 0
                loss_times += 1
            if win_times == W:
                wins += 1
                break
            if loss_times == L:
                break
    return Fraction(wins, times)

from fractions import Fraction
from random import Random

def solve() -> Fraction:
    W = 3
    L = 10
    F13 = Fraction(1, 13)

    dp: list[list[Fraction]] = [[Fraction(0) for j in range(W)] for i in range(L + 1)]

    for i in range(L - 1, -1, -1):
        dp[i][0] = Fraction(13, 14) * dp[i + 1][0] + Fraction(1, 14)
        dp[i][2] = F13 + F13 * dp[i][0] + F13 * dp[i + 1][1]
        dp[i][1] = F13 * dp[i][2] + F13 * dp[i][0] + F13 * dp[i + 1][1]

    return dp[0][0]

但是！ 能不能求在任意 W 和 L 下的胜率呢？

在任意 W 和 L 下，有如下递推式：

$$ \left\{ \begin{aligned} & P(l, W) = \frac{1}{3} + \frac{1}{3} P(l, 0) + \frac{1}{3} P(l + 1, 0)\\ & …\\ & P(l, w) = \frac{1}{3} P(l, w + 1) + \frac{1}{3} P(l, 0) + \frac{1}{3} P(l + 1, 0)\\ & …\\ & P(l, 1) = \frac{1}{3} P(l, 2) + \frac{1}{3} P(l, 0) + \frac{1}{3} P(l + 1, 0)\\ & P(l, 0) = \frac{1}{3} P(l, 1) + \frac{1}{3} P(l, 0) + \frac{1}{3} P(l + 1, 0) \end{aligned} \right. $$

代入：

$$ \begin{aligned} P(l, 0) & = \frac{1}{3} P(l, 1) + \frac{1}{3} P(l, 0) + \frac{1}{3} P(l + 1, 0)\\ & = \frac{1}{3} (\frac{1}{3} P(l, 2) + \frac{1}{3} P(l, 0) + \frac{1}{3}P(l + 1, 0)) + \frac{1}{3} P(l, 0) + \frac{1}{3} P(l + 1, 0)\\ & = \frac{1}{3^2} P(l, 2) + (\frac{1}{3} + \frac{1}{3^2}) P(l, 0) + (\frac{1}{3} + \frac{1}{3^2}) P(l + 1, 0)\\ & = …\\ & = \frac{1}{3^W} + (\frac{1}{3} + \frac{1}{3^2} + … + \frac{1}{3^W}) P(l, 0) + (\frac{1}{3} + \frac{1}{3^2} + … + \frac{1}{3^W}) P(l + 1, 0)\\ & = \frac{1}{3^W} + S_W * P(l, 0) + S_W * P(l + 1, 0) \end{aligned} $$

其中 $S_W$ 是首项为 $\frac{1}{3}$，公比为 $\frac{1}{3}$，项数为 $W$ 的等比数列的和。求解，得

$$ P(l, 0) = \frac{1}{1 - S_W} (\frac{1}{3^W} + S_W * P(l + 1, 0)) $$

剩下的 $P(l, w)$ 都可以通过 $P(l, 0)$ 推出

代码如下：

from fractions import Fraction
from random import Random

def solve(L: int, W: int) -> Fraction:
    F13 = Fraction(1, 3)
    SW = F13 * (1 - F13**W) / (1 - F13)
    SW_REV = Fraction(1, (1 - SW))

    dp: list[list[Fraction]] = [
        [*(Fraction(0) for j in range(W)), Fraction(1)] for i in range(L + 1)
    ]

    for i in range(L - 1, -1, -1):
        dp[i][0] = SW_REV * (Fraction(1, 3**W) + SW * dp[i + 1][0])
        for j in range(W - 1, 0, -1):
            dp[i][j] = F13 * dp[i][j + 1] + F13 * dp[i][0] + F13 * dp[i + 1][1]

    return dp[0][0]

if __name__ == "__main__":
    L, W = 28, 4
    print(float(solve(L, W)), float(solve_number(1000, L, W)))

如果有更巧妙的办法，也欢迎大家讨论哇~

一觉醒来突然想到一个问题，既然都知道 $P(l, 0) = \frac{1}{1 - S_W} (\frac{1}{3^W} + S_W * P(l + 1, 0))$ 了，令

$$ \begin{aligned} S_{W2} & = \frac{S_W}{1 - S_W}\\ S_{W3} & = \frac{1}{1 - S_W} * \frac{1}{3^W} \end{aligned} $$

于是

$$ \begin{aligned} P(0, 0) & = S_{W3} + S_{W2} P(1, 0)\\ & = S_{W3} + S_{W2} S_{W3} + (S_{W2})^2 P(2, 0)\\ & = S_{W3} + S_{W2} S_{W3} + (S_{W2})^2 S_{W3} + (S_{W2})^3 P(3, 0)\\ & = …\\ & = S_{W3}(1 + S_{W2} + (S_{W2})^2 + … + (S_{W2})^{L - 1}) + (S_{W2})^{L} P(L, 0)\\ & = S_{W3} * \frac{1 - (S_{W2})^{L}}{1 - S_{W2}} + (S_{W2})^L P(L, 0)\\ & = S_{W3} * \frac{1 - (S_{W2})^{L}}{1 - S_{W2}} \end{aligned} $$

其中 $P(L, 0) = 0$ （因为甲都已经输了）

代码如下：

from fractions import Fraction
from random import Random

def solve(L: int, W: int) -> Fraction:
    F13 = Fraction(1, 3)
    SW = F13 * (1 - F13**W) / (1 - F13)
    SW2 = SW / (1 - SW)
    SW3 = 1 / ((1 - SW) * 3**W)
    SW2_SUM = Fraction(1 - SW2**L, 1 - SW2)

    return SW3 * SW2_SUM

if __name__ == "__main__":
    L, W = 28, 4
    print(float(solve(L, W)), float(solve_number(1000, L, W)))