贝叶斯原理在硬币实验中的应用

贝叶斯原理在硬币实验中的应用

  相关资料：《Pattern Recognition and Machine Learning》 by Christopher M. Bishop

  以“贝叶斯原理在硬币实验中的应用”为背景，对其进行代码实现，以助于对相关理论部分的理解

原理

掷硬币问题的过拟合

掷硬币问题描述

  非常经典的掷硬币情景：将一枚硬币连续掷N次，每次的结果都会是正面朝上或反面朝上，将这些结果全部记录下来。现在要做的，是仅根据这些结果，来计算硬币正面朝上的概率

最大似然估计法的缺陷

  撇开生活经验不谈，假定硬币正面朝上的概率应该为$\mu$，且记硬币朝上时，观测值$x=1$，反之，观测值$x=0$。则每次观测值的分布可以记录为如下形式（Bernoulli分布）：

$$Bern(x|\mu )={\mu }^x(1-\mu {)}^{1-x}$$

  $|$符号右边的$\mu$表示这是一个控制分布的参数

  现在，假定对此硬币进行了N次观测（掷了N次），得到数据集（观测值的集合） D = { x 1 ,   … ,   x N } \mathcal D=\{x_1,\ \dots,\ x_N\} D={x1​, …, xN​}，则根据各次观测之间相互独立的特性，可以计算出数据集处于各种状态（即对于各种可能的N个观测值的组合）的概率：

$$p(\mathcal{D}|\mu )=\underset{n=1}{\overset{N}{\Pi }}p(x_n|\mu )$$

  两边取对数，得到：

ln ⁡ p ( D ∣ μ ) = ∑ n = 1 N ln ⁡ p ( x n ∣ μ ) = ∑ n = 1 N { x n ln ⁡ μ + ( 1 − x n ) ln ⁡ ( 1 − μ ) } \ln p(\mathcal D|\mu)=\sum\limits_{n=1}^N\ln p(x_n|\mu)=\sum\limits_{n=1}^N\{x_n\ln\mu + (1-x_n)\ln(1-\mu)\} lnp(D∣μ)=n=1∑N​lnp(xn​∣μ)=n=1∑N​{xn​lnμ+(1−xn​)ln(1−μ)}

  试图找到上述$\ln p(\mathcal{D}|\mu )$的最大值，即求出满足下列等式的${\mu }_{ML}$：

$$\frac{\partial \ln p(\mathcal{D}|\mu )}{\partial \mu }=0$$

  得到：${\mu }_{ML}=\frac{1}{N}\sum _{n=1}^N{x}_n$. 对于上述提及的掷硬币问题，假定有m次是正面朝上，则最后的结果会是${\mu }_{ML}=\frac{m}{N}$。这便是使用极大似然估计法，计算估测的硬币正面朝上的概率

  这个估算出来的概率，在N无限大的时候是准确的，但在N非常小的时候就会出现问题。例如：N=3，连续三次硬币都是正面朝上，则${\mu }_{ML}=1$，这显然是不合理的。这是一种过拟合的问题

使用贝叶斯思想进行改进

  贝叶斯思想可以表示为：

$$posterior\propto likelihood\times prior$$

  即后验概率正比于似然度与先验概率的乘积。在这里，参数$\mu$对应的先验概率分布为$p(\mu )$

  在掷硬币问题中，受Bernoulli分布的启发，将先验概率的形式设为：

$$Beta(\mu |a,b)=\frac{\Gamma (a+b)}{\Gamma (a)\Gamma (b)}{\mu }^{a-1}(1-\mu {)}^{b-1}$$

  其中，$\Gamma$为gamma函数，其定义为：

$$\Gamma (x)={\int }_0^{\infty }{u}^{x-1}{e}^{-u}du$$

  这使得Beta函数满足归一化条件：

$${\int }_0^{1}Beta(\mu |a,b)d\mu =1$$

  此外，似然度选用binomial分布：

$$Bim(m|N,\mu )=C_N^m{\mu }^m(1-\mu {)}^{N-m}$$

  则后验概率也将满足Beta函数的形式：

$$p(\mu |m,l,a,b)=\frac{\Gamma (m+a+l+b)}{\Gamma (m+a)\Gamma (l+b)}{\mu }^{m+a-1}(1-\mu {)}^{l+b-1}$$

  其中，$l=N-m$

  现在，对于掷硬币过程中的任意时刻，都有着当前的先验概率。在此基础上再进行投掷，则可以根据投掷的情况确定似然度，与之相乘从而得到后验概率，随后再拿后验概率来替换先验概率，即可使概率分布不断地被迭代更新

  注意到初始时刻，没有任何记录，因此在$Beta(\mu |a,b)$中，$a=b=1$，对应的函数为常函数，其值为1

实现

Gamma函数的计算

  由于Gamma函数在先验概率和后验概率表达式中频繁出现，因此计算Gamma函数的值是不可避免的。因为Gamma函数是一个广义积分，使用解析的方式求解会非常地困难。所以在这里，借助了python的scipy模块来进行数值积分计算

  例如，对$\Gamma (1)$的求值相关代码如下：

from scipy import integrate
import numpy as np
import math

def gammaFunc(u, x):
    return u ** (x - 1) * math.exp(-u)

