抛硬币出现连续正面的概率计算

Mon, 4th November 2019Edit on GithubcoinprobabilityFibonacci

摘要

mathe于2008年7月引用百度知道中一个抛硬币的概率问题
抛硬币100次,出现10次以上连续正面的概率是多少?
我们将它推广一下,抛硬币n次,其中出现过t次连续正面的概率是多少?
shshsh_0510发现这个问题广义t阶Fibonacci数列 有关系, 其中不出现t次连续正面的概率为 Fn+2(t)2n\frac{F_{n+2}^{(t)}}{2^n}
或者说对于本题,就是要计算1F102(10)21001-\frac{F_{102}^{(10)}}{2^{100}}
最后我们得出广义t阶Fibonacci数列Fn(t)=rt1(1rt)2t(t+1)rtrtnF_n^{(t)}=\lfloor\frac{r_t^{-1}(1-r_t)}{2t-(t+1)r_t} r_t^{n}\rceil, 其中rtr_t是方程 xt+12xt+1=0x^{t+1}-2x^t+1=0 的唯一一个大于1的实数根,x\lfloor x\rceil表示x四舍五入到最接近的整数后的结果。

具体内容

在mathe转载了内容 而shshsh_0510发现了wolfram网站上给出了问题和广义t阶Fibonacci数列的关系以后 mathe先通过矩阵方法验证 了这种关系。

广义Fibonacci数列

我们需要计算硬币抛n次过程中出现t次连续正面的概率。 我们对抛硬币过程中出现的状态进行分类: 其中s0s_0表示还没有出现t次连续正面的情况,而且当前最后一次硬币出现状态为反面或初试状态。 其中s1s_1表示还没有出现t次连续正面的情况,而且当前最后一次硬币出现状态为正面但是最后才连续1次正面 其中s2s_2表示还没有出现t次连续正面的情况,而且当前最后一次硬币出现状态为正面而且最后正好连续2次正面 ... 其中st1s_{t-1}表示还没有出现t次连续正面的情况,而且当前最后一次硬币出现状态为正面而且最后正好连续t-1次正面 其中sts_t表示已经出现t次连续正面的情况。 所以开始抛硬币之前处在状态s0s_0,然后抛硬币过程中,如果当前状态为sks_k,其中k<tk\lt t,那么抛一次硬币有12\frac12的概率 转移到s0s_0,还有12\frac12的概率转移到sk+1s_{k+1};而如果当前状态为sts_t,那么状态永远不会再改变。 所以我们可以得到状态转移矩阵M为 [121212120120000012000001200000121]\begin{bmatrix}\frac12&\frac12&\frac12&\dots&\frac12&0\\ \frac12&0&0&\dots&0&0\\0&\frac12&0&\dots&0&0\\0&0&\frac12&\dots&0&0\\ \dots&\dots&\dots&\dots&\dots&\dots\\ 0&0&0&\dots&\frac12&1\end{bmatrix} 其中我们知道对于初试状态,状态只可能为s0s_0,对应概率分布为b=[10000]b=\begin{bmatrix}1\\ 0\\0\\ \dots\\0\\0\end{bmatrix} 其中b是一个长度为t+1的列向量,b的第h个分量表示状态在sh1s_{h-1}的概率 而经过k步以后,状态分布的概率变化为MkbM^k b 而我们需要的值为经过n步以后到达状态sts_t的概率,也就是MnbM^n b最后一个分量,或者我们记 u=(0,0,0,...,0,1)u=(0,0,0,...,0,1) 那么所求的结果就是an=uMnba_n = u M^n b 为此我们可以先对矩阵M进行分析,由于矩阵2M更加简单,我们改为分析矩阵2M 而矩阵2M的特征多项式很简单,为 g(x)=(2x)(1+x+...+xt1xt)g(x)=(2-x)(1+x+...+x^{t-1}-x^t) =(xt+12xt+1)(2x)x1=\frac{(x^{t+1}-2x^t+1)(2-x)}{x-1} 所以矩阵M的特征多项式为f(x)=g(2x)f(x)=g(2x) 所以我们知道如果计算结果为ana_n, 我们知道数列一个特征根为1,而由于n趋向无穷时an1a_n\to 1, 可以假设an=1+u1(x12)n+u2(x22)n+...+ut(xt2)na_n=1+u_1 (\frac{x_1}2)^n+u_2 (\frac{x_2}2)^n+...+u_t (\frac{x_t}2)^n 其中x1,x2,...,xtx_1,x_2,...,x_t时方程xt+12xt+1x^{t+1}-2x^t+1除去1以后的n个根 我们可以假设bn=2n(1an)b_n=2^n (1-a_n) 那么我们可以知道xt+12xt+1x1\frac{x^{t+1}-2x^t+1}{x-1}是数列bnb_n的特征多项式 也就是说,数列bnb_n满足递推式bn+t+1=2bn+tbnb_{n+t+1}=2b_{n+t}-b_n或者bn+t=bn+t1+bn+t2+...+bnb_{n+t}=b_{n+t-1}+b_{n+t-2}+...+b_n,这个应该就是所谓的广义Fibonacci数列

