对于许多人来说，贝叶斯统计仍然有些陌生。由于贝叶斯统计中存在一些主观先验，如果没有测试数据的支持，仍然很难理解他的理论。本文汇集了作者最近在普林斯顿举办的研讨会上的演讲幻灯片，阐明了为什么贝叶斯方法不仅逻辑合理，而且易于使用。这里将以三种不同的方式实现相同的推理问题。数据我们的示例是在具有倾斜背景的嘈杂数据中寻找峰值的问题，这可能出现在粒子物理学和其他多组分事件过程中。开源：%matplotlibinline%configInlineBackend.figure_format='svg'importmatplotlib.pyplotaspltimportnumpyasnpdefsignal(theta,x):l,m,s,a,b=thetapeak=l*np.exp(-(m-x)**;2/(2*s**2))背景=a+b*xreturnpeak+背景defplot_results(x,y,y_err,samples=None,predictions=None):fig=plt.figure()ax=fig.gca()ax.errrable(x,y,yerr=y_err,fmt=".k",capsize=0,label="Data")x0=np.linspace(-0.2,1.2,100)ax.plot(x0,100;signal(theta,x0),"r",label="True",order=0)ifsamplesisnotNone:inds=np.随机的。randint(len(samples),size=50)fori,indenumenurate(inds):theta_=samples[ind]ifi==0:label='后验'else:label=Noneax.plot(x0,signal(theta_,x0),"C0",alpha=0.1,order=-1,label=label)elifpredictionsisnotNone:。fori,predenumerate(pr版本）：ifi==0:label='Posterior'else:label=Noneax.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#随机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)有了数据，我们可以引入三种马尔可夫链蒙特卡洛方法MarkovChainMonteCarloemcee是用纯python实现的，它只需要评估后验的对数作为参数θ的函数。在这里使用对数是有用的，因为它使指数分布族的分析评估更容易，并且因为它更能处理通常很好地出现的非常小的数字。importemceedeflog_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):return0return-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，进度=真）；结果如下：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,rasterized=True)ax.set_xlim(0,1000)ax.set_ylabel(labels[i])轴[-1].set_xlabel("步数");现在我们需要细化链，因为我们的样本是相关的这里有一个方法来计算每个参数的自相关性，我们可以合并所有样本：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("Remainingsamples:",mcmc_samples.shape)#ResultAutocorrelationtime:[122.5162686675.87228105137.19550954.6357251379.0331587]Remainingsamples:(4260,5)也由这个-Dankeyeman的创建者提供一个有用的包角来可视化样本：importcornercorner.corner(mcmc_samples,labels=labels,truths=theta);虽然后验样本是推理的主要基础，但参数轮廓本身很难解释。但是使用样本生成新数据要简单得多，因为这种可视化我们对数据空间有了更多的了解。以下是来自50个随机样本的模型评估：plot_results(x,y,y_err,samples=mcmc_samples)HamiltonianMonteCarloHamiltonianMonteCarlo梯度在高维设置中提供更多指导。为了实现一般推理，我们需要一个框架来计算任意概率模型的梯度。这里的关键部分是自动微分。我们需要的是一个计算框架，可以跟踪参数的各种运行路径。为了简单起见，我们使用的框架是jax。因为一般情况下，numpy中实现的函数可以用jax中类比代替，jax可以自动计算函数的梯度。还需要能够计算概率分布的梯度。这在几种概率编程语言中是可能的，这里我们选择了NumPyro。让我们看看如何进行自动推导：样本('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()#结果如下：sample:100%|██████████|4000/4000[00:03<00:00，1022.99it/s，17stepsofsize2.08e-01。符合。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.00Numberofdivergences:0还是使用corner可视化Numpyro的mcmc结构：因为我们已经实现了整个概率模型（与emcee相反，我们只实现后验），所以后验预测可以直接从样本下面，我们将噪声设置为零以获得纯模型的无噪声表示：fromnumpyro.inferimportPredictive#从后验进行预测hmc_samples=mcmc.get_samples()predictive=Predictive(model,hmc_samples)#needtosetnoisetozero#因为整个模型包含噪声贡献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之间的映射），它可以形成似然或后验的近似代理。与产生无噪声模型的传统模拟情况的一个重要区别是，需要将噪声添加到模拟中，并且噪声模型应尽可能匹配观察到的噪声。否则我们无法区分由于噪声引起的数据变化和由于参数变化引起的数据变化。importtorchfromsbiimportutilsasutilslow=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))return让我们看看噪声模拟器的输出：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,信号(theta,x0),"k",label="truth")现在，我们使用sbi从这些模拟模拟中训练神经后验估计器(NPE)。fromsbi.inference.baseimportinferthis_simulator=lambdatheta:simulator(theta,torch.tensor(x),torch.tensor(y_err))posterior=infer(this_simulator,prior,method='SNPE',num_simulations=10000)NPE使用条件归一化Streamline了解如何根据一些数据生成后验分布：运行10000次模拟。：0%||0/10000[00:00
