在我们开启第二章之前，先去回顾一下第一章的主要内容
我们从对统计建模，概率、条件概率、随机变量以及概率分布的讨论，延申至贝叶斯理论的知识。我们紧接着用一个硬币的问题来介绍基础的贝叶斯模型和数据分析。我们用经典的骰子例子介绍贝叶斯统计中概率分布以及不确定性。我们尝试展现什么是先验并放入相同的数据环境（相同的似然，相同的Beta问题）中进行校准。
最后我们去讨论解释贝叶斯分析的结果。我们假设在有限的数据集下，或者是试验、观测值或是仿真，真正的分布通常是不被知道的。为了在仅有的观测数据下了解真实分布，我们建立了概率模型。其中包括两个基本要素：先验概率和似然函数。之后我们讨论了贝叶斯推断以及得到一个后验分布，这个分布总结了关于这个问题的所有信息。同时也可以利用函数调用后验分布进行预测，得到厚颜预测分布。但是后验分布的作用取决于模型和数据的质量，最终我们总结了贝叶斯数据分析的主要环节。
2 概率相关编程
从这一章开始我们将进入大规模Python编程环节
2.1 概率编程
2.1.1 硬币问题
在第一章中我们用PyMC来解决了一个硬币的问题。我们现在来用相同的人工数据来回顾一下，我们给这个硬币的偏向概率为，被称为真实偏向：
np.random.seed(123)
trials = 4
theta_real = 0.35 # unknown value in a real experiment
data = pz.Binomial(n=1, p=theta_real).rvs(trials)我们得到了一组人工数据，回忆之前我们用伯努利分布作为似然函数，包含参数n=1和p=。用Beta分布作为先验概率，包含参数
。（这个先验和等概率先验[0,1]在这个例子中是等效的）
把这两个放到PyMC中：
with pm.Model() as our_first_model:
θ = pm.Beta('θ', alpha=1., beta=1.)
y = pm.Bernoulli('y', p=θ, observed=data)
idata = pm.sample(1000)第一行表示新建立一个模型，可以看作一个语法糖，第二行是建立一个先验函数，第三行是似然行数，形式有点类似于第二行的先验函数，其中data数据可以是列表，元组、数组或者是pandas的Dataframe，最后一行就是在校准好的概率分布上进行1000次的取样，这里可能需要等上一会。
Auto-assigning NUTS sampler...
Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (4 chains in 4 jobs)
NUTS: [θ]
Sampling 4 chains for 1_000 tune and 1_000 draw iterations
(4_000 + 4_000 draws total) took 1 seconds.没有错误下这应该是返回的结果。第一行表示自动选择的一个抽样器，叫做NUTs，这个不用管，第二行是初始化NUTs中的东西，不管。第三行是表示将会有四个平行抽样的过程，我们把一个平行的过程叫做Chain（在一条链子上，这几条链子不交叉），这个4个chain是它自动选的，也可以在sample函数中利用参数chains=？这个函数来规定有几个chain。第四行表示抽样器中自动选中的模型中的需要校准的参数（就是这个案例中的偏向概率），最终也是在用数据校准参数，得到我们想要的后验分布。但是不一定一直会这种情况，PyMC会根据参数来调整抽样器，不同的抽样器对不同类型的参数校准效果不同，会选择最好的哪个抽样器，可以在sample函数中利用参数step=?来人为的选择抽样器。
最后一行是一个进度条，这个案例中我们要求1000次抽样，但是PyMC进行了8000次，有4条平行链中1000次随机抽样来调试抽样器，这些随机抽样默认被丢弃，PyMC用来增强模型的可信度。之后再有1000次抽样（4x1000）用来得到我们最终想要的后验分布，总共4000+4000次抽样过程。我们可以用sample函数中参数tune=
az.plot_trace(idata)?更改调参过程的迭代轮数，用draw更改这个过程的抽样数量。
2.2 后验分布分析
我们利用代码可以将刚刚校准好的后验分布可视化出来，用ArviZ的plot_trace函数。
az.plot_trace(idata)
我们得到两个子图，左侧子图时核密度估计（Kernel Density Estimation，KDE）图像，四条曲线表示这四个平行抽样的过程，非常的相似。右侧子图，可以得到每一次抽样中得到的确切值。
利用参数可以得到所有平行抽样中的一个单个KDE并且得到一个排列图
az.plot_trace(idata, kind="rank_bars", combined=True, rank_kwargs={"colors": "k"})
（如果你没有设定样式，你的得到的图应该是彩色的，不过图像信息没有丢失）右侧的子图排序图表示我们对选取的参数的信心，在这个图中，有四条chain的直方图，由于随机采样，均匀性会出现一些小偏差，但均匀性的大偏差是链正在探索后验的不同区域的信号。用代码得到对于参数的相关值：
az.summary(idata, kind="stats").round(2)| mean平均值 | sd标准差 | hdi_3%HDI前边有3% | hdi_97%HDI后边有97% | |
|---|---|---|---|---|
| θ | 0.34 | 0.18 | 0.04 | 0.67 |
这里的HDI是选取的最小长度在x轴上的取值，（这属于数据统计中的内容，具体内容可参考相关知识）我们也可以直接画一个后验分布：
az.plot_posterior(idata, figsize=(12, 4))
这个函数通常如果是离散的会画直方图，如果是连续的会画KDE。我们也可以在函数中加入point_estimate=?来看中位数和众数，图中横线表示HDI在94%的取值区间，这个置信度可以利用函数中hdi_prob=?来设置。
2.3 基于后验的决定
SD密度比
在案例中，我们得到一个还可以的后验分布，均值很接近我们最开始设置的0.35，但是我们想知道到底有多大的信心在这个图中可以支持我们说这个硬币的偏向概率就是0.35，或者说在真实世界中说这个硬币就是非常公平的half-half硬币。一个用来评价在给定数据下后验的支持性的方法是用那个取值下的后验概率密度比上先验概率密度，被叫做Savage-Dickey Density Ratio。
az.plot_bf(idata, var_name="θ",
prior=np.random.uniform(0, 1, 10000), ref_val=0.5, figsize=(12, 4))
黑色线表示先验，浅灰表示后验，比如在0.5这个点的支持性就用后验比上先验密度，上边这个数值中BF_10就表示先验比上后验，相反01表示后验比上先验，在这个案例中我们可以认为，在后验分布下得知偏向概率为0.5是先验概率下的1.32倍。这个倍数如下可以被认为：1到3.2倍数一般，3.2到10还可以，10到100很强，大于100非常强。
ROPE实际等效区域
当然在实际中，我们可能不会很武断说这个硬币的偏向概率就为0.5，而是说在一个区间0.45到0.55之间，我们叫这个区间为实际等效区域（Region Of Practical Equivalence, ROPE），一旦ROPE被定下，我们可以用这个来和HDI比较：
- HDI包含ROPE，这个硬币确实如此（公平）
- HDI部分包含ROPE，这个硬币可能如此或不然（公平或不公平，不够确信）
- HDI和ROPE没交集，不包含，这个硬币不然（不公平）
az.plot_posterior(idata, rope=[0.45, .55], figsize=(12, 4))
另一个工具可以帮助我们比较后验分布并决定的是下边内容，它来看左右两侧的概率密度大小区别
az.plot_posterior(idata, ref_val=0.5, figsize=(12, 4))
损失函数
我们定义损失函数来定量分析真实值和我们预测值之间的差距，通常有以下几种常见的损失函数：
- 绝对值损失函数
- 平方损失函数
- 0-1损失函数
在真实情况中我们不会知道真实参数值，但可以用后验的参数取值，来看从0到1之间每个值与其的损失值，可以得到不同损失函数的最小值在哪个地方，
grid = np.linspace(0, 1, 200)
θ_pos = idata.posterior['θ']
lossf_a = [np.mean(abs(i - θ_pos)) for i in grid]
lossf_b = [np.mean((i - θ_pos)**2) for i in grid]
_, ax = plt.subplots(figsize=(12, 3))
for lossf, c in zip([lossf_a, lossf_b], ['C0', 'C1']):
mini = np.argmin(lossf)
ax.plot(grid, lossf, c)
ax.plot(grid[mini], lossf[mini], 'o', color=c)
ax.annotate('{:.2f}'.format(grid[mini]),
(grid[mini], lossf[mini] + 0.03), color=c)
ax.set_yticks([])
ax.set_xlabel(r'$\hat \theta$')
巧合与否，这里都有：平方损失的最小值是点估计的均值，绝对值损失是中位数，0-1损失是众数。这表示如果我们用了一种方法来表示点估计（比如用后验的均值表示估计值大小），同时也就无意识地选择了对应地损失函数。在明面上选择损失函数就可以根据问题来定义，而不是被方法左右你的损失函数地定义。当我们需要不对称的损失的时候，比如大于一个值后，损失的变多或变少，（根据问题定义），我们可以定义一个复杂的损失函数
lossf = []
for i in grid:
if i < 0.5:
f = 1/np.median(θ_pos / np.abs(i**2 - θ_pos))
else:
f = np.mean((i - θ_pos)**2 + np.exp(-i)) - 0.25
lossf.append(f)
mini = np.argmin(lossf)
_, ax = plt.subplots(figsize=(12, 3))
ax.plot(grid, lossf)
ax.plot(grid[mini], lossf[mini], 'o')
ax.annotate('{:.2f}'.format(grid[mini]),
(grid[mini] + 0.01, lossf[mini] + 0.1))
ax.set_yticks([])
ax.set_xlabel(r'$\hat \theta$')
2.4 高斯
我们上边一直在用Beta和伯努利，在问题中我们要用更加复杂的模型。在统计学中，如果一个事件发生的结果被不少于200个参数影响，那么在数据量足够的情况下，这个结果的分布一定呈现高斯分布（中心极限定理，CLT）
高斯推断
这里举例子为蛋白质的化学转化，这个数据中包括每个蛋白质的转化数值（不重要），做一下boxplot（箱型图）
data = np.loadtxt("data/chemical_shifts.csv")
_, ax = plt.subplots(figsize=(12, 3))
ax.boxplot(data, vert=False)
在左侧这个两条竖线之间包含的部分，就为符合高斯分布的数据点，这个箱子的左右两端表示四分位点，中间的分割线表示中位数，但是有两个外边的点在这之外，我们先带着他俩，先建立一个模型，比如
我们定义是一个等概率函数，上下限要包含这个数据的最大小数值，
是一个一半的高斯函数，在缺乏更多信息的情况下，我们可以选择与数据规模相比较大的值。
with pm.Model() as model_g:
μ = pm.Uniform('μ', lower=40, upper=70)
σ = pm.HalfNormal('σ', sigma=5)
Y = pm.Normal('Y', mu=μ, sigma=σ, observed=data)
idata_g = pm.sample(random_seed=4591)
az.plot_trace(idata_g)
az.plot_pair(idata_g, kind='kde', marginals=True)
在这个模型中，有两个参数，那么trace的时候做出的图是两个参数的图。用plot_pair函数可以做出二维的KDE图像。
2.5 后验预测
当得到我们想要的后验分布之后，我们可以来生成预测值，生成的新分布叫做后验预测分布，我们可以看作是给定模型和数据的未来数据（猜想的未来数据），用PyMC来生成分布，并用ArviZ可视化得到的分布内容。
pm.sample_posterior_predictive(
idata_g, model=model_g, extend_inferencedata=True, random_seed=4591
)
az.plot_ppc(idata_g, num_pp_samples=100, figsize=(12, 4))
黑色线是观测数据的KDE，虚线是后验预测的KDE均值，灰色线是100个抽样中的每一个KDE，这个灰色线越杂乱，表示我们对预测数据的越不自信。你可能会发现模型预测的整体比观测值要靠右，可能是因为有两个右侧的异常点引起的。
标签：plot,idata,函数,3rd,Python,lossf,Analysis,后验,我们 From： https://blog.csdn.net/weixin_60302467/article/details/143462246