共计 6629 个字符，预计需要花费 17 分钟才能阅读完成。
对许多人来说，贝叶斯统计依然有些生疏。因为贝叶斯统计中会有一些主观的先验，在没有测试数据的反对下理解他的实践还是有一些艰难的。本文整顿的是作者最近在普林斯顿的一个研讨会上做的演讲幻灯片，这样能够说明为什么贝叶斯办法不仅在逻辑上是正当的，而且应用起来也很简略。这里将以三种不同的形式实现雷同的推理问题。
数据
咱们的例子是在具备歪斜背景的噪声数据中找到峰值的问题，这可能呈现在粒子物理学和其余多重量事件过程中。
首先生成数据：
%matplotlibinline
%configInlineBackend.figure_format='svg'
importmatplotlib.pyplotasplt
importnumpyasnp
defsignal(theta, x):
l, m, s, a, b=theta
peak=l*np.exp(-(m-x)**2/ (2*s**2))
background =a+b*x
returnpeak+background
defplot_results(x, y, y_err, samples=None, predictions=None):
fig=plt.figure()
ax=fig.gca()
ax.errorbar(x, y, yerr=y_err, fmt=".k", capsize=0, label="Data")
x0=np.linspace(-0.2, 1.2, 100)
ax.plot(x0, signal(theta, x0), "r", label="Truth", zorder=0)
ifsamplesisnotNone:
inds=np.random.randint(len(samples), size=50)
fori,indinenumerate(inds):
theta_=samples[ind]
ifi==0:
label='Posterior'
else:
label=None
ax.plot(x0, signal(theta_, x0), "C0", alpha=0.1, zorder=-1, label=label)
elifpredictionsisnotNone:
fori, predinenumerate(predictions):
ifi==0:
label='Posterior'
else:
label=None
ax.plot(x0, pred, "C0", alpha=0.1, zorder=-1, label=label)
ax.legend(frameon=False)
ax.set_xlabel("x")
ax.set_ylabel("y")
fig.tight_layout()
plt.close();
returnfig
# random x locations
N=40
np.random.seed(0)
x=np.random.rand(N)
# evaluate the true model at the given x values
theta= [1, 0.5, 0.1, -0.1, 0.4]
y=signal(theta, x)
# add heteroscedastic Gaussian uncertainties only in y direction
y_err=np.random.uniform(0.05, 0.25, size=N)
y=y+np.random.normal(0, y_err)
# plot
plot_results(x, y, y_err)有了数据咱们能够介绍三种办法了
马尔可夫链蒙特卡罗 Markov Chain Monte Carlo
emcee 是用纯 python 实现的，它只须要评估后验的对数作为参数 θ 的函数。这里应用对数很有用，因为它使指数分布族的剖析评估更容易，并且因为它更好地解决通常呈现的十分小的数字。
importemcee
deflog_likelihood(theta, x, y, yerr):
y_model=signal(theta, x)
chi2= (y-y_model)**2/ (yerr**2)
returnnp.sum(-chi2/2)
deflog_prior(theta):
ifall(theta>-2) and (theta[2] >0) andall(theta<2):
return0
return-np.inf
deflog_posterior(theta, x, y, yerr):
lp=log_prior(theta)
ifnp.isfinite(lp):
lp+=log_likelihood(theta, x, y, yerr)
returnlp
# create a small ball around the MLE the initialize each walker
nwalkers, ndim=30, 5
theta_guess= [0.5, 0.6, 0.2, -0.2, 0.1]
pos=theta_guess+1e-4*np.random.randn(nwalkers, ndim)
# run emcee
sampler=emcee.EnsembleSampler(nwalkers, ndim, log_posterior, args=(x, y, y_err))
sampler.run_mcmc(pos, 10000, progress=True);后果如下：
100%|██████████| 10000/10000 [00:05<00:00, 1856.57it/s]咱们应该始终查看生成的链，确定 burn-in period，并且须要人肉察看平稳性：
fig, axes=plt.subplots(ndim, sharex=True)
mcmc_samples=sampler.get_chain()
labels= ["l", "m", "s", "a", "b"]
foriinrange(ndim):
ax=axes[i]
ax.plot(mcmc_samples[:, :, i], "k", alpha=0.3, rasterized=True)
ax.set_xlim(0, 1000)
ax.set_ylabel(labels[i])
axes[-1].set_xlabel("step number");当初咱们须要细化链因为咱们的样本是相干的。这里有一个办法来计算每个参数的自相干，咱们能够将所有的样本联合起来:
tau=sampler.get_autocorr_time()
print("Autocorrelation time:", tau)
mcmc_samples=sampler.get_chain(discard=300, thin=np.int32(np.max(tau)/2), flat=True)
print("Remaining samples:", mcmc_samples.shape)
#后果
Autocorrelationtime: [122.51626866 75.87228105137.195509 54.63572513 79.0331587]
Remainingsamples: (4260, 5)emcee 的创建者 Dan Foreman-Mackey 还提供了这一有用的包 corner 来可视化样本：
importcorner
corner.corner(mcmc_samples, labels=labels, truths=theta);尽管后验样本是推理的次要根据，但参数轮廓自身却很难解释。然而应用样本来生成新数据则要简略得多，因为这个可视化咱们对数据空间有更多的了解。以下是来自 50 个随机样本的模型评估:
plot_results(x, y, y_err, samples=mcmc_samples)哈密尔顿蒙特卡洛 Hamiltonian Monte Carlo
梯度在高维设置中提供了更多领导。为了实现个别推理，咱们须要一个框架来计算任意概率模型的梯度。这里要害的本局部是主动微分，咱们须要的是能够跟踪参数的各种操作门路的计算框架。为了简略起见，咱们应用的框架是 jax。因为个别状况下在 numpy 中实现的函数都能够在 jax 中的进行类比的替换，而 jax 能够主动计算函数的梯度。
另外还须要计算概率分布梯度的能力。有几种概率编程语言中能够实现，这里咱们抉择了 NumPyro。让咱们看看如何进行主动推理：
importjax.numpyasjnp
importjax.randomasrandom
importnumpyro
importnumpyro.distributionsasdist
fromnumpyro.inferimportMCMC, NUTS
defmodel(x, y=None, y_err=0.1):
# define parameters (incl. prior ranges)
l=numpyro.sample('l', dist.Uniform(-2, 2))
m=numpyro.sample('m', dist.Uniform(-2, 2))
s=numpyro.sample('s', dist.Uniform(0, 2))
a=numpyro.sample('a', dist.Uniform(-2, 2))
b=numpyro.sample('b', dist.Uniform(-2, 2))
# implement the model
# needs jax numpy for differentiability here
peak=l*jnp.exp(-(m-x)**2/ (2*s**2))
background =a+b*x
y_model=peak+background
# notice that we clamp the outcome of this sampling to the observation y
numpyro.sample('obs', dist.Normal(y_model, y_err), obs=y)
# need to split the key for jax's random implementation
rng_key=random.PRNGKey(0)
rng_key, rng_key_=random.split(rng_key)
# run HMC with NUTS
kernel=NUTS(model, target_accept_prob=0.9)
mcmc=MCMC(kernel, num_warmup=1000, num_samples=3000)
mcmc.run(rng_key_, x=x, y=y, y_err=y_err)
mcmc.print_summary()
#后果如下：sample: 100%|██████████|4000/4000 [00:03<00:00, 1022.99it/s, 17stepsofsize2.08e-01. acc. prob=0.94]
mean std median 5.0% 95.0% n_eff r_hat
a -0.13 0.05 -0.13 -0.22 -0.05 1151.15 1.00
b 0.46 0.07 0.46 0.36 0.57 1237.44 1.00
l 0.98 0.05 0.98 0.89 1.06 1874.34 1.00
m 0.50 0.01 0.50 0.49 0.51 1546.56 1.01
s 0.11 0.01 0.11 0.09 0.12 1446.08 1.00
Numberofdivergences: 0还是应用 corner 可视化 Numpyro 的 mcmc 构造:
因为咱们曾经实现了整个概率模型 (与 emcee 相同，咱们只实现后验)，所以能够间接从样本中创立后验预测。上面，咱们将噪声设置为零，以失去纯模型的无噪声示意:
fromnumpyro.inferimportPredictive
# make predictions from posterior
hmc_samples=mcmc.get_samples()
predictive=Predictive(model, hmc_samples)
# need to set noise to zero
# since the full model contains noise contribution
predictions=predictive(rng_key_, x=x0, y_err=0)['obs']
# select 50 predictions to show
inds=random.randint(rng_key_, (50,) , 0, mcmc.num_samples)
predictions=predictions[inds]
plot_results(x, y, y_err, predictions=predictions)基于仿真的推理 Simulation-based Inference
在某些状况下，咱们不能或不想计算可能性。所以咱们只能一个失去一个仿真器（即学习输出之间的映射 θ 和仿真器的输入 D），这个仿真器能够造成似然或后验的近似代替。与产生无噪声模型的传统模仿案例的一个重要区别是，须要在模仿中增加噪声并且噪声模型应尽可能与观测噪声匹配。否则咱们无奈辨别因为噪声引起的数据变动和参数变动引起的数据变动。
importtorch
fromsbiimportutilsasutils
low=torch.zeros(ndim)
low[3] =-1
high=1*torch.ones(ndim)
high[0] =2
prior=utils.BoxUniform(low=low, high=high)
defsimulator(theta, x, y_err):
# signal model
l, m, s, a, b=theta
peak=l*torch.exp(-(m-x)**2/ (2*s**2))
background =a+b*x
y_model=peak+background
# add noise consistent with observations
y=y_model+y_err*torch.randn(len(x))
returny让咱们来看看噪声仿真器的输入:
plt.errorbar(x, this_simulator(torch.tensor(theta)), yerr=y_err, fmt=".r", capsize=0)
plt.errorbar(x, y, yerr=y_err, fmt=".k", capsize=0)
plt.plot(x0, signal(theta, x0), "k", label="truth")当初，咱们应用 sbi 从这些模仿仿真中训练神经后验预计 (NPE)。
fromsbi.inference.baseimportinfer
this_simulator=lambdatheta: simulator(theta, torch.tensor(x), torch.tensor(y_err))
posterior=infer(this_simulator, prior, method='SNPE', num_simulations=10000)NPE 应用条件归一化流来学习如何在给定一些数据的状况下生成后验散布：
Running 10000 simulations.: 0%| | 0/10000 [00:00<?, ?it/s]
Neural network successfully converged after 172 epochs.在推理时，以理论数据 y 为条件简略地评估这个神经后验：
sbi_samples=posterior.sample((10000,), x=torch.tensor(y))
sbi_samples=sbi_samples.detach().numpy()能够看到速度十分快简直不须要什么工夫。
Drawing 10000 posterior samples: 0%| | 0/10000 [00:00<?, ?it/s]而后咱们再次可视化后验样本：
corner.corner(sbi_samples, labels=labels, truths=theta); plot_results(x, y, y_err, samples=sbi_samples)能够看到仿真 SBI 的的后果不如 MCMC 和 HMC 的后果。然而它们能够通过对更多模仿进行训练以及通过调整网络的架构来改良（尽管并不确定改完后就会有进步）。
然而咱们能够看到即便在没有拟然性的状况下，SBI 也能够进行近似贝叶斯推理。
https://avoid.overfit.cn/post/7d210cd0e4424371a7d931b6ee247fc7
作者：Peter Melchior
