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

Question

我已经编写了一些代码，执行 Monte Carlo 模拟并生成与时间相关的信号强度曲线。这种曲线的形状取决于各种参数，其中两个是我的合作者想通过我正在模拟的实验的“真实版本”来确定的。

我们准备将她的实验数据与我的模拟曲线进行比较，但现在我卡住了，因为我还没有能够执行任何拟合（到目前为止，我已经用模拟的嘈杂数据替换了实验数据进行测试）。我尝试使用 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所建议的，我上传了一个样本图。我真的不认为提供样本数据会有很大帮助，因为每个x值（时间）只有一个y值（信号强度）...

- canavanin

---

如果这个问题还没有解决，如果您能展示一下您正在使用的数据和拟合程序的示例，那么会很有帮助。 - ev-br

看起来是一个有趣的应用程序，你应该添加一些图表或示例数据，这样你的问题会得到更多的关注。 - Andrea Zonca

@Andrea Zonca 我已经按照你的建议上传了一个图表。不过我省略了样本数据，因为每个y值只有一个x值，这可能无法提供太多信息...谢谢！ - canavanin

您上传的情节不存在（404）。 - xApple

Answer 1

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

请查看此页面第5.4节中的教程，看看它是否有帮助。

另外，您的基础函数可能很难拟合（这个函数是什么？），但听起来您可能还有其他问题。

此外 - 尽管StackOverflow很棒，但直接向Scipy-User邮件列表发布一些示例代码和更多细节，您可能会得到更专业的回答。

- Russ

---

非常感谢您的回复。针对您的问题：目前我们并不知道是否存在底层函数。我们是通过多个步骤（模拟小分子与蛋白质的结合以及该复合物的解离）来获取模拟曲线的。 - canavanin

所以...听起来你并不是正在进行曲线拟合！？听起来你是在尝试获得一个"有多接近"你的模拟数据和实际数据的图形吗？另一方面...你也试图从比较中提取出两个参数...这确实听起来像是曲线拟合。你能否想出一个表示模拟中多个步骤结果的表达式呢？如果可以，那就使用LM。它需要这样做。如果不能，我不知道该说什么。这可能更适合在其他论坛讨论，因为它不太是一个编程问题。也许是统计或数学网站？此外，你需要更多的细节。 - Russ

是的，也许我应该寻求其他论坛的帮助，但这对我来说完全是新的领域，所以我不知道接下来该怎么做（也不知道应该在我的问题中包含哪些信息...）。但还是谢谢！ - canavanin

我会尽可能详细地说明你在模拟数据时的步骤，包括蒙特卡洛模拟的流程以及输入参数的数量和引入步骤。然后，请指出你试图从一组真实数据中估计哪些参数。祝好运！ - Russ

Answer 2

如果您不知道全局期望的函数形式，但可以预测在当前系统状态下“下一个”点的期望值，那么您可以考虑使用卡尔曼滤波器（是的，在拟合上下文中，“滤波器”听起来有些傻，但这个名称是历史遗留问题，现在很难更改）。

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

1. 定义表示空间
2. 将数据（模拟或实验）表达为表示空间中的数据。
3. 定义更新过程，从给定的表示中获取“下一个”表示
4. 从由拟合器返回的系列表示中提取所需的曲线

似乎至少有一个现有的Python软件包支持此功能（请注意，此处的接口与我所熟悉的接口不同，因此我无法提供太多关于使用它的建议）。

- dmckee --- ex-moderator kitten

---

非常感谢您的回复，我会查看您所建议的内容！ - canavanin

Answer 3

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

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

你想做的是将你的PDF和数据点转换为卡方（chi^2）或-ln(似然度)或你选择的其他指标，并使用PyMinuit来最小化该指标。它可以被配置为非常详细，以便您可以找出哪里出了问题（如果出了问题）。

- Piti Ongmongkolkul

---

非常感谢你告诉我这个包，我已经试用了一下，它似乎真的能胜任这项工作。我还需要再熟悉一下它，但看起来你把我们的项目带回正轨了！我真的非常感激你。 - canavanin

Answer 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 来管理你的计算（它专门为这种情况设计）。

最后，绘制结果，查看最小值（可能有多个）。这是一种“愚蠢”的方法，但在许多其他方法卡住的情况下也会起作用。

如果这个计算量太大，你可以分两步进行：先进行粗略的网格搜索，然后在粗略空间的最小值附近进行细粒度搜索。

- luispedro