有一种赌博游戏，投掷一枚硬币，如果正面向上，赌徒将获得参赌赌资的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，则可以进行单次的不会进行连续累计的游戏。