Project Euler 852

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

Project Euler 852

题目

Coins in a Box

This game has a box of \(N\) unfair coins and \(N\) fair coins. Fair coins have probability \(50\%\) of landing heads while unfair coins have probability \(75\%\) of landing heads.

The player begins with a score of \(0\) which may become negative during play.

At each round the player randomly picks a coin from the box and guesses its type: fair or unfair. Before guessing they may toss the coin any number of times; however, each toss subtracts \(1\) from their score. The decision to stop tossing and make a guess can be made at any time. After guessing the player’s score is increased by \(20\) if they are right and decreased by \(50\) if they are wrong. Then the coin type is revealed to the player and the coin is discarded.

After \(2N\) rounds the box will be empty and the game is over. Let \(S(N)\) be the expected score of the player at the end of the game assuming that they play optimally in order to maximize their expected score.

You are given \(S(1) = 20.558591\) rounded to \(6\) digits after the decimal point.

Find \(S(50)\). Give your answer rounded to \(6\) digits after the decimal point.

解决方案

令\(N=50\)。关键观察：玩家在一轮内的所有掷币与何时停止只影响这一轮的得分，不影响之后盒子里剩余硬币的组成（因为无论怎么掷，最后都会揭示类型并丢弃这枚硬币）。因此整局最优期望可以写成对每个“轮开始时剩余硬币数量状态”的加权求和。

设某一轮开始时盒子里还有 \(f\) 枚公平、\(u\) 枚偏置（\(f+u\ge1\)）。定义\(E(f,u)\)表示在该轮开始处于状态\((f,u)\)时，这一轮的最优期望得分。

那么整局期望满足\(\displaystyle S(N)=\sum_{(f,u)\in \mathcal{A}_N\setminus\{(0,0)\}} P_{f,u}E(f,u),\)其中\(\mathcal{A}_N=\{(f,u)\mid 0\le f\le N,0\le u\le N,f+u\ge 1\}.\)表示可达状态集合，\(P_{f,u}\) 是某一轮开始时盒子处于 \((f,u)\)的概率。这个 \(P\) 完全由不放回随机抽样决定，与策略无关，后面可以用简单 DP 计算。

下面先求单轮最优值 \(E(f,u)\)。

这一轮开始时，先验：\(\Pr(\text{fair})=\dfrac{f}{f+u}, \Pr(\text{unfair})=\dfrac{u}{f+u}.\)

若已经掷了 \(n\) 次，其中 \(k\) 次为正面（反面 \(n-k\)），则在给定“只关心正面次数 \(k\)”的条件下：若硬币公平：\(\Pr(k\mid \text{fair})=\dbinom{n}{k}\left(\dfrac12\right)^n\)；若硬币偏置：\(\Pr(k\mid \text{unfair})=\dbinom{n}{k}\left(\dfrac34\right)^k\left(\dfrac14\right)^{n-k}=\dbinom{n}{k}\dfrac{3^k}{4^n}.\)乘上先验得到未归一化权重

\[ w_f=\frac{f}{f+u}\binom{n}{k}\left(\frac12\right)^n, w_u=\frac{u}{f+u}\binom{n}{k}\frac{3^k}{4^n}. \]

计算后验时 \(\dbinom{n}{k}\) 会相互抵消，且整体比例不重要。为了数值稳定，把两者同乘 \(2^{2n}\)（不改变比值）：\(W_f=f\cdot 2^n, W_u=u\cdot 3^k.\)于是

\[ p_f=\Pr(\text{fair}\mid n,k)=\frac{W_f}{W_f+W_u}, p_u=\Pr(\text{unfair}\mid n,k)=\frac{W_u}{W_f+W_u}=1-p_f. \]

定义价值函数 \(V_{f,u}(n,k)\)：在本轮已经掷了 \(n\) 次得到 \(k\) 次正面后，不计此前这 \(n\) 次掷币扣分（这些扣分由递推里的 \(-1\) 处理），从现在开始在最优策略下本轮还能获得的期望得分。显然\(E(f,u)=V_{f,u}(0,0).\)我们有以下两个选择：

  1. 现在就猜。此时最优就是猜更可能的类型，正确概率 \(\max(p_f,p_u)\)，错误概率 \(\min(p_f,p_u)\)。不再付掷币成本，因此立即猜的期望得分为\(G(n,k)=20\max(p_f,p_u)-50\min(p_f,p_u)=70\max(p_f,p_u)-50.\)
  2. 再掷一次。下一掷为正面的条件概率是混合分布：\(p_H=\Pr(H\mid n,k)=\dfrac12 p_f+\dfrac34 p_u,p_T=\Pr(T\mid n,k)=\dfrac12 p_f+\dfrac14 p_u=1-p_H.\)若再掷一次，需要立刻付出 \(1\) 分代价。掷出正面后状态变为 \((n+1,k+1)\)，掷出反面后变为 \((n+1,k)\)。因此至少再掷一次的期望为\(C(n,k)=p_H\bigl(V_{f,u}(n+1,k+1)-1\bigr)+p_T\bigl(V_{f,u}(n+1,k)-1\bigr).\)

