python - python中使用leastsq同时拟合数据

Question

我已经很长时间没有编程了，而且从来都不擅长编程，但这是我正在努力解决的一项重要任务。我正在尝试拟合两组数据(x - 时间，y1 和 y2 - 应从文本文件读取的不同值列)。对于每个数据集(y1 和 y2)，我都有一个适合它们的函数。在这两个函数中，我有几个需要安装的参数。对于某些时间值，“y”的数据不存在，因此任务是在“y”缺失时以某种方式对其进行编程，并在没有这些值的情况下拟合数据。参数应该适合两个函数，因此应该同时拟合。这就是我无法解决的问题。

要在我的例子中定义函数，需要使用拉普拉斯逆变换，所以我使用了并且没有任何问题。

import numpy as np
import pylab as plt
from scipy.optimize import leastsq
from cmath import *

# Here is the Laplace functions 

def Fp(s, td, m0, kon, koff):
   gs=s+kon-kon*koff/(s+koff)
   sr=np.sqrt(gs*td)
   return (m0/(s*s))*sr/sinh(sr)

def Fd(s, td, m0, kon, koff):
   gs=s+kon-kon*koff/(s+koff)
   sr=np.sqrt(gs*td)
   fu=koff/(koff+kon)
   fs=fu+koff*(1-fu)/(s+koff)
   return (m0/s)*fs*2*tanh(0.5*sr)/sr

# Define the trig functions cot(phi) and csc(phi)

def cot(phi):
   return 1.0/tan(phi)
def csc(phi):
   return 1.0/sin(phi)

# inverse Laplace transform for two functions which are going to be fitted next

def Qpt(t, td, m0, kon, koff): 
    shift = 0.1;
    ans = 0.0;
    N=30
    h = 2*pi/N;
    c1 = 0.5017
    c2 = 0.6407
    c3 = 0.6122
    c4 = 0. + 0.2645j

    for k in range(0,N):
         theta = -pi + (k+1./2)*h;
         z = shift + N/t*(c1*theta*cot(c2*theta) - c3 + c4*theta); 
         dz = N/t*(-c1*c2*theta*(csc(c2*theta)**2)+c1*cot(c2*theta)+c4);
         ans += exp(z*t)*Fp(z, td, m0, kon, koff)*dz;
    return ((h/(2j*pi))*ans).real

def Qdt(t,td, m0, kon, koff): 
    shift = 0.1;
    ans = 0.0;
    N=30
    h = 2*pi/N;
    c1 = 0.5017
    c2 = 0.6407
    c3 = 0.6122
    c4 = 0. + 0.2645j
    for k in range(0,N):
        theta = -pi + (k+1./2)*h;
        z = shift + N/t*(c1*theta*cot(c2*theta) - c3 + c4*theta); 
        dz = N/t*(-c1*c2*theta*(csc(c2*theta)**2)+c1*cot(c2*theta)+c4);
        ans += exp(z*t)*Fd(z, td, m0, kon, koff)*dz;
    return ((h/(2j*pi))*ans).real

#now we have Qp and Qd as theoretical functions

我编译了它并询问程序几个值，Qp 和 Qd 定义正确。关于上面部分的唯一问题:是否有可能以某种方式同时定义两个函数而不进行两次变换？

然后我添加了同时拟合此功能的部分，它给了我一个错误:

TypeError: only length-1 arrays can be converted to Python scalars 

所以这是我的合适部分:

# FITTING PART

def residuals(pars, t, qd, qp):
    td = np.array(pars[0])
    m0 = np.array(pars[1])
    kon = np.array(pars[2])
    koff = np.array(pars[3])
    diff1 = Qdt(t,td, m0, kon, koff) - qd
    diff2 = Qpt(t,td, m0, kon, koff) - qp
    return np.concatenate((diff1[np.where(qd!=-1)], diff2[np.where(qp!=-1)]))

# for both functions with all the values
t = np.array([0.5, 2, 5, 10, 15, 20, 30, 40, 60, 90, 120, 180])
qd = np.array([0.145043746,0.273566338,0.437829373,0.637962531,-1,0.898107567,-1,1.186340492,1.359184345,-1,1.480552058,1.548143954])
qp = np.array([-1,-1,0.002701867,0.006485195,0.014034067,-1,0.06650739,-1,0.309055933,0.645945584,1.000811933,-1])


# initial values 
par_init = np.array([1, 1, 1, 1])

best, cov, info, mesg, ier = leastsq(residuals, par_init, args=(t, qd, qp), full_output=True)

print(" best-fit parameters: ", best)


#for each function separately to plot them and fitted functions as well
xqd= [0.5, 2, 5, 10, 20, 40, 60, 120, 180]
xqp= [5, 10, 15, 30, 60, 90, 120]
yqd= [0.145043746, 0.273566338, 0.437829373, 0.637962531, 0.898107567, 1.186340492, 1.359184345, 1.480552058, 1.548143954]
yqp= [0.002701867, 0.006485195, 0.014034067, 0.06650739, 0.309055933, 0.645945584, 1.000811933]

tt=np.linspace(0,185,100)
qd_fit=Qdt(tt,best[0], best[1], best[2], best[3])
qp_fit=Qdp(tt,best[0], best[1], best[2], best[3])

plt.plot(xqd,yqd,'bD:',xqp,yqp,'r^:', tt,qd_fit,'b',tt,qp_fit,'r')

plt.grid()
plt.show()

我将不胜感激任何帮助和建议!我迫切需要纠正这个错误!

提前致谢!

Answer 1

要同时将多个模型函数拟合到不同的数据集，您需要让残差函数执行以下操作:

1. 将所有可变参数放入第一个数组中——leastsq() 的第二个参数
2. 接收所有要拟合的数据(此处为 x、y1 和 y2)——必须与 lesssq() 的“args”参数匹配
3. 解压“当前参数值”数组(残差第一个参数)
4. 使用适当参数的值计算每个数据集的子残差
5. 将结果连接到单个一维数组作为返回值

一个简单的残差可能如下所示:

import numpy as np
from scipy.optimize import leastsq

def residual_two_functions(pars, x, y1, y2):
    off1   = pars[0]
    slope1 = pars[1]
    off2   = pars[2]
    slope2 = pars[3]
    diff1 = y1 - (off1 + slope1 * x)
    diff2 = y2 - (off2 + slope2 * x)
    return np.concatenate((diff1, diff2))

# create two tests data sets
NPTS = 201
x  = np.linspace(0, 10, NPTS)
y1 = -0.7 +  1.7*x + np.random.normal(scale=0.01, size=NPTS)
y2 =  5.2 + 50.1*x + np.random.normal(scale=0.02, size=NPTS)

# initial values 
par_init = np.array([0, 1, 10, 100])

best, cov, info, message, ier = leastsq(residual_two_functions,
                                        par_init, args=(x, y1, y2),
                                        full_output=True)

print(" Best-Fit Parameters: ",  best)