Evaluating crystallographic likelihood functions using numerical quadratures

Zwart, P.H.; Perryman, E.D.

doi:10.1107/S2059798320008372

研究论文

结构
生物学

编号：2059-7983

第76卷| 第8部分| 2020年8月| 第736-750页

https://doi.org/10.107/S2059798320008372

打开

访问

用数值求积计算晶体似然函数

彼得鲁斯·兹瓦特 ^a、，^b条 ^*和埃利奥特·佩里曼 ^a、，^b、，^c（c）

^一美国加州伯克利回旋加速器路1号劳伦斯伯克利国家实验室计算研究部能源研究应用高等数学中心，邮编：94720，^b条美国加州伯克利回旋加速器路1号劳伦斯伯克利国家实验室分子生物物理和综合生物成像部，邮编94720^c（c）田纳西大学诺克斯维尔分校，诺克斯维尔，田纳西州37916，美国
^*通信电子邮件：电话：phzwart@lbl.gov

英国剑桥大学R.J.Read编辑(收到日期：2020年1月13日； 2020年6月23日接受；在线2020年7月27日)

晶体学应用中基于强度的似然函数有可能提高从边缘衍射数据导出的结构的质量。然而，由于能够有效计算这些目标函数，它们的使用变得复杂。在这里，开发了一种数值求积，它允许在晶体学应用中快速评估基于强度的似然函数。通过使用一系列可变变换，包括非线性域压缩操作，构造了一个精确、稳健和高效的求积。该方法具有灵活性，可以相对容易地合并不同的噪声模型。

关键词：最大似然;精细化;数值积分.

类似文章

1.简介

