沿着一篇 MBAR 论文检查自己的计算

先看图上的点到底来自哪里

冻结阅读版本为 Shirts 与 Chodera 的 $\text{arXiv:0801.1426v3}$，对应已发表论文 J. Chem. Phys. 129, 124105。第五节的图 1 比较一条恒力轨迹与多条轨迹合并得到的 DNA 发夹自由能。测得的伸长包含 DNA 发夹、两段双链 DNA 连接链与连接部位。横轴是这套连接体系的伸长，纵轴是相对于任意零点的平均力势。误差条描述所用统计处理下的估计波动。

原实验有十六条恒力记录，温度 $296.15\,\mathrm K$，每条五秒、每 $0.1\,\mathrm{ms}$ 采样一次。力从 $12.35$ 到 $14.41\,\mathrm{pN}$，重建目标为 $14.19\,\mathrm{pN}$。分析使用五十个不等宽伸长区间，因此式 (14) 要除箱宽。作者还说明力测量约有百分之十的相对误差。

读图时先问：单条轨迹在哪段伸长上样本少？其他力条件能否补上那一段？多窗口误差条较小，说明合并数据对这份平衡估计有帮助；它不自动消除力标定误差，也不提供无偏反应速率。

来源：Michael R. Shirts and John D. Chodera

确定这一次究竟复算什么

本课程复算论文式 (11) 的自由能自洽方程、目标权重与式 (14) 的密度换算。输入是自建二维势的带种子模拟轨迹。没有下载或重新分析原图的 DNA 发夹记录，也没有把二维山谷曲线标成原论文实验数据。

原文在 arXiv 保留作者授权的非独占传播许可；这不等于原始实验数据获得开放再分发许可。页面没有复制原图或原始轨迹。课程示意图与教学数据由本项目生成，数据文件同时保存参数、种子、软件版本和生成脚本。

这样的有限复算仍然有用：模型配分函数可以另行积分，方程可以交给另一套成熟求解器。若结果不符，就能把问题缩小到采样、求解或投影的某一步。原实验的复现则需要另外取得完整记录、校准说明与使用许可。

来源：Michael R. Shirts and John D. Chodera

先用独立求解器查方程

本次固定计算使用七个窗口，每窗两条独立副本，每副本保存一千帧，共一万四千个构型。浏览器默认每副本五百帧，共七千帧；这里用的是另存的基准记录。完全相同的输入分别交给浏览器模型的 MBAR 求解程序与 PyMBAR 4.0.3。比较前，先减去各自第一个窗口的自由能，统一零点。PyMBAR 官方文档要求独立样本；这里只比较同一份记录的方程解，相关帧的误差另用分块和独立种子评估。

七个约化窗口自由能的最大差约为 $1.26\times10^{-7}$。这个结果检验两套算法对同一份数据是否解出相同方程；它没有增加这份轨迹的覆盖，也没有把有限样本估计变成真值。

独立的二维 Simpson 积分通过网格加密检查，每个投影箱的概率变化约 $2.09\times10^{-12}$。使用四百八十格中点求积时，与该参考的每箱最大差约 $4.78\times10^{-6}$。这些误差对应本次网格和参数；换一个设置，要重新做网格加密。

来源：PyMBAR 4.0.3

再让独立种子挑战采样

默认模型左右对称，右半区的积分概率精确为 $0.5$。八个独立种子的短采样计算给出平均 $0.481$，种子间标准差约 $0.129$。把积分步长减半的另一批计算，平均约 $0.436$、标准差约 $0.200$。这些短轨迹的波动明显大于求解器之间的差。

步长减半后平均值没有自动更准。随机路径的采样差异仍很大，八个种子也不足以把很小的步长偏差与这种波动可靠分开。应增加独立交换和重复规模，再比较均值、误差与步长趋势。

在每副本一千帧的基准中，故意让隐藏坐标受困时，窗口投影仍可连通，右半区概率的绝对误差却达到约 $0.457$。没有重叠的预设则拒绝报告全局自由能。前者提醒我们检查完整构型的平衡；后者保留数据真正缺失的相对高度。

从低维检查走向研究问题

已发表的方法给出估计方程及其适用条件；本次可复算数值比较检查二维模型上的程序是否算出相同权重。自动选择集体变量、高维分子充分采样，以及从加速轨迹恢复可靠动力学，仍是需要另外验证的问题。

研究新坐标时，可先把它放回二维基准：同一投影位置的 committor 是否集中？从不同隐藏谷出发能否给出相同自由能？训练种子之外是否仍有改善？低维模型容易看见失败，也便于解释失败。

提出新方法的下一步，应把模型、观测量和失败标准提前固定。保留全部种子，包括失败的种子；比较同等计算量下的误差分布，再决定是否值得进入高维计算。最后应报告整组重复计算的误差分布。