数值计算结果

有了上面结论,我们就可以比较轻松数值计算题目中要求的结果了。

也就是为了对于给定的t,如果我们要计算对应通项公式,可以通过解方程xt+12xt+1=0x^{t+1}-2x^t+1=0的所有解,然后就可以有通项公式了。
谁来数值计算一下F102(10)F_{102}^{(10)}F1002(10)F_{1002}^{(10)}? 其中F0(10)=0,F1(10)=1,F2(10)=1,F3(10)=2,F4(10)=4,F5(10)=8,F6(10)=16,F7(10)=32,F8(10)=64,F9(10)=128,F10(10)=256,F11(10)=512,F12(10)=1023F_0^{(10)}=0,F_1^{(10)}=1,F_2^{(10)}=1,F_3^{(10)}=2,F_4^{(10)}=4,F_5^{(10)}=8,F_6^{(10)}=16,F_7^{(10)}=32,F_8^{(10)}=64,F_9^{(10)}=128,F_{10}^{(10)}=256,F_{11}^{(10)}=512,F_{12}^{(10)}=1023Fn+11(10)=2Fn+10(10)Fn(10)F_{n+11}^{(10)}=2F_{n+10}^{(10)}-F_n^{(10)} 而抛100次,存在连续正面10次以上的概率为1F102(10)21001-\frac{F_{102}^{(10)}}{2^{100}}

无心人很快给出了数值结果
F102(10)=1211700015849788251502892752696F_{102}^{(10)} = 1211700015849788251502892752696 至于那个分式 = 0.9558627713595783 于是投100次硬币出现连续10次以上的概率为1-0.9558627713595783=0.0441372286404217。

后面mathe继续用后面分析的公式解进行数值计算给出了精度更高的结果
而对于t=10,使用计算器就可以算出rt=1.9990186327101011386634092391292r_t=1.9990186327101011386634092391292 所以b(102)=1211700015849788251502892752685.8,b(1002)=6.5849587983847754109895686668823e+300b(102) =1211700015849788251502892752685.8,b(1002) =6.5849587983847754109895686668823e+300 而投100次不出现连续10个正面的概率为b(102)2100=0.95586277135957831225932140641244\frac{b(102)}{2^{100}}=0.95586277135957831225932140641244, 而投1000次不出现连续10个正面的概率为b(1002)21000=0.61455024758751836408966755676318\frac{b(1002)}{2^{1000}}=0.61455024758751836408966755676318,

公式解分析

不过mathe显然不满足于对前面的数值解,此后花费了很多精力去分析其公式解。由于已经知道其特征多项式为xt+12xt+1x1\frac{x^{t+1}-2x^t+1}{x-1}, 所以其公式解必然和这个多项式的根有关系。
他先做了如下的理论推导