根据实验观测结果估计模型参数在自然科学中起着核心作用，并且基于相似性的方法的使用已被证明能够对“最佳猜测”值及其相关置信区间进行稳健估计（Rossi，2018 ).最大似然这一估计可追溯到19世纪高斯（Gauss，1809）的零星使用 , 1816 , 1823 )和黑根（1867年 )，并由费希尔（1915）进一步发展 )威尔克斯（1938）)，以及奈曼和皮尔逊（奈曼和斯科特，1948年 ; 皮尔逊，1970年 ). 在晶体学界，Beu等。(1962 )是第一个明确使用最大似然估计，将其应用于格参数精细化在粉末衍射中。在对这项工作的后期反应中，威尔逊（1980 )声明`使用最大似然是不必要的，并且可能会遭到一些反对”，然后重新定义了Beu的工作等。(1962)需要注意的是，在随机变量正态性假设下，最小二乘估计方法等价于似然形式。在对大分子分析产生重大影响后，结构科学中使用非正态分布的基于极大似然的方法开始流行。对于这些类型的样品，结构解决方案和精细化由于启动模型非常不完整或质量低，问题往往很成问题，使得标准最小二乘法表现不佳。在20世纪80年代和90年代，基于likelihood的方法成为主流，最终形成了常规确定和完善以前被认为有问题的结构的能力（Lunin&Urzhumtsev，1984 ; 里德，1986年 ; Bricogne和Gilmore，1990年 ; de La Fortelle&Bricogne，1997年 ; Pannu&Read，1996年 ; 穆尔舒多夫等。, 1997 ). 这一成功的关键因素是开发交叉验证技术，以减少控制似然函数行为的超参数估计中的偏差（Lunin和Skovoroda，1995）; Pannu&Read，1996年). 21世纪初，里德及其同事将似然形式主义扩展到了分子置换设置，从而大大提高了从边缘启动模型求解结构的能力（麦考伊等。, 2005 ; 斯托罗尼等。, 2004 ; Read，2001年 ). 首次使用近似似然方法从反常或导数数据中检测重原子源于Terwilliger&Eisenberg（1983）)，他使用了原始的Patterson相关函数下部结构解决方案。Bricogne（1997）展示了这种方法 )等价于基于Ric的似然函数的二阶近似。最近的一个发展是在子结构的位置上使用了更精细的似然公式（Bunkóczi等。, 2015 )显示了重原子定位能力的显著提高。在密度修正中，似然公式的使用显著增加了其收敛半径（Terwilliger，2000 ; 考坦，2000年 ; 斯科巴克等。, 2010 ).

如上述示例所示，通过将基于相似性的方法应用于各种晶体学问题，已经取得了令人瞩目的进展。在所有描述的场景中，通过推导特定于问题的似然函数并将其应用于具有挑战性的结构确定问题，取得了关键进展。在大多数情况下，对实验误差的彻底处理只起到次要作用，从而使近似值在中等或低噪声环境中工作良好。处理晶体学似然函数中随机噪声的主要挑战是如何有效卷积类Rice分布函数，以模拟结构系数从一个带有误差的不完全模型出发，用适当的分布对实验噪声进行建模。在这份手稿中，我们开发了正交方法来克服这些困难。我们通过使用一系列变量的变化来实现这一点，这些变量易于使用标准方法进行直接的数值积分。导出的方法在模型中有直接的应用精细化和分子替换，而一般方法也可以推广到其他晶体学场景。在本文的其余部分中，我们将对基于似然的方法进行一般介绍，提供数值积分技术的相关背景，开发一种自适应求积方法，将其应用于Rice型似然函数并验证其结果。

1.1.最大似然形式主义

模型参数的估计θ给定一些数据集 $[｛\cal X｝]$ = {x个₁…，x个_j个…，x个_N个}通过可能性形式化以以下方式执行。给定概率密度函数（PDF）（f）(x个_j个|θ)单个观察结果x个_j个给定模型参数θ，在观测值独立的假设下，整个数据集的联合概率等于单个PDF的乘积，

$[f（｛\cal X｝|\boldtheta）=\textstyle\prod\limits_{j=1｝^{N} （f）（x_j|\boldtheta）。\等式（1）]$

数据的概率 $[{\cal X}]$ 给定模型参数θ已知为给定数据的模型参数的可能性：

$[L（黑体θ|{\cal X}）=f（黑体X}|\boldtheta）.\eqno（2）]$

一个自然的选择最佳估计通过找到θ使似然函数最大化。此选项称为最大似然估计值（MLE）。似然函数本身 $[L（黑体字|{\cal X}）]$ 可以被视为概率分布，允许获得MLE的置信极限估计（Rossi，2018). MLE的测定通常通过优化对数似然进行：

$[\ln L（\boldtheta|{\cal X}）=\textstyle\sum\limits_{j=1}^{N}\ln f（X_{j}|\boldtheta）.\eqno（3）]$

通常，似然函数所需的分布必须通过一个称为边缘化的过程来获得。在此集成过程中，集成了一个所谓的有害参数，

$[f（x|\boldtheta）=\textstyle\int\limits_{-\infty}^{\infty}f（x，y|\bolttheta）\，{\rm d}y，\eqno（4）]$

其中，在条件独立的假设下，

$[f（x，y|\boldtheta）=f（x|\bolttheta）f（y|x，\boldtheta）.\eqno（5）]$

根据所涉及的分布的数学形式，这种边缘化可能从一个微不足道的分析练习到一个具有数值挑战性的问题。在晶体学背景下的似然函数中，这种边缘化需要考虑实验噪声的影响，其有效计算是本次交流的重点。

1.2. 动机

晶体应用中最常见的似然函数指定了真实的给定计算值的结构系数振幅结构系数源自一个有误差的模型（Sim，1959 ; Srinivasan&Parthasarathy，1976年 ; 卢扎提，1952年 ; 伍尔夫森，1956年 ; Lunin和Urzhumtsev，1984年):

$[f_{rma}（f|f_{rmC}，\alpha，\beta）={2F}\over{varepsilon\beta}}\exp\left（-{f^{2}+\alpha^{2} F类_{\rm C}^{2}}\ over{\varepsilon\beta}}\ right）I_{0}\ left（{{2\alpha FF_{\rmC}}\ ever{varepsilon\beta}}\ rift），\eqno（6）]$

$[\eqaligno{f_{rmc}（f|f_{rmC}，\alpha，\beta）&=\left（{{2}\over{varepsilon\pi\beta}}\right）^{1/2}\exp\left^{2} F类_{\rm C}^{2}}\ over{2\varepsilon\beta}}\right）\cosh\left（{{2\alpha FF_{\rm-C}}\-over{2\ varepsiron\beta{}}\ right）。\cr&&（7）}]$

（f）_一和（f）_c（c）是无中心反射和中心反射的分布（所谓的Rice分布），∊是对称增强因子，F类是真实的结构系数振幅F类_C类是当前模型结构系数振幅，而α和β是似然分布参数（Lunin和Urzhumtsev，1984). 对于精细化在给定实验数据的原子模型中，基于模型的结构因子振幅的可能性是必需的，并且可以通过对未知的、无误差的结构因子幅度的边缘化来获得。继Pannu和Read之后（1996年)假设实验强度分布之间存在条件独立性我_o个和振幅F类，我们获得

$[\eqaligno{L（F{\rm C}|I{\rmo}）&=F（8）}]$

哪里（f）(F类|F类_C类，α，β)由表达式（6）给出或（7）和（f）(我_o个|σ_我²，F类)等于具有平均值的正态分布F类²和方差σ_我²该积分等效于Pannu&Read（1996）推导的MLI目标函数). 由于该积分没有快速收敛的级数近似或简单的闭合形式分析表达式，因此已经开发了各种方法，正如Read&McCoy（2016 )，包括一种矩量法类型的方法，以找到基于强度的似然函数的合理分析近似值。

在这项工作中，我们研究了使用数值积分方法来获得积分（8)同时还考虑了估计标准偏差中的不确定性。上述方法将Rice函数与高斯分布进行卷积，基本上假设平均值的标准偏差是准确已知的。鉴于标准偏差和平均值都是从相同的实验数据中得出的，这个假设显然是次优的，尤其是当数据冗余度较低时。为了考虑观察到的标准偏差中可能存在的错误，我们将使用t吨-分布而非正态分布，当真实标准偏差由实验数据的估计值近似时，这是分布的选择（Student，1908 ). 这项工作的目的是推导一种有效的方法来获得目标函数，在处理非常边缘的数据时，例如从时间分辨串行晶体学或自由电子激光数据获得的数据，可以提供额外的性能提升，其中实验误差通常大于使用基于标准旋转的方法获得的误差或具有非标准误差模型的误差（布鲁斯特等。, 2019 ). 此外，高质量的数据集很少受到衍射几何的限制，这表明如果使用适当的目标函数，可以随时获得更多的边缘数据，从而可能提高最终模型的质量。在本手稿的其余部分中，我们开发并比较了一些数值积分方案，这些方案旨在快速评估基于强度的似然函数及其导数，其中考虑了平均强度及其相关标准偏差中存在的实验误差。

2.方法

为了评估各种数值积分方案和近似方法，首先将方程重铸为归一化结构因子振幅E类和归一化强度Z轴_o个框架，使用σ_A类所涉及分布的公式，假设对1空间组，这样的话∊=1（阅读，1997 ). 无误差结构因子振幅的联合概率分布E类和实验强度Z轴_o个，给定计算的标准化结构系数 E类_C类，模型质量参数σ_A类，观测值的估计标准偏差σ_Z轴和有效的自由度 ν，读取

$[\eqaligno{f_{rma}（E，Z_{rmo}|E_{rmC}，\sigma_{rmA}，\sigma_{Z}^{2}，\ nu）&={2E}\在{1-\sigma{rmA}^{2}}\exp\左（-{{E^{2{+\sigma_{rm a}^{2} 电子_{\rm C}^{2}}\ over{1-\sigma_{\rmA}^{2}}\ right）\cr&\\quad{times}\I_{0}\ left（{{2\sigma{\rmA}EE_{\orm C}}\-over{1-\sigma_}\rm A}^{2]}\ right）f（Z_{\rmao}|E，\sigma{Z}^{2\nu）&（9）}]$

无偏反射和

$[\eqaligno{f_{\rm-c}（E，Z_{\rmao}|E_{\rm-c}，\sigma_{\RMA}，\sigma_{Z}^{2}，\nu）&=left[{2}\ over{\pi^{2} E类_｛\rm C｝^｛2｝｝｝\over｛2-2\西格玛_｛\rm A｝^｛2｝｝｝｝\right）\cr和quad｛\times｝\\cosh\left（｛\alpha EE_｛\rm C｝｝\over｛1-\西格玛_｛\rm A｝^｛2｝｝｝\right）f（Z_｛\rm o｝|E，\西格玛_｛Z｝^｛2｝，\nu）&（10）｝]$

用于中心反射。当观测到的平均强度分布Z轴_o个由t吨-分布（学生，1908)位置参数等于E类²，我们有

$[\eqaligno{f（Z_{\rmo}|E，\sigma_{Z}^{2}，\nu）&=\Gamma\left（{{\nu+1}\over{2}}\right）\left^{2}_{Z} }}\right）^{1/2}\left[\Gamma\left（{{\nu}\over{2}}\right）\right]^{-1}\cr&\\quad{times}\\left[1+{{left（Z_{o} -E类^{2} \右）^{2}}\在{\nu\sigma_{Z}^{2{}}\右]^{-（\nu+1）/2}上，&（11）}]$

哪里ν是有效的吗自由度与有效样本量有关N个_效率，

$[\nu=N_｛\rm eff｝-1。\等式（12）]$

有效样本量可以被视为观测强度的冗余，也可以在数据处理过程中使用韦尔奇-萨特思韦特方程进行估算（韦尔奇，1947 )考虑数据处理中实施的加权协议（布鲁斯特等。, 2019). 这个t吨-分布是指在给定一组观测值的样本均值和样本方差的情况下，选择的分布（Student，1908). 正态分布的使用基本上假设方差没有不确定性σ_Z轴²，但仅限于观察到的平均值Z轴_o个. Thet吨-分布类似于正态分布，但尾部较重，因此预计将产生对观测强度和模型强度之间较大偏差惩罚较小的似然函数。什么时候？ν趋于无穷大，则上述分布收敛于正态分布，

$[f（Z_{\rmo}|E^{2}，\sigma_{\RMo}^{2{）={1}\over{\left（{2\pi\sigma-{Z}^{2]}\right）^{1/2}}\exp\left[-{{（Z_\\rmo}-E^{2neneneep）。\等式（13）]$

上述联合概率分布需要被边缘化E类在里面 $[{\bb R}^{+}]$ 要获得利息分配：

$[f_{\cdot}（Z_{\rmo}|E_{\rm C}，\sigma_{\rma A}，\sigma_{Z}^{2}，\nu）=\textstyle\int\limits_{0}^{\infty}f_{\ cdot}]$

2.1. 差异通货膨胀

避免执行上述积分的常见方法是增大Rice函数的方差（1−σ_A类²)通过“观测到的结构系数振幅”的方差，得出（1−σ_A类²+σ_E类²)（格林，1979年 ). 这种方法避免了执行集成的需要，但在许多不同的方面都是次优的。因为我们不观察振幅，所以需要根据观测到的强度数据估计振幅及其方差。进行强度-振幅转换的一种常见方法是通过贝叶斯估计（French&Wilson，1978）)在原子均匀分布的假设下单位单元格。尽管在大多数情况下都使用这种所谓的威尔逊先验，但当对正半线上的结构因子振幅的可能值使用一个常量、不适当的先验时，可以得到稍微不同的结果（Sivia&David，1994 ). 这导致强度-振幅转换，不依赖于平均强度的准确估计，可能会因伪对称、衍射各向异性或孪晶的影响而变得复杂：

$[E_{\rmo}=\left\{{1}\over{2}}\left[Z_{\rmo}+（{Z_{rmo}^{2}+2\sigma^{2}_{Z} }）^{1/2}\right]\right\}^{{1/2}}，\eqno（15）]$

$[\displaystyle\sigma_{E}^{2}={{\sigma_{Z}^{2%}\在{4\左上方（{Z_{rmo}^{2]+2\sigma_2}}\右）^{1/2}}.\eqno（16）]$

更多详细信息见附录E类虽然此过程允许直接进行强度-振幅转换，即使强度为负，并且随后可用于膨胀Rice函数的方差，但它不能替代完全积分。鉴于方差膨胀方法的简单性及其在许多晶体学应用中的广泛应用，我们将使用该方法作为基准，使用基于威尔逊先验（表示为法语-威尔逊）和概述的统一非信息先验（标记为西维亚）的转换方案。

2.2. 数值积分方法

对于表达式（8）等不当积分，存在几种传统的数值积分近似标准方法包括具有截断积分范围的基于梯形的方法、使用拉普拉斯方法、基于蒙特卡罗的方法或基于正交多项式的方法（Davis&Rabinowitz，1984 ). 尽管直接使用梯形积分方案很有吸引力，但除非使用非常精细的采样网格，否则分布参数的某些组合的被积函数的形状将导致相当大的机会错过函数的大部分质量。在使用拉普拉斯近似时，被积函数由适当缩放和平移的高斯函数近似，被积矩阵可能会明显偏离高斯函数，也会导致性能不佳。图1总结了这些挑战，其中针对不同的参数选择可视化了许多典型的被积函数形状。下面将概述一些数值积分和近似方法，包括如何基本事实作为比较的基础。这里，由于拉普拉斯近似的简单性和梯形规则，我们将把自己局限于拉普拉斯逼近，因为它们在应用于实线上的解析函数时具有出色的收敛性，并且与经典高斯求积关系密切（Trefethen&Weideman，2014）). 不考虑使用（准）蒙特卡罗方法，因为这些方法通常被用作高维积分的“最后手段”（Cools，2002）).

图1
不同参数设置的偏心分布和中心分布的积分形状显示了计算边际似然时出现的各种函数形状。当实验误差相对于强度相对较大时，函数的高质量区域跨越积分域的适当部分E类≤ 6 (一). 当实验数据的误差相对较小时，被积函数的大部分质量集中在较小的区域(b条，c（c）). 如果是t吨-基于分布的噪声模型，与正常噪声模型相比，分布的尾部被抬高。这些形状的多样性使得标准求积或拉普拉斯近似的统一应用效率低下且次优。

2.3. 变量的变化与拉普拉斯近似

分析和数值积分通常会因被积函数变量的变化而大大简化（Davis&Rabinowitz，1984). 变分定理与函数的积分有关小时(u个)变量发生变化时u个=ψ(x个),

$[{\textstyle\int\limits_{a}^{b}}h（u）\，{\rm d}u={\texttyle\int\limits_{\psi^{-1}（a）}^{\psi ^{-1{（b）}h[\psi（x）]$

由于被积函数的形状因变量的改变而改变，所以使用所谓的拉普拉斯近似很有吸引力。在拉普拉斯近似中，被积函数由具有适当选择的平均值和长度标度的标度平方指数函数近似（Peng，2018 ). 拉普拉斯近似可由被积函数对数的截断泰勒展开式导出：

$[\eqalignno{{textstyle\int\limits_{a}^{b}}f（x）\，{rm d}x&={textstyle\nint\limiss_{a}^{b}}\exp[g（x）]\，{orm d}x\cr&={\textstyle\int\ limits_a}^{b2}}\exp\left[{textstyle\sum\limits{n=0}^{infty}}{g^{（n）}（x）{0}）}\在{n！}}上（x-x{0}）^{n}\right]\，{rmd}x\cr&\simeq{\textstyle\int\limits_{-\infty}^{\infty}}\exp[g（x{0{）]\exp\left[\textstyle｛｛1｝\over｛2｝｝g^｛\prime\prime｝（x_｛0｝）（x-x_｛0｝）^｛2｝\right]\，｛\rm d｝x，&（18）｝]$

哪里克(x个)=英寸[（f）(x个)]和x个₀是最大值的位置（f）(x个)，暗示着克′′(x个₀)=0。注意，在方程式（18）的最后一步中)假设是（f）(x个)不接近时变为0x个₀足够快，可以集成[一，b条]产生与积分相同的结果 $[\bb R}]$ 。虽然此近似值不适用于所有可能的克(x个)事实证明，它是贝叶斯分析中边缘化分布的一个成功工具（Kass&Steffey，1989 )和晶体学应用（Murshudov等。, 2011 ).

因此，上述表达式得出

$[{\textstyle\int\limits_{a}^{b}}f（x）\，{\rmd}x\simeqf（x{0}）\left[{2\pi}\ over{-g^{prime\prime}（x{0:}）}\right]^{1/2}。\等式（19）]$

这种近似的有效性取决于x个₀（应包含在原始积分域内）克′′(x个₀)以及高阶导数的速度克(x个)消失在周围x个₀上述变量转换策略有助于提高近似表达式的性能（8).

2.4. 正交法

尽管变量变化法与拉普拉斯近似法相结合有可能产生精确的积分近似，根据差分图或一阶或高阶优化方法的需要，获得对数似然导数的合理估计，使用拉普拉斯方法似乎不那么简单。困难在于需要获得被积函数最大值位置的导数，因为该值是计算导数的变量的函数。此外，引入t吨-基于噪声的模型在分布中引入了重尾，而高斯近似的性能较差。因此，使用求积方法很有意义。在用求积进行数值积分时，感兴趣的积分由选定函数值的加权和近似。因此，上述拉普拉斯近似可被视为一个单点求积，其中函数值的位置位于被积函数的最大值处，并且相关权重由被积函数局部高斯拟合导出。扩展求积方法提供了一种通过增加采样点数量来提高积分精度的简单方法，但也避免了使用拉普拉斯近似时遇到的计算被积函数最大值位置导数的问题。然而，假设正交法需要大量术语才能获得足够的精度（Read&McCoy，2016)这可能会使它们成为实际晶体学应用中不具吸引力的目标。为了避免或至少改善这些问题，我们设计了结合拉普拉斯近似和基本数值求积优点的求积（附录A类).

图2描述了我们集成方法的高级概述通过将幂变换与被积函数的双曲变换相结合，我们将积分域从[0，∞]变换到[0，1]。而第一次电力变换（附录C类)允许被积函数具有更像高斯的特征，第二个变分操作非线性地将低质量区域压缩到相对较小的线段上，同时将被积函数的高质量区域近似线性地变换到新积分域的中间（附录A类). 该双重变换函数随后可以使用等距梯形积分方案进行积分。与拉普拉斯近似一样，第二个变量转换操作需要了解幂变换被积函数的最大值，这可以使用标准优化方法获得（附录B类和C类). 在最后一步中，可以通过应用逆变换在原始变量中重写在双重变换变量的域上表示的结果求积。随后的进一步简化允许我们将整个积分重铸为加权Rice函数的总和，其中噪声观测和其他误差的影响是隐藏的在E类在 $[{\bb R}^{+}]$ 以及相关重量（附录D类),

$[\eqaligno{Q（E_{\rm C}|Z_{\rmo}，\sigma_{\rma A}，\sigma_{Z}^{2}，\nu）&=\ln L_{\cdot}（E_{\ rm C{|Z{\rmO}，\ sigma{\rmA}，\tigma_\Z}^{2}，_nu）\cr&=\ln\left[\textstyle\sum\limits_{j=1}^{N} w个_{j} \，f_{\cdot}（E_{j}|E_{\rm C}，\sigma_{\rmA}）\右]，&（20）}]$

哪里E类_j个是正交采样点和w个_j个是关联的权重。采样点和权重取决于E类_C类，Z轴_o个，σ_A类，σ_Z轴和ν。使用的正交采样可以是N个-点功率变换的双曲正交或基于（功率变换的）拉普拉斯近似的单点正交。更多详细信息见附录A类–D类 .

图2
所开发的数值积分程序描述为一系列步骤。一般的想法是使用一系列变量变换，从而在[0,1]上生成一个光滑函数，该函数可以通过梯形积分方案轻松积分。一旦建立了求积点，积分就可以写成加权Rice函数的和。详见正文。

2.5. 衍生品

实际使用基于相似hood的目标函数需要计算其导数，以便将其用于基于梯度的优化方法。来自表达式（20)，关于Y（Y）∈ {E类_C类，σ_A类，ν}可以如下获得：

$[\eqaligno{Q^{prime}_{Y}（E_{\rm C}|Z_{\rmo}，\sigma{\rmA}，\ sigma_{Z}^{2}，\nu）&={{\rmd}\ over{\rmd}Y}}Q（E_{\rm C}|Z_{\rmo}，\sigma_{\rma A}，\fu）\cr&=\exp\left[-Q（E_{\rm C}|Z_{\rmo}，\sigma_{\rma A}，\sigma_{Z}^{2}，\nu）\right]\cr&\\quad{times}\{\textstyle\sum\limits_{j=1}^{N}}{{\rmd}w_{j}\，f_{\cdot}（E_{j}|E_{\rm C}，\sigma_{\rmA}）}\在{{\rmd}Y}}上&(21)}]$

大米成分的衍生物（f）_·(E类_j个|E类_C类，σ_A类)关于E类_C类附录中列出了B类.

3.结果和讨论

评估所提出的积分方法的第一步是建立我们希望近似的积分的基本事实。为此，构造了一个等间距的非幂变换梯形求积，将E类=0至E类=6，使用50 使用所有分布参数组合的000个采样点，如表1所示，在强度上的高斯误差的假设下。将此积分结果与使用1500点双曲求积得到的结果进行比较，表明这两种积分方法给出了类似的结果。因此，我们将地面实况作为使用1500个或更多样本点通过双曲线求积获得的结果。对于偏心分布和中心分布，设置功率转换变量γ到2提供了良好的结果，如表2所示和3，其中报告了对数被积函数中相对误差的平均值和标准偏差（以百分比表示）。对许多不同的近似方案进行了测试，并使用对数被积函数中的平均相对误差对结果进行了比较。因为方差-通货膨胀近似值实际上并不执行边际化，而是执行更多特别的将低脂测量纳入修正，其相对于对数似然的相对误差并不是对其性能的公平衡量，也无法洞察其优缺点。相反，我们将比较对数似然目标函数相对于E类_C类对于所有近似，因为与方差膨胀方法相比，该度量与计算全积分时出现的不同规范化无关。此外，假设对数似然函数形式的梯度是用于完成或重建结构的三维差分或梯度图的傅里叶系数，将各种近似的梯度与从完全似然函数中获得的梯度进行比较，可以对不同近似的强度提供有价值的见解。梯度的使用当然是基于能够估计σ_A类，在这种情况下，可以使用固定分辨率外壳中的简单行搜索来执行。这些测试及其结果的详细信息见下文。

表1
比较积分方法的参数界限

参数	起点	终点	采样点
E类_C类	0.1	6	20
σ_A类	0	0.95	10
Z轴_o个	−5.0	50	20
Z轴/σ_Z轴	0.5	10	20

表2
积分结果：偏心分布

在整个参数范围内，相对对数似然的平均误差和标准偏差以百分比的形式报告。

方法	γ= 1	γ= 2	γ= 3
拉普拉斯近似	−0.142/0.874	0.294/0.971	0.485/1.135
正交(N个= 3)	0.191/0.778	0.152/0.831	0.281/1.058
正交(N个= 5)	0.130/0.377	0.126/0.481	0.196/0.627
正交(N个= 7)	0.085/0.218	0.074/0.309年	0.116/0.428

表3
集成结果：中心分布

在整个参数范围内，相对对数似然的平均误差和标准偏差以百分比的形式报告。正交结果γ=1不存在，因为双曲求积方案不保证函数在原点处为零。

方法	γ= 1	γ= 2	γ= 3
拉普拉斯近似	−1.766/8.170	0.357/1.729	0.738/1.841
正交(N个= 3)	—	0.300/1.617	0.725/1.850
正交(N个= 5)	—	0.391/0.990	0.438/1.183
正交(N个= 7)	—	0.269/0.750	0.311/0.943

3.1. 比较集成方法

图3显示了使用许多不同近似值的积分结果的比较对于根据附录中概述的程序生成的数据集F类。对于显示的结果σ_A类设置为0.70，并选择固定的错误率，以便Z轴/σ_Z轴〉 = 1.0. 冗余设置为4，导致ν= 3. 对于Z轴_o个，E类_C类对，一个使用1500点双曲求积计算的似然函数被视为一个t吨-分布和假设正态分布的误差模型。将这些值与拉普拉斯近似值（一点求积）以及两种误差模型的七点和49点求积进行比较。虽然拉普拉斯近似在正态误差模型中表现相对较好，但它无法正确处理t吨-分布，使用求积可以获得更好的结果。使用由七个或更多采样点组成的正交得到了令人满意的结果。原则上，可以开发通用的启发式算法，以根据误差模型的超参数定制求积的特定精度。正如预期，t吨-低的分布ν由于存在较重的尾部，与源自正态分布的值相比，值需要较大的正交才能获得可比较的误差。

图3
使用拉普拉斯近似和基于求积的方法对正态和正态分布的似然函数的相对平均误差t吨-基于噪声模型，用于偏心(一)和中心反射(b条). 虚线水平线设置为1%作为视觉参考。

3.2. 比较似然函数

为了更好地直观地了解目标函数的行为，我们直接为几个输入参数绘制它们。图4描述了似然函数L（左）(E类_C类|Z轴_o个)对于偏心反射，仅使用French–Wilson协议来估计振幅，同时不膨胀方差，并使用French-Wilson和Sivia方法的方差膨胀方法，以及使用高斯误差模型和t吨-分布变量。所示的所有函数均已在0≤范围内进行数值标准化E类_C类≤ 12. 当比较弱强度和负强度的曲线时，使用估算值的技术之间存在显著差异E类_o个基于非信息先验（法语–Wilson&Sivia）与完全积分得出的先验（图4一, 4b条和4c（c）). 如果观测值的相关标准偏差较低，则近似值之间的差异较小。正常误差模型和t吨-分布表现在似然函数近似的尾部行为中，而最大值的位置似乎相对不变（图4d日). 假设正态误差模型和t吨-类型错误模型在估计σ_A类基于相应的似然逼近。使用附录中概述的协议构建了一个有错误的合成数据集F类。使用固定的错误级别选择错误，以便预期的〈Z轴_o个/σ_Z轴〉为0.5时ν→ ∞ （见附录F类). 产生的结果Z轴_o个，σ_Z轴和E类_C类值用于确定σ_A类通过黄金分割驱动的似然最大化程序（基弗，1953 ). 由此得出的估计σ_A类以及不同冗余值的相关估计标准偏差(ν+1）如图5所示。对于大值ν的估计值σ_A类对于两个误差模型都是等效的，在较低的冗余值下，正态误差模型系统地低估了σ_A类。当使用French–Wilson协议时σ_A类估计被低估得更多（图6).

图4
对于不同的输入参数，在许多不同近似下的归一化似然函数的形状表明，使用负强度或高噪声值的点估计会导致与理想似然函数之间的显著偏差。The difference between at吨-基于噪声模型和正态噪声模型虽然较小，但显著影响了似然函数的尾部行为。图中给出了法国-威尔逊和西维亚方法的振幅和标准偏差估计值。

图5
基于lilihood的行为σ_A类-当数据具有t吨-采用基于似然的方法处理基于正态噪声的噪声模型：在估计σ_A类低冗余。

图6
使用不同近似方案计算的梯度的比较t吨-基于噪波ν= 3. (一)–(c（c）)和(d日)–(（f）)用降低的噪声水平描述梯度近似的行为。坡度使用最大似然估计σ_A类使用相应的近似值。正常模型和t吨-基于模型的性能明显优于French–Wilson和Sivia方法，而当t吨-采用基于模型的方法。

3.3. 比较对数似然梯度

通过直接比较一组选定参数组合的梯度，可以进一步了解似然函数近似的行为。数值试验表明，使用幂变换函数的1500点双曲求积（带γ对于偏心分布和中心分布都设置为2）与用50计算的有限差分梯度无法区分 000点梯形法。为了研究各种近似的质量精细化场景，我们使用附录中概述的随机抽样方法构建一个合成数据集F类.冗余4(ν=3）用于这些测试中。梯度使用49点积分计算，使用值为σ_A类根据似然函数的相应近似估计。这些比较的结果如图6所示并总结在表4中和5。坡度的质量由相关系数至真实值。结果表明，对于Z轴₀/σ_Z轴〉较大，所有梯度计算方法都收敛于那些使用基于强度的完全似然函数获得的方法，具有实验误差和Studentt吨噪声模型，但对于高噪声级和中等噪声级，方差膨胀法的表现明显不佳。而正常和t吨-基于相关系数，风格噪声模型似乎很小，在高噪声和低冗余设置下，在单个反射中可以看到显著的偏差。这些异常梯度可能会对结构完工梯度图的质量产生负面影响。

表4
通过计算使用1500点求积计算的梯度与正确σ的相关性，比较模拟数据的似然梯度_A类值（0.70）和使用四种不同近似方法获得的值，如正文所述，以及最大似然σ的估计_A类给出了似然函数的近似

报告的条目是以下各项的估计值σ_A类梯度相关。

方法	〈Z轴/σ_Z轴〉 = 0.25	〈Z轴/σ_Z轴〉 = 0.5	〈Z轴/σ_Z轴〉 = 1.5
西维亚	0.35/68.8%	0.52/82.1%	0.64/93.5%
法语–威尔逊	0.44/80.5%	0.53/87.7%	0.64/94.3%
正常	0.68/96.9%	0.66/97.4%	0.67/97.9%
学生的t吨	0.73/100%	0.70/100%	0.71/100%

表5
通过计算使用1500点求积计算的梯度与正确σ的相关性，比较模拟数据的似然梯度_A类值（0.90）和使用四种不同近似方法获得的值，如正文所述，以及最大似然σ的估计_A类给出了似然函数的近似

报告的条目是以下各项的估计值σ_A类梯度相关。

方法	〈Z轴/σ_Z轴〉 = 0.25	〈Z轴/σ_Z轴〉 = 0.5	〈Z轴/σ_Z轴？=1.5
西维亚	0.50/63.1%	0.64/73.6%	0.85/93.2%
法语–威尔逊	0.64/60.3%	0.65/74.9%	0.79/88.4%
正常	0.91/97.0%	0.88/96.9%	0.87/97.6%
学生的t吨	0.94/98.7%	0.91/99.9%	0.89/99.9%

4.结论

开发并比较了有效确定基于强度的似然函数及其梯度的数值方法。虽然拉普拉斯近似在正态噪声模型下对似然函数本身的估计表现得相当好，但我们的结果表明，使用数值求积可以显著提高似然及其相关梯度。鉴于对数似然函数的导数是基于梯度的关键成分精细化方法和用于计算结构完工的差分图，该方法可以提高现有方法的收敛性精细化和建模方法。虽然在实际情况下，最佳正交阶数或噪声模型应该是什么尚不清楚，但我们的结果表明，正常噪声的采样点可能低于15个，而正常噪声的则可能低于49个t吨-类型错误。在算法上，最昂贵的运算是寻找被积函数最大值的迭代过程。提出的基于牛顿的方法通常能在50个函数求值范围内很好地收敛，即使在没有预定好的行搜索起点的情况下也是如此。双曲求积的构造不需要任何迭代优化，也不需要相关梯度和函数值的后续计算。考虑到精细化或其他最大似然在结晶学应用中，使用所提出的方法计算目标函数可能只对工作流的总运行时间产生最小的影响，同时为基于强度的完全似然提供了一个快速收敛的近似值，该似然考虑了平均强度及其方差估计中的实验误差。虽然只有完全集成到晶体学软件包中才能确定使用概述的方法可以获得实际好处的情况，但这里的测试表明，有可能进行重大改进。此外，与现有近似方法相比，所提出的求积方法易于适应不同的误差模型选择，这是一个很大的优点，例如，可以在精细化和分阶段目标（Sharma等。, 2017 ).

附录A

双曲线求积

给定一个函数克(x个)，使用x个≥0且克(x个)≥0，我们寻求计算其在正半线上的积分：

$[G=\textstyle\int\limits_{0}^｛\fnty｝G（x）\，｛\rm d｝x.\eqno（22）]$

设置

$[h（x）=ln g（x）.\eqno（23）]$

定义的上确界克(x个)由x个₀这样的话小时′(x个₀) = 0. 对于我们感兴趣的函数类，克（0）等于0，例如由于正文中概述的功率变换，以及 $[\lim\limits_{x\to\infty}g（x）]$ 也是0。根据移位和重缩放的逻辑函数定义以下变量变化，

$[t={{exp（kx）-1}\在{\exp（kx）+\exp的（kx{0}）}}上，\eqno（24）]$

哪里k个是一个正数。请注意t吨(x个=0）=0和 $[\lim\limits_{x\to\infty}t（x）]$ = 1. 逆函数为

$[x（t）=x_{0}-{{1}\over{k}}\ln\left[{{\exp（x_{0}k)（1-t）}\在{1+t\exp（kx{0}）}}\右]\eqno（25）]上$

并具有关于t吨等于

$[x^{\prime}（t）={{\exp（-kt）[\exp。\方程式（26）]$

价值观x个₀确定双曲线压缩格式的近似“拐点”和常数k个控制中点周围的坡度。安N个-现在可以通过均匀采样来构造点求积t吨介于0和1之间，

$[t_{j}={{j}\在{N+1}}上，\eqno（27）]$

对于1≤j个≤N个考虑到这两者克（0）和 $[\lim\limits_{x\to\infty}g（x）]$ 为零，积分G公司现在可以通过梯形积分规则计算，

$[G={{1}\在{N+1}}\textstyle\sum\limits_{j=1}上^{N} 克[x（t{j}）]x^{\素数}（t{j}）。\等式（28）]$

如果k个被选为

$[k=\left[｛｛-2h^｛\prime\prime｝（x_｛0｝）｝\over｛\pi｝｝\right]^｛1/2｝\eqno（29）]$

那么上面的求积N个=1表示拉普拉斯近似x个₀|小时′′(x个₀)|^1/2很大，因为|x个₀−x个（1/2）|变为零。如果在幂变换变量的分布上构造双曲求积，则这些导出的权重可以乘以该变换的雅可比矩阵，从而可以在原始变量集中进行最终的数值计算。

附录B

分配和衍生品

B1.水稻功能，无着丝粒

无偏Rice分布的对数及其导数E类和E类_C类如下所示。

$[\eqaligno{h_{\rma}，{\rmRice}}（E，E_{\rmC}，\sigma{\rmA}）&=\ln f_{\rma}（E|E_{\rmaC}、\sigma{\rmA}）\cr&=\lin 2+\ln E-\ln左}\ over{1-\sigma_{\rm a}^{2}}\cr&\\quad+\\ln eI_{0}\ left（{{2\sigma{\rmA}EE_{\rma C}}\ over the{1-\sigma_}\rm a}^{2}}\ right），&（30）}]$

$[\eqaligno{h^{prime E}{{\rma}，{\rmRice}}（E，E_{\rmC}，\sigma_{\rmA}）&={\rmd}\ over{\rmd}E}}\ln f_{\rma}（E|E_{\rmaC}、\sigma{\rmA}）\cr&={1}\over{E}}-{2（E-\sigma{\rm a}E_\rm{1-\sigma_{\rm a}^{2}}}\cr&\\quad+\{2E_{rm C}\sigma a}}上\左（\displaystyle{2EE_{\rm C}\sigma_{\rma A}}\ over{1-\sigma A}^{2}}}\ right）}\ over{I_{0}\ left（\disposystyle}{2EE_{\rm-C}\simma{\rm-A}}\ ever{1-\ sigma{\rmA}^{2}}}\right）{}-1\ right]，&（31）}]$

$[eqaligno{h^{prime\prime E}{{rma}，{rmRice}}（E，E_{rmC}，\sigma_{rmA}）&={{rmd}^{2}}\over{rmd{E^{2{}}\ln f_{rma{（E|E_{rmc}，\sigma{rmA{）\cr&=-{1}\over{E^2}}}-{{2}\在{1-\sigma_{\rmA}^{2}}+\左（{2E_{\rmC}\sigma _{\rmaA}}在{1-\ sigma _{\rmA}^{2}}\右）}在{E（1-\sigma_{\rm A}^{2}）}\left[{I_{1}\left（\displaystyle{2EE_{\rma C}\sigma A}}\over{1-\simma_{\orm A}^{2{}}\right）}\ over{I_0}\left}^{2}}\right）}}\ right]\cr&\\quad-\\left（\显示样式{{2E_{rm C}\sigma_{rm A}}\ over{1-\sigma{\rm A{^{2{}}\右）^{2neneneep\left[{I_{1}\ left（\显示样式{2EE_{\rm C}\sigma_{\orm A}}\ over{1-\sigma A}^{2}}\ right）}\ over{I_0}\ left[\显示样式}{2EE_{\rm C}\ sigma_{\rmA}}\ ever{1-\sigma_{\rma A}^{2}}\ right]}，&（32）}]$

$[\eqaligno{h^{\prime E_{\rm C}}{{\rma}，{\rmRice}}（E，E_{\rma C}，\sigma_{\rmA}）&={\rmd}在{\rmd}E_{\orm C}{}}\ln f_{\rmaa}（E|E_{\RMC}，\sigma_{\RMA}\rm a}E_{\rm C}\左（\displaystyle{2EE_{\rm C}\sigma_{\rma A}}\ over{1-\sigma A}^{2}}}\ right）}\ over{I_{0}\ left（\disposystyle}{2EE_{\rmC}\simma{\rmA}}\ ever{1-\sigma_{\orm A}^{2}{}}}\right）{}-1\ right]&(33)}]$

B2.大米功能，中心

中心Rice分布的对数及其导数E类和E类_C类如下所示。

$[eqaligno{h_{\rmc}，{\rmRice}}（E，E_{\rmC}，\sigma{\rmA}）&=ln f_{\rmac}（E|E_{\rmaC}、\sigma{\rmA}）\cr&={1}\over{2}}[\ln2-\ln\pi-\ln（1-\sigma-{\rm A}^2}）}^{2} E类_{\rm C}^{2}}\ over{2（1-\sigma_{\rmA}^{2}）}\cr&\\quad+\\ln\cosh\left$

$[\eqaligno{h^{primeE}{{rmc}，{rmRice}}（E，E_{rmC}，\sigma_{rmA}）&={{rmd}\over{rmd{E}}\lnf_{rmc{^{2}-1}}+{{\sigma_{\rm A}E_{\rmC}}\在{1-\sigma A}^{2}}\tanh\left上（{{\sigma_}\rmA}EE_{\rmaC}}\over{1-\σA}^}}}\right），&（35）}]$

$[\eqaligno{h^{prime\prime E}{{rm c}，{rm Rice}}（E，E_{rm c}，\sigma_{rm A}）&={{rmd}E^{2}}\ln f_{rmc}（E|E_{rmC}，\sigma_{rmA}）^{2}-1}}+左（{{\sigma_{\rm A}E_{\rmC}}\在{1-\sigma A}^{2}}\右上）^{2{\cr&\\quad-\\left C}}\在{1-\sigma_{\rm A}^{2}}\右）\right]^{2{，&（36）}]上$

$[\eqaligno{h^{\prime E_{\rm C}}{{\rmc}，{\rmRice}}^{2} E类_{\rm C}}\在{\sigma_{\rmA}上^{2}-1}}+{{\sigma_{\rm A}E}\在{1-\sigma A}^{2}}\tanh\left上（{{\σA}EE_{\rma C}}\在}1-\simma_{\orm A}^}}\right上）。&(37)}]$

B3.学生的t吨-分配

的对数t吨-正文及其派生词中规定的关于E类如下所示。

$[\eqalinno｛h｛\rm o｝（Z_｛\rm o｝，E，\ sigma_｛Z｝^｛2｝，\nu）&=\ln f（Z_｛\rm o｝| E，\ sigma_｛\rm o｝^｛2｝，\nu^{2}_{\rm o}\right）-\ln\Gamma\left（{\nu}\over{2}\rift）\cr&\quad-\{{\nu+1}\over{2}}\ln\left[1+{{\rmo}-E^{2}\ right）^{2{}\over-{nu\sigma_{Z}^{2neneneep}\right]，&（38）}]$

$[\eqaligno{h^{prime E}{\rmo}（Z_{\rmo}，E，\sigma_{Z}^{2}，\nu）&={\rmd}\over{\rmd}E}}\lnf（Z_\rmo{|E，\sigma_\rm o}^2}}^{2}\nu+（Z_{\rmo}-E^{2{）^{2neneneep}，&（39）}]$

$[eqaligno{h^{prime\prime\prime E}_{rmo}{{{[\sigma_{Z}^{2}\nu（Z_{\rmo}-3E^{2{）+（Z_}\rmo{-E^{2]）^{2｝（E^{2]+Z_{rmo}）]}\在{[\sigma_}Z}^2}\nu+（Z{\rmo}-E^{2neneneep）^{2]]^{2◄上。\cr&&（40）}]$

什么时候？ν很大t吨-分布可以近似为正态分布：

$[h_{\rmo}（Z_{\RMo}，E，\sigma_{Z}^{2}）=-{{1}\over{2}}\ln2\pi+ln\sigma_{Z}（Z）-{{左（Z_{\rmo}-E^{2}\右）^{2{}\在{2\sigma_{Z}^{2neneneep}}上，\eqno（41）]$

$[h^{\prime E}_{\rmo}（Z_{\rmao}，E，\sigma_{Z}^{2}）={2E\左（Z__{\RMo}-E^{2{右）}\上{\sigma{Z}，\eqno（42）]$

$[h^{\prime\prime E}_{\rmo}（Z_{\rmao}，E，\sigma_{Z}^{2}）={6E^{2{}\ over{\sigma{Z}^2}}.\eqno（43）]$

附录C

发现x个₀

如正文所述，通过双曲求积进行数值积分在很大程度上可以通过变量的变化来辅助

$[E=x^{{\gamma}}\eqno（44）]$

雅可比矩阵（见表达式17)

$[{{\rmd}E（x）}\在{{\rmad}x}}=\gammax^{\gamma-1}.\eqno（45）]$

被积函数最大值的位置，x个₀，使用牛顿寻根算法的直接应用，使用下面概述的一阶和二阶导数。设置

$[h_{t，{\rma}}（E）=h_{\rma}（E，\cdots）+h_{\rmo}（E.，\cdot）\eqno（46）]$

和

$[h_{t，{\rm c}}（E）=h_{\rmc}（E，\cdots）+h_{\rmo}（E.，\cdot）.\eqno（47）]$

使用速记小时_t吨(E类，……）为了表示这两个函数中的任何一个，我们首先应用幂变换（表达式44）如上所述。因为我们需要找到这个函数的最大值，我们需要关于x个也。所得函数及其关于x个转换变量后的操作由

$[\eqaligno{\hat{h}（小时）_{t} （x，\gamma）&=\ln\gamma+（\gamma-1）\lnx+h{t}{t}（x，\gamma）&={{（1-\gamma）}\在{x^{2}}}+\gamma^{2} x个^{{2\gamma-2}}h^{prime\prime}_{t}（E=x^{gamma}，\cdots）\cr&\\quad+\\gamma（\gamma-1）x^{\gamma-2{h^{prime}_t}(48)}]$

的分析形式小时_t吨(⋯),小时′_t吨（…）和小时′′_t吨附录中给出了（…）B类牛顿搜索的衰减起始值可以通过对（例如）15个等距值的集合执行单个基于牛顿的更新来找到x个在0和之间采样x个 = 6^1/γ.得到的更新采样点的被积加权平均值通常在十次迭代内细化到上确界。

附录D

可能性综合

使用上述方法，完全似然函数可以表示为加权Rice函数的总和（表达式30和34)，其中E类基于从功率变换变量导出的正交进行采样E类=x个^γ使用上述双曲线抽样方案。考虑到功率变换和双曲求积的组合，求积的采样节点等于

$[E_{j}=x_{j}^{gamma}，\eqno（49）]$

$[x{j}=x_{0}-{{1}\over{k}}\ln\left[{{\exp（x_{0}k)（1-t{j}）}\over{1+t{j{\exp（kx{0}）{\right]，\eqno（50）]$

哪里t吨_j个，x个₀和k个按照附录中的规定进行定义和计算A类和C类且1≤j个≤N个。现在可以设置正交权重以吸收双曲线采样、功率变换和作用于观测强度及其相关标准偏差的误差模型：

$[\eqaligno{w{j}&=\gammax{j}^{\gamma-1}\次{{exp（-kt{j}）[\exp（kx{0}）+\exp将{{1}\乘以{N+1}}。&(51)}]$

这样就产生了一个加权Rice函数的总和，该函数近似于完全似然函数，

$[L_{\cdot}（E_{\rm C}|Z_{\rmo}，\sigma{\rmA}，\sigma{Z}^{2}）=\textstyle\sum\limits_{j=1}^{N} w个_{j} \，f_{\cdot}（E_{j}|E_{\rm C}，\sigma_{\rmA}），\eqno（52）]$

哪里（f）_·（？）是无着丝粒或有着丝粒Rice函数。

当用幂变换拉普拉斯近似代替求积方法来近似似然函数时，我们得到了一个加权Rice函数

$[L_{\cdot}（E_{\rm C}|Z_{\rmo}，\sigma{\rmA}，\sigma{Z}^{2}）=w_{0}，f_{\cdot}$

重量由

$[w_{0}=\gamma x_{0{^{\gamma-1}\次f（Z_{\rmo}|E_{j}，\sigma_{Z}^{2}，\ nu）\次\左[{{2\pi}\在{hat{h}^{prime\prime}（x_{0}）}}\右]^{1/2}上，\eqno（54）]$

哪里E类₀=x个₀^γ和 $[\hat{h}^{prime\prime}{t}（x{0}）]$ 在表达式（48）中定义.

附录E

结构因子振幅估计

为了使用膨胀方差修正作为全数值积分的近似值，我们需要能够估计反射振幅及其与观测强度的标准偏差及其标准偏差。虽然此过程通常使用标准的法式-威尔逊估算程序执行，但也可以采用另一种方法，遵循Sivia&David（1994）开发的方法). 在……之前穿制服E类，因此

$[f（E）=\cases{1&$E\geq 0$\cr 0&$E\，\lt\，\，0$}，\eqno（55）]$

导致条件分布

$[f（E|Z_{\rmo}，\sigma_{Z}^{2}）\propto E\exp\left[-{{（E^{2} -Z轴_{\rmo}）^{2}}\在{2\sigma上^{2}_{Z} }}\右]。\等式（56）]$

此分布的正态近似值可以通过矩方法获得，或者，如本文所述，通过最大值后验的使用等于上述分布模式的平均值和基于模式位置处对数似然二阶导数估计的标准偏差进行近似：

$[E_{\rmo}=\sup\limits_{E}\ln\left[f（E|Z_{\rm o}，\sigma_{Z}^{2}）\right]，\eqno（57）]$

$[\sigma_{E}^{2}=-\left\{{{{\rmd}^{2}}\ over{\rm d}E^{2{}}\ln\left[f（E|Z_{rmo}，\sigma{Z}^{2}）\right]\right\}^{{-1}}{E=E_{\rmo}}}}。\等号（58）]$

很容易得到解析表达式，从而得出

$[E_{\rmo}=\left\{{1}\over{2}}[Z_{\rmo}+（{Z_{rmo}^{2}+2\sigma^{2}_{Z} }）^{1/2}]\right\}^{{1/2}}，\eqno（59）]$

$[\sigma_{E}^{2}={{\sigma _{Z}^{2}}\在{4\左上方（{Z_{rmo}^{2]+2\sigma-{Z}^{2{}}\右）^{1/2}}.\eqno（60）]$

此近似值的质量主要取决于Z轴_o个和σ_Z轴²。请注意，上的标准错误传播E类_o个=Z轴_o个^1/2产量

$[\sigma_{E}^{2}={{\sigma _{Z}^{2]}\在{4Z_{rmo}}}上，\eqno（61）]$

可以看出它收敛于表达式（59）对于以下情况Z轴_o个明显大于σ_Z轴.

附录F

模拟合成数据

数值试验和基准的数据是通过使用以下程序从基础分布中取样获得的。对于无中心反射，通过从具有零平均值和指定方差的正态分布中连续提取来生成“真”复杂结构因子和“扰动”复杂结构系数：

$[\eqaligno{{\bfE}{\rm True}&=（X{\Re}，X{\Im}）\semif（X_{\cdot}）=N（0,1/2），\cr\boldDelta{\rmC}&=&=\sigma_{\rm A}{\bf E}_{\orm True}+\boldDelta_{\rma C}，&（62）}]$

哪里N个(μ，σ²)表示正态分布。对于中心反射，使用以下步骤：

$[\eqaligno{{bfE}{\rm True}&=（X{\Re}，0）\simf（X{\ Re}）=N（0,1），\cr\boldDelta{\rmC}&=}_{\rm True}+\boldDelta_{\rma C}&(63)}]$

以以下方式添加噪音。给定目标方差σ²_目标，我们生成ν+1个正态随机变量用于计算样本均值和样本方差：

$[\eqaligno{Z_{rm True}&=|{bf E}{rm真}|^{2}，\cr Z_{j}&=X\simf（X）=N[Z_{rma真}，（nu+1）\sigma^{2}_{\rm目标}]，\cr Z_{\rmo}&=[1/（\nu+1）]\textstyle\sum Z_{j}，\cr\sigma_{Z}^{2}&=\{1/[\nu（\nu+1）]^{2{}\}\textstyle\sum（Z_{j} -Z轴_{\rm o}）^{2}。&(64)}]$

所述测试中采用了两种误差模型，即固定误差比或a固定误差水平.在固定误差比法中，σ_目标每个模拟强度都不同，并选择为Z轴_真的/τ，其中τ等于所需 $[{\bb E}（Z_{\rm真}/\sigma_{\rma目标}）]$ 级别。对于固定误差水平方法，σ_目标固定为1/τ适用于所有强度。假设无中心反射与中心反射的比率为9:1，对完整的数据集进行了模拟。

致谢

上述算法是在一组Python3例程中实现的，可根据要求提供。这项工作的一些部分是为了部分满足伯克利实验室本科生研究（BLUR）项目的要求而准备的，该项目由伯克利实验室的劳动力发展和教育部管理。本文内容仅由作者负责，不一定代表NIH的官方观点。

资金筹措信息

这项研究部分得到了高级科学计算研究和基础能源科学项目的支持，这些项目由美国能源部科学办公室（DOE）根据合同DE-AC02-05CH11231提供支持。进一步的支持来自美国国立卫生研究院（NIH）国家普通医学研究所，获得5R21GM129649-02奖。

工具书类

Beu，K.E.，Musil，F.J.&Whitney，D.R.（1962年）。阿克塔·克里斯特。 15, 1292–1301. 交叉参考 IUCr日志谷歌学者
 Brewster，A.S.、Bhowmick，A.、Bolotovsky，R.、Mendez，D.、Zwart，P.H.和Sauter，N.K.（2019年）。阿克塔·克里斯特。D类75, 959–968. 交叉参考 IUCr日志谷歌学者
 Bricogne，G.（1997年）。CCP4研究周末会议记录。阶段化的最新进展由K.S.Wilson、G.Davies、A.W.Ashton和S.Bailey编辑，第159–178页。沃灵顿：达斯伯里实验室。谷歌学者
 Bricogne，G.和Gilmore，C.J.（1990年）。阿克塔·克里斯特。A类46, 284–297. 交叉参考中国科学院科学网 IUCr日志谷歌学者
 Bunkóczi，G.，McCoy，A.J.，Echols，N.，Grosse-Kunstleve，R.W.，Adams，P.D.，Holton，J.M.，Read，R.J.&Terwilliger，T.C.（2015）。自然方法，12, 127–130. 科学网公共医学谷歌学者
 Cools，R.（2002）。J.计算。申请。数学。 149, 1–12. 交叉参考谷歌学者
 Cowtan，K.（2000年）。阿克塔·克里斯特。D类56, 1612–1621. 科学网交叉参考中国科学院 IUCr日志谷歌学者
 Davis，P.J.和Rabinowitz，P.（1984）。数值积分方法第二版，纽约：学术出版社。谷歌学者
 Fisher，R.A.（1915年）。生物特征，10, 507–521. 谷歌学者
 French，S.&Wilson，K.（1978年）。阿克塔·克里斯特。A类34, 517–525. 交叉参考中国科学院 IUCr日志科学网谷歌学者
 Gauss，C.F.（1809）。Solem Ambientium松果线虫Corporum Coelestium理论汉堡：珀斯和贝瑟。谷歌学者
 Gauss，C.F.（1816）。Z.天文。韦旺特·维斯。 1, 1816. 谷歌学者
 Gauss，C.F.（1823）。微小恶臭的理论组合观察哥廷根：亨利科斯·迪特里希。谷歌学者
 Green，E.A.（1979年）。阿克塔·克里斯特。A类35, 351–359. 交叉参考中国科学院 IUCr日志科学网谷歌学者
 Hagen，G.（1867年）。Wahrscheinlichkeits-Rechnung公司。柏林：安永会计师事务所。谷歌学者
 Kass，R.E.和Steffey，D.（1989）。美国统计协会。 84, 717–726. 交叉参考谷歌学者
 Kiefer，J.（1953年）。程序。美国数学。Soc公司。 4, 502–506. 交叉参考谷歌学者
 La Fortelle，E.de&Bricogne，G.（1997）。方法酶制剂。 276, 472–494. 公共医学科学网谷歌学者
 Lunin，V.Y.和Skovoroda，T.P.（1995）。阿克塔·克里斯特。A类51, 880–887. 交叉参考中国科学院科学网 IUCr日志谷歌学者
 Lunin，V.Y.和Urzhumtsev，A.G.（1984）。阿克塔·克里斯特。A类40, 269–277. 交叉参考中国科学院科学网 IUCr日志谷歌学者
 Luzzati，V.（1952年）。阿克塔·克里斯特。 5, 802–810. 交叉参考 IUCr日志科学网谷歌学者
 McCoy，A.J.、Grosse-Kunstleve，R.W.、Storoni，L.C.和Read，R.J.（2005）。阿克塔·克里斯特。D类61，458–464页科学网交叉参考中国科学院 IUCr日志谷歌学者
 Murshudov，G.N.、Skubák，P.、Lebedev，A.A.、Pannu，N.S.、Steiner，R.A.、Nicholls，R.A、Winn，M.D.、Long，F.&Vagin，A.（2011）。阿克塔·克里斯特。D类67, 355–367. 科学网交叉参考中国科学院 IUCr日志谷歌学者
 Murshudov，G.N.、Vagin，A.A.和Dodson，E.J.（1997）。阿克塔·克里斯特。D类53, 240–255. 交叉参考中国科学院科学网 IUCr日志谷歌学者
 Neyman，J.和Scott，E.L.（1948年）。计量经济学，16, 1–32. 交叉参考谷歌学者
 Pannu，N.S.和Read，R.J.（1996年）。阿克塔·克里斯特。A类52, 659–668. 交叉参考中国科学院科学网 IUCr日志谷歌学者
 皮尔逊，E.S.（1970年）。统计与概率史研究由E.S.Pearson和M.G.Kendall编辑，第411-413页。伦敦：查尔斯·格里芬。谷歌学者
 彭瑞德（2018）。高级统计计算.https://leanpub.com/advstatcomp. 谷歌学者
 Read，R.J.（1986年）。阿克塔·克里斯特。A类42, 140–149. 交叉参考中国科学院科学网 IUCr日志谷歌学者
 Read，R.J.（1997）。方法酶制剂。 277, 110–128. 交叉参考公共医学中国科学院科学网谷歌学者
 Read，R.J.（2001）。阿克塔·克里斯特。D类57, 1373–1382. 科学网交叉参考中国科学院 IUCr日志谷歌学者
 Read，R.J.和McCoy，A.J.（2016）。阿克塔·克里斯特。D类72, 375–387. 科学网交叉参考 IUCr日志谷歌学者
 Rossi，R.J.（2018）。数理统计：基于可能性的推理导论。霍博肯：约翰·威利父子公司。谷歌学者
 Sharma，A.、Johansson，L.、Dunevall，E.、Wahlgren，W.Y.、Neutze，R.和Katona，G.（2017）。阿克塔·克里斯特。A类73, 93–101. 科学网交叉参考 IUCr日志谷歌学者
 Sim，G.A.（1959年）。阿克塔·克里斯特。 12，813–815交叉参考 IUCr日志科学网谷歌学者
 Sivia，D.S.和David，W.I.F.（1994）。阿克塔·克里斯特。A类50, 703–714. 交叉参考中国科学院科学网 IUCr日志谷歌学者
 Skubák，P.、Waterreus，W.-J和Pannu，N.S.（2010年）。阿克塔·克里斯特。D类66, 783–788. 科学网交叉参考 IUCr日志谷歌学者
 Srinivasan，R.和Parthasarathy，S.（1976年）。X射线晶体学中的一些统计应用第1版，牛津：佩加蒙出版社。谷歌学者
 Storoni，L.C.、McCoy，A.J.和Read，R.J.（2004）。阿克塔·克里斯特。D类60, 432–438. 科学网交叉参考中国科学院 IUCr日志谷歌学者
 学生（1908）。生物特征，6, 1–25. 谷歌学者
 Terwilliger，T.C.（2000）。阿克塔·克里斯特。D类56，965–972页科学网交叉参考中国科学院 IUCr日志谷歌学者
 Terwilliger，T.C.和Eisenberg，D.（1983年）。阿克塔·克里斯特。A类39, 813–817. 交叉参考中国科学院科学网 IUCr日志谷歌学者
 Trefethen，L.N.和Weideman，J.A.C.（2014）。SIAM版本。 56, 385–458. 交叉参考谷歌学者
 Welch，B.L.（1947）。生物特征，34, 28–35. 中国科学院公共医学科学网谷歌学者
 Wilks，S.S.（1938年）。安。数学。斯达。 9, 60–62. 交叉参考谷歌学者
 Wilson，A.J.C.（1980）。阿克塔·克里斯特。A类36，937–944页交叉参考中国科学院 IUCr日志谷歌学者
 Woolfson，M.M.（1956年）。阿克塔·克里斯特。 9, 804–810. 交叉参考中国科学院 IUCr日志科学网谷歌学者

这是一篇根据知识共享署名（CC-BY）许可它允许在任何介质中不受限制地使用、分发和复制，前提是引用了原始作者和来源。