有一种赌博游戏,投掷一枚硬币,如果正面向上,赌徒将获得参赌赌资的1+a倍,如果反面向上,赌徒将获得参赌赌资的1-b倍,赌局连续进行。现在利用蒙特卡洛模拟方法来说明ab取值对赌局结果的影响。

为方便起见,设置赌徒的初始资金为1 a=0.2。

a=0.2

假设连续进行M次赌博,b取值为0.2,M取值为1000。

def dubo(M,b):
    rec=np.ones(M)
    for dujv in range(1,M):
        if np.random.rand()>0.5:
            rec[dujv]=rec[dujv-1]*(1+a)
        else:
            rec[dujv]=rec[dujv-1]*(1-b)
    return rec
res=dubo(1000,0.2)
def showline(sz):
    plt.figure(figsize=(16,9),frameon=False,dpi=120)
    plt.plot(sz,'.-')
    plt.xlabel('赌局次数')
    plt.ylabel('赌徒持有赌资')
    plt.show()
showline(res)

通过这个赌资随赌局变化图上看出,虽然赌资最大值达到8,赌徒进行1000次赌博之后赌资趋近于0,在大约前400次赌博之后即可认为该赌徒已经彻底输光了。

如果减小b的值,看一下情况会不会有所变化。
取b=0.15

res2=dubo(1000,0.15)
showline(res2)

当b取0.15时,赌徒的赌资迅速增加,等到进行后面几轮时,赌资已经达到了惊人的$10^{6}$数量级,这时候可以认为赌场已经实质性破产了。

从b=0.2时赌徒输光到b=0.15赌场破产,在其中会有一个过渡。现在开始寻找这个临界值。

由于每次赌博的赌资变化有着不确定性,为了能够反映当赌局次数趋向于正无穷的情况,可以进行多次模拟并取平均。为了简化问题,只对每次模拟最终赌资是否大于初始赌资,如果模拟结束时最终赌资大于初始赌资,则赌徒获胜。通过统计不同b取值的情况下赌徒的胜率,找出b的临界点。

为研究问题方便,每个b的取值进行5000次模拟,每次模拟进行1000轮赌局。

Mt=5000
Lt=1000
def dubot(b):
    r=np.zeros_like(b)
    for bdog in range(0,len(b)):
        ui=np.zeros(1)
        bb=b[bdog]
        for cdog in range(0,Mt):
            uii=np.ones(1)
            for ldog in range(0,Lt):
                if np.random.rand()>0.5:
                    uii=uii*(1+a)
                else:
                    uii=uii*(1-bb)
            if (uii>1):
                ui=ui+1
        r[bdog]=ui/Mt
    return r

bbb2=np.linspace(0.1,0.2,21)
rrr2=dubot(bbb2)

可以很轻易的画出赌徒胜率与b的关系图

plt.figure(figsize=(15,9),frameon=False,dpi=120)
plt.plot(bbb2,rrr2,'x-')
plt.xlabel('b')
plt.ylabel('赌徒胜率')
# plt.yticks(np.arange(0,1,0.1))
# plt.xticks(np.arange(0.15,0.2,.005))
plt.show()

从图中可以看出来大约在0.16到0.18之间存在着一个胜率为0.5的临界点,可以继续放大图像查看。

bbb=np.linspace(0.15,0.2,101)
rrr=dubot(bbb)
plt.figure(figsize=(15,9),frameon=False,dpi=120)
plt.plot(bbb,rrr,'.-')
plt.xlabel('b')
plt.ylabel('赌徒胜率')
plt.yticks(np.arange(0,1,0.1))
plt.xticks(np.arange(0.15,0.2,.005))
plt.show()

从图中可以看出随着b的增大,赌徒的胜率逐渐减少以至于趋向于0,当b接近0.167时胜率接近于0.5。

从理论上来说,b=0.167作为临界相变点可以从细致平衡上面来理解。此时b的值满足关系$(1-b)(1+a)=1$即在第一次输钱之后,第二次再赢钱后的结果将回到第一次输钱之前。如果b值大于临界值,在经过足够多的赌博轮数之后可以认为是必输的,至少游戏结束时比初始钱数大的可能性几乎为0。

当b值大于临界值的情况下,寻求进行赌局轮数的最优解。

现在取b=0.19的情况,进行五万次模拟,每次模拟2000轮。

在b等于0.19的情况下,如果仅仅考虑数学期望,这是大于1的,按道理应该是应该一直玩下去的,但实际上虽然期望值大于1
但达到这个结果却并不容易

MM=50000
NN=2000
quu=np.ones([MM,NN])
for ui in range(1,NN):
    quu[:,ui]= quu[:,ui-1]*(np.random.randint(2,size=MM)*.39+0.81)

首先看一下如果进行n次游戏后的赌资的平均值。

ssl2=np.mean(quu,axis=0)
plt.figure(figsize=(15,9),frameon=False,dpi=120)
plt.plot(ssl2,'.-')
plt.xlabel('次数')
plt.ylabel('赌资平均数')
plt.show()

通过对这50000次的赌博进行平均,可以发现平均值都是大于1,完全可以一直玩下去。

看一下中位数,马上就会发现问题,所有这五万次的模拟表明中位数全部小于1,即会有大于一半的可能性会输钱。

ssl=np.median(quu,axis=0)
plt.figure(figsize=(15,9),frameon=False,dpi=120)
plt.plot(ssl,'.-')
plt.xlabel('次数')
plt.ylabel('赌资中位数')
plt.show()

为了更直观的看出分布情况,现在截取第125次、250次、500次、1000次结束后的直方分布图。

plt.figure(figsize=(15,9),frameon=False,dpi=120)
sj=plt.hist(quu[:,125],log=True,bins=np.linspace(0,10,101))
plt.xlabel('赌资')
plt.ylabel('频数')
plt.show()

第125次

plt.figure(figsize=(15,9),frameon=False,dpi=120)
sj=plt.hist(quu[:,250],log=True,bins=np.linspace(0,10,101))
plt.xlabel('赌资')
plt.ylabel('频数')
plt.show()

第250次

plt.figure(figsize=(15,9),frameon=False,dpi=120)
sj=plt.hist(quu[:,500],log=True,bins=np.linspace(0,10,101))
plt.xlabel('赌资')
plt.ylabel('频数')
plt.show()

第500次

plt.figure(figsize=(15,9),frameon=False,dpi=120)
sj=plt.hist(quu[:,1000],log=True,bins=np.linspace(0,10,101))
plt.xlabel('赌资')
plt.ylabel('频数')
plt.show()

第1000次

从中可以明显看出来随着次数的增加,分布会逐渐向0靠拢,即全部赔光的可能越来越多。

到现在可以得出结论,尽管单次赌博的期望值是大于1的,但从统计上来看随着次数的增加,赢钱的可能性越来越小,虽然一旦赢钱所得到的数字会很大以至于将整体的期望值提高到1以上。归根结底原因在于,当b大于临界点的数值情况下,减少一次后需要增加一次以上才能回到初始值,而增加与减少是等概率的,就是这个原因造成了虽然期望值是大于1的,但实际操作起来钱却越来越少。

所以得出的结论是如果b大于临界点的数值,就不要进行连续游戏。如果期望值大于1,则可以进行单次的不会进行连续累计的游戏。