3

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

我们已经准备好将她的实验数据与我的模拟曲线进行比较，但现在我被卡住了，因为我还不能进行任何拟合（到目前为止，我已经用模拟噪声数据替换了实验数据进行测试）。我尝试使用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 值（信号强度）......

4

4 回答 4

3

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

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

您的基础功能也可能难以适应（功能是什么？），但听起来您可能还有其他问题。

此外 - 与 StackOverflow 一样棒，通过直接发布到Scipy-User 邮件列表以及一些示例代码和更多详细信息，您可能会获得更多关于 scipy 问题的知识回复。

于 2011-09-16T13:23:06.733 回答
2

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

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

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

于 2011-09-26T07:36:24.397 回答
2

如果您不知道全局的预期函数形式，但可以在给定系统当前状态的情况下预测“下一个”点的预期值，您可能会考虑使用卡尔曼滤波器（是的，“滤波器”在合适的上下文，但名称是历史的，现在不能轻易更改）。

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

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

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

于 2011-09-21T15:07:41.670 回答
1

由于您只有两个参数，因此您应该只进行网格搜索。

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

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

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

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

于 2011-09-27T02:29:31.593 回答