变参数复依赖函数的scipy.optimize.curve_fit

Question

我对于使用Python来适应我的数据还比较陌生，所以请原谅我缺乏编程技巧。然而，我一直无法找到解决我目前的曲线拟合尝试所造成的错误的方法。我认为这些误差是由于我的模型函数对两个变量中的一个(即Kd)的复杂依赖所致。我想知道是什么具体导致了这个问题，以及我如何调整我的定义或适合的软件包，以避免它。以下是最低限度的工作示例：

# Import libraries
import scipy as scipy
from scipy import stats
import numpy as np
from scipy.optimize import curve_fit

np.set_printoptions(precision=4)
ConcSyringeTotal = 9.5 ## total monomer concentration in the syringe [M]tot, in mM
Vinj = 10 ## volume injected in each injection, in uL
Vinit = 1250 ## volume of solvent initially in the sample cell, in uL
Vcell = 1000 ## cell volume, in uL (only the heat change within this volume is measured)
Injections = np.arange(2.00,26.00,1.00)
print Injections
Energy = np.array([136.953, 105.119, 84.414, 69.373, 60.898, 52.813, 46.187, 39.653, 33.894, 29.975, 27.315, 24.200, 21.643, 19.080, 16.158, 13.454, 13.218, 11.568, 10.742, 9.547, 8.693, 7.334, 6.111, 4.741])
print Energy

def DimerDissociation(injection, Kd, DHd): ## a dimer dissociation model for an ITC dilution experiment
    ## returns the heat flow (y-data, in ucal) per injection (x-data, unitless)
    ## fit for the dissociation constant (Kd, in mM = mmol/L, umol/mL, nmol/uL) 
    ## and the enthalpy of dissociation (DHd, in ucal/nmol = kcal/mol)
    ##
    ## concentration (in mM) of the free monomer in the cell after equilibration of the i-th injection
    VolumeAdded = 6+(injection-1)*Vinj ## in uL
    VolumeTotal = Vinit + VolumeAdded ## in uL
    CellTotal = ConcSyringeTotal*VolumeAdded ## Total in the cell after the i-th injection, in nmol
    ConcCellTotal = CellTotal/VolumeTotal ## Total concentration in the cell after the i-th injection, in mM
    ConcCellMonomer_roots =  np.roots([1, Kd/2, -Kd*ConcCellTotal/2]) 
    ConcCellMonomer_real = ConcCellMonomer_roots.real[abs(ConcCellMonomer_roots.imag)<1e-5]
    ConcCellMonomer_positive = ConcCellMonomer_real[ConcCellMonomer_real>0]
    ConcCellMonomer = ConcCellMonomer_positive[ConcCellMonomer_positive<ConcCellTotal]
    ##
    ## concentration (in mM) of the free monomer in the syringe
    ConcSyringeMonomer_roots =  np.roots([1, Kd/2, -Kd*ConcSyringeTotal/2]) 
    ConcSyringeMonomer_real = ConcSyringeMonomer_roots.real[abs(ConcSyringeMonomer_roots.imag)<1e-5]
    ConcSyringeMonomer_positive = ConcSyringeMonomer_real[ConcSyringeMonomer_real>0]
    ConcSyringeMonomer = ConcSyringeMonomer_positive[ConcSyringeMonomer_positive<ConcSyringeTotal]
    ## nmol of the free monomer injected from the syringe
    SyringeMonomerInjected = Vinj*ConcSyringeMonomer[0]
    ##
    ## concentration (in mM) of the free monomer in the cell before the i-th injection
    VolumeAddedPre = 6+(injection-2)*Vinj
    VolumeTotalPre = Vinit + VolumeAddedPre
    CellTotalPre = ConcSyringeTotal*VolumeAddedPre
    ConcCellTotalPre = CellTotalPre/VolumeTotalPre
    ConcCellMonomerPre_roots =  np.roots([1, Kd/2, -Kd*ConcCellTotalPre/2]) 
    ConcCellMonomerPre_real = ConcCellMonomerPre_roots.real[abs(ConcCellMonomerPre_roots.imag)<1e-5]
    ConcCellMonomerPre_positive = ConcCellMonomerPre_real[ConcCellMonomerPre_real>0]
    ConcCellMonomerPre = ConcCellMonomerPre_positive[ConcCellMonomerPre_positive<ConcCellTotalPre]
    ## nmol of the free monomer in the cell before the i-th injection
    CellMonomerPre = VolumeTotalPre*ConcCellMonomerPre[0]
    ##
    ## concentration of the free monomer before equilibration of the i-th injection, in mM
    ConcCellMonomerBefore = (CellMonomerPre+SyringeMonomerInjected)/VolumeAdded
    ## concentration of the free monomer after equilibration of the i-th injection, in mM
    ConcCellMonomerAfter = ConcCellMonomer[0]
    ## change in concentration of the free monomer over the equilibration of the i-th injection, in mM
    ConcCellMonomerChange = ConcCellMonomerAfter - ConcCellMonomerBefore
    ##
    return Vcell*DHd*ConcCellMonomerChange
