全局拟合使用曲线曲线图

Question

我有一个关于使用scipy.optimize.curve_fit进行全局拟合的快速问题。根据我的理解，在设置脚本的局部拟合和全局拟合之间的唯一区别是连接函数的不同。以下面的脚本为例：

input_data = [protein, ligand]
titration_data=input('Load titration data')

def fun(_, kd):
    a = protein
    b = protein + ligand
    c = ligand
    return np.array((b + kd - np.sqrt(((b + kd)**2) - 4*a*c))/(2*a))

kD=[]
for values in titration_data:
    intensity=[values]
    intensity_array=np.array(intensity)
    x = ligand
    y = intensity_array.flatten()
    popt, pcov = curve_fit(fun, x, y)

输入数据是6x2矩阵，滴定数据也是8x6矩阵。将每一行滴定数据分别拟合到模型中，得到kd值。这是局部匹配，现在我要将其更改为全局匹配。基于我对全局拟合的理解，我尝试了以下脚本：

^{pr2}$

根据我的理解，这应该给我一个单一的kd输出，同时拟合所有的数据。但是，我在尝试实现这一点时遇到了一条错误消息：

^{3}$

这里的问题是glob_fun实际上是一个函数数组（据我的理解，对于全局拟合来说应该是这样）。然而，它似乎并没有使用该函数的输出（基于它为kD选择的任何内容），而是使用了数组本身的一个函数。因此，你不能减去一个函数的错误（或者至少，这是我对错误的理解）。在

编辑： 我添加了数据，这样误差和函数是可重复的。在

import numpy as np
from scipy.optimize import curve_fit

concentration= np.array([[0.6 , 0.59642147, 0.5859375 , 0.56603774, 0.53003534,0.41899441],
[0.06 , 0.11928429, 0.29296875, 0.62264151, 1.21908127,3.05865922]])
protein = concentration[0,:]
ligand = concentration[1,:]

input_data = [protein, ligand]
titration_data=np.array([[0, 0, 0.29888413, 0.45540198, 0.72436899,1],
 [0,0,0.11930228, 0.35815982, 0.59396978, 1],
 [0,0,0.30214337, 0.46685577, 0.79007708, 1],
 [0,0,0.27204954, 0.56702549, 0.84013344, 1],
 [0,0,0.266836,   0.43993175, 0.74044123, 1],
 [0,0,0.28179148, 0.42406587, 0.77048624, 1],
 [0,0,0.2281092,  0.50336244, 0.79089151, 0.87029517],
 [0,0,0.18317694, 0.55478412, 0.78448465, 1]]).flatten()

glob=[]
for values in titration_data:
    def fun(_, kd):
        a = protein
        b = protein + ligand
        c = ligand
        return np.array((b + kd - np.sqrt(((b + kd)**2) - 4*a*c))/(2*a))
        print (fun)
    glob.append(fun)

def glob_fun(_,kd):
  return np.array(glob).flatten()

x = ligand
y = titration_data
popt, pcov = curve_fit(glob_fun, x, y)

Answer 1

您已经成功地执行了对单个数据集的拟合。现在，您需要同时对多个数据集执行同一函数的全局拟合。数据集在一个多维数组中，以前执行的成功的单次拟合的每个数据集沿内轴运行。但是，^{}需要

a length M array

对于其参数ydata。据我所知，这意味着您将无法使用[[0], [1]]，例如：

>>> from scipy.optimize import curve_fit
>>> curve_fit(lambda x, a: x, [[0], [1]], [[0], [1]])
ValueError: object too deep for desired array
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/home/user/.local/lib/python3.6/site-packages/scipy/optimize/minpack.py", line 744, in curve_fit
    res = leastsq(func, p0, Dfun=jac, full_output=1, **kwargs)
  File "/home/user/.local/lib/python3.6/site-packages/scipy/optimize/minpack.py", line 394, in leastsq
    gtol, maxfev, epsfcn, factor, diag)
minpack.error: Result from function call is not a proper array of floats.

正如您已经发现的，一个解决方案是将数组展平，这样来自每个单独拟合的每个数据集都一个接一个地串在一起。我想，这不是真正的“全局拟合”，而是“串联拟合”。在

为了演示如何使用curve_fit，我编写了以下最小的示例：

• 首先，我们创建一些带有随机噪声的x形状(m,)和形状(n, m)的y。（示例数据正在打印中，如果您想查看它的话。）
• 然后，y中的每一行y_i使用函数f进行局部拟合。（这不是全局拟合所必需的，但很高兴在绘图中看到结果线以便进行比较。）
• 最后，对整个y进行全局拟合：我们将不得不使用一个函数lambda x, a, b: np.tile(f(x, a, b), len(y))，它将f应用于x，并重复结果len(y)次（因为y中有n或{}行来适应，每个数据集对应一个）。接下来，相同的a和{}被用于y中的每一行，我们得到了一个全局拟合。（与单独的a和{}相比，对于每个数据集，每个单一拟合都是如此。）

^{pr2}$

结果例如：

x = [0 1 2 3 4]
y = [[ 0.17209542  1.02497865  1.84162787  3.0763016   3.76940871]
 [-0.05657471  0.96686915  2.20283785  3.09199915  3.78047165]
 [-0.53504594  1.21865205  2.35021432  3.02407509  4.22551247]]

标记的点是输入数据y，彩色线是对这些点（相同颜色）的单一拟合，黑色虚线是对所有组合点的全局拟合。在

将此示例应用于您的代码应该可以得到如下结果：

import numpy as np
from scipy.optimize import curve_fit

concentration = np.array(
    [
        [0.6, 0.59642147, 0.5859375, 0.56603774, 0.53003534, 0.41899441],
        [0.06, 0.11928429, 0.29296875, 0.62264151, 1.21908127, 3.05865922],
    ]
)

protein = concentration[0, :]
ligand = concentration[1, :]

titration_data = np.array(
    [
        [0, 0, 0.29888413, 0.45540198, 0.72436899, 1],
        [0, 0, 0.11930228, 0.35815982, 0.59396978, 1],
        [0, 0, 0.30214337, 0.46685577, 0.79007708, 1],
        [0, 0, 0.27204954, 0.56702549, 0.84013344, 1],
        [0, 0, 0.266836, 0.43993175, 0.74044123, 1],
        [0, 0, 0.28179148, 0.42406587, 0.77048624, 1],
        [0, 0, 0.2281092, 0.50336244, 0.79089151, 0.87029517],
        [0, 0, 0.18317694, 0.55478412, 0.78448465, 1],
    ]
)

def fun(_, kd):
    a = protein
    b = protein + ligand
    c = ligand
    return np.array((b + kd - np.sqrt(((b + kd) ** 2) - 4 * a * c)) / (2 * a))

def glob_fun(_, kd):
    return np.tile(fun(_, kd), len(titration_data))

x = ligand
y = titration_data
popt, pcov = curve_fit(glob_fun, x, y.ravel())