软件X：最大熵重构矩问题的Python软件

软件X 10（2019）100353原始软件出版物PyMaxEnt：一个用于最大熵矩重建的Python软件托尼·萨阿德1，阿，Giovanna Ruai犹他大学化学工程系，盐湖城，UT 84102，美国ar t i cl e i nf o文章历史记录：收到2019年收到修订版2019年10月21日接受2019年10月21日保留字：最大熵重构逆矩问题粒度分布a b st ra ctPyMaxEnt是一个软件，它实现了最大熵的原理，以重建给定有限数量的已知矩的函数分布。该软件支持连续和离散重建，并且通过单个函数调用非常容易使用。在本文中，我们开始验证和验证软件对几个测试，从离散概率分布的重建有偏见的骰子一直到多模态高斯和贝塔分布。PyMaxEnt用Python编写，为社区提供了一个健壮且易于使用的实现©2019作者由爱思唯尔公司出版这是CC BY许可下的开放获取文章（http://creativecommons.org/licenses/by/4.0/）中找到。代码元数据当前代码版本v1.0用于此代码版本的代码/存储库的永久链接https://github.com/ElsevierSoftwareX/SOFTX_2019_241Code Ocean compute capsule NA法律代码许可证MIT许可证使用git的代码版本控制系统使用Python 3.x的软件代码语言、工具和服务编译要求、操作环境依赖性使用numpy和scipy安装Python 3如果可用，链接到开发人员文档/手册http://github.com/saadgroup/pymaxent/问题支持电子邮件tony. chemeng.utah.edu1. 动机和意义来自实验和数值观测的数据分析通常需要后处理来重建原始物理量[1]。例如，核反应堆中的中子通量分布不容易获得，但必须根据数值或实验测量的中子矩重建[2，3]。这种原始量的重建通常受到噪声结果和不充分信息的阻碍，从而使重建问题欠确定[1]。这种情况的一个经典例子是逆矩问题[4，5]。*通讯作者。电子邮件地址：tony. chemeng.utah.edu（T. Saad）。1 助理教授https://doi.org/10.1016/j.softx.2019.100353如下所述：给定有限数量的已知时刻-例如，平均属性例如，在颗粒系统的模拟中，人们经常求解颗粒密度分布的平均尺寸和高阶矩[6]。给定这些矩，目标是重建未知的粒子密度分布。根据定义，逆矩问题是欠定的，因为有无限多个实值函数可以产生所需的矩。为了获得唯一的一个这样的约束集是通过调用最大熵原理获得的[5，7]。2352-7110/©2019作者。 由Elsevier B.V.出版。这是一篇开放获取的文章，使用CC BY许可证（http://creativecommons.org/licenses/by/4.0/）。可在ScienceDirect上获得目录列表SoftwareX期刊主页：www.elsevier.com/locate/softx2T. Saad和G. Ruai / SoftwareX 10（2019）100353⎡˜˜˜≡λkK˜˜˜˜˜≤≤+δλk =0w及k=0=k=0我δλδpi=1⎢∂λ01∂λ11λk⎢˜µ1˜µ2···˜µk+1⎥在考虑所有已知信息的情况下，具有最大信息熵的系统是最有可能的µk=这个泛函是一个最大值，当函数的导数x i p i.（十二）E01K∑=∑（∑）∑+我1.1. 最 大熵，则雅可比矩阵由下式给出：最大熵方法基于以下概念：使信息熵最大化的分布是这在统计上是最有可能发生的[5]。信息熵被定义为信息产生的平均速率∂µ0∂λ0⎢∂˜µ∂µ0∂λ1∂˜µ···∂˜µ0⎤∂˜µ⎥∂λ⎡˜µ0˜µ1···˜µk⎤一个系统。信息熵也可以定义为：Jµ2⎢∂λ∂µ2∂λ···µ2=µ2µ3···µk+2。（九）一个衡量系统有多少是未知的。 越高- 是的... 没关系.+.. ⎥⎦信息熵越小，我们对过程的了解就越少好吧... ⎥⎦克雷奇克1···kk状态，因为它是定义了最少信息量的系统。在数学上，分布p（x）的信息熵S由以下积分给出：S=−p（x）lnp（x）dx，（1）其中，k是分布的支持。给定p（x）的有限个矩，我们寻找：一个分布p（x），它使S在这些已知的矩下最大化。这可以表述为一个变分问题，拉格朗日乘数。我们的目的是找到p（x），使（1）中给出的信息熵S最大化，注意，µk是基于重构分布的第k 牛顿求解器将基于找到λk，使得||µk−µk||低于某个特定的公差。1.3. 初始猜测牛顿法可能对初始猜测很敏感。结果发现，如果使用以下初始猜测，则在我们所有的实验中都可以实现收敛{−ln<$2πi=0，∫xkp （ x ） dx=µk;k=0 ， 1 ， ... ， N ，（2）其中（N1）是已知时刻的个数。注意已知μk为0k N。然后，通过引入拉格朗日乘子λk，我们定义了熵泛函这些猜测基于具有零均值和单位方差（μ=0，σ=1）的高斯分布。1.4. 离散情况类似的分析可以应用于离散的情况。用于一组H<$S+∑λk（<$xk p（x）dx−µk），（3）我们首先将离散信息熵定义为：或者更准确地说，NS= −∑pilnpi，（11）NH[−p（x）lnp（x）+∑λkxkp（x）]dx−∑ λkµk。（四）我其中k阶矩由下式给出：k=0k=0∑kH关于p（x）和λk为零：δHδHδp（x）我与连续情形类似，我们定义熵泛函0. 这些导数中的第一个导数没有为我们提供任何额外的信息，因为它返回在（2）中定义的约束。然而，这些导数中的第二个估计为：N作为H S+∑λ（∑xk p−µ）lnp（x）= − 1 + ∑λkx k。（五）KiIkK iN（）k=0（5）的一般解是=− ∑pi ln pi + ∑λk ∑xkpi −µk、（十三）p（x）=e−1+∑Nλkxke∑Nλkxk，（6）我或者更准确地说，k=0i其中我们设置（λ0−1）→ λ0。1.2. 溶液Nik=0Nλkxk pi−pilnpi−k=0λkµk。（十四）为了找到最大熵解，我们必须求解拉格朗日乘子λk。为此，我们必须解决如果关于pi和λk的泛函导数为非线性方程零：δHK=0和δH=0第一个方程是-∫∫Ω eλ0+λ1x+···+λkxkdxµ0，λ+λx+···+λxk.Ω我···克∂λ0···克∂λ1克λkλ初始值=0否则。（十=dx=µ1，λkxk−1− lnpi= 0，（101ΩΩHT. Saad和G. Ruai / SoftwareX 10（2019）1003533Ω01Kk=0我我k=0我转动由式（12）给出的力矩约束。第二个函数导数返回（七δHδpiN∫xkeλ+λx+···+λxkdx=k=0它对pi的解很简单可以使用全局收敛的牛顿求解器来计算λk。p=e∑Nλkxk−1=e∑Nλkxk，（16）基于最大熵解的矩，其中我们设置λ0−1→λ0。与连续情形类似，克雷蒂安Ωxkeλ0+λ1x+···+λkxkdx，（8）全局收敛的牛顿求解器可用于求解拉格朗日乘子λk。我k.雅可比矩阵也很容易计算。如果我们把4T. Saad和G. Ruai / SoftwareX 10（2019）100353因此，∑∑μ=（px）=（x）=3。5. 要运行此测试，请执行以下操作：2. 软件描述Fig. 1. PyMaxEnt功能概述。2.3. 示例代码段分析该软件是用Python编写的，因为它的流行性和易用性，以及通过SciPy提供的强大的多维非线性求解器。PyMaxEnt软件为连续和离散最大熵重建提供了一个单一的接口，以及下面描述的一些有用的功能。2.1. 软件构架PyMaxEnt 的 实 现 非 常 简 单 ， 并 且 在 一 个 名 为pymaxent.pyPython文件中形成。为了正确的操作，三个独立的助手例程实现连续和离散的情况下。这些是一个数值积分器，一个残差计算器，和牛顿求解器，如图所示。1.一、牛顿求解器基于SciPy从用户的角度来看，对主例程进行单个函数调用，2.2. 软件功能要使用PyMaxEnt，只需调用函数reconstruct。函数签名为sol，bndas = reconstruct（moments，rndvar，bnds）这里，moments是一个必需的已知矩列表或数组，rndvar是一个可选参数，包含随机变量的离散值，bnds是一个元组[a，b]，包含结果分布的预期边界。当提供rndvar时，重建假设离散分布。代码返回两个量：（1）解的函数形式，即（6），可以绘制，（2）包含拉格朗日乘子的numpy数组。在离散情况下，函数解只是一个包含对应于（16）的函数值的numpy值数组。最后，PyMaxEnt提供了一个名为moments的助手例程，用于计算函数的前k个矩，并用于验证目的。下面是一些使用PyMaxEnt软件的例子。从一个例子开始，用两个矩和随机变量x的六个可能值来重建基本离散分布，从pymaxent导入* mu =[1，3.5]x =[1，2，3，4，5，6]sol，rndas = reconstruct（mu，rndvar=x）类似地，对于连续分布，传递输入矩的列表。然而，在这种情况下，必须指定边界（bnds）以指示这是连续重建。这从pymaxent导入* mu =[1，0，0.04]sol，bndas = reconstruct（mu，bnds=[-1，1]）# plot重建的解决方案x = np.linspace（-1，1）plot（x，sol（x））3. 说明性实例为了测试和验证PyMaxEnt，我们开始研究一些概率分布的重建，从离散到连续的复杂分布，具有无限的支持。3.1. 验证3.1.1. 平衡模对于一个平衡的骰子，面值在1到6之间，一面着地的概率是1/6。要求的数值是16PyMaxEnt，我们只需设置以下内容frompymaxentimport reconstruct mu= [1，3.5]#设置矩#指定随机变量x =[1，2，3，4，5，6]sol，Rndas = reconstruct（mu=mu，rndvar=x）打印（溶胶）软件返回的sol值为[0.166667 0.166667 0.166667 0.166667 0.166667]其是平衡模具的期望均匀概率分布。T. Saad和G. Ruai / SoftwareX 10（2019）1003535----==µµ1，精确√√√12π++图二. 最 大熵重建的双峰高斯分布与μ00。25，σ0 1/14，µ10。5和σ11/ 20。(a)实例重建和（b）作为输入矩数量的函数的重构中的面积误差图3.第三章。三 峰高斯分布µ014、µ12/4、µ23/4、σ01/14和σ1σ21/20的最大熵重建。(a)示例重建和（b）作为输入矩数量的函数的重建中的面积误差3.1.2. 偏置模具考虑具有以下概率的偏置裸片：p（2）=p（3）=1/ 12，p（4）=2/ 12，p（5）=3/ 12和p（6）=4/ 12。期望值为µ1= 4。417.为了运行此测试，我们使用frompymaxentimport reconstruct mu= [1，4.417]#设置矩#指定随机变量x =[1，2，3，4，5，6]sol，Rndas = reconstruct（mu=mu，rndvar=x）打印（溶胶）如[8]中所述然后，可以通过输入曲线的面积（对于概率密度函数，在所有情况下都是1）来归一化该差请注意，对于所有情况，计算的矩与精确矩的误差最大为1 .一、49× 10−8，这是牛顿解算器的默认公差。3.2.1. 多模态高斯分布我们考虑由下式给出的双峰高斯分布：f（x）=1Exp[−（x−µ0）2]+1Exp[−（x−µ1）2]0软件返回的sol值为[0.0614560 0.0859580 0.120229 0.168163 0.235209 0.328985]2σ0√2π2σ22σ1√2π2σ12.（十七）这对应于非均匀概率分布，偏置裸片，裸片落在较小值上的可能性较高换句话说，落在数字1上的概率由重构解p（1）0给出。0614560，而落在2上的概率由p（2）0给出。0859580等。当然，正如预期的那样，在计算的概率中存在误差，因为非线性求解器作用于分布的矩而不是分布本身。例如，在这个图图2a示出了双峰高斯与不同高度的峰和分别使用5和10个矩的最大熵相关。显然，已知时刻的数量越多，反射就越精确图在图2b中，示出了在重建中产生的面积误差与增加数量的矩的关系。一个明显的下降，在错误是可见的高达10个时刻，错误中的行为可见。这是因为，在在这种情况下，矩的相对误差为µ0，精确−µ0，计算0，精确10−13 和µ1，精确−µ1，计算 -10-73.2. 预测性和误差如此大量的矩、圆误差在数值积分和非线性求解器中占主导地位。重建的性质包含一个多项式的指数，这反过来又导致了预期的增加rounding错误。对于三峰高斯分布，我们使用为了测试PyMaxEnt的预测性，我们选择了一些非平凡的分布，为这些分布生成了矩，（x−µ0）22σ2（x−µ1）22σ2（x−µ2）22σ2然后，当我们改变输入到软件中的矩的数量时，查看重构分布的形状相对于所选分布形状的误差f（x）=3σe0e3σ1 2π1e3σ2 2π2.（十八）根据两条曲线之间的差异面积计算重建分布与精确分布之间的差异重建的分布如图所示。其中5和13阶矩用于重建精确分布。在这−−−≈1106T. Saad和G. Ruai / SoftwareX 10（2019）100353=-=-=[客户端]图四、 具有α3和β9的β分布的最大熵重建。(a)示例重建和（b）作为函数的重建中的面积误差输入矩的数量图五. 对数正态分布的最大熵重建，μ0。2和σ0。25. (a)实例重建和（b）重建中的面积误差作为输入矩的数量的函数在这种情况下，需要高阶矩来捕获该分布的凸性的变化。图图3b显示了随着我们增加矩的数量，同样，我们看到圆形误差在第12时刻占主导地位。3.2.2. Beta分布贝塔分布可能是一个具有挑战性的分布重建。它是由4. 影响本文提出的方法已用于核反应堆中子注量率分布的重建[2，3]。此外，它还被杨百翰大学和犹他大学的同事和合作者用于未发表的工作，用于重建二氧化碳矿化反应器[6]和煤燃烧模拟中新的研究问题可以在以下背景下进行：f（x） Γ（α+β）xΓ（α）Γ（β）α−1（1−x）β−1，（十九）高性能计算用于反应颗粒传输，如煤燃烧、气化、流化床和核反应堆。传输单个粒子（或粒子云其中，Γ是广义Gamma函数，α和β是决定β分布形状的控制参数这种情况的结果如图所示。四、最大熵重建能够捕捉到这种分布的关键特征，特别是对于在域边界处具有奇异性的情况。3.3. 半无限支集当在无限或半无限支撑上积分时，数值积分器可能会遇到麻烦这是观察到的，特别是与分布，慢慢衰减到无穷大，如对数正态分布。然而，如果在有限支撑上进行积分，仍然能够恢复对分布的合理近似。作为一个例子，使用对数正态分布，在区间0， 10上进行数值积分。由此产生的重建如图所示。 五、是一个非常昂贵的命题，人们经常诉诸的方法，解决这些粒子分布的时刻。重建母颗粒尺寸分布的能力使得能够详细了解这些系统，例如基于其尺寸和成分了解煤颗粒如何在煤燃烧器中移动。该代码还将用于我们团队正在开展的一个新项目，该项目涉及大气中有害颗粒物（PM2.5）的建模以及如何使用无人机准确据我们所知，我们不知道有类似于PyMaxEnt的公开软件。我们相信，将PyMaxEnt提供给社区是朝着正确方向迈出的一步。5. 结论在这篇文章中，我们开始演示如何使用Python软件，利用最大熵原理从给定的时刻软件，T. Saad和G. Ruai / SoftwareX 10（2019）1003537PyMaxEnt被证明可以重建具有有限，半无限和无限支持的离散和连续分布。由于舍入误差的存在，使用大量矩时必须小心。牛顿解算器在所有测试中均表现出稳健性，而不会改变求根算法的初始猜测。如有必要，用户可以很容易地改变非线性求解器的初始猜测。未来的工作将包括支持更高精度的浮点数，以减少舍入错误的数量。竞合利益作者声明，他们没有已知的竞争性财务利益或个人关系，可能会影响本文报告的工作致谢作者非常感谢来自犹他大学研究生办公室项目（UROP）和化学工程系的引用[1]von der Linden W.最大熵数据分析。应用物理A1995;60（2）：155-65.[2]放大图片作者：Crawford Douglas S，Saad Tony，用MCNP 5、Attila-7.1.0和GODIVA实 验验 证 和 确认 重 建中 子 通 量的 最 大熵 方 法 。 Ann Nucl Energy2013;53：188[3]克劳福德·道格拉斯·斯图尔特，给特里·亚瑟打电话.用MCNP 5和Attila-7.1.0验证一维含能中子扩散方程的解析能量矩2012年。[4]詹姆斯·亚历山大将军塔玛金·雅各布·大卫问题的时刻，卷。1.一、美国数学学会一九四三年[5]Mead Lawrence R ， Papanicolaou Nikos. 矩 问 题 中 的 最 大 熵 。 J MathPhys1984;25（8）：2404-17.[6]Abboud Alex W，Schroeder Ben B，Saad Tony，Smith Sean T，Harris DerekD，Lignell David O.大涡模拟与一维湍流对湍流降水的数值比较。AIChEJ2015;61（10）：3185-97。[7]Bandyopadhyay K，Bhattacharya Arun K， Biswas Parthapratim ，DraboldDA. 最大熵和矩问题：一个稳定的算法。PhysRev E2005;71（5）：057701。[8]Jekel Charles F，Venter Gerhard，Venter Martin P，Stander Nielen，HaftkaRaphael T. 用 逆 分 析 从 磁 滞 回 线 识 别 材 料 参 数 的 相 似 性 度 量 。 Int J MaterForm2019;12（3）：355-78.