DimerDissociation_opt, DimerDissociation_cov = curve_fit(DimerDissociation, Injections, Energy, p0=[0.4,10])
DimerDissociation_stdev = np.sqrt(np.diag(DimerDissociation_cov))
print "optimized parameters:", DimerDissociation_opt
print "covariance matrix:", DimerDissociation_cov
print "standard deviation of fit parameters:", DimerDissociation_stdev

以及相关的错误：

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-38-b5ef2361feed> in <module>()
     52     ##
     53     return Vcell*DHd*ConcCellMonomerChange
---> 54 DimerDissociation_opt, DimerDissociation_cov = curve_fit(DimerDissociation, Injections, Energy, p0=[0.4,10])
     55 DimerDissociation_stdev = np.sqrt(np.diag(DimerDissociation_cov))
     56 print "optimized parameters:", DimerDissociation_opt

//anaconda/lib/python2.7/site-packages/scipy/optimize/minpack.pyc in curve_fit(f, xdata, ydata, p0, sigma, absolute_sigma, **kw)
    553     # Remove full_output from kw, otherwise we're passing it in twice.
    554     return_full = kw.pop('full_output', False)
--> 555     res = leastsq(func, p0, args=args, full_output=1, **kw)
    556     (popt, pcov, infodict, errmsg, ier) = res
    557 

//anaconda/lib/python2.7/site-packages/scipy/optimize/minpack.pyc in leastsq(func, x0, args, Dfun, full_output, col_deriv, ftol, xtol, gtol, maxfev, epsfcn, factor, diag)
    367     if not isinstance(args, tuple):
    368         args = (args,)
--> 369     shape, dtype = _check_func('leastsq', 'func', func, x0, args, n)
    370     m = shape[0]
    371     if n > m:

//anaconda/lib/python2.7/site-packages/scipy/optimize/minpack.pyc in _check_func(checker, argname, thefunc, x0, args, numinputs, output_shape)
     18 def _check_func(checker, argname, thefunc, x0, args, numinputs,
     19                 output_shape=None):
---> 20     res = atleast_1d(thefunc(*((x0[:numinputs],) + args)))
     21     if (output_shape is not None) and (shape(res) != output_shape):
     22         if (output_shape[0] != 1):

//anaconda/lib/python2.7/site-packages/scipy/optimize/minpack.pyc in _general_function(params, xdata, ydata, function)
    443 
    444 def _general_function(params, xdata, ydata, function):
--> 445     return function(xdata, *params) - ydata
    446 
    447 

<ipython-input-38-b5ef2361feed> in DimerDissociation(injection, Kd, DHd)
     19     CellTotal = ConcSyringeTotal*VolumeAdded ## Total in the cell after the i-th injection, in nmol
     20     ConcCellTotal = CellTotal/VolumeTotal ## Total concentration in the cell after the i-th injection, in mM
---> 21     ConcCellMonomer_roots =  np.roots([1, Kd/2, -Kd*ConcCellTotal/2])
     22     ConcCellMonomer_real = ConcCellMonomer_roots.real[abs(ConcCellMonomer_roots.imag)<1e-5]
     23     ConcCellMonomer_positive = ConcCellMonomer_real[ConcCellMonomer_real>0]

//anaconda/lib/python2.7/site-packages/numpy/lib/polynomial.pyc in roots(p)
    199     """
    200     # If input is scalar, this makes it an array
--> 201     p = atleast_1d(p)
    202     if len(p.shape) != 1:
    203         raise ValueError("Input must be a rank-1 array.")

//anaconda/lib/python2.7/site-packages/numpy/core/shape_base.pyc in atleast_1d(*arys)
     47     res = []
     48     for ary in arys:
---> 49         ary = asanyarray(ary)
     50         if len(ary.shape) == 0 :
     51             result = ary.reshape(1)

//anaconda/lib/python2.7/site-packages/numpy/core/numeric.pyc in asanyarray(a, dtype, order)
    512 
    513     """
--> 514     return array(a, dtype, copy=False, order=order, subok=True)
    515 
    516 def ascontiguousarray(a, dtype=None):

ValueError: setting an array element with a sequence.

Answer 1

问题是，numpy.curve_fit将xdata作为数组传递给目标函数。这意味着injection上的所有DimerDissociation操作实际上都是数组操作。因此，ConcCellTotal也是一个数组(通过在代码中的第27行插入print type(ConcCellTotal)来检查这一点)。这意味着对np.roots的调用看起来像np.roots([scalar, scalar, array])，这是错误的来源。

当我处理这些东西时，我总是会掉头，但我认为优化器的目标函数应该是完全向量化的；每次调用时，它都需要返回一个数组，每个注入值都有一个能量值。

我通过显式地将ConcCellMonomer_roots作为数组修正了上面的错误，我还提交了一些关于变量状态的幼稚报告：

