抛硬币1000次，至少连续10次正面朝上的概率 详细解答

1. 起因

看论文的时候，论文里简单提了一下。大概意思是，抛硬币1000次，至少连续10次正面朝上的概率比较大。我无聊就算了一下（后来就想拍死这个无聊的自己T^T）。

2. 问题陈述

一开始没什么思路，查了很多资料，答案五花八门。但是有一个观点让我很惊喜，他强调了一下问题本身。
这个问题到底什么意思？

1. 角度1
   抛硬币1000次，如果至少连续10次正面朝上，则称事件A发生。（连续11次、12次，或者有好多次“连续10次”都不管，只要至少有一次连续10次就看作事件A发生）
2. 角度2
   抛硬币1000次，共有$2^{1000}$种情况。可以看作一个基数为$2^{1000}$的集合，记为$S$。集合$S$中的每一个元素都是长度为1000，仅包含字符0，1的字符串。集合$A\subset S$，如果“1111111111”（10个1）是某一字符串的子串，那么该字符串在集合$A$中。实际上集合$S$就是样本空间，集合$A$就是随机事件A。
   我们想知道$|A|$。

3. 问题简化版本

好了，现在我们对问题已经比较清楚了。先看一个简化版本，这样有助于思考：连续抛10次硬币，至少连续两次正面朝上的概率是多少？
但好像不太好算欸：

1. 思考1 分解事件A

事件$A$=事件$B_2$(恰好连续两次正面朝上)+事件$B_3$(恰好连续三次正面朝上)+···+事件$B_{10}$(恰好连续十次正面朝上)

• 事件$B_2$=事件$B_{2,1}$(恰好一次连续两次正面朝上)+事件$B_{2,2}$(恰好两次连续两次正面朝上)+事件$B_{2,3}$(恰好三次连续两次正面朝上)（四次是不可能的）
• 事件$B_3$=事件$B_{3,1}$(恰好一次连续三次正面朝上)+···
• ···
  （我选择放弃）

2. 思考2 反事件
   $P(A)=1-P(\overline{A})$ 但事件$\overline{A}$也很难求

3.1 遇事不决 先Brutal一下

$2^{10}$好像也挺大的（但比$2^{1000}$小得多），就直接暴力搜索一下吧：

# 暴力法直接搜索
# 投掷硬币m次 其中n次连续正面朝上的情况种数/概率
def brutal(m=10, n=2):
	# 表示连续n次的子串，也就是关键字
    s = '1'*n
    count = 0
    # 直接枚举样本空间
    for i in range(2**m):
        if s in bin(i)[2:]:
            count += 1
    # 返回符合条件的种数，以及概率P
    return count, count/(2**m)

结果出来了，count=880, p=0.859375。这个结果可能还不够清晰，我们再看个更简单的例子：只抛四次硬币。
这时就比较直观了：(O表示正面，X表示反面)

-	-	-	-
XXXX	XOXX	OXXX	OOXX
XXXO	XOXO	OXXO	OOXO
XXOX	XOOX	OXOX	OOOX
XXOO	XOOO	OXOO	OOOO

发现brutal函数在这个输入下的输出是正确的

3.2 换个思路

我们先定义一个数列 $A_n^k$ 表示投掷硬币n次，连续k次正面朝上的情况种数（数列下标还是喜欢n啊）
现在我们想求$A_{10}^2$：

要抛10次硬币，第一次要么正面，要么反面：（@表示任意）
O@@@@@@@@@ (情况1)
X@@@@@@@@@ (情况2)

• 对于(情况1)，再考虑第二次抛硬币：
  OO@@@@@@@@@ (情况1.1)
  OX@@@@@@@@@ (情况1.2)
  • 对于(情况1.1)，已经满足条件了，剩下8次就随便吧 $2^8$
  • 对于(情况1.2)，既然已经抛出反面了，又得重新开始了。$A_8^2$
• 对于(情况2)，不就是投掷硬币9次的情况种数吗？$A_9^2$

以上，我们可以得出$A_{10}^2=2^8+A_8^2+A_9^2$
推广一下：$A_n^2=2^{n-2}+A_{n-1}^2+A_{n-2}^2$
再推广一下：$A_n^k=2^{n-k}+{\sum }_{i=1}^k{A}_{n-i}^k(1)$

用本节的推导方式可以验证$k=10$时，公式(1)是正确的

4. 用递推公式算算通项吧

