使用Lmfit拟合多个数据集，每个数据集具有共享参数和一个固定参数

Question

我是python的初学者，我遇到了一个问题。我想使用lmfit模块将一个函数拟合到来自.csv文件的3组(x，y)数据，以及一些共享参数(a，b，d)和一个固定的单独参数(c)。我的数据格式如下: x1 y1 x2 y2 x3 y3 ..。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。

。。。。。。。。。。。。

这里是一个指数函数的例子，以及我尝试拟合到每个数据集和绘图的尝试。我在定义全局残差和为c参数的每个数据集定义固定值时遇到了问题。这是我到目前为止所得到的：

 import numpy as np
import pandas as pd
from lmfit import minimize, Parameters, report_fit
import matplotlib.pyplot as plt


#Load data
df=pd.read_csv('mydata.csv')
df[['x1', 'y1', 'x2', 'y2', 'x3','y3']][:6]

# set up the data
xs = np.array(df[df.columns[0::2]]) 
ys = np.array(df[df.columns[1::2]]) 

#define function
def biexp(xs,a,b,c,d):
    return a*np.exp(-xs*b*c)+(1-a)*np.exp(-xs*d*c)
 
 #Define a function that calculates biexp for data set i
def biexp_dataset(params, i, xs):
    a = params['a_%i' % (i+1)].value
    b = params['b_%i' % (i+1)].value
    c = params['c_%i' % (i+1)].value
    d = params['d_%i' % (i+1)].value
    return biexp(xs,a,b,c,d)
    
#Define the real function to minimize, which calculates the total residual for all fits, modeled by biexp function
"""Where I'm stuck"""
def objective(params, xs,ys):
    nys, nxs= df.shape[1]/2
    resid = 0.0*ys[:]
    # make residual per data set
    for i in range(nys):
       resid = ys[:, i] - biexp_dataset(params,i,xs)
    # now flatten this to a 1D array, as minimize() needs
    return resid.flatten()

#Define parameters
fit_params = Parameters()
for i in range(df.shape[1]/2): 
    fit_params.add( 'a_%i' % (i+1), value=1, min=0,  max=9,vary=True)
    fit_params.add( 'b_%i' % (i+1), value=0.3, min=0.0, vary=True)
    fit_params.add( 'd_%i' % (i+1), value=0.5, min=0.0,  vary=True)
"""How can I define c for each data sets?"""
"""Iwant for x1,y1 => c1=2; x2,y2=>c2=0.4; x3,y3=>c3=5, for example"""

#  constrain all values of a,b and d to have the same value
for i in (2, 3):
    fit_params['a_%i' % i].expr='a_1'
    fit_params['b_%i' % i].expr='b_1'
    fit_params['d_%i' % i].expr='d_1'

      
# run the global fit to all the data sets
result=minimize(objective, fit_params, args=(xs, ys))
report_fit(result.params)

# plot the data sets and fits
plt.figure()
for i in range(df.shape[1]/2):
    y_fit = biexp_dataset(fit_params, i, xs)
    plt.plot(xs[:,i], ys[:, i], 'o', xs[:,i], y_fit, '-')

plt.show()

所以基本上我要问的是:如何为每个数据集定义一个固定的c参数值？c1=0.4；例如，c2=1;c3=5

我在目标定义中做错了什么，因为当我运行代码时它指出: TypeError:不能将序列乘以'float‘类型的非整型

任何帮助都将不胜感激..。

更新

因此，我使用我想要拟合的真实函数编写了一些新代码(并且我已经使用originLab中的全局拟合选项拟合了它)。因此，我的数据集共享相同的D1、D2、tau1和tau2。然而，该模型作为一个新的固定变量，即实验时间(t_exps)。现在我在代码中没有任何错误，但我得到的不是6个数据集的6条拟合曲线，而是36条曲线拟合。我迷路了..。我的代码出了什么问题：

import numpy as np
import pandas as pd
from lmfit import minimize, Parameters, report_fit
import matplotlib.pyplot as plt

#set up the data
df=pd.read_excel('BTMS_hex.xlsx')

df = df.dropna() #drop rows that have invalid values (extra line in data file)
# set up the data
xs = np.array(df[df.columns[0::2]]).astype(float)
ys = np.log(np.array(df[df.columns[1::2]]).astype(float))
t_exps = [0.04857, 0.0989, 0.14923, 0.1993, 0.49957, 0.99957]


#Create function
def karger(xs,D1,D2,tau1,tau2,texp):
    term1=D1+D2+(1/(xs/texp))*((1/tau1)+(1/(tau2)))
    term2=np.sqrt(((D1-D2+(1/(xs/texp))*((1/tau1)+(1/(tau2))))**2)+(4/(((xs/texp)**2)*tau1*tau2)))
    DA=0.5*(term1-term2)
    DB=0.5*(term1+term2)
    P1=tau1/(tau1+tau2)
    P2=1-P1
    CB=(P1*D1+P2*D2-DA)/(DB-DA)
    return np.log((1-CB)*np.exp(-xs*DA)+CB*np.exp(-xs*DB))

def karger_dataset(params, xs,texp):
    D1 = params['D1'].value
    D2 = params['D2'].value
    tau1 = params['tau1'].value
    tau2 = params['tau2'].value 
    return karger(xs, D1,D2,tau1,tau2,texp)
    
def objective(params, xs,ys):
    resid = np.array([])
    for i in range(xs.shape[1]):
        y_pred = karger_dataset(params, xs[:,i],t_exps[i])
        resid=np.concatenate((resid,(ys[:, i] - y_pred)))     
    return resid

#Define parameters
fit_params = Parameters()
fit_params.add( 'D1', value=2E-1,min=0, vary=False)
fit_params.add( 'D2', value=1E-2, min=0, vary=True)
fit_params.add( 'tau1', value=8.2, min=0, vary=True)
fit_params.add( 'tau2', value=0.5, min=0,  vary=True)

# run the global fit to all the data sets
result=minimize(objective, fit_params, args=(xs, ys))
report_fit(result.params)

print result.residual

# plot the data sets and fits
plt.figure()
for i in range(df.shape[1]/2):
    y_fit = karger_dataset(fit_params, xs,t_exps[i])
    plt.plot(xs[:,i], ys[:, i], 'o', xs[:,i], y_fit, '-')

"plt.yscale('log')"
plt.axis([0,1600, -10,0.1])
plt.savefig('books_read.png')

Answer 1

为

所以基本上我要问的是：，我如何为每个数据集定义一个固定的c参数值？，

，is，set，for，set，for？，我该如何为每个数据集定义一个固定的c参数值？c1=0.4；例如，c2=1;c3=5

为什么不直接做呢？

fit_params.add('c1', value=0.4, vary=False)
fit_params.add('c2', value=1.0, vary=False)
fit_params.add('c3', value=5.0, vary=False)

在定义其他参数的for i in range(df.shape[1]/2):循环之后？

为

我在目标定义中做错了什么，因为当我运行代码时，它指出: TypeError:不能将序列乘以‘

’类型的非整型。

如果没有看到完整的回溯(这可能会告诉您发生异常的代码行)，我会猜测在某个地方有一个列表，而不是一个numpy数组乘以一个浮点数。

我不太理解为什么你会得到这个错误，而不是“我找不到参数'c1'”的错误，但这表明它可能在你开始运行你的fit之前就发生了。完整的回溯将澄清这一点。

Answer 2

这里有一个相关的问题和答案，可能对复习有用：

https://stats.stackexchange.com/questions/242729/chi-squared-of-simultaneous-least-squares-fitted-data/