def DimerDissociation(injection, Kd, DHd): 
    print 'Called DimerDissociation'
    VolumeAdded = 6.0+(injection-1.0)*Vinj ## in uL
    VolumeTotal = Vinit + VolumeAdded ## in uL
    CellTotal = ConcSyringeTotal*VolumeAdded ## Total in the cell after the i-th injection, in nmol
    ConcCellTotal = CellTotal/VolumeTotal ## Total concentration in the cell after the i-th injection, in mM
    print 'total\t',np.shape(ConcCellTotal)
    ConcCellMonomer_roots =  np.asarray([np.roots([1.0, Kd/2.0, -Kd*i/2.0]) for i in ConcCellTotal])
    print 'roots\t',np.shape(ConcCellMonomer_roots)
    ConcCellMonomer_real = ConcCellMonomer_roots.real[abs(ConcCellMonomer_roots.imag)<1e-5]
    print 'real\t',np.shape(ConcCellMonomer_real)
    ConcCellMonomer_positive = ConcCellMonomer_real[ConcCellMonomer_real>0]
    print 'positive\t',np.shape(ConcCellMonomer_positive)
    ConcCellMonomer = ConcCellMonomer_positive[ConcCellMonomer_positive<ConcCellTotal]
    print 'monomer\t',np.shape(ConcCellMonomer)

我还使用ConcCellMonomerPre_roots对np.asarray进行了相应的更正。对于这些编辑，我让优化器迭代几次，直到ConcCellMonomer_roots包含一些假想的值。一旦发生这种情况，ConCellMonomer_real就不再是ConcCellTotal的形状了，因此ConcCellMonomer_positive[ConcCellMonomer_positive<ConcCellTotal]行会抛出一个广播错误。对DimerDissociation的调用提供了以下输出：

Called DimerDissociation
total   (24,)
roots   (24, 2)
real    (48,)
positive(24,)
monomer (24,)

直到最后一次迭代：

Called DimerDissociation
total   (24,)
roots   (24, 2)
real    (4,)
positive(4,)
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "C:\Anaconda\lib\site-packages\spyderlib\widgets\externalshell\sitecustomize.py", line 540, in runfile
    execfile(filename, namespace)
  File "C:/Users/Devin/Documents/Python Scripts/SO.py", line 66, in <module>
    DimerDissociation_opt, DimerDissociation_cov = curve_fit(DimerDissociation, Injections, Energy, p0=[0.4,10])
  File "C:\Anaconda\lib\site-packages\scipy\optimize\minpack.py", line 533, in curve_fit
    res = leastsq(func, p0, args=args, full_output=1, **kw)
  File "C:\Anaconda\lib\site-packages\scipy\optimize\minpack.py", line 378, in leastsq
    gtol, maxfev, epsfcn, factor, diag)
  File "C:\Anaconda\lib\site-packages\scipy\optimize\minpack.py", line 444, in _general_function
    return function(xdata, *params) - ydata
  File "C:/Users/Devin/Documents/Python Scripts/SO.py", line 35, in DimerDissociation
    ConcCellMonomer = ConcCellMonomer_positive[ConcCellMonomer_positive<ConcCellTotal]
ValueError: operands could not be broadcast together with shapes (4) (24) 

希望这能让你走上正轨，尽管我不是这里的专家，其他人可能会有更好的想法。

Answer 2

我无法重现你的错误。我注意到的第一个问题是np.roots的使用。roots(p)返回由p (特别是p[0] + p[1] * x + p[2] * x**2 + ... )中的系数指定的多项式的根。第三个系数，-Kd*ConcCellTotal/2是injections的一个函数，它是一个数组。np.roots没有文档化的签名，可以将数组作为p的成员传递。

你能编辑和澄清吗？

-Ravi

一个演示curve_fit如何工作的玩具示例：

import numpy as np
from scipy.optimize import curve_fit

x_in = np.array([-3.0,-2.0,-1.0,0.0,1.0,2.0,3.0])

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

y_in = f(x_in,3,2)
parameters_fit,cov = curve_fit(f,x_in,y_in)
y_out = parameters_fit[0]*x_in+parameters_fit[1]
print parameters_fit
print y_in
print y_out

y_in = f(x_in,10,75)
parameters_fit,cov = curve_fit(f,x_in,y_in)
y_out = parameters_fit[0]*x_in+parameters_fit[1]
print parameters_fit
print y_in
print y_out

目标函数f以x值数组和一个或多个参数作为参数.curve_fit以目标函数、x值x_in数组和y值y_in数组作为参数。然后给出了参数a和b的一些值，并在x_in上对目标函数进行了评估，给出了一个数组y_out。它计算、y_in、和y_out之间的均方根误差，然后调整其、a、和b的值，直到均方根误差最小化。

关键在于a和b的初始值是如何被选择的(如果没有像OP那样提供它们)，以及它们是如何被调整的。这很复杂，但对我们scipy.optimize用户来说，完全理解并不是绝对必要的。