4.1 线性代数解法

（这里巨坑，算了好久还算错了。最后用线性代数的解法才解出来）
递推公式：（隐含$k=2$）
$$A_n=2^{n-2}+A_{n-1}+A_{n-2}$$配一下：
$$A_n-2^n=A_{n-1}-2^{n-1}+A_{n-2}-2^{n-2}$$令$B_n=A_n-2^n$，则特征方程：
$$\begin{gathered} x^2=x+1 \\ {x}_1=\frac{1+\sqrt{5}}{2},x_2=\frac{1-\sqrt{5}}{2} \end{gathered}$$又$B_1=A_1-2=-2,B_2=A_2-2^2=-3$ ($A_1,A_2$的值应该比较显然)，则$$\{\begin{array}{l}{c}_1{x}_1+c_2{x}_2=-2\\ c_1{x}_1^2+c_2{x}_2^2=-3\end{array}$$解得$$\{\begin{array}{l}{c}_1=-\frac{1}{2}(\frac{3}{\sqrt{5}}+1)\\ c_2=\frac{1}{2}(\frac{3}{\sqrt{5}}-1)\end{array}$$即$B_n=-\frac{1}{2}(\frac{3}{\sqrt{5}}+1)(\frac{1+\sqrt{5}}{2}{)}^n+\frac{1}{2}(\frac{3}{\sqrt{5}}-1)(\frac{1-\sqrt{5}}{2}{)}^n$
则$A_n=2^n-\frac{1}{2}(\frac{3}{\sqrt{5}}+1)(\frac{1+\sqrt{5}}{2}{)}^n+\frac{1}{2}(\frac{3}{\sqrt{5}}-1)(\frac{1-\sqrt{5}}{2}{)}^n(2)$

4.2 Mathematica法

命令如下：
RSolve[{a[n] == 2^(n - 2) + a[n - 1] + a[n - 2], a[1] == 0, a[2] == 1}, a[n], n]

将斐波那契数列和卢卡斯数列展开就可以得到公式(2)
（第一次碰到卢卡斯数列，刚开始一脸懵）
检查数列的项，来进行验证。

5. 回到原问题

递归计算

公式(1)已经给出了我们所需要的通项公式，但如果按第4节的解法，是解不出来的。线性代数法里会出现一个一元十次方程（菜鸡不会解这个特征方程），Mathematica则差点把我电脑烧坏（谨慎尝试，尽快放弃计算）
这个通项公式的确很难解，但我们可以绕开它：

# 递推公式法
# 投掷硬币m次 其中n次连续正面朝上的情况种数/概率
def recurrence(m=10, n=2):
	# 通项公式的最开始几项
    if m < n:
        return 0
    if m == n:
        return 1
    l = [0]*n
    l[-1] = 1
    multi = 1
    # 不断向前递推
    for _ in range(n, m):
        multi *= 2
        tmp = sum(l)+multi
        del l[0]
        l.append(tmp)
    return l[-1], l[-1]/(2**m)

这个算法，比3.1节中的算法快很多 $O(2^n)\to O(n)$

结果

$A_{1000}^{10}$大约$10^{301}$，这里就不列出来了。P=0.38544975241248164。
我们还是验证一下吧

1. 在n很小的时候，其实很好算的（把10次正面朝上看作一个整体），然后验证recurrence函数的输出是否正确
2. 更有说服力的验证，就是随机事件模拟了

# 随机事件模拟
# 投掷硬币m次 其中n次连续正面朝上的概率
def simulate(m=10, n=2, size=10000000):
    count = 0
    for i in range(size):
    	# 一次模拟出m次投掷的结果
        sample = np.random.randint(0, 2, (m))
        for j in range(m-n):
        	# 滑动长度为n的窗口 判断是否满足条件
            if np.sum(sample[j:j+n]) == n:
                count += 1
                break
        # 进度条君
        if (i+1) % (size//100) == 0:
            print('\r %d%%' % ((i+1)//(size//100)), end='')
    print()
    print(count/size)
    return count/size

在size=1000000时，结果为0.384777，已经比较接近了（注意：费时操作）

6. 尾声

虽然没能得到最终的通项公式，但还是较完美地解决了这个问题。概率论真的太难了。当时我算了几个小时，还算错了，随意用Mathematica跑了一下，竟然算得十分漂亮。辛辛苦苦，还输给了一个数学软件T^T.