Dirichlet过程混合模型在Pyro中的应用

贝叶斯非参数模型是参数数量随提供的数据量自由增长的模型；因此，不是训练几个复杂度不同的模型并对它们进行比较，而是能够设计一个其复杂度随着观察到的更多数据而增长的模型。实践中贝叶斯非参数的典型例子是Dirichlet过程混合模型(DPMM)。DPMM允许从业者在其数据的几何结构中的不同簇的数量未知的情况下构建混合模型-换句话说，允许随着观察到的数据的增多而增长簇的数量。这一特性使得DPMM对于探索性数据分析非常有用，在探索性数据分析中，我们几乎不知道有问题的数据的各个方面；本演示文稿旨在演示这一事实。

Dirichlet过程是离散概率分布上的一族概率分布族。形式上，Dirichlet过程(DP)由某种基本概率分布\(G_0：\ω\to\mathbb{R}\)和一个通常表示为\(\α)的正的实数标度参数来指定。参数为\(G_0：\omega\to\mathbb{R}\)和\(\α\)的Dirichlet过程的样本\(G\)本身是\(\Ω\)上的分布。对于\(\Omega\)的任何不相交分区\(\Omega_1，.，\Omega_k\)和任何样本\(G\sim DP(G_0，\alpha)\)，我们有：

本质上，这是对样本空间\(\Ω\)进行离散划分，然后使用基础分布\(G_0\)在其上构建离散分布。虽然Dirichlet过程在公式上相当抽象，但在各种图形模型中作为先验是非常有用的。在下面的方案中，这一事实变得更容易看到。

想象一下，一家餐厅有无限张桌子(由正整数索引)，一次只接受一位顾客。第(N)位客户根据以下概率选择座位：

以概率\(\frac{n_t}{\alpha+n-1}\)坐在\(t\)桌，其中\(n_t\)是在桌\(t\)的人数。

如果我们从欧米茄上的基分布(G_0)中关联一个抽签给每个表\(t\)，然后将非正规化概率质量\(n_t\)关联到该抽签上，则得到的欧米茄上的分布等价于狄利克雷过程的抽签(DP(G_0，α)\)。

此外，我们可以很容易地将其扩展到定义非参数混合模型的生成过程：每个表\(T)至少有一个客户席位与一组聚类参数\(\theta_t\)相关联，这些参数本身来自于某个基分布\(G0\)。对于每个新的观测值，首先根据上述概率将该观测值分配到一个表中；然后，从该表的聚类参数参数化的分布中提取该观测值。如果观测被分配到新的表，则从\(G_0\)中绘制一组新的群集参数，然后从由这些群集参数参数化的分布中绘制观测。

虽然Dirichlet过程混合模型的这种表述是直观的，但在概率规划框架中也很难对其进行推理。这激发了DPMM的另一种表述，经验表明，这种表述更有利于推理(例如，Blee和Jordan，2004年)。

对于每个观测值\(n\in\{1，.，N\}\)，绘制\(z_n\sim\pi(\beta_{1：\infty})\)，然后绘制\(x_n\sim f(\theta_{z_n})\)。

这里，可以更容易地看到Dirichlet过程混合模型的无限本质。此外，所有的\(\beta_i\)都是独立的，因此在概率编程框架中执行推理要容易得多。

从Torch.Distributions导入constraintsimport torch.nn.function as fimport熊猫作为pdimport numpy作为np导入matplotlib.pylot作为pltimport matplotlib.cm作为cmfrom tqdm导入tqdm从pyro.分发导入*从pyro.optim导入Adam从pyro.infer导入svi，trace_elbo，Predictiveassert pyro.__version__.startswith('；1'；)pyro.enable_validation(True)#可以帮助调试pyro.set_。

