本文通过贝叶斯方法分析抛币赌局，使用Python进行模拟。作业中涉及三种先验：一无所知、硬币均匀、硬币有偏。通过对100次和1000次投掷硬币的数据分析，计算后验密度、赢钱概率和边缘似然。结果表明，随着投掷次数增加，后验密度更加集中，赢钱概率更接近真实概率，边缘似然帮助选择合适的先验。
【机器学习作业分享3】贝叶斯方法体验
介绍
本篇文章是第三次机器学习作业的求解分享。在课上学习了贝叶斯方法，本次作业主要是运用这一方法，对抛币赌局进行相关计算。使用的语言为python，编辑器为Jupyter Notebook。批评指正，互相交流。
抛币赌局
抛币赌局是世纪之交街边一种常见娱乐活动。与今天的赌博原理本质上一致。玩家支付1元作为本金参加游戏，如果玩家赢得游戏，则庄家不仅返还本金，还要支付额外一元作为奖金给玩家。想法，如果玩家输掉游戏，则庄家不返还本金。玩法如下：
连续抛掷10次硬币（或者同时抛掷10枚硬币），若有六面/次以上（不含六面/次）硬币朝上，则视为庄家赢，相反则玩家获胜。
这个赌局的关键是庄家的硬币是否是一个均匀的硬币，即其正面朝上的概率是否为0.5。
作业要求
① 以
r
=
0.7
r= 0.7
r=0.7的二项分布模拟抛币赌局，分别生成100次和1000次数据。
② 定义三种先验分布，一无所知、硬币均匀、硬币有偏
③ 以①生成的数据作为输入，分别计算②三种先验条件下的后验密度、赢钱概率和边缘似然。
注意：①赢钱概率包括用期望值点估计计算出的概率 和用参数后验期望计算出的概率
② 三种先验的控形参数可以直接引用，也可自己定义
作业求解
首先，将必要的函数库导入
# 函数库的导入与设置
import numpy as np
import math
import scipy.stats as sp
import matplotlib.pyplot as plt
from scipy.special import gamma
plt.rcParams['font.sans-serif']=['SimHei']
plt.rcParams['axes.unicode_minus'] = False
本次实验的理论基础为贝叶斯定理：
p
(
r
∣
y
N
)
=
P
(
y
N
∣
r
)
p
(
r
)
P
(
y
N
)
p(r|y_N) = \frac{P(y_N|r)p(r)}{P(y_N)}
p(r∣yN)=P(yN)P(yN∣r)p(r)
其中
P
(
y
N
∣
r
)
P(y_N|r)
P(yN∣r)为似然值
p
(
r
)
p(r)
p(r)为先验概率密度
P
(
y
N
)
P(y_N)
P(yN)为边缘密度
任务一：生成数据
首先模拟抛币100次和1000次。在之后我们将用这两次模拟的结果来作为 r r r的估计值
# 分别生成100次和1000次数据
g_data_100 = np.random.binomial(100,0.7,1) # 100个硬币同时扔
print("扔100次硬币，得到r的期望值为：" + str(g_data_100[0]/100))
g_data_1000 = np.random.binomial(1000,0.7,1)
print("扔1000次硬币，得到r的期望值为：" + str(g_data_1000[0]/1000))
g_N = [100, 1000]
g_y_N = [g_data_100[0], g_data_1000[0]]
g_r_hat =[g_data_100[0]/100, g_data_1000[0]/1000]
生成 r r r的离散数据，主要用于后文的绘图。
# 生成r的离散数据，用于绘图
g_r_vals = np.linspace(0,1,100)[:,None]
#g_r_vals
任务二：定义三种先验
三种先验分别为：一无所知、硬币均匀和硬币有偏。先验的表达式如下所示：
p
(
r
)
=
Γ
(
α
+
β
)
Γ
(
α
)
Γ
(
β
)
r
α
−
1
(
1
−
r
)
β
−
1
p(r) = \frac{\Gamma(\alpha+\beta)}{\Gamma(\alpha)\Gamma(\beta)}r^{\alpha-1}(1-r)^{\beta-1}
p(r)=Γ(α)Γ(β)Γ(α+β)rα−1(1−r)β−1
这三种情况分别是通过改变控形参数
α
\alpha
α和
β
\beta
β来实现的
我们定义这三种情况的
α
\alpha
α
β
\beta
β分别为：
什么也不知道：
α
=
1
\alpha=1
α=1，
β
=
1
\beta=1
β=1
硬币均匀：
α
=
50
\alpha=50
α=50，
β
=
50
\beta=50
β=50
硬币有偏：
α
=
5
\alpha=5
α=5，
β
=
1
\beta=1
β=1
# 定义α和β
g_a = [1, 50, 5]
g_b = [1, 50, 1]
三种先验下， r r r的概率密度展示：
# 生成所有先验的数据
g_p_r = [] # 先验
for i in range(len(g_a)): # 不同种情况
# g_p_r_i = []
g_p_r.append(sp.beta.pdf(g_r_vals,g_a[i],g_b[i])) # β函数
plt.plot(g_r_vals,g_p_r[0],label="一无所知")
plt.plot(g_r_vals,g_p_r[1],label="硬币均匀")
plt.plot(g_r_vals,g_p_r[2],label="硬币有偏")
plt.legend()
绘制出三种先验的概率密度图，如下：
任务三：计算后验密度、赢钱概率和边缘似然。
首先，来计算后验密度
经过推导，后验密度为：
p
(
r
∣
y
N
)
=
Γ
(
α
+
β
+
N
)
Γ
(
α
+
y
N
)
Γ
(
β
+
N
−
y
N
)
r
α
+
y
N
−
1
(
1
−
r
)
β
+
N
−
y
N
−
1
p(r|y_N) = \frac{\Gamma(\alpha+\beta+N)}{\Gamma(\alpha+y_N)\Gamma(\beta+N-y_N)}r^{\alpha+y_{N}-1}(1-r)^{\beta+N-y_{N}-1}
p(r∣yN)=Γ(α+yN)Γ(β+N−yN)Γ(α+β+N)rα+yN−1(1−r)β+N−yN−1
其中，不妨令：
δ
=
y
N
+
α
\delta = y_N + \alpha
δ=yN+α,
γ
=
N
−
y
N
+
β
\gamma = N - y_N + \beta
γ=N−yN+β
下面计算，在两种投币次数条件下，三种先验所对应的不同的后验值
g_p_r_after = [] # 所有的后验
for i in range(len(g_N)):
g_p_r_after_i = []
for j in range(len(g_a)):
g_p_r_after_i.append(sp.beta.pdf(g_r_vals,g_a[j] +g_y_N[i] ,g_b[j] + g_N[i] - g_y_N[i]))
g_p_r_after.append(g_p_r_after_i)
选择颜色
b = plt.colormaps['Accent']( ### 这个地方可以研究colormaps来调整颜色变化
np.linspace(0.15, 0.85, 10))
color1 = plt.colormaps['Pastel1']( ### 这个地方可以研究colormaps来调整颜色变化
np.linspace(0.1, 0.99, 10))
colors=['#1f77b4','#ff7f0e','#2ca02c','#d62728','#9467bd','#8c564b','#e377c2','#7f7f7f','#bcbd22','#17becf']
绘制先验与后验的对比图
fig, axes = plt.subplots(1,2, figsize = (20, 8))
axes[0].set_title("投掷100次硬币的后验密度变化")
axes[0].plot(g_r_vals,g_p_r[0],color=colors[0],marker=" ",linestyle="-",linewidth=3,alpha=0.2,
label="投掷前 AND 一无所知")
axes[0].plot(g_r_vals,g_p_r[1],color=colors[1],marker=" ",linestyle="-",linewidth=3,alpha=0.2,
label="投掷前 AND 硬币均匀")
axes[0].plot(g_r_vals,g_p_r[2],color=colors[2],marker=" ",linestyle="-",linewidth=3,alpha=0.2,
label="投掷前 AND 硬币有偏")
axes[0].plot(g_r_vals,g_p_r_after[0][0],color=colors[0],marker=" ",linestyle="-",linewidth=3,alpha=1,
label="投掷后 AND 一无所知")
axes[0].plot(g_r_vals,g_p_r_after[0][1],color=colors[1],marker=" ",linestyle="-",linewidth=3,alpha=1,
label="投掷后 AND 硬币均")
axes[0].plot(g_r_vals,g_p_r_after[0][2],color=colors[2],marker=" ",linestyle="-",linewidth=3,alpha=1,
label="投掷后 AND 硬币有偏")
axes[0].legend()
axes[0].grid()
###########################################################################################################
axes[1].set_title("投掷1000次硬币的后验密度变化")
axes[1].plot(g_r_vals,g_p_r[0],color=colors[0],marker=" ",linestyle="-",linewidth=3,alpha=0.2,
label="投掷前 AND 一无所知")
axes[1].plot(g_r_vals,g_p_r[1],color=colors[1],marker=" ",linestyle="-",linewidth=3,alpha=0.2,
label="投掷前 AND 硬币均匀")
axes[1].plot(g_r_vals,g_p_r[2],color=colors[2],marker=" ",linestyle="-",linewidth=3,alpha=0.2,
label="投掷前 AND 硬币有偏")
axes[1].plot(g_r_vals,g_p_r_after[1][0],color=colors[0],marker=" ",linestyle="-",linewidth=7,alpha=1,
label="投掷后 AND 一无所知")
axes[1].plot(g_r_vals,g_p_r_after[1][1],color=colors[1],marker=" ",linestyle="-",linewidth=3,alpha=1,
label="投掷后 AND 硬币均匀")
axes[1].plot(g_r_vals,g_p_r_after[1][2],color=colors[2],marker=" ",linestyle="-",linewidth=2,alpha=1,
label="投掷后 AND 硬币有偏")
axes[1].legend()
axes[1].grid()
上图分别展示了100次和1000次投掷硬币后后验密度的变化。可以看到，无论是哪一种先验，在进行很多次实验之后，根据观测值的变化，后验都会趋于均值在0.7左右的正态分布。
对比左右两张图，可以看到，投掷的次数越多，对于先验的改变就越大。同时，在先验中，控形参数越大，就越不容易被改变
带入一开始100次投币和1000投币的 r r r的期望值，计算后验密度得：
print("对硬币一无所知条件下，投掷100次之后得" + "r ="+str(g_r_hat[0]) + "的后验密度为：" + str(sp.beta.pdf(g_r_hat[0], g_a[0] +g_y_N[0] ,g_b[0] + g_N[0] - g_y_N[0])))
print("硬币均匀条件下，投掷100次之后得" + "r="+str(g_r_hat[0])+"的后验密度为：" + str(sp.beta.pdf(g_r_hat[0], g_a[1] +g_y_N[0] ,g_b[1] + g_N[0] - g_y_N[0])))
print("硬币有偏条件下，投掷100次之后得" + "r="+str(g_r_hat[0])+"的后验密度为：" + str(sp.beta.pdf(g_r_hat[0], g_a[2] +g_y_N[0] ,g_b[2] + g_N[0] - g_y_N[0])))
print()
print("对硬币一无所知条件下，投掷1000次之后得" + "r="+str(g_r_hat[1])+"的后验密度为：" + str(sp.beta.pdf(g_r_hat[1], g_a[0] +g_y_N[1] ,g_b[0] + g_N[1] - g_y_N[1])))
print("硬币均匀条件下，投掷1000次之后得" + "r="+str(g_r_hat[1])+"的后验密度为：" + str(sp.beta.pdf(g_r_hat[1], g_a[1] +g_y_N[1] ,g_b[1] + g_N[1] - g_y_N[1])))
print("硬币有偏条件下，投掷1000次之后得" + "r="+str(g_r_hat[1])+"的后验密度为：" + str(sp.beta.pdf(g_r_hat[1], g_a[2] +g_y_N[1] ,g_b[2] + g_N[1] - g_y_N[1])))
计算结果如下：
之后，计算赢钱概率
赢钱概率如下计算：
P
(
Y
≤
6
)
=
1
−
P
(
Y
≥
7
)
P(Y\leq6) = 1-P(Y\geq7)
P(Y≤6)=1−P(Y≥7)
而
P
(
Y
=
y
)
P(Y=y)
P(Y=y)为：
P
（
Y
=
y
）
=
(
n
k
)
r
y
(
1
−
r
)
N
−
y
P（Y = y） = \binom{n}{k}r^{y}(1-r)^{N-y}
P（Y=y）=(kn)ry(1−r)N−y
计算的方法有两种，第一种是用
r
r
r期望值来代替实际的
r
r
r值，第二种是计算赢钱概率关于
r
r
r的期望
E
p
(
r
∣
y
N
)
P
(
Y
n
e
w
=
y
n
e
w
∣
r
)
=
(
N
n
e
w
y
n
e
w
)
Γ
(
δ
+
γ
)
Γ
(
δ
)
Γ
(
γ
)
Γ
(
δ
+
y
n
e
w
)
Γ
(
γ
+
N
n
e
w
−
y
n
e
w
)
Γ
(
γ
+
δ
+
N
n
e
w
)
E_{p(r|y_N)}{P(Y_{new}=y_{new}|r)} = \binom {N_{new}} {y_{new}} \frac{\Gamma(\delta+\gamma)}{\Gamma(\delta)\Gamma(\gamma)}\frac{\Gamma(\delta+y_{new})\Gamma(\gamma+N_{new}-y_{new})}{\Gamma(\gamma+\delta+N_{new})}
Ep(r∣yN)P(Ynew=ynew∣r)=(ynewNnew)Γ(δ)Γ(γ)Γ(δ+γ)Γ(γ+δ+Nnew)Γ(δ+ynew)Γ(γ+Nnew−ynew)
运用r的估计值计算赢钱概率
首先我们用第一种来计算。
r
r
r的期望如下所示：
E
p
(
r
∣
y
N
)
R
=
δ
δ
+
γ
E_{p(r|y_N)}{R} = \frac{\delta}{\delta+\gamma}
Ep(r∣yN)R=δ+γδ
δ
=
α
+
y
N
\delta = \alpha + y_N
δ=α+yN,
\quad
γ
=
β
+
N
−
y
N
\gamma = \beta + N-y_N
γ=β+N−yN
# 投掷了100次和1000次之后 控形参数的变化
g_delta = []
g_gamma = []
for i in range(len(g_N)):
g_delta.append(np.array(g_a) + g_y_N[i])
g_gamma.append(np.array(g_b) + g_N[i] - g_y_N[i])
g_delta = np.array(g_delta) # array化
g_gamma = np.array(g_gamma)# array 化
g_E_r = g_delta/(g_delta + g_gamma) # 计算E 2*3 维的 两种投掷方式，各三种先验
赢钱概率计算：
g_P_win = []
for i in range(len(g_N)):
P_i = []
for j in range(len(g_E_r[0])):
P_j = 0
for k in range(7,11): # 累计概率
C = math.factorial(10)//(math.factorial(k)*math.factorial(10 - k)) # 组合数
tmp = C * g_E_r[i][j]**k * (1 - g_E_r[i][j])**(10-k) # 概率
P_j = tmp + P_j
P_i.append(1 - P_j)
g_P_win.append(P_i)
结果展示：
# 赢钱的概率为：
print("采用点估计方法计算获胜概率：")
print("投掷" + str(g_N[0]) +"次" + ",在一无所知的先验下获胜的概率为：" + str(round(g_P_win[0][0],3)) )
print("投掷" + str(g_N[0]) +"次" + ",在硬币均匀先验下获胜的概率为：" + str(round(g_P_win[0][1],3)) )
print("投掷" + str(g_N[0]) +"次" + ",在硬币有偏先验下获胜的概率为：" + str(round(g_P_win[0][2],3)) )
print()
print("投掷" + str(g_N[1]) +"次" + ",在一无所知的先验下获胜的概率为：" + str(round(g_P_win[1][0],3)) )
print("投掷" + str(g_N[1]) +"次" + ",在硬币均匀先验下获胜的概率为：" + str(round(g_P_win[1][1],3)) )
print("投掷" + str(g_N[1]) +"次" + ",在硬币有偏先验下获胜的概率为：" + str(round(g_P_win[1][2],3)) )
下面采用对后验积分的方法来计算获胜概率
在这个地方要注意，如果我们直接用
Γ
\Gamma
Γ函数来求的化，会导致溢出。在这里，我们将式子转化为下式来求
E
p
(
r
∣
y
N
)
P
(
Y
n
e
w
=
y
n
e
w
∣
r
)
=
(
N
n
e
w
y
n
e
w
)
Γ
(
δ
+
γ
)
Γ
(
δ
)
Γ
(
γ
)
Γ
(
δ
+
y
n
e
w
)
Γ
(
γ
+
N
n
e
w
−
y
n
e
w
)
Γ
(
γ
+
δ
+
N
n
e
w
)
=
(
N
n
e
w
y
n
e
w
)
β
(
r
,
δ
,
γ
)
β
(
r
,
δ
+
y
n
e
w
,
γ
+
N
n
e
w
−
y
n
e
w
)
⋅
r
y
n
e
w
⋅
(
1
−
r
N
n
e
w
−
y
n
e
w
)
\begin{aligned} E_{p(r|y_N)}{P(Y_{new}=y_{new}|r)} =& \binom {N_{new}} {y_{new}} \frac{\Gamma(\delta+\gamma)}{\Gamma(\delta)\Gamma(\gamma)}\frac{\Gamma(\delta+y_{new})\Gamma(\gamma+N_{new}-y_{new})}{\Gamma(\gamma+\delta+N_{new})}\\ =& \binom {N_{new}} {y_{new}}\frac{\beta(r,\delta,\gamma)}{\beta(r,\delta+y_{new},\gamma+N_{new}-y_{new})}\cdot r^{y_{new}}\cdot (1-r^{N_{new}-y_{new}})\end{aligned}
Ep(r∣yN)P(Ynew=ynew∣r)==(ynewNnew)Γ(δ)Γ(γ)Γ(δ+γ)Γ(γ+δ+Nnew)Γ(δ+ynew)Γ(γ+Nnew−ynew)(ynewNnew)β(r,δ+ynew,γ+Nnew−ynew)β(r,δ,γ)⋅rynew⋅(1−rNnew−ynew)
在实际计算的时候，我们取
r
=
0.5
r=0.5
r=0.5
g_N_new = 10 # 研究十次的概率
g_P_win_E = []
r = 0.5
for i in range(len(g_delta)):
p_i = []
for j in range(len(g_delta[0])):
p_i_j = 0
for k in range(7,11):
tmp = 0
C = math.factorial(10)//(math.factorial(k)*math.factorial(10 - k)) # 组合数
tmp_1 = sp.beta.pdf(r,g_delta[i][j],g_gamma[i][j])/sp.beta.pdf(r,g_delta[i][j]+k, g_gamma[i][j]+10-k)
tmp_2 = (r**k) * (1-r)**(10-k)
tmp = C * tmp_1*tmp_2
p_i_j =p_i_j + tmp
p_i.append(1 - p_i_j)
g_P_win_E.append(p_i)
第二种方法赢钱的概率计算：
# 赢钱的概率为：
print("采用期望积分估计方法计算获胜概率：")
print("投掷" + str(g_N[0]) +"次" + ",在一无所知的先验下获胜的概率为：" + str(round(g_P_win_E[0][0],3)) )
print("投掷" + str(g_N[0]) +"次" + ",在硬币均匀先验下获胜的概率为：" + str(round(g_P_win_E[0][1],3)) )
print("投掷" + str(g_N[0]) +"次" + ",在硬币有偏先验下获胜的概率为：" + str(round(g_P_win_E[0][2],3)) )
print()
print("投掷" + str(g_N[1]) +"次" + ",在一无所知的先验下获胜的概率为：" + str(round(g_P_win_E[1][0],3)) )
print("投掷" + str(g_N[1]) +"次" + ",在硬币均匀先验下获胜的概率为：" + str(round(g_P_win_E[1][1],3)) )
print("投掷" + str(g_N[1]) +"次" + ",在硬币有偏先验下获胜的概率为：" + str(round(g_P_win_E[1][2],3)) )
比较两种方法的计算结果，可以看到，投掷次数为100时，二者计算的赢钱概率有一定差距；但当投掷次数为1000时，赢钱的概率已经没有差距了。这体现了大数定理的作用
最后，计算边缘似然
边缘似然的计算公式如下如下所示：
p
(
y
N
∣
α
,
β
)
=
(
N
n
e
w
y
n
e
w
)
Γ
(
α
+
β
)
Γ
(
α
)
Γ
(
β
)
Γ
(
α
+
y
n
e
w
)
Γ
(
β
+
N
n
e
w
−
y
n
e
w
)
Γ
(
β
+
α
+
N
n
e
w
)
=
(
N
n
e
w
y
n
e
w
)
β
(
r
,
α
,
β
)
β
(
r
,
α
+
y
n
e
w
,
β
+
N
n
e
w
−
y
n
e
w
)
⋅
r
y
n
e
w
⋅
(
1
−
r
N
n
e
w
−
y
n
e
w
)
\begin{aligned} p(y_N|\alpha,\beta) =& \binom {N_{new}} {y_{new}} \frac{\Gamma(\alpha+\beta)}{\Gamma(\alpha)\Gamma(\beta)}\frac{\Gamma(\alpha+y_{new})\Gamma(\beta+N_{new}-y_{new})}{\Gamma(\beta+\alpha+N_{new})} \\ =&\binom {N_{new}} {y_{new}}\frac{\beta(r,\alpha,\beta)}{\beta(r,\alpha+y_{new},\beta+N_{new}-y_{new})}\cdot r^{y_{new}}\cdot (1-r^{N_{new}-y_{new}}) \end{aligned}
p(yN∣α,β)==(ynewNnew)Γ(α)Γ(β)Γ(α+β)Γ(β+α+Nnew)Γ(α+ynew)Γ(β+Nnew−ynew)(ynewNnew)β(r,α+ynew,β+Nnew−ynew)β(r,α,β)⋅rynew⋅(1−rNnew−ynew)
代码部分与之前的类似，只是略做修改。
g_N_new = 10 # 研究十次的边缘似然
g_p_edge = []
r = 0.5
for i in range(len(g_N)):
p_i = []
for j in range(len(g_a)):
C = math.factorial(g_N[i])//(math.factorial(g_y_N[i])*math.factorial(g_N[i] - g_y_N[i])) # 组合数
tmp_1 = sp.beta.pdf(r,g_a[j],g_b[j])/sp.beta.pdf(r,g_a[j]+g_y_N[i], g_b[j]+g_N[i] - g_y_N[i])
tmp_2 = (r**g_y_N[i]) * (1-r)**(g_N[i]-g_y_N[i])
tmp = C * tmp_1*tmp_2
p_i.append(tmp)
g_p_edge.append(p_i)
展示边缘似然：
print("边缘似然：")
print("投掷" + str(g_N[0]) +"次" + ",在一无所知的先验下边缘似然为：" + str(round(g_p_edge[0][0],6)) )
print("投掷" + str(g_N[0]) +"次" + ",在硬币均匀先验下边缘似然为：" + str(round(g_p_edge[0][1],6)) )
print("投掷" + str(g_N[0]) +"次" + ",在硬币有偏先验下边缘似然为：" + str(round(g_p_edge[0][2],6)) )
print()
print("投掷" + str(g_N[1]) +"次" + ",在一无所知的先验下边缘似然为：" + str(round(g_p_edge[1][0],6)) )
print("投掷" + str(g_N[1]) +"次" + ",在硬币均匀先验下边缘似然为：" + str(round(g_p_edge[1][1],6)) )
print("投掷" + str(g_N[1]) +"次" + ",在硬币有偏先验下边缘似然为：" + str(round(g_p_edge[1][2],6)) )
我们可依据边缘似然的计算结果来选择先验。观察计算结果，当投掷1000次时，硬币有偏先验下边缘似然最高，故我们选择硬币有偏这一先验是合理的。
扫一扫
2348
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
