三门问题（基于贝叶斯分析和python的实现）

三门问题

背景

三门问题（Monty Hall problem）亦称为蒙提霍尔问题，出自美国的电视游戏节目Let’s Make a Deal。问题名字来自该节目的主持人蒙提·霍尔（Monty Hall）。游戏规则如下：

• 现在有三扇门，只有一扇门有汽车，其余两扇门的都是山羊。
• 汽车事前是等可能地被放置于三扇门的其中一扇后面。
• 参赛者在三扇门中挑选一扇。他在挑选前并不知道任意一扇门后面是什麽。
• 主持人知道每扇门后面有什么。
• 如果参赛者挑了一扇有山羊的门，主持人必须挑另一扇有山羊的门。
• 如果参赛者挑了一扇有汽车的门，主持人等可能地在另外两扇有山羊的门中挑一扇门。
• 参赛者会被问是否保持他的原来选择，还是转而选择剩下的那一扇门。

问题是：转换选择可以增加参赛者拿到汽车的机会吗？

这个问题的首次出现，可能是在1889年约瑟夫·贝特朗所著的 Calcul des probabilités 一书中。 在这本书中，这条问题被称为“贝特朗箱子悖论”（Bertrand’s Box Paradox）,也有人把此叫做金盒问题，叙述如下：

现有三个箱子，每个箱子都有两个空档，中间用木板隔开。第一个箱子里面是两块金条，第二个箱子里面是两块银条，第三个箱子里面是一块金条和一块银条。随机抽取一个箱子，然后随机打开一个空档。如果里面是金条，那另外一个空档里是金条的概率是多少？

这是一个有趣的概率学问题，接下来本篇将用两种方法计算保持初始选择不变和改变初始选择下获得车子的概率

贝叶斯分析

针对这个题目，我们假设有A、B、C三个门，参赛者选择了A门，主持人打开了B门，然后要参赛者在A门和C门之间抉择换还是不换。现在我们分两种情况来求解得到汽车的概率：

1.保持选择（参赛者选择A门）
主持人打开门的这个动作对参赛者的选择没有造成任何影响，所以参赛者获得汽车🚗的概率是

$$P(A)=\frac{1}{3}$$