采集pyro-ppl下载https://files.pythonhosted.org/packages/75/c8/3a52883ab27c5503385a983e4451a1e2944c2baab66c0f7ffc1f84bca329/pyro_ppl-1.1.0-py3-none-any.whl(429kB)|█|430kB 3.4MB/s采集pyro-api>；=0.1.1下载https://files.pythonhosted.org/packages/c2/bc/6cdbd1929e32fff62a33592633c2cc0393c7f7739131ccc9c9c4e28ac8dd/pyro_api-0.1.1-py3-none-any.whlRequirement已满足：Torch&>；=1.3.0 in/usr/local/lib/python3.6/dist-Packages(来自pyro-ppl)(1.3.1)采集tqdm>；=4.36下载https://files.pythonhosted.org/packages/7f/32/5144caf0478b1f26bd9d97f510a47336cf4ac0f96c6bc3b5af20d4173920/tqdm-4.40.2-py2.py3-none-any.whl(55KB)|█|61KB 4.2MB/s已满足需求：OPT-EINSUM>；=2.3.2 in/usr/local/lib/python3.6/dist-Packages(来自pyro-ppl)(3.1.0)已满足要求：numpy>；=1.7 in/usr/local/lib/python3.6/dist-Packages(from pyro-ppl)(1.17.4)安装收集的软件包：pyro-api，tqdm，pyro-ppl找到现有安装：tqdm 4.28.1卸载tqdm-4.28.1：成功卸载tqdm-4.28.1成功安装pyro-api-0.1.1 pyro-ppl-1.1.0 tqdm-4.40.2。

我们首先在由四个2D高斯混合生成的合成数据集上演示Dirichlet过程混合模型的能力：

