利用中心极限定理求解圣彼得堡悖论问题的近似曲线

此文为《概率论》课程小项目。

关于圣彼得堡悖论的一些思考

记 \(N\) 为 游戏的轮数，则 \(N \sim Ge(\frac{1}{2}),P(N=k)=2^{-k},k=1,2,3,...\)

奖金 \(X=2^N\)，\(E(X)=E(2^N)=\sum_{k=1}^{+\infty} 2^k\times 2^{-k}=\sum_{k=1}^{+\infty} 1 = +\infty\)

理论上如果能玩无限轮游戏，便可获得无限的奖金。

下面作模拟：

import random
import matplotlib.pyplot as plt

MaxN = 10000000

def getAward():
    award = 1
    while(1):
        award *= 2
        if (random.random() <= 0.5):
            break
    return award

sumAward = 0

x_v = []
y_v = []
for i in range(0, MaxN):
    sumAward += getAward()
    x_v.append(i + 1)
    y_v.append(sumAward / (i + 1))

plt.xlabel("Count of rounds")
plt.ylabel("Average award")
plt.plot(x_v, y_v)

[<matplotlib.lines.Line2D at 0x1f843cec1c0>]

​

该算法的期望复杂度为\(O(\sum_{k=1}^{N}\sum_{j=1}^{+\infty})(\frac{1}{2})^j=O(N)\)

通过模拟发现现实中，在有限的游戏轮数的情况下，奖金并不是无限。

随着轮数的增加，平均奖金在逐渐升高。

下面通过引入置信度的方法来大致计算出上面这条曲线的函数。

考虑求条件期望 \(E(X|X\le M)\)，其中 \(M\) 是使 \(P\{ X \le M\} \ge 1-p\) 的最大值。

\(1-p\) 即为置信度。

若令 \(p=10^{-3}\)，则 \(M=2^{10}=1024\)，则 \(E(X|X \le 1024)=10/(1-p) \approx 10\)

考虑 \(n\) 轮的情况。

记 \(X_i\) 表示第 \(i\) 轮的奖金。

则 \(X=X_1+X_2+...+X_n=\mu_n\)。

我们想要知道 \(X/n\) 的相关性质（由于 \(E(X_i)=+\infty\) ，所以不能用辛钦大数定律）

还是考虑计算条件期望 \(E(X/n|X \le M)\)，其中 \(M\) 是使 \(P\{ X \le M\} \ge 1-p\) 的最大值，设 \(p=2^{-m}\)。

考虑如何求 \(M\)。

记 \(Y_i\) 为 \(X_i\) 去掉 \(>M\) 之后部分的新的随机变量。

此时 \(\mu_n = \sum_{i=1}^n Y_i\)
则 \(P\{Y_i=k\}=2^{-k}/(1-2^{-m})\)

\(E(Y_i)=m/(1-2^{-m})\)

\(D(Y_i)=(2^{m+1}-2)/(1-2^{-m})\)

\(n\rightarrow \infty\) 时，对 \(\mu_n\) 成立中心极限定理：

\(P\{\mu_n \le M\}=P\{\frac{\mu_n-nE(Y_i)}{\sqrt{n\times D(Y_i)}} \le \frac{M-nE(Y_i)}{\sqrt{n\times D(Y_i)}} \} = \Phi(\frac{M-nE(Y_i)}{\sqrt{n\times D(Y_i)}})\)

不放令 \(m=30,n=10^7\)，解 \(M\) 如下：

from scipy.stats import norm

n = 10000000
m = 30
p = pow(0.5, m)
Ey = m/(1-p)
Dy = (pow(2, m + 1) - 2) / (1 - p)

x = norm.ppf(1 - p)

M = x * pow(n * Dy, 0.5) + n * Ey

print("M = ", M)

M =  1180628405.2149527

解出 \(M\) 后，\(E(X/n|X \le M)\) 就好算了。

可以近似认为就是 \(log M\)

import math

print(math.log(M)/math.log(2))

30.136907811683777

考虑对比两条曲线：

y_v_2 = []

for i in range(0, MaxN):
    n = i + 1
    M = x * pow(n * Dy, 0.5) + n * Ey
    y_v_2.append(math.log(M)/math.log(2))

plt.plot(x_v, y_v, label = "real")
plt.plot(x_v, y_v_2, label = "estimate", color = 'red')
plt.legend()
plt.show()
    

​

​

不难发现还挺准的。