对于许多人来说，贝叶斯统计仍然有些陌生。由于贝叶斯统计中存在一些主观先验，如果没有测试数据的支持，仍然很难理解他的理论。本文汇集了作者最近在普林斯顿举办的研讨会上的演讲幻灯片，阐明了为什么贝叶斯方法不仅逻辑合理，而且易于使用。这里将以三种不同的方式实现相同的推理问题。数据我们的示例是在具有倾斜背景的嘈杂数据中寻找峰值的问题，这可能出现在粒子物理学和其他多组分事件过程中。首创数据：%matplotlibinline%configInlineBackend.figure_format='svg'importmatplotlib.pyplotaspltimportnumpyasnpdefsignal(theta,x):l,m,s,a,b=thetapeak=l*np.exp(-(m-x)**2/(2*s**2))background=a+b*xreturnpeak+backgrounddefplot_results(x,y,y_err,samples=None,predictions=None):图=plt.figure()ax=fig.gca()ax.errorbar(x,y,yerr=y_err,fmt=".k",capsize=0,label="数据")x0=np.linspace(-0.2,1.2,100)ax.plot(x0,signal(theta,x0),"r",label="Truth",zorder=0)如果样本不是None:inds=np.random.randint(len(samples),size=50)fori,indinenumerate(inds):theta_=samples[ind]ifi==0:label='Posterior'else:label=Noneax.plot(x0,signal(theta_,x0),"C0",alpha=0.1,zorder=-1,label=label)elif预测不是None:fori,predinenumerate(predictions):ifi==0:label='Posterior'else:label=Noneax.plot(x0,pred,"C0",alpha=0.1,zorder=-1,label=label)ax.legend(frameon=假)ax.set_xlabel("x")ax.set_ylabel("y")fig.tight_layout()plt.close();returnfig#随机x位置N=40np.random.seed(0)x=np.random.rand(N)#在给定的x值下评估真实模型theta=[1,0.5,0.1,-0.1,0.4]y=signal(theta,x)#仅在y方向添加异方差高斯不确定性y_err=np.random.uniform(0.05,0.25,size=N)y=y+np.random.normal(0,y_err)#plotplot_results(x,y,y_err)有了数据我们可以引入三种方法马尔可夫链蒙特CarloMarkovChainMonteCarloemcee是用纯python实现的，它只需要评估后验的对数作为参数θ的函数。对数的使用在这里很有用，因为它使指数分布族的分析评估更容易，并且因为它可以更好地处理经常出现的非常小的数字。导入emceedeflog_likelihood(theta,x,y,yerr):y_model=signal(theta,x)chi2=(y-y_model)**2/(yerr**2)returnnp.sum(-chi2/2)deflog_prior（theta）：如果所有（theta>-2）和（theta[2]>0）和所有（theta<2）：返回0return-np.infdeflog_posterior（theta，x，y，yerr）：lp=log_prior(theta)ifnp.isfinite(lp):lp+=log_likelihood(theta,x,y,yerr)returnlp#在MLE周围创建一个小球，初始化每个walkernwalkers,ndim=30,5theta_guess=[0.5,0.6,0.2,-0.2,0.1]pos=theta_guess+1e-4*np.random.randn(nwalkers,ndim)#运行emceesampler=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]我们应该经常检查生成的链，确定老化周期，需要人肉来观察平稳性：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,光栅化=True)ax.set_xlim(0,1000)ax.set_ylabel(labels[i])axes[-1].set_xlabel("stepnumber");现在我们需要细化链，因为我们的样本是相关的这里有一个计算每个参数自相关的方法，我们可以组合所有样本：tau=sampler.get_autocorr_time()print("Autocorrelationtime:",tau)mcmc_samples=sampler.get_chain(discard=300,thin=np.int32(np.max(tau)/2),flat=True)print("剩余样本：",mcmc_samples.shape)#结果自相关时间：[122.5162686675.87228105137.19550954.6357251379.0331587]Remaining(4s)emceecreatorDanForeman-Mackey还提供了这个有用的包角来可视化样本：importcornercorner.corner(mcmc_samples,labels=labels,truths=theta);虽然后验样本是推理的主要基础，但参数轮廓本身很难解释。但是使用样本生成新数据要简单得多，因为这种可视化我们对数据空间有了更多的了解。这是50个随机样本的模型评估：plot_results(x,y,y_err,samples=mcmc_samples)HamiltonianMonteCarloHamiltonianMonteCarlo梯度在高维设置中提供更多指导。为了实现一般推理，我们需要一个框架来计算任意概率模型的梯度。这里的关键部分是自动微分。我们需要的是一个计算框架，可以跟踪参数的各种运行路径。为了简单起见，我们使用的框架是jax。因为一般情况下，numpy中实现的函数可以用jax中类比代替，jax可以自动计算函数的梯度。还需要能够计算概率分布的梯度。这在几种概率编程语言中是可能的，这里我们选择了NumPyro。让我们看看如何进行自动推理：importjax.numpyasjnpimportjax.randomasrandomimportnumpyroimportnumpyro.distributionsasdistfromnumpyro.inferimportMCMC,NUTSdefmodel(x,y=None,y_err=0.1):#定义参数（包括先验范围）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))#实现模型#这里需要jaxnumpy来区分peak=l*jnp.exp(-(m-x)**2/(2*s**2))background=a+b*xy_model=peak+background#注意我们将这个采样的结果限制在观察ynumpyro.sample('obs',dist.Normal(y_model,y_err),obs=y)#需要为jax拆分密钥随机实现rng_key=random.PRNGKey(0)rng_key,rng_key_=random.split(rng_key)#使用NUTS运行HMCkernel=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()#结果如下:样本:100%|██████████|4000/4000[00:03<00:00,1022.99it/s,大小为2.08e-01的17个步长。符合。prob=0.94]meanstdmedian5.0%95.0%n_effr_hata-0.130.05-0.13-0.22-0.051151.151.00b0.460.070.460.360.571237.441.00l0.980.050.980.891.061874.341.00m0.500.010.500.490.511546.561.01s0.110.010.110.090.121446.081.00发散数：0仍然使用角来可视化Numpyro的mcmc结构：由于我们已经实现了整个概率模型（相对于emcee，我们只实现后验），我们可以直接从中创建后验预测样本下面，我们将噪声设置为零以获得纯模型的无噪声表示：fromnumpyro.inferimportPredictive#makepredictionsfromposteriorhmc_samples=mcmc.get_samples()predictive=Predictive(model,hmc_samples)#needtoset噪声为零#因为整个模型包含噪声贡献predictions=predictive(rng_key_,x=x0,y_err=0)['obs']#选择50个预测来显示inds=random.randint(rng_key_,(50,),0,mcmc.num_samples)predictions=predictions[inds]plot_results(x,y,y_err,predictions=predictions)基于模拟的推理在某些情况下，我们不能或不想计算可能性。所以我们只能得到一个模拟器（即学习模拟器的输入θ和输出D之间的映射），它可以形成似然或后验的近似代理。与产生无噪声模型的传统模拟情况的一个重要区别是，需要将噪声添加到模拟中，并且噪声模型应尽可能匹配观察到的噪声。否则我们无法区分由于噪声引起的数据变化和由于参数变化引起的数据变化。从sbi导入utils作为utilslow=torch.zeros(ndim)low[3]=-1high=1*torch.ones(ndim)high[0]=2prior=utils.BoxUniform(low=low,high=high)defsimulator(theta,x,y_err):#信号模型l,m,s,a,b=thetapeak=l*torch.exp(-(m-x)**2/(2*s**2))background=a+b*xy_model=peak+background#添加与观测一致的噪声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)。从sbi.inference.baseimportinferthis_simulator=lambdatheta:simulator(theta,torch.tensor(x),torch.tensor(y_err))posterior=infer(this_simulator,prior,method='SNPE',num_simulations=10000)NPE使用条件归一化流来学习如何在给定一些数据的情况下生成后验分布：运行10000次模拟。：0%||0/10000[00:00