data=torch.cat((MultivariateNormal(-8*torch.ones(2)，torch.ye(2)).sample([50])，MultivariateNormal(8*torch.ones(2)，torch.ye(2)).sample([50])，MultivariateNormal(torch.tensor([1.5，2])，torch.ye(2)).sample([50])，MultivariateNormal(torch.tensor(2)，torch.tensor(。0]，数据[：，1])plt.title("；来自4个高斯混合()plt.show()N=data.form[0]的数据样本。

在此示例中，群集参数\(\θ_i\)是描述具有单位协方差的多变量高斯的平均值的二维向量。因此，Dirichlet过程的基分布(G_0)也是多元高斯分布(即共轭先验分布)，尽管这种选择在计算上没有那么有用，因为我们不是使用Pyro进行坐标上升变分推断，而是使用Pyro进行黑箱变分推断。

首先，让我们定义“折断”函数来生成我们的权重，给定我们的样本\(\beta\)：

def Mix_Weight(Beta)：beta1m_umprod=(1-beta).umprod(-1)返回F.pad(beta，(0，1)，value=1)*F.pad(beta1m_umprod，(1，0)，value=1)。

接下来，让我们定义一下我们的模型。参考本教程第一部分中介绍的折杆模型的定义可能会有所帮助。

注意，所有\(\beta_i\)样本在条件上是独立的，因此我们使用大小为T-1的pyro.plate对它们建模；我们对集群参数\(\mui\)的所有样本执行相同的操作。然后，我们使用下面的采样\(\beta\)值(第9行)构造一个分类分布，其参数是混合权重，并对该分类中的每个数据点的聚类分配\(z_n\)进行采样。最后，我们从多元高斯分布中抽样，其平均值正好是对应于我们为点\(x_n\)画出的赋值\(z_n\)的簇参数。这可以在下面的Pyro代码中看到：

def model(Data)：with pyro.plate("；beta_plate"；，T-1)：beta=pyro.sample("；beta"；，Beta(1，alpha))with pyro.plate("；mu_plate"；，T)：MU=pyro.sample("；mu"；，MultivariateNormal(torch.zeros(2)，5*torch.ye(2)，带pyro.plate("；beta_plate"；，T-1)：beta=pyro.sample("；beta"；，5*torch.ye(2))。，N)：Z=pyro.sample("；z"；，范畴(Mix_Weight(Beta)pyro.sample("；obs"；，MultivariateNormal(mu[z]，torch.ye(2))，obs=data)。

我们在变分推理过程中优化的变分族\(q(\β，\θ，z)\)由下式给出：

注意，由于我们不能对模型假设的无限团簇进行计算建模，所以我们在(T)团簇上截断了我们的变分族。这不会影响我们的模型；相反，它是在推理阶段进行的简化，以便于处理。

该向导完全按照上面我们的变分族\(q(\β，\θ，z)\)的定义来构造。对于我们模型中的每个样本，我们都有条件独立的Beta分布，对于每个簇参数，都有T-1个条件独立的多元高斯分布，对于每个簇分配，都有N个条件独立的分类分布。我们的模型中的每个β样本都有条件独立的Beta分布，每个簇参数都有条件独立的多元高斯分布，每个簇分配都有N个条件独立的分类分布。

因此，我们的变分参数(pyro.param)是参数化我们的变分Beta分布的第二个参数的\(T-1\)许多正标量(第一个形状参数固定在\(1\)，如在模型定义中一样)，以及参数化我们的变分多元高斯分布的\(T\)许多二维向量(我们不参数化高斯的协方差矩阵，尽管这应该在分析真实世界的数据集时进行以获得更大的灵活性)，而我们的变分参数(pyro.param)是参数化我们的变分Beta分布的第二个参数的\(T-1\)许多正标量(第一个形状参数固定在\(1\)，如在模型定义中一样)，以及

定义向导(Data)：kappa=pyro.param('；kappa'；，lambda：Uniform(0，2).sample([T-1])，Constraints.Constraints.Positive)tau=pyro.param('；tau'；，lambda：MultivariateNormal(torch.zeros(2)，3*torch.ye(2)).sample([T]))phi=pyro.param('；tau'；，lambda：MultivariateNormal(torch.zeros(2)，3*torch.ye(2)).sample([T]))phi=pyro.param。，lambda：Dirichlet(1/T*torch.ones(T)).sample([N])，Constraint=constraints.simplex)with pyro.plate("；beta_plate"；，T-1)：q_beta=pyro.sample("；beta"；，Beta(torch.ones(T-1)，kappa))with pyro.plate("；mu_plate"；，T)：q。，MultivariateNormal(tau，torch.ye(2))，带pyro.plate("；data"；，N)：Z=pyro.sample("；z"；，分类(Phi))。

在执行推理时，我们将数据集中最大聚类数的“猜测”设置为\(T=6\)。我们定义了优化算法(pyro.Optim.Adam)和Pyro SVI对象，并训练模型进行1000次迭代。

在进行推理之后，我们构造均值(我们的变分近似中每个因素的期望值)的贝叶斯估计，并在原始数据集上用红色绘制它们。请注意，我们还根据学习到的变分分布删除了分配给它们的低于特定权重的任何群集，然后重新规格化这些权重，使它们的总和为1：

t=6optim=Adam({"；lr"；：0.05})SVI=SVI(model，Guide，Optim，Loss=Trace_Elbo())Loss=[]def Train(Num_Iterations)：pyro.clear_param_store()for j in tqdm(range(Num_Iterations))：Loss=svi.step(数据)Lost。append(Loss)def truncate(alpha，center，Weight)：Threshold=。TRUE_CENTERS=CENTERS[Weights&>Threshold]TRUE_WEIGHTS=Weights[Weights&>Threshold]/torch.sum(Weights[Weights&>Threshold])返回TRUE_CENTERS，TRUE_Weight tsalpha=0.1Train(1000)#我们使用中心和权重的后验平均值tau和phi对模型参数进行点估计Bayes_Centers_01，Bayes_Weight_01=Truncate(alpha，pyro.param("；tau。).Detach()，dim=0)α=1.5Train(1000)#我们使用中心和权重的τ和φ的后验均值对我们的模型参数进行点估计Bayes_Centers_15，Bayes_Weight_15=Truncate(alpha，pyro.param("；tau"；).Detach()，torch.means("；phi"；pyro.param("；phi"；).DETACH()，dim=0)plt.Figure(FigSize=(15，5))plt.subploy(1，2，1)plt.Scatter(Data[：，0]，Data[：，1]，color="；blue"；)plt.scatter(Bayes_Centers_01[：，0]，Bayes_Centers_01[：，1]，color="；红色"；)plt.subploy(1，2，2)plt.Scatter(Data[：，0]，Data[：，1]，color="；)plt.sublot(1，2)plt.dister(Data[：，0]，Data[：，1]，color="。蓝色"；)plt.Scatter(贝叶斯_中心_15[：，0]，贝叶斯_中心_15[：，1]，color="；red"；)plt.tight_layout()plt.show()。

100%|█|1000/1000[00：06<；00：00,153.93it/s]100%|█|1000/1000[00：06<；00：00,152.39it/s]。

上面的曲线图显示了缩放超参数\(\α\)的影响。较大的\(\alpha\)会产生更厚尾的权重分布，而较小的\(\alpha\)会在较少的星团上放置更多质量。特别地，中间的群集看起来可以生成单个高斯(尽管实际上它是由两个不同的高斯生成的)，因此\(\α\)的设置允许实践者进一步编码他们关于数据包含多少群集的先验信念。

如前所述，当探索潜在几何结构完全未知的数据集时，Dirichlet过程混合模型确实大放异彩。为了证明这一点，我们根据过去300年的太阳黑子计数数据(由比利时皇家天文台提供)拟合了DPMM：

df=pd.read_csv('；http://www.sidc.be/silso/DATA/SN_y_tot_V2.0.csv'；，9月='；；'；，名称=['；时间'；，'；太阳黑子。年份'；]，用法=[0，1])Data=torch.tensor(df['；sunspot.year'；].values，dtype=torch.flat32)N=data.Shape[0]plt.hist(df['；太阳黑子.年份&39；].value，bin=40)plt.title("；年数与太阳黑子计数"；)plt.xlabel("；太阳黑子计数"；)plt.ylabel("；年数"；)plt.show()

在本例中，集群参数\(\theta_i\)是速率参数，因为我们构建的是泊松分布的比例混合。同样，选择\(G_0\)作为共轭先验，在本例中是伽马分布，尽管这对于通过Pyro进行推理仍然没有严格意义。以下是该模型的实现：

def model(Data)：with pyro.plate("；beta_plate"；，T-1)：beta=pyro.sample("；beta"；，beta(1，alpha))with pyro.plate("；lambda_plate"；，T)：lmbda=pyro.sample("；lambda"；，Gamma(3，0.05))with pyro.plate("；lambda"；，Gamma(3，0.05))with pyro.plate("；lambda"；，Gamma(3，0.05))with pyro.plate("；lambda"；，T)。，分类(Mix_Weight(Beta))pyro.sample("；obs"；，Poisson(lmbda[z])，obs=data)def Guide(Data)：kappa=pyro.param('；kappa'；，lambda：均匀(0，2).sample([T-1])，Constraint=constraints.Positive)tau_0=pyro.param('；tau_0&#。tau_1'；，lambda：LogNormal(-1，1).sample([T])，Constraint=Constraints.Positive)Phi=pyro.param('；Phi'；，lambda：Dirichlet(1/T*torch.ones(T)).Sample([N])，Constraint=Constraints.simplex)带焦板("；beta_plate"；，T-1)：Q_beta。，Beta(torch.ones(T-1)，kappa)，带烟盘("；lambda_plate"；，T)：q_lambda=pyro.sample("；lambda"；，Gamma(tau_0，tau_1))带烟盘("；data"；，N)：Z=pyro.sample("；z"；，classical(phi。：0.0 5})svi=svi(模型，导轨，最优，损耗=跟踪_ELBO())损耗=[]列车(N_ITER)样本=火炬.arange(0，30 0).type(火炬.浮动)tau0_Optimal=pyro.param("；tau_0"；).detach()tau1_optimal=pyro.param("；tau_1"；).detach()kappa_optimal=pyro.param("；kappa"；).DETACH()#对于群集参数和权重，我们使用τ和kappa的后验均值对我们的潜在变量进行点估计Bayes_Rates=(τ0_最优/τ1_最优)Bayes_Weight=Mix_Weight(1./(1.+kappa_Optimal))def Mixine_of_Poisson(权重，速率，样本)：Return(权重*Poisson(rates).log_prob(samples.unsqueeze(-1)).exp()).sum(-1)likelihood=Mixture_of_Poisson(Bayes_Weight，Bayes_Rates，样本)plt.title("；太阳黑子计数的年数vs.)plt.hist(data，bin=60，Density=True，lw=0，alpha=0.75)；plt.lot(样本，可能性，标签=#34；估计的混合密度)plt.Legend()plt.show()。

100%|█|1500/1500[00：09<；00：00,160.10it/s]。

上图是集群参数的贝叶斯估计器的混合密度，按其相应的权重加权。与高斯例子一样，我们通过分别计算λ和β的后验均值，得到了每个聚类参数的贝叶斯估计及其相应的权重。

下面是使用PYRO进行推理期间SVI迭代上损失函数(负Trace_ELBO)行为的曲线图，以及ELBO‘时间序列’与迭代次数的自相关曲线图。我们可以看到，大约500次迭代，损失停止显著减少，所以我们可以假设需要大约500次迭代才能实现收敛。自相关图在滞后500左右达到非常接近于0的自相关，进一步证实了这一假设。请注意，这些都是试探性的，并不一定意味着收敛。

ELBO_PLOT=plt.Figure(FigSize=(15，5))ELBO_AX=ELBO_PLOT.ADD_SUBPLOT(1，2，1)ELBO_AX.SET_TITLE("；太阳黑子Data"；)elbo_ax.set_ylabel("；ELBO"；)elbo_ax.set_xlabel("；Iteration号"；)elbo_ax.lot(np.arange(N_Iter)，亏损)autocorr_ax=elbo_plot.add_sublot(1，2，2)autocorr_ax.acorr(np.asarray(亏损)，deTrend=lambda x：x-x.means()，maxlags=750，usevline=false，mark='；，'；)autocorr_ax.set_xlim(0,500)autocorr_，lw=1)autocorr_ax.set_title("；太阳黑子Data"；)autocorr_ax.set_xlabel("；Lag"；)autocorr_ax.set_ylabel("；Autocorrelation"；)elbo_plot.tight_layout()elbo_plot.show()上火红BBVI的ELBO与LAG的自相关。

由于我们计算了适合于太阳黑子长期数据的DPMM的近似后验，我们可以利用一些内在的度量，如对数预测、后验离散指数和后验预测检验。

由于Dirichlet过程混合模型的后验预测分布本身是具有解析近似的比例-混合分布(Blee和Jordan，2004)，这使得它特别服从上述度量：

特别是，为了计算对数预测，我们首先使用数据的训练子样本对我们的模型进行变分推断后计算后验预测分布(如上定义)。因此，对数预测是在测试子样本中的每个点评估的预测密度的对数值：

\[\log p(x_{new}|X)=\log\mathbb{E}_{\beta，\theta|X}\Left[p(x_{new}|\beta，\theta)\right]\]。

由于训练样本和测试样本都来自同一数据集，因此我们预计该模型会为测试样本分配高概率，尽管在推理过程中没有看到它们。这给出了一个度量，用来选择我们的超参数\(T\)、\(\alpha\)和\(G_0\)的值：我们希望选择最大化此值的值。

我们在下面使用不同的\(\alpha\)值执行此过程，以查看最佳设置是什么。

#保留原始数据的10%来测试upondf_test=df.sample(frac=0.1)data=torch.tensor(df.drop(df_test.index)['；sunspot.year'；].values，dtype=torch.loat)data_test=torch.tensor(df_test['；sunpot.year'；].values，dtype=torch.Float)N=data.form[0]N_test=data_test.shape[0]alphas=[0.05，0.1，0.5，0.75，0.9，1.，1.25，1.5，2，2.5，3]alpha中的val：alpha=val T=20 SVI=SVI(model，Guide，Optim，Loss=Trace_Elbo()系列(500)S=100#蒙特卡罗样本数。在后验预测类中构建：后验=预测性(GUIDE，Num_Samples=S，Return_Sites=["；beta&34；，"；lambda&34；])(数据)post_pred_weights=Mix_Weight(后验["；beta&34；])post_pred_cluster=后验["；lambda"；]#LOG_PROB Shape=N_TEST x S LOG_PROB=(POST_PROD_Weight ts.log()+Poisson(post_pred_clusters).log_prob(data.reshape(-1，1，1)).logsumexp(-1)Mean_LOG_PROB=LOG_PROSE.logsumexp(-1)-np.log(S)LOG_PORTER_PRODUCTION=Mean_LOG_PROT.SUM(-1)LOG_Predictives.append(LOG_PORTER_PROFITVES.SUM。

..