我们现在看看能否给出t阶Fibonacci数列b(n)=Fn(t)b(n)=F_n^{(t)}的通项公式
我们采用初试值b(t+2)=b(t+3)=...=b(0)=0,b(1)=1,b(2)=1b(-t+2)=b(-t+3)=...=b(0)=0,b(1)=1,b(2)=1
为了方便起见,我们可以定义b(t+1)=1b(-t+1)=1,在这样定义以后b(1)的值也可以通过递推公式来递推得到了。
我们定义h(y)=yt+12y+1=k(y)(y1)h(y)=y^{t+1}-2y+1=k(y)(y-1),那么我们知道h(1y)yt+1y1\frac{h(\frac 1y)y^{t+1}}{y-1}是数列的特征方程,所以h(y)h(y)的所有根是数列特征值的倒数。
对于复平面上的圆z=d|z|=d其中23d<1\sqrt{\frac 23}\le d\lt 1,我们有
zt+1=dt+12d12z1|z^{t+1}|=d^{t+1}\le 2d-1\le |2z-1| 所以根据儒歇定理,我们知道h(y)h(y)在圆zd|z|\le d内部解的数目同方程2z+1=0-2z+1=0的数目相同。也就是在z<1|z|\lt 1内方程h(z)=0h(z)=0只有一个解。而如果z=1|z|=1满足h(z)=0h(z)=0,那么容易得到z=1,而且根的重数为1。
由此我们得到h(z)有一个根在单位圆内部,t-1个根在单位圆外部,此外还有一个根1。
假设单位圆内部根为z1z_1,单位圆外部的根为z2,z3,...,ztz_2,z_3,...,z_t,于是我们可以假设数列b(n)b(n)的通项公式为
u1z1n+u2z2n+...+utztnu_1z_1^{-n}+u_2z_2^{-n}+...+u_{t}z_t^{-n}
其中u1,u2,...,utu_1,u_2,...,u_t为待定系数。
采用类似随机游走中的概率问题的方法可以得到其中 us=1k(zs)u_s=\frac 1{k'(z_s)}
而我们知道
k(z)=(zz1)(zz2)(zz3)...(zzt)k(z)=(z-z_1)(z-z_2)(z-z_3)...(z-z_t)
k(zs)=hs(zszh)k'(z_s)=\prod_{h\ne s} (z_s-z_h)
h(zs)=(z1)hs(zszh)h'(z_s)=(z-1)\prod_{h\ne s} (z_s-z_h)
所以k(zs)=h(zs)zs1k'(z-s)=\frac{h'(z_s)}{z_s-1}
us=zs1h(zs)u_s=\frac{z_s-1}{h'(z_s)}
h(zs)=(t+1)zst2=(t+1)(2zs1)2=2t(t+1)zs1h'(z_s)=(t+1)z_s^t-2=(t+1)(2-z_s^{-1})-2=2t-(t+1)z_s^{-1}
所以
us=zs12t(t+1)zs1u_s=\frac{z_s-1}{2t-(t+1)z_s^{-1}}

rt=z11r_t=z_1^{-1}
而递推式中其他各式绝对值都远远小于,所以我们得到,对于n充分大
b(n)=rt1(1rt)2t(t+1)rtrtnb(n)=\lfloor\frac{r_t^{-1}(1-r_t)}{2t-(t+1)r_t} r_t^{n}\rceil
而他怀疑这个表达式对所有的n都成立,不过这是还没有验证,但是如果仅仅数值计算,不需要精确值,由于r|r|接近2,而其他项都递减,所以用上面公式精度已经非常高了。

继续对b(n)=rt1(1rt)2t(t+1)rtrtnb(n)=\frac{r_t^{-1}(1-r_t)}{2t-(t+1)r_t} r_t^{n}做误差分析可以得出
除了上面一项以外,实际上还有t-1项uszsn,s=2,3,...,tu_s z_s^{-n}, s=2,3,...,t.
其中对于n1n\ge 1,有uszsnuszs1+zs12t(t+1)<2t1|u_s z_s^{-n}|\le |u_s*z_s|\le \frac{1+|z_s^{-1}|}{|2t-(t+1)|}\lt \frac2{t-1} 所以t-1项之和的绝对值小于2(t1)t1=2\frac{2(t-1)}{t-1}=2
也就是说使用这个公式b(n)b(n)的误差必然小于2. 而则换成计算概率,第n项的计算误差小于12n1\frac1{2^{n-1}},对于比较大的n仅使用第一项计算精度非常高,实际计算过程中可能r采用的精度都没有这么高。

后面他使用数值计算进行验算发现:

而对于公式b(n)=u1rtnb(n)=\lfloor u_1 r_t^n\rceil,计算前面若干项我们可以发现猜想均成立:

n b(n) u1rtnu_1r_t^n
1 1 0.50222005921650437539314900083299
2 1 1.0039472560945626042296717505261
3 2 2.0069092711912102894563100627129
4 4 4.0118490272698787619203081600742
5 8 8.0197609571323822998698197638132
6 16 16.031651583188626895454838284701
7 32 32.047570227910457178387829251568
8 64 64.063690018678506437470207581223
9 128 128.06451002750246161567818042825
10 256 256.00334173386700758850659443979
11 512 511.75545016205159804636690872342
12 1023 1023.0086802648866917173406684459
13 2045 2045.0134132736788208304516651412
14 4088 4088.0199172761664313714470202193
15 8172 8172.0279855250629839809737322779

然后他发现这个问题中特征方程xt+12xt+1x^{t+1}-2x^t+1随机游走中的概率问题中进t退1的特殊情情况,最后可以得出我们使用上面的近似公式那么b(n)的误差为q(n)q(n+1)q(n)-q(n+1),其中q(n)是那个问题中5#中定义的当A在B后面n步时能够返回原点的概率.所以我们只要能够证明0<q(n)<120<q(n)<\frac12那么就证明了只使用一项然后四舍五入的方法是正确的。而且也可以解释上面例子中除了第一项以外,后面所有项的误差都是正的情况(因为q(n)应该是单调的) 而且显然对于固定的n,q(n)关于t是递减的,所以如果我们能够对于t=2的情况计算出q(1)12q(1)\le\frac12,那么对于所有这一类题目,都能够直接使用那个近似计算方法,然后四舍五入,而结果就会是精确的。

以此为基础,他又开始对上面猜测正确性的分析
为了证明猜想公式的正确性,即b(n)=r1(r1)(t+1)r2trnb(n)=\lfloor\frac{r^{-1}(r-1)}{(t+1)r-2t}r^n\rceil
我们先证明
r1(r1)(t+1)r2t<12\frac{r^{-1}(r-1)}{(t+1)r-2t}\lt\frac12
其中1<r<21\lt r\lt 2而且rt+12rt+1=0r^{t+1}-2r^t+1=0,
首先我们容易得出(t+1)r2t>0(t+1)r-2t>0
这个是因为方程h(x)=xt+12xt+1h(x)=x^{t+1}-2x^t+1中,显然h(2tt+1)<0,h(2)>0h(\frac{2t}{t+1})\lt 0,h(2)\gt 0,而且前面已经证明r是h(x)唯一一个绝对值大于1的解,所以2tt+1<r<2\frac {2t}{t+1}\lt r\lt 2
所以我们需要证明2(1r1)<(t+1)r2t2(1-r^{-1})\lt(t+1)r-2t也就是(2r)r<2t+1(2-r)r\lt \frac2{t+1}
由于我们知道rt(2r)=1r^t(2-r)=1,所以只要证明rt1>t+12r^{t-1}\gt\frac{t+1}2也就是r>t+12t1r\gt \sqrt[t-1]{\frac{t+1}2}(其中t2t\ge2
而其中右边容易得出在t=2t=2时有最大值1.5,而r>2tt+1r\gt\frac{2t}{t+1},所以在t4t\ge 4时都已经有不等式成立。
而对于t=3时,我们知道r3=1.83>1.5r_3=1.83\dots\gt 1.5,对于t=2,r2=1.618>1.5r_2=1.618\dots\gt 1.5
所以我们证明了r1(r1)(t+1)r2t<12\frac{r^{-1}(r-1)}{(t+1)r-2t}\lt\frac 12
而更加容易的,我们可以证明r1(t+1)r2t>12\frac{r-1}{(t+1)r-2t}\gt\frac12 (这个等价于r<2r\lt 2),
而又有前面的不等式得到r1(t+1)r2t<r2<1\frac{r-1}{(t+1)r-2t}\lt\frac r2\lt1

现在我们定义序列e(n)=b(n)r1(1r)2t(t+1)rrne(n)=b(n)-\frac{r^{-1}(1-r)}{2t-(t+1)r}r^n也就是这个误差,其中对于所有的nt+2n\ge-t+2全部定义。
而其中b(n)b(n)n0n\le 0时定义为b(t+2)=b(t+3)=...=b(0)=0b(-t+2)=b(-t+3)=...=b(0)=0(这样所有这些项都可以用同一个递推公式推导了)
显然e(n)e(n)也满足递推公式e(n+t)=e(n+t1)+e(n+t2)+...+e(n)e(n+t)=e(n+t-1)+e(n+t-2)+...+e(n),这个递推式对于所有的nt+2n\ge -t+2都成立。
而且我们知道e(0)=b(0)r1(1r)2t(t+1)r=r1(1r)2t(t+1)re(0)=b(0)-\frac{r^{-1}(1-r)}{2t-(t+1)r}=-\frac{r^{-1}(1-r)}{2t-(t+1)r},所以12<e(0)<0-\frac12\lt e(0)\lt 0
同样,我们知道12<e(1)<0,12<e(2)<0,,12<e(t+2)<0-\frac 12\lt e(-1)\lt 0,-\frac 12\lt e(-2)\lt 0,\dots,-\frac 12\lt e(-t+2)\lt 0
e(1)=1r1(t+1)r2te(1)=1-\frac{r-1}{(t+1)r-2t},所以0<e(1)<120\lt e(1)\lt \frac12

我们现在定义一个模型,有一个人在实数轴的整点上行走(开始坐标大于1),他每次分别以均等的概率随机后退t格或前进1格,知道到达坐标位置[t+2,1][-t+2,1]之间停止;如果他最后停留在坐标x,那么我们说他最终的得分为e(x)e(x),那么请问如果他的起始坐标为x,那么他最终得分的期望值是多少呢?(容易证明他最终会停止下来的概率为1) 非常有意思,对于任何一个位置n,他最终得分的期望值为e(n)e(n),由于每个最终得分都在12-\frac1212\frac12之间,我们可以知道所有的e(n)e(n)也必然只能在12-\frac1212\frac12之间。也就是这种方法可以证明上面的猜想。

这个是因为假设开始在位置n的得分为f(n)f(n),我们知道对于任意n2n\ge 2f(n)=f(n+1)+f(nt)2f(n)=\frac{f(n+1)+f(n-t)}2或者说f(n+1)=2f(n)f(nt)f(n+1)=2f(n)-f(n-t),而且对于n1n\le 1有边界条件f(n)=e(n)f(n)=e(n) 根据f(n)f(n)的递推公式我们知道它对应的特征方程为xt+12xt+1x^{t+1}-2x^t+1,这个方程有一个根的绝对值大于1,还有一个根为1,其余t-1个根绝对值小于1,根据 随机游走中的概率问题中的结论,由于数列f(n)f(n)有界,那么绝对值大于1的根所对应的项的系数必然为0。也就是f(n)f(n)的通项公式可以由其余t个特征根决定,也就是前t项已经唯一决定了数列f(n)f(n),而由于e(n)e(n)也满足同样递推式,必然得出f(n)=e(n)f(n)=e(n)

由此我们可以证明了前面猜想的正确性,即
b(n)=r1(r1)(t+1)r2trnb(n)=\lfloor\frac{r^{-1}(r-1)}{(t+1)r-2t}r^n\rceil
而抛硬币n次出现连续t次正面的概率为
p(n)=1r1(t+1)r2trn+12n;r>1;rt+12rt+1=0p(n)=1-\frac{\lfloor\frac{r-1}{(t+1)r-2t}r^{n+1}\rceil}{2^n}; r\gt 1; r^{t+1}-2r^t+1=0

其它应用

广义Fibonacci数列还会以其它一些形式出现,比如倪举鹏提出过的寿命10年的兔子问题:
第一年有1对兔子,以后每一年活着的兔子都能繁殖出一倍兔子来(活着的兔子乘以2)。兔子寿命都是10年。第n年有多少对兔子.
实际上计算的就是9阶广义Fibonacci数列

Github