S, err = integrate.quad(gammaFunc, 0, np.inf, args = (1))
print("积分结果：%f" % S)

  gammaFunc(u, x)中的第一个参数u是积分变量，而第二个参数x则是任意指定的变量。任意指定的变量是在调用integrate.quad()方法的时候，由args参数进行传递。例如，这里的args = (1)，就代表着令gammaFunc的第二个参数x=1

  运行后，得到结果：

积分结果：1.000000

  现在，查看x取不同整数时的计算结果：

from scipy import integrate
import numpy as np
import math
import matplotlib.pyplot as plt

def gammaFunc(u, x):
    return u ** (x - 1) * math.exp(-u)

numList = np.arange(1, 51, 1)
resultList = []
for eachNum in numList:
    S, err = integrate.quad(gammaFunc, 0, np.inf, args = (eachNum))
    resultList.append(S)
fig, ax = plt.subplots(1, 1, figsize = (16, 9))
ax.stem(numList, resultList)
ax.set_yscale("log")
ax.set_ylim(bottom = 0.5)
fig.savefig("gamma_int.png")
plt.show()

  观察到$\Gamma (x)$的值是随着x的增加而急剧增加的（图中纵轴使用了对数尺度）

beta函数的计算

  在能够计算$\Gamma (x)$后，就可以对不同a, b下的beta函数进行计算了：

from scipy import integrate
import numpy as np
import math
import matplotlib.pyplot as plt

def gammaFunc(u, x):
    return u ** (x - 1) * math.exp(-u)

def gammaFuncV(x):
    S, err = integrate.quad(gammaFunc, 0, np.inf, args = (x))
    return S

def betaFunc(mu, a, b):
    beta = gammaFuncV(a + b)/(gammaFuncV(a) * gammaFuncV(b)) * (mu ** (a - 1)) * ((1 - mu) ** (b - 1))
    return beta

fig, ax = plt.subplots(6, 6, figsize = (36, 36))
i = 0
for a in range(1, 18, 3):
    j = 0
    for b in range(1, 18, 3):
        mu = np.arange(0, 1, 0.01)
        beta = betaFunc(mu, a, b)
        ax[i][j].plot(mu, beta, 'b-')
        ax[i][j].set_title("a = %d, b = %d" % (a, b))
        ax[i][j].set_xlim([0, 1])
        j += 1
    i += 1
fig.savefig("beta.png")

  由此大致作出图像，以观察不同的a、b值对beta函数曲线的影响：

  注意到a和b值越相近时，beta函数曲线的极值越靠近中间

掷硬币实验的贝叶斯方法模拟

  标题可能不准确，因为不知道该怎么起标题了

  因为先验概率和后验概率的函数形式都是beta函数，且后验概率是在先验概率的基础上改变a、b值。所以在利用贝叶斯原理计算后验概率时，只需在先验概率的基础上更改a、b值即可，而不必专门去乘似然度

  在硬币实验中，当硬币正面朝上时，似然度函数为$\mu$，反之为$1-\mu$. 即正面朝上时，更新a=a+1，反之，更新b=b+1

  由此即可进行模拟实验：

import random
from scipy import integrate
import numpy as np
import math
import matplotlib.pyplot as plt

observeTS = [0, 1, 2, 3, 5, 7, 12, 18, 30] # 在掷硬币多少次后对当前的beta函数进行观察

def gammaFunc(u, x):
    return u ** (x - 1) * math.exp(-u)

def gammaFuncV(x):
    S, err = integrate.quad(gammaFunc, 0, np.inf, args = (x))
    return S

def betaFunc(mu, a, b):
    beta = gammaFuncV(a + b)/(gammaFuncV(a) * gammaFuncV(b)) * (mu ** (a - 1)) * ((1 - mu) ** (b - 1))
    return beta

# 掷硬币，返回0或1
def throwCoin():
    return 0 if(random.random() * 2 < 1) else 1

# idx转二维下标
def idx2xy(idx):
    x = idx // 3
    y = idx % 3
    return [x, y]

# 初始化beta函数的横坐标
mu = np.arange(0, 1.01, 0.01)
# 初始化画布及子图下标
fig, ax = plt.subplots(3, 3, figsize = (9, 9))
idx = 0
# 初始化beta函数的参数a、b
a = 1
b = 1
# 初始化掷硬币的次数
N = 0
# 初始化beta函数
beta = betaFunc(mu, a, b)
if(N in observeTS):
    [x, y] = idx2xy(idx)
    ax[x][y].plot(mu, beta, 'b-')
    ax[x][y].set_title("N = %d" % N)
    idx += 1
for dummyIterIdx in range(31):
    observeValue = throwCoin()
    if(observeValue == 1):
        a += 1
        beta = betaFunc(mu, a, b)
    else:
        b += 1
        beta = betaFunc(mu, a, b)
    N += 1
    if(N in observeTS):
        [x, y] = idx2xy(idx)
        ax[x][y].plot(mu, beta, 'b-')
        ax[x][y].set_title("N = %d" % N)
        idx += 1
fig.savefig("coin_exp.png")

  当N（掷硬币的次数）逐渐增大时，beta函数的最大值也会逐渐趋近于中心。这里计算到N=30的时候，已经快达到这个数值计算所能承受的极限了（再大一些就会有OverflowError: (34, ‘Result too large’)这个报错出现了——数值太大了），所以就没有再往后模拟了

  不过按照这个情况来看，感觉这个beta函数在N趋近于无穷大的时候，会趋近于$\delta$函数（即在某个点处无穷大，在其他地方全为0）