使用 scipy.curve_fit 进行全局拟合

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 值。这是局部拟合，现在我想将其更改为全局拟合。根据我对全局拟合的理解，我尝试了以下脚本：

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

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)

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

popt, pcov = curve_fit(glob_fun, x, y)
return func(xdata, *params) - ydata
TypeError: unsupported operand type(s) for -: 'function' and 'float'

这里的问题是 glob_fun 实际上是一个函数数组（据我了解，它应该是全局拟合）。但是，似乎不是使用该函数的输出（基于它为 kD 选择的任何内容），而是将其最小化为 ydata，而是使用数组本身的函数之一。因此，您不能减去函数的错误（或者至少，这是我对错误的理解）。

编辑： 我已经添加了数据，因此错误和功能是可重现的。

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

您已经成功地对单个数据集进行了拟合。现在，您希望同时对多个数据集执行相同函数的全局拟合。数据集位于多维数组中，其中来自先前执行的成功单次拟合的每个数据集沿内轴运行。但是，scipy.optimize.curve_fit 预计

一个长度为 M 的数组

因为它的论点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,) 和y，形状为(n, m)，带有随机噪声。 （示例数据正在打印中，如果您想看一下。）
• 然后，y 中的每一行 y_i 正在本地拟合，使用函数 f。 （这对于全局拟合不是必需的，但很高兴在图中看到结果线以进行比较。）
• 最后，全局拟合整个y：我们必须使用函数lambda x, a, b: np.tile(f(x, a, b), len(y)) 代替lambda x, a, b: np.tile(f(x, a, b), len(y))，它将f 应用于x 并将结果重复len(y) 次（因为y 中有n 或len(y) 行适合，每个数据集一个）使用np.tile。随后，相同的a 和b 用于y 中的每一行，我们得到一个全局拟合。 （与单独的 a 和 b 相比，每个单独拟合每个数据集。）

import matplotlib.pyplot as plt
import numpy as np
from scipy.optimize import curve_fit

m = 5
n = 3
x = np.arange(m)
y = np.array([x + np.random.normal(0, 0.2, len(x)) for _ in range(n)])
print("x =", x)
print("y =", y)

def f(x, a, b):
    return a * x + b

# single fits to each dataset
for y_i in y:
    popt, pcov = curve_fit(f, x, y_i)
    plt.plot(x, y_i, linestyle="", marker="x")
    plt.plot(x, f(x, *popt), color=plt.gca().lines[-1].get_color())

# global fit to concatenated dataset
popt, pcov = curve_fit(lambda x, a, b: np.tile(f(x, a, b), len(y)), x, y.ravel())
plt.plot(x, f(x, *popt), linestyle="--", color="black")

plt.show()

例如导致：

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())