Project Euler 689

题目

Binary Series

For \(0 \le x \lt 1\), define \(d_i(x)\) to be the \(i\text{th}\) digit after the binary point of the binary representation of \(x\).

For example \(d_2(0.25) = 1\), \(d_i(0.25) = 0\) for \(i \ne 2\).

Let \(f(x) = \displaystyle{\sum_{i=1}^{\infty}\frac{d_i(x)}{i^2}}\).

Let \(p(a)\) be probability that \(f(x) \gt a\), given that \(x\) is uniformly distributed between 0 and 1.

Find \(p(0.5)\). Give your answer rounded to \(8\) digits after the decimal point.

解决方案

取一个截断位置 \(k\)，将 \(X\) 分解为\(X=S_k+R_k,\)其中\(\displaystyle{S_k=\sum_{i=1}^{k}\frac{B_i}{i^2},R_k=\sum_{i=k+1}^{\infty}\frac{B_i}{i^2}.}\)

这样，前 \(k\) 项 \(S_k\) 只含有限个独立 Bernoulli 变量，可以精确枚举分布；而尾项 \(R_k\) 是大量极小独立项的加权和，可以用中心极限定理做正态近似。

因此，目标概率变为\(p(a)=\mathbb P(S_k+R_k>a).\)

对任意 \(S_k=s\)，有条件概率\(\mathbb P(X>a\mid S_k=s)=\mathbb P(R_k>a-s).\)因此\(p(a)=\mathbb E\big[\mathbb P(R_k>a-S_k\mid S_k)\big]=\mathbb E\big[g(a-S_k)\big],\)

其中定义尾项的生存函数\(g(t)=\mathbb P(R_k>t).\)也就是说，只要我们能对所有可能的 \(S_k\) 取值计算 \(g(a-S_k)\)，再取平均即可。

由于 \(B_i\sim \text{B}\left(1,\dfrac12\right)\)，有\(\mathbb E[B_i]=\dfrac12, \mathrm{Var}(B_i)=\dfrac14.\)于是尾项均值为\(\displaystyle{\mu_k=\mathbb E[R_k]=\sum_{i=k+1}^{\infty}\frac{\mathbb E[B_i]}{i^2}=\frac12\sum_{i=k+1}^{\infty}\frac1{i^2}.}\)

利用巴塞尔级数\(\displaystyle{\sum_{i=1}^{\infty}\frac1{i^2}=\zeta(2)=\frac{\pi^2}{6},}\)可得

\[ \mu_k=\frac12\left(\frac{\pi^2}{6}-\sum_{i=1}^{k}\frac1{i^2}\right). \]

同理，尾项方差为\(\displaystyle{\sigma_k^2=\mathrm{Var}(R_k)=\sum_{i=k+1}^{\infty}\mathrm{Var}!\left(\frac{B_i}{i^2}\right)=\sum_{i=k+1}^{\infty}\frac{\mathrm{Var}(B_i)}{i^4}=\frac14\sum_{i=k+1}^{\infty}\frac1{i^4}.}\)又有级数\(\displaystyle{\sum_{i=1}^{\infty}\frac1{i^4}=\zeta(4)=\frac{\pi^4}{90},}\)因此

\[ \sigma_k^2 =\frac14\left(\frac{\pi^4}{90}-\sum_{i=1}^{k}\frac1{i^4}\right). \]

当 \(k\) 足够大时，尾项 \(R_k\) 是许多“很小的独立随机量”之和，其分布会非常接近正态分布：\(R_k\approx \mathcal N(\mu_k,\sigma_k^2).\)

于是对任意实数 \(t\)，有\(g(t)=\mathbb P(R_k>t)\approx 1-\Phi\left(\dfrac{t-\mu_k}{\sigma_k}\right),\)其中 \(\Phi\) 是标准正态分布的累积分布函数（CDF）。

进一步利用互补误差函数 \(\mathrm{erfc}\) 与正态分布的关系：\(1-\Phi(z)=\dfrac{1}{2}\mathrm{erfc}\left(\dfrac{z}{\sqrt2}\right),\)

代入 \(z=\dfrac{t-\mu_k}{\sigma_k}\)，得到\(\displaystyle{g(t)\approx\frac12\mathrm{erfc}\left(\frac{t-\mu_k}{\sqrt2\sigma_k}\right)=\frac12\mathrm{erfc}\left(\frac{t-\mu_k}{\sqrt{2\sigma_k^2}}\right).}\)

因此

\[ \mathbb P(X>a\mid S_k=s)=\mathbb P(R_k>a-s)\approx\frac12\mathrm{erfc}\left(\frac{a-s-\mu_k}{\sqrt{2\sigma_k^2}}\right). \]

由于 \(S_k\) 由 \(k\) 个独立 Bernoulli 位决定，其所有可能取值为\(\displaystyle{S_k \in \left\{\sum_{i=1}^{k}\frac{\varepsilon_i}{i^2}\ :\varepsilon_i\in\{0,1\}\right\},}\)共有 \(2^k\) 个等概率情形，每个概率为 \(\frac{1}{2^k}\)。

于是有

\[ p(a)=\mathbb E[g(a-S_k)] \approx \frac{1}{2^k}\sum_{s\in \mathcal S_k} \frac12\mathrm{erfc}\left(\frac{a-s-\mu_k}{\sqrt{2\sigma_k^2}}\right), \]

其中 \(\mathcal S_k\) 表示所有 \(S_k\) 可能取值的集合（重数计入即按所有位模式枚举）。

算法上，可从集合 \(\{0\}\) 开始迭代构造所有子集和：对每个 \(i=1,2,\dots,k\)，把所有已有的和再加上 \(\frac{1}{i^2}\) 追加进去。因此最终\(S_k\)一共有\(2^k\)个元素。

当 \(k\) 增大时：

• 枚举规模 \(2^k\) 指数增长；
• 但尾项 \(R_k\) 的方差 \(\sigma_k^2\) 快速衰减，正态近似精度显著提升。

实践中取 \(k=22\) 时，已经足够满足精度。

代码

import numpy as np
from math import pi
from scipy.special import erfc

a = 0.5
def p_binary_series(a,k = 22):
    """
    p(a) = P( sum_{i>=1} B_i/i^2 > a ), B_i ~ Bernoulli(1/2) iid
    Exact enumeration for first k bits + normal approximation for the tail.
    """
    i = np.arange(1, k + 1, dtype=np.float64)

    # Tail mean/variance:
    # mu = 1/2*(pi^2/6 - sum_{1..k} 1/i^2)
    # var = 1/4*(pi^4/90 - sum_{1..k} 1/i^4)
    s2 = np.sum(1.0 / (i * i))
    s4 = np.sum(1.0 / (i**4))
    mu = 0.5 * (pi * pi / 6.0 - s2)
    var = 0.25 * (pi**4 / 90.0 - s4)

    # Enumerate all subset sums of {1/i^2}_{i=1..k}
    sums = np.array([0.0], dtype=np.float64)
    for t in range(1, k + 1):
        w = 1.0 / (t * t)
        sums = np.concatenate([sums, sums + w])

    # Tail prob under normal approximation:
    # P(R > a - s) = 0.5 * erfc((a - s - mu) / sqrt(2*var))
    z = (a - sums - mu) / np.sqrt(2.0 * var)
    tail_prob = 0.5 * erfc(z)

    return float(tail_prob.mean())


ans = p_binary_series(a, 22)
print(f"{ans:.8f}")