用Python拟合模拟和实验数据点
我写了一些代码,用来做蒙特卡洛模拟,生成信号强度和时间的曲线。这种曲线的形状取决于几个参数,其中有两个参数是我的合作者想通过“真实实验”来确定的,我正在模拟这个实验。
我们准备将她的实验数据与我模拟的曲线进行比较,但现在我遇到了麻烦,因为我还没有能够进行任何拟合(目前我只是用模拟的噪声数据来测试)。我尝试使用了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 个回答
如果你不知道全局的预期功能形式,但可以根据系统当前的状态预测“下一个”点的预期值,那么你可以考虑使用卡尔曼滤波器(虽然“滤波器”这个名字听起来有点奇怪,但这是历史原因,不能轻易改变)。
底层的数学可能看起来有点吓人,但重要的是你不需要完全理解它。你通常需要做到以下几点:
- 定义一个表示空间
- 在这个表示空间中表达你的数据(可以是模拟数据或实验数据)。
- 定义一个更新过程,从给定的表示中获取“下一个”表示
- 从拟合器返回的一系列表示中提取出你想要的曲线
似乎至少有一个现有的Python库可以支持这个功能(注意,这里的接口和我熟悉的不同,我不能提供太多使用建议)。
如果你想用一个任意的函数来拟合你的数据,通常会用到 Levenberg–Marquardt 算法,而 scipy.optimize.leastsq
就是使用这个算法,所以你很可能在用对的函数。
你可以看看 这页 第5.4节的教程,看看是否对你有帮助。
也有可能你的基础函数比较难拟合(你用的是什么函数?),但听起来你可能还有其他问题。
另外,虽然 StackOverflow 很棒,但如果你在 scipy 方面有问题,直接去 Scipy-User 邮件列表 发帖,附上示例代码和更多细节,可能会得到更专业的回答。
这不是一个直接的答案,但你可以试试PyMinuit。
http://code.google.com/p/pyminuit/
你想做的事情是把你的pdf和数据点转换成chi^2或者-ln(似然)或者你选择的其他指标,然后用PyMinuit来最小化这个指标。它可以设置得非常详细,这样你就能找到问题出在哪里(如果真的出错了的话)。