综上得到状态转移方程：\(V_{f,u}(n,k)=\max(G(n,k),C(n,k)).\)这已经足以计算所有 \(E(f,u)\)。

注意到：如果决定“至少再掷一次”，那么从现在起无论后面信息多么充分，最终最多也只能做到必然猜对得 \(20\) 分，但至少要再扣 \(1\) 分，因此若继续掷至少一次，则最终期望 \(\le 20-1=19\).因此只要当前立即猜的期望 \(G(n,k)\ge 19\)，继续掷不可能更优，必有\(V_{f,u}(n,k)=G(n,k) (G(n,k)\ge 19).\)

实际计算中再设一个最大深度 \(n\le D=4N+2\) 的截断：当 \(n\) 达到该上限且仍未触发 \(G(n,k)\ge19\) 时，用上界 \(19\) 作为返回值即可（因为此时最优值也不可能超过 \(19\)）。由于要在 \(D\) 次掷币后仍保持强不确定性，本身在真实路径上概率极小，对要求的 \(10^{-6}\) 精度影响可忽略。

现在计算整局 \(S(N)\)。令 \(P_{f,u}\) 为“某一轮开始时剩余 \((f,u)\)”的概率。初始\(P_{N,N}=1.\)

在状态 \((f,u)\) 开始这一轮时，本轮期望贡献为 \(E(f,u)\)。随后随机抽到并丢弃一枚：抽到公平的概率 \(\dfrac{f}{f+u}\)，转移到 \((f-1,u)\)；抽到偏置的概率 \(\dfrac{u}{f+u}\)，转移到 \((f,u-1)\)。因此一边累加答案一边转移即可。

代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
from functools import lru_cache

from tools import gcd

N = 50
MAXD = 4 * N + 2
SHIFT = MAXD + 2

pow2 = [1.0] * (MAXD + 2)
for i in range(1, MAXD + 2):
    pow2[i] = pow2[i - 1] * 2.0

pow3 = [1.0] * (MAXD + 2)
for i in range(1, MAXD + 2):
    pow3[i] = pow3[i - 1] * 3.0


@lru_cache(None)
def E_ratio(a, b):
    if a == 0 or b == 0:
        return 20.0

    dp = [0.0] * (SHIFT * SHIFT)
    base = (MAXD + 1) * SHIFT
    for k in range(SHIFT):
        dp[base + k] = 19.0

    for n in range(MAXD, -1, -1):
        off = n * SHIFT
        off1 = (n + 1) * SHIFT
        fairW = a * pow2[n]
        for k in range(n, -1, -1):
            unfairW = b * pow3[k]
            total = fairW + unfairW
            pf = fairW / total

            if pf >= 0.5:
                guess = 70.0 * pf - 50.0
            else:
                guess = 70.0 * (1.0 - pf) - 50.0

            if guess >= 19.0:
                dp[off + k] = guess
                continue

            pu = 1.0 - pf
            pH = 0.5 * pf + 0.75 * pu
            cand = pH * (dp[off1 + k + 1] - 1.0) + (1.0 - pH) * (dp[off1 + k] - 1.0)
            dp[off + k] = max(cand, guess)

    return dp[0]


def E(f, u):
    if f == 0 or u == 0:
        return 20.0
    g = gcd(f, u)
    return E_ratio(f // g, u // g)


def S(N):
    prob = [[0.0] * (N + 1) for _ in range(N + 1)]
    prob[N][N] = 1.0
    ans = 0.0

    for f in range(N, -1, -1):
        for u in range(N, -1, -1):
            if f == 0 and u == 0:
                continue
            p = prob[f][u]
            if p == 0.0:
                continue
            ans += p * E(f, u)
            tot = f + u
            if f > 0:
                prob[f - 1][u] += p * (f / tot)
            if u > 0:
                prob[f][u - 1] += p * (u / tot)

    return ans


val = S(N)
print(f"{val:.6f}")