2.改变选择（参赛者选择C门）
在主持人打开了B门，参赛者改变了ta的选择选取了C门，求C门为汽车🚗的概率。
在此情况下，我们可以用贝叶斯公式来求解概率：
$$P(C车|打开B门）=\frac{P(C车）P(打开B门|C车）}{P(打开B门）}$$

• P(C车）:
  $$P(C车）=\frac{1}{3}$$
• P(打开B门|C车）：
  已知C门为车的条件下，主持人只能选择打开B门了，故打开B门的概率为1
• P(打开B门）：
  参赛者选择了A门，主持人打开了B门，$$P(打开B门）=\frac{2}{3}\ast \frac{1}{2}+\frac{1}{3}\ast \frac{1}{2}=\frac{1}{2}$$
  $$P(C车|打开B门）=\frac{\frac{1}{3}\ast 1}{\frac{1}{2}}=\frac{2}{3}$$

因此，参赛者选择了A门，主持人打开了B门，参赛者继续选择A门获得车子的概率为1/3，而改变选择选择C门获得车子的概率为2/3

python模拟

• 固定实验次数（10000）

import random 

strategy = ['stay','switch']

def three_door(strategy,N):
    win = 0
    for i in range(N):
        doors = [1,2,3]
        random.shuffle(doors)
        #假定1为车，2，3为羊
        first_choice = random.choice(doors)
        if first_choice == 1:
            doors.remove(random.choice((2,3)))
        else:
            doors = [1,first_choice]
        if strategy == 'stay':
            second_choice = first_choice
        elif strategy == 'switch':
            doors.remove(first_choice)
            second_choice = doors[0]
        if second_choice == 1:
            win += 1
    pb = win/N
    print('采用'+strategy+'方法，获得车子的概率为：'+str(pb)+'(模拟次数为'+str(N)+')')

three_door('stay',10000)
three_door('switch',10000)

采用stay方法，获得车子的概率为：0.3351(模拟次数为10000)
采用switch方法，获得车子的概率为：0.6642(模拟次数为10000)

• 实验次数与概率的关系图

import random 
import matplotlib.pyplot as plt

#保持选择
def stay(N):
    win = 0
    for i in range(N):
        doors = ['car','goat','goat'] 
        random.shuffle(doors)               #随机doors里面三个物品的排列顺序
        first_choice = random.choice(doors)  #随机选择一个门
        if first_choice == 'car':
            win  += 1
    return win/N               #返回获得车子的概率
#改变选择
def switch(N):     #N为实验次数
    win = 0
    for i in range(N):
        doors = ['car','goat','goat']
        random.shuffle(doors)
        first_choice = random.choice(doors)           
        doors.remove(first_choice)                 #改变选择，所以移除第一次的选择
        doors.remove('goat')                       
        if doors == ['car']:               
            win += 1
    return win/N

x=np.linspace(1,10000,10000)  #设置x变量的范围
y=[]               #设置一个空列表
z=[]
m=[]
n=[]

for i in range(1,10001):
    y.append(switch(i))
for j in range(1,10001):
    z.append(stay(j))
for k in range(1,10001):
    m.append(0.66666666)
    n.append(0.33333333)
plt.xlabel("Numbers of experiment") 
plt.ylabel("Probability of gaining car")
plt.plot(x,y,label='switch')
plt.plot(x,z,label='stay')
plt.plot(x,m,color = 'k',linestyle=':',zorder=2)
plt.plot(x,n,color = 'k',linestyle=':',zorder=2)
plt.grid(linestyle='-.')
plt.legend()
plt.show()        

由上图可以看出，N越大，保持选择的概率越趋向于0.33333，改变选择的概率越趋向于0.6666666

拓展（概率之悖论篇）

三门问题为一个经典的概率悖论，但是概率学上还有许多其他有趣的悖论。

1.生日悖论，也叫生日问题，是指如果一个房间里有23个或23个以上的人，那么至少有两个人的生日相同的概率要大于50%（如下图）；而对于60或者更多的人，这种概率要大于99%。

import random
import numpy as np
import matplotlib.pyplot as plt

#检测birthday中是否有同一天的生日
def same(bir):
    for i in bir:
        num=0
        for j in bir:
            if i==j:
                num=num+1
        if num > 1:
            return 1
    return 0

x=np.linspace(1,100,100)  #设置x变量的范围
y=[]               #设置一个空列表
for i in range(1,101):
    count=0
    for j in range(10000):  #做10000次随机实验
        birthday = []     #设置一个空列表
        for z in range(i):   #每个i产生一个随机数n
            n = random.randint(1, 365)
            birthday.append(n)   #将所有的n放到birthday列表中
        if same(birthday)==1:    #若birthday相同则计一次数
            count += 1
    bit=count/10000             #求概率
    y.append(bit)   
#绘图 
plt.xlabel("Number of people") 
plt.ylabel("Probability of a pair")
plt.plot(x,y,color='r')
plt.text(23, 0.5, (23,0.5), ha='left', rotation=15, wrap=True)
plt.scatter(23, 0.5, marker='o', s=120, zorder=2)
plt.grid(linestyle='-.')
plt.show()

2.假阳性悖论

假设人群中有1%的人罹患某疾病，而其他人是健康的。我们随机选出任一个体，假设检验动作实施在未患病的人身上时，有1%的机率其结果为假阳性；实施在患病的人身上时，有1%的机率其结果为假阴性。

计算可知：如果患病，而被检出为阳性的概率为99%；而如果被测定为阳性者，实际上患病的概率仅为50%。也就是说有一半实际上是假阳性。

数学家John Allen Paulos 在他的《数学盲》一书中指出医生、律师以及其他受过很好教育的非统计学家经常会犯这样的错误。这种错误可以通过用实数而不是概率来描述数据的方法来避免。

3.酒鬼问题

已知某酒鬼有90%的日子都会出去喝酒，喝酒只去固定三家酒吧。今天警察找了其中两家酒吧都没有找到酒鬼。问:酒鬼在第三家酒吧的概率？

可能你会认为依然是90%，其实答案是75%，并很好理解：在警察开始找之前，酒鬼在三家酒吧以及家里的概率分别是30%、30%、30%、10%，而在排除了前两家酒吧以后，酒鬼不在第三家酒吧就在家里。因此，在第三家酒吧的概率为30%/（30%+10%）=75%。

4.本福特定律

本福特定律，也称为本福德法则，说明一堆从实际生活得出的数据中，以1为首位数字的数的出现机率约为总数的三成，接近期望值1/9的3倍。推广来说，越大的数，以它为首几位的数出现的机率就越低。

应用：1972年，Hal Varian提出这个定律来用作检查支持某些公共计划的经济数据有否欺瞒之处。1992年，Mark J. Nigrini便在其博士论文"The Detection of Income Tax Evasion Through an Analysis of Digital Frequencies."（Ph.D. thesis. Cincinnati, OH: University of Cincinnati, 1992.）提出以它检查是否有伪帐。
推而广之，它能用于在会计、金融甚至选举中出现的数据。该定律被华盛顿邮报上的一篇文章引用，该文章以此为基础声称2009年伊朗总统大选中有造假[1]。
若所用的数据有指定数值范围；或不是以概率分布出现的数据，如正态分布的数据；这个定律则不准确。

import matplotlib.pyplot as plt
import math

def firstDigital(x):
    while x >= 10:
        x //= 10
    return x

if __name__ == "__main__":
    n = 1
    frequency = [0] * 9  #记录第一位数字中1-9数字出现次数
    times=[1000] * 9 
    pb = [0] * 9
    for i in range(1,1001):
        n *= i
        m = firstDigital(n) - 1
        frequency[m] += 1
        pb[m] = frequency[m]/times[m]
    print('1-9数字出现的概率为'+str(pb)+'')
    plt.xlabel("number") 
    plt.ylabel("Probability")
    plt.plot(pb,"r-",linewidth=2)
    plt.plot(pb, "go", markersize=8)
    plt.grid(True)
    plt.show()

1-9数字出现的概率为[0.293, 0.176, 0.124, 0.102, 0.069, 0.087, 0.051, 0.051, 0.047]

参考网页

• 三门问题
• 趣味概率之悖论篇
• 生日悖论
• 本福特定律