St. Petersburg Paradox II

终于开始写后续了… 强化学习开头讲了个老虎机问题,又把这个陈年老题给想起来…

先简单回顾下问题

圣彼得堡悖论: 掷硬币,若第一次掷出正面,你就赚2元.若第一次掷出反面,那就要再掷一次,若第二次掷的是正面,你便赚4元.若第二次掷出反面,那就要掷第三次,若第三次掷的是正面,你便赚8元,如此类推,即可能掷一次游戏便结束,也可能反复掷没完没了.问题是:你最多肯付多少钱参加这个游戏?

Solution

\[\begin{eqnarray*} E_{now}(1) &=& 2 \\ E_{now}(N) &=& E_{now}(N-1) + 1 \\ E_{win}(1) &=& 1 \\ E_{win}(N) &=& E_{now}(N-1)/2 + E_{win}(N-1) + 1/2 \\ \end{eqnarray*}\]

解得:

\[\begin{eqnarray*} E_{now}(N) &=& n + 1 \\ E_{win}(N) &=& (n - 1)(n + 4) / 4 + 1 \\ \end{eqnarray*}\]

Updates

上一篇post已经解出了投n次硬币之后获益的期望,假设每开一局需要消耗c,那么我们要解$\arg_x \min E_{win}(x) > c(x+1)/2$,获得最少的投掷次数已使获利期望为正.

可惜之后写了个模拟程序,发现大多数情况下,投掷x次获益为负.这是因为整个获益的分布方差太大导致的.理论上来讲,这时候只需要求一下中位数就好了,这样50%的概率获利.

然而递推公式并不能求得中位数,目前想到的方法依然需要把整个分布展开.也就是需要2^n项全部计算.

既然不能计算median,那么另一个方案是计算几何平均,对获益取log,然后算数平均,然后再exp回来.但是递推依然不支持,所以还没有验证.

还有一个是算方差,然后套Chebyshev’s Inequality, \(P(|X-E(X)| \ge b) \le \frac{Var(X)}{b^2}\). 但我怀疑这个界十分松散,根本不管用,毕竟Var(X)轻松彪上十几位数字.

纯理论上的方法既然不管用,那就来点工程暴力学了. 下面一个简单的思路是,通过某种方式把极大的获利在统计时去掉.

比如我要投掷100次,那么一百次全为反面,会导致当前收益为$2^{100}$,虽然平均之后他只为平均值贡献了1…坑爹…

如何截断统计100次全为反面的统计,比较简单的一个方法方式就是,规定这一百次投掷中,第50次投掷强制是正面,这样就把51~100次连续反面的可能性全部截断了. 事实证明这个方式非常有效,直接将平均值减半了.

比如平均值是2000,那么截断之后平均值变成了1000.简单讲50%的概率1000,50%的概率3000.

Implementation

注意到这个递推公式里面有两个常量 $E_{now}(1)$和$E_{win}(1)$. 其实我们可以把他们扩展成变量,即从任意中间状态推导至$k$次投币之后. 改进之后的公式:

\[\begin{eqnarray*} E_{now}(0) &=& n \\ E_{now}(k) &=& E_{now}(k-1) + 1 \\ E_{win}(0) &=& w \\ E_{win}(k) &=& E_{now}(k-1)/2 + E_{win}(k-1) + 1/2 \\ \end{eqnarray*}\]

解得:

\[\begin{eqnarray*} E_{now}(k) &=& k + n \\ E_{win}(k) &=& k(k + 1 + 2n) / 4 + w \\ \end{eqnarray*}\]

这个扩展之后,不禁方便我们加入一个强制,而且可以加入任意多位置的强制

 
def solve(now, win, k):
    return k + now, win + k * (k + 1 + 2 * now) / 4

def solve1(now, win):
    return now + 1, win + (1 + now) / 2

def solveN(n):
    return solve(1, 0, n)

def win(now, win):
    return 2 * now, win

def lose(now, win):
    return 2, win + now + 1

def mustWinSeq(n, seq):
    Enow = 1.0
    Ewin = 0.0
    cur = 0
    for i in seq:
        if i - 1 > cur:
            Enow, Ewin = solve(Enow, Ewin, i - 1 - cur)
        Enow, Ewin = win(Enow, Ewin)
        cur = i
    Enow, Ewin = solve(Enow, Ewin, n - cur)
    return Enow, Ewin

def mustLoseSeq(n, seq):
    Enow = 1.0
    Ewin = 0.0
    cur = 0
    for i in seq:
        if i - 1 > cur:
            Enow, Ewin = solve(Enow, Ewin, i - 1 - cur)
        Enow, Ewin = lose(Enow, Ewin)
        cur = i
    Enow, Ewin = solve(Enow, Ewin, n - cur)
    return Enow, Ewin

def measureImprove(n, seq):
    retNow, retWin = solveN(n)
    retNow, retWin = lose(retNow, retWin)
    
    Enow, Ewin = mustWinSeq(n, seq)
    Enow, Ewin = lose(Enow, Ewin)
    print("solveN:", retWin)
    print("mustWinSeq:", Ewin)

    EnowL, EwinL = mustLoseSeq(n, seq)
    EnowL, EwinL = lose(EnowL, EwinL)
    print("mustLoseSeq:", EwinL)
    
    p = 0.5**len(seq)
    improve = retWin - (retWin - p*Ewin) / (1-p)
    retWin -= improve
    return retWin, improve, 1-p

    
measureImprove(1000,[500])
# solveN: 251752.0
# mustWinSeq: 376750.5
# mustLoseSeq: 126753.5
# (126753.5, 124998.5, 0.5)

可以看到投掷1000次,获益期望本来是251752, 现在是1/2的概率获益126753.5, 1/2的概率获益376750.5