使用 Python 拟合模拟和实验数据点答案

【问题标题】：Fit simulated and experimental data points with Python使用 Python 拟合模拟和实验数据点
【发布时间】：2011-11-18 16:37:45
【问题描述】：

我编写了一些代码来执行蒙特卡罗模拟并生成信号强度与时间的曲线。这种曲线的形状取决于各种参数，我的合作者希望通过我正在模拟的实验的“现实版本”来确定其中两个参数。

我们已准备好将她的实验数据与我的模拟曲线进行比较，但现在我被卡住了，因为我还不能进行任何拟合（到目前为止，我已经用模拟噪声数据替换了实验数据以进行测试）。我尝试使用scipy.optimize.leastsq，它以代码2退出（根据文档，这意味着拟合成功），但它主要只是返回值（不完全相同，但接近）我作为初始猜测输入，无论它们与真实值的距离有多近或多远。如果它确实报告了不同的值，那么得到的曲线仍然与真实曲线有很大的不同。

另一个观察结果是infodict['nfev'] 总是包含

The relative error between two consecutive iterates is at most 0.000000

在使用我的模拟噪声数据时，我玩弄了两个参数的真实值处于相同数量级（因为我认为使用的步长可能只会对其中一个产生明显影响）数量级，我改变了步长（参数epsfcn），但无济于事。

有谁知道我可能做错了什么，或者我可以使用什么拟合函数来代替leastsq 的？如果是这样：提前非常感谢！

编辑

按照 Russ 的建议，我现在将提供有关如何执行模拟的一些细节：我们正在研究小分子与大分子的结合。这发生的概率取决于它们的相互亲和力（亲和力是要从实验数据中提取的值之一）。一旦发生结合，我们还模拟复合体再次分解需要多长时间（解离时间常数是我们感兴趣的第二个参数）。还有许多其他参数，但它们仅在计算预期信号强度时才相关，因此与实际模拟无关。

我们从给定数量的小分子开始，每个小分子的状态都会模拟多个时间步长。在每个时间步，我们使用亲和力值来确定该分子是否与大分子结合，以防它未结合。如果它已经绑定，我们使用解离时间常数和它已经绑定的时间量来确定它是否在此步骤中解离。

在这两种情况下，参数（亲和力、解离时间常数）都用于计算概率，然后将其与随机数（0 到 1 之间）进行比较，并在此比较中取决于小分子的状态（绑定/未绑定）更改。

编辑 2

没有明确定义的函数可以确定所得曲线的形状，而且，即使形状明显可重现，但每个单独的数据点都存在随机性。 因此，我现在尝试使用optimize.fmin 而不是leastsq，但它不会收敛并在执行最大迭代次数后简单地退出。

编辑 3

根据 Andrea 的建议，我上传了 sample plot。我真的不认为提供样本数据会有很大帮助，它只是每个 x 值（时间）的一个 y 值（信号强度）......

【问题讨论】：

如果问题仍然存在，如果您可以展示您正在使用的数据和拟合过程的示例，将会有所帮助。
它看起来像一个有趣的应用程序，你应该放一些图或一些示例数据，这样你的问题会得到更多的关注。
@Andrea Zonca 正如你所建议的，我已经上传了一个情节。不过，我省略了样本数据，因为每个 y 值只有一个 x 值，这可能不会告诉任何人太多...谢谢！
你上传的情节是404。

标签： python scipy scientific-computing curve-fitting least-squares

【解决方案1】：

为了用任意函数拟合您的数据，您通常需要 Levenberg–Marquardt 算法，这就是 scipy.optimize.leastsq 使用的，因此您很可能使用正确的函数。

查看this page 第 5.4 节中的教程，看看是否有帮助。

也有可能你的底层函数难以适应（函数是什么？），但听起来你可能有其他问题。

此外 - 与 StackOverflow 一样出色，通过直接向 Scipy-User mailing list 发布一些示例代码和更多详细信息，您可能会获得更多知识渊博的 scipy 问题回复。

【讨论】：

非常感谢您的回复。回答您的问题：没有基础功能（或者至少我们不知道）。我们通过多个步骤获得了模拟曲线（模拟小分子与蛋白质的结合以及该复合物的解离）。
所以...听起来你不是曲线拟合！？听起来您正在尝试了解您的模拟数据与实际数据的“接近程度”吗？ OTOH...您还试图从比较中提取两个参数...听起来确实像曲线拟合。你能想出一个表示模拟中多个步骤的结果的表达式吗？如果是这样，那就扔LM。它需要它。如果没有，我不知道该说什么。这可能更适合不同的论坛，因为它不是一个编程问题。也许是统计数据或数学网站？您还需要更多细节。
是的，也许我应该到其他论坛寻求帮助，但这对我来说是全新的领域，所以我真的不知道下一步该做什么（也不知道我应该在我的题...）。不过谢谢！
我会尽可能详细地说明您将如何获得模拟数字。即：您的蒙特卡洛流程是什么样的，您有多少输入参数以及何时在您的步骤中引入它们。然后指出您尝试从一组真实数据中获取swag 的哪些参数。祝你好运！

【解决方案2】：

不完全是答案，但也可以尝试 PyMinuit。

http://code.google.com/p/pyminuit/

您想要做的是将您的 pdf 和数据点转换为 chi^2 或 -ln(likelihood) 或您选择的指标，并使用 PyMinuit 最小化该指标。可以将其配置为非常详细，这样您就可以找出哪里出了问题（如果确实出了问题）。

【讨论】：

非常感谢您告诉我有关此软件包的信息，我玩过它，它似乎确实可以完成工作。我认为我仍然需要更多地熟悉它，但看起来你让我们的项目重回正轨！我真的非常感谢你。

【解决方案3】：

如果您不知道全局的预期函数形式，但可以在给定系统当前状态的情况下预测“下一个”点的预期值，您可以考虑使用 @987654321 @（是的，“过滤器”在合适的上下文中听起来很愚蠢，但名称是历史性的，现在不能轻易更改）。

基础数学看起来有点吓人，但重要的一点是您不必理解它。您通常需要能够

定义表示空间
在表示空间中表达您的数据（模拟或实验）。
定义一个从给定的获取“下一个”表示的更新过程
从拟合器返回的序列表示中提取所需曲线

似乎至少有one existing python package 支持这一点（注意这里的界面与我习惯的不同，我无法提供太多关于使用它的建议）。

【讨论】：

非常感谢您的回复，我会看看您的建议！

【解决方案4】：

因为你只有两个参数，你应该做一个网格搜索。

results = {}
for p0 in parameter_space_for_first_parameter:
     for p1 in parameter_space_for_second_parameter:
           results[p0,p1] = compare(p0,p1)

如果您负担得起计算，compare 应该执行多次运行（使用不同的初始化）并计算均值和标准差。你可以尝试使用我的包jug 来管理你的计算（它就是为这类事情而设计的）。

最后，绘制结果，查看最小值（可能有几个）。这是一种“愚蠢”的方法，但它适用于其他方法卡住的许多情况。

如果计算量太大，您可以分两次执行：粗粒度网格，然后是接近粗粒度空间最小值的细粒度网格。

【讨论】：