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

我非常感谢任何帮助和建议！我迫切需要解决这个错误！

提前谢谢你们！

数值分析数据预处理缺失值处理数据拟合参数优化最小二乘法逆拉普拉斯变换同时拟合

在Python中使用leastsq进行数据拟合

1 个回答

撰写回答