- 贝叶斯简介
先验：体现了对参数的了解，设定参数的分布
似然：引入观测数据，反映在给定参数下，观测数据的可信度
后验：基于观测数据，获得参数的分布和可信度
贝叶斯定理：后验正比于似然乘以先验
- 抛硬币案例
以抛硬币为例
参数为硬币朝上的概率A（表示抛一次硬币时正面朝上的可能性），定义为P(A)，此为先验（此处指的是概率P的概率）
似然为已知该参数A，观测值B（N次抛硬币正面朝上的次数）对应的概率，也就是已知硬币朝上的概率，实验中出现这个观测值的概率
后验则为通过一系列观测值B，获得参数A的分布
首先选择似然分布，此处选择了二项分布，表示N次抛硬币实验y次正面朝上的概率分布
（似然是指的观测值的可能分布）
（此时二项分布表达式内既包含了参数A，也包含了观测值B，也包含了试验次数N）
实际上就是已知硬币朝上的概率P（参数A），实验了N次，求实验中发生y次硬币朝上（参数B）的概率
其次为参数A选择先验分布，其只与参数A直接相关
选用了beta分布，因为其输出为0到1，与P一致
计算后验，将似然去乘以先验
由于只考虑了其正比的比例，所以可以把与A和B不相关的系数都去掉，得到的后验分布表达式和先验分布很相似
也是一个beta分布，此时就可以使先验beta分布的系数和后验beta分布的系数构成一个等式，后验的分布就可以使用先验的系数
- 注意点
贝叶斯分析的结果是后验分布而不是某个值，描述了根据给定数据和模型得到不同数值的可能性
后验最可能的值是根据后验分布的形态决定的（也就是后验分布的峰值）
后验分布的离散程度与我们对参数的不确定性有关，分布越离散，不确定性越大
给定足够多的数据时，两个不同先验的贝叶斯模型会趋近于相同的结果，这个过程也称之为收敛
不同后验收敛到相同分布取决于数据和模型
- 先验的选择
使用信息量较小的先验，比如只使用先验来限定其取值范围等等，称为正则化先验
可以使用多个先验来验证先验的正确性
- 最大后验密度区间（HPD）
包含一定比例概率密度的最小区间，一般有95%HPD或者98%HPD
如果HPD区间为[2,5]，其含义是参数位于[2,5]的概率是0.95
在贝叶斯学派中，可以告诉我们参数取值的概率，这在频率学中是不可能的
因为频率学认为参数是一个固定值，置信区间只能包含或者不包含参数的真实值
def naive_hpd(post): sns.kdeplot(post) HPD = np.percentile(post, [2.5, 97.5]) print(HPD) plt.plot(HPD, [0, 0], label='HPD {:.2f} {:.2f}'.format(*HPD), linewidth=8, color='k') plt.legend(fontsize=16); plt.xlabel(r"$\theta$", fontsize=14) plt.gca().axes.get_yaxis().set_ticks([]) np.random.seed(1) post = stats.beta.rvs(5, 11, size=1000) #根据beta分布（5，11），生成数量为1000的array naive_hpd(post) plt.xlim(0, 1) plt.show()
实际上是找到2.5，97.5概率对应的值即可
- 后验预测检查
一旦获得后验分布，就可以根据后验获得未来的数据来进行预测
通过比较未来数据和已有数据，可以进行一致性检查
- 抛硬币代码实现
theta_real = 0.35 #参数真值 trials = [0, 1, 2, 3, 4, 8, 16, 32, 50, 150] # 实验次数 data = [0, 1, 1, 1, 1, 4, 6, 9, 13, 48] # 每次实验的正面朝上的次数 beta_params = [(1, 1), (0.5, 0.5), (20, 20)] # beta先验分布的参数对 dist = stats.beta x = np.linspace(0, 1, 100) # 0到1之间均匀采样100个数，用于计算后验分布的概率密度 print(x) for idx, N in enumerate(trials): # 一次循环表示一次实验过程 print(idx) print(N) if idx == 0: plt.subplot(4,3, 2) else: plt.subplot(4,3, idx+3) y = data[idx] #每一次实验对应的实验数据 print(y) for (a_prior, b_prior), c in zip(beta_params, ('b', 'r', 'g')): # 使用不同的先验参数，绘制三张图，分别对应蓝色，红色，绿色 p_theta_given_y = dist.pdf(x, a_prior + y, b_prior + N - y) # 将beta先验分布参数更换为后验分布的参数（此处就相当于引入了观测值y及N），并求解x对应的概率密度，此处已经获得了后验分布的参数 print(p_theta_given_y) plt.plot(x, p_theta_given_y, c) plt.fill_between(x, 0, p_theta_given_y, color=c, alpha=0.6) plt.axvline(theta_real, ymax=0.3, color='k') # 用竖直黑线绘制出参数的真值 plt.plot(0, 0, label="{:d} experiments\n{:d} heads".format(N, y), alpha=0) plt.xlim(0,1) plt.ylim(0,12) plt.xlabel(r"$\theta$") plt.legend() plt.gca().axes.get_yaxis().set_visible(False) plt.tight_layout() plt.show()
观测值B包含了实验次数（experiments）以及朝上次数（heads）
将值代入后验分布的参数中，并均匀采样参数A的可能值导入后验分布中，即得到了参数A每个可能值的概率密度分布
核心为
p_theta_given_y = dist.pdf(x, a_prior + y, b_prior + N - y)
此处就实现了使用观测值来更改先验，使其变成后验，从而得到参数的可能值分布
（参数的后验由先验和观测值共同组成）
- 自动的后验
实际上通过分析方法获得后验分布在实际中是很难成立的
通过其他非分析解的方法，也可以对后验进行计算
非马尔可夫方法：
网格计算，二次近似，变分方法（用简单的分布来近似后验）
马尔科夫方法：
马尔科夫链蒙特卡洛（MCMC）
蒙特卡洛：一种随机采样的方法
马尔科夫链：每个状态转移到其他状态的概率只与当前状态有关，此状态链就成为马尔科夫链
为了找到细节平衡的马尔科夫链，可以使用Metropolis-Hasting算法采样到一条链
该算法拥有自己的方式（即拒绝和接受的标准，采样是随机的，但是接受和拒绝的标准是人为的）（接收新的值的概率正比于新旧值之比）来确定采样方向，从而得到一条采样轨迹或者采样链
采样值就是后验的近似，出现次数最多的值就是后验中最可能的值
混合蒙特卡洛方法：避免采样过慢且结果自相关，需要计算梯度
NUTS算法是其中一种HMC算法，不需要手动设置参数
- 使用网格算法硬币案例
def posterior_grid(grid_points=100, heads=6, tosses=9): """ A grid implementation for the coin-flip problem """ # 参数构成的网格（0和1之前构成的网格） grid = np.linspace(0, 1, grid_points) print(grid) # 定义先验，将5重复网格数次 prior = np.repeat(5, grid_points) print(prior) #prior = (grid <= 0.4).astype(int) # truncated #prior = abs(grid - 0.5) # "M" prior # 计算每个网格点的似然值,引入可观测值和网格点 likelihood = stats.binom.pmf(heads, tosses, grid) print(likelihood) # 计算似然乘以先验，似然相当于就是更新先验的一个参数因子 unstd_posterior = likelihood * prior print(unstd_posterior) # 标准化后验，使其值为1 posterior = unstd_posterior / unstd_posterior.sum() print(posterior) return grid, posterior points = 15 h, n = 1, 4 # 四次实验只有一次正面朝上 grid, posterior = posterior_grid(points, h, n) plt.plot(grid, posterior, 'o-') plt.plot(0, 0, label='heads = {}\ntosses = {}'.format(h, n), alpha=0) plt.xlabel(r'$\theta$', fontsize=14) plt.legend(loc=0, fontsize=14) plt.show()
网格：使用0到1的均布数据
[0. 0.07142857 0.14285714 0.21428571 0.28571429 0.35714286
0.42857143 0.5 0.57142857 0.64285714 0.71428571 0.78571429
0.85714286 0.92857143 1. ]
先验：定义全为5的均匀分布
[5 5 5 5 5 5 5 5 5 5 5 5 5 5 5]
似然：引入观测值（此处是引入试验次数和朝上次数）设定一个似然分布，然后导入网格，即可获得每个网格值对应的概率分布
[0. 0.22875885 0.35985006 0.41576426 0.41649313 0.37952936
0.31986672 0.25 0.17992503 0.11713869 0.0666389 0.03092461
0.00999584 0.0013536 0. ]
（似然可以看做是更新先验的一个因素）
后验：先验乘以似然来更新先验获得后验
[0. 1.14379425 1.79925031 2.07882132 2.08246564 1.89764681
1.59933361 1.25 0.89962516 0.58569346 0.3331945 0.15462307
0.04997918 0.00676801 0. ]
将后验归一化
[0. 0.08239883 0.12961782 0.14975809 0.15002063 0.1367063
0.11521584 0.09004988 0.06480891 0.0421933 0.0240033 0.01113903
0.0036005 0.00048757 0. ]
考虑网格数为15，实验4次，朝上1次
考虑网格数为1000，实验200次，朝上85次
- 求圆的面积蒙特卡洛代码实现
蒙特卡洛主要用于表征随机性
在变成为2R的正方形内随机撒N个点
在正方形内画一个半径为R的圆，计算在圆内点的个数
（使用关于x，y的不等式来判断其是否在圆内）
得出pi的估计值
圆内点数除以总点数等于圆面积与正方形面积之比（pi*R*R/4*R*R）
(使用点的个数来代表面积)
N = 10000 x, y = np.random.uniform(-1, 1, size=(2, N)) # 随机撒点，为中心为原点的正方形 inside = (x**2 + y**2) <= 1 # 计算位于圆内的点的个数 pi = inside.sum()*4/N # 换算出pi error = abs((pi - np.pi)/pi)* 100 outside = np.invert(inside) plt.plot(x[inside], y[inside], 'b.') plt.plot(x[outside], y[outside], 'r.') plt.plot(0, 0, label='$\hat \pi$ = {:4.3f}\nerror = {:4.3f}%'.format(pi, error), alpha=0) plt.axis('square') plt.legend(frameon=True, framealpha=0.9, fontsize=16); plt.show()
获得pi值
- Metropolis-Hasting算法代码实现
def metropolis(func, steps=10000): """A very simple Metropolis implementation""" samples = np.zeros(steps) # 指定一维的0矩阵，作为采样链的初始值 old_x = func.mean() # 计算beta分布的平均值，作为采样的起始位置 old_prob = func.pdf(old_x) # 计算采样起始值的概率值，应该具有比较大的概率值，用于后续的比较 for i in range(steps): # 一次循环即采样一次 new_x = old_x + np.random.normal(0, 1) # 以beta分布平均值作为起始值，将正态分布作为更新的偏差值 new_prob = func.pdf(new_x) # 计算采样值对应的概率 acceptance = new_prob/old_prob # 计算接受概率，即新旧测量值概率之比 if acceptance >= np.random.random(): #随机取值作为筛选条件，更新采样值 samples[i] = new_x old_x = new_x old_prob = new_prob else: samples[i] = old_x return samples np.random.seed(345) func = stats.beta(0.4, 2) # 定义了一个指定参数的beta分布 samples = metropolis(func=func) # 对beta分布进行采样，返回一条采样链 print(samples) x = np.linspace(0.01, .99, 100) y = func.pdf(x) plt.xlim(0, 1) plt.plot(x, y, 'r-', lw=3, label='True distribution') # 真实的参数分布及其对应概率密度函数，也就是最先定义的beta分布 plt.hist(samples, bins=100, density=True, edgecolor='k', label='Estimated distribution') # 根据采样链，获得其概率密度分布 plt.xlabel('$x$', fontsize=14) plt.ylabel('$pdf(x)$', fontsize=14) plt.legend(fontsize=14) plt.show()
该算法的作用主要是
在给定一个分布（比如此处的beta分布）后，通过采样获得一条采样链，该采样链对应的概率密度分布与原分布基本一致
此处就实现了一个由分布到采样数据之间的转换过程
- 基于pymc3的代码实现
np.random.seed(123) n_experiments = 4 theta_real = 0.35 # 生成数据使用的真值，是需要估计的值 data = stats.bernoulli.rvs(p=theta_real, size=n_experiments) # 按照二项分布及真值生成数据，作为了观测数据 print(data) if __name__ == '__main__': with pm.Model() as our_first_model: # 指定参数的先验（也就是硬币朝上的概率） # 按照先验theta来生成采样链 theta = pm.Beta('theta', alpha=1, beta=1) # 定义似然的分布（也就是已知硬币朝上的概率时，在实验中硬币朝上的概率分布，也就是观测值的可能性），此处传入了观测值 # 似然不会进行采样，因为其是一个观测变量 # 似然表示的是观测值的可能性，其基于参数和观测值，作为先验的更新因子来获得后验 y = pm.Bernoulli('y', p=theta, observed=data) #y = pm.Binomial('theta',n=n_experimentos, p=theta, observed=sum(datos)) start = pm.find_MAP() step = pm.Metropolis() trace = pm.sample(1000, step=step, start=start) # 从最大后验值处进行采样，此处获得是参数theta的后验采样trace burnin = 100 # no burnin # 该参数用于处理老化，因为采样的开始部分可能效果较差 chain = trace[burnin:] # trace就是获得的后验链，由之前的采样链计算得到 pm.traceplot(chain) # 每一条采样链的后验概率密度及其采样图像 pm.forestplot(chain) # 每一条采样链的森林图，表示了其一定的HDI区间 pm.autocorrplot(chain) # 每一条采样链的自相关图像 pm.plot_posterior(chain) # 后验的hpd曲线 plt.show() print(pm.summary(chain)) # 获得后验的详细参数
指定参数先验
指定似然分布
对先验进行采样
计算出后验
每条采样链的森林图
后验的所有信息
