Python、Lmfit、globalfit

Question

我正在尝试使用1个共享参数和6个其他参数进行全局拟合。dataComplete是我从文本文件解析的数据，剂量=x值。在这里，bValue随着浓度的变化而变化，而其他的保持不变。在进行全局拟合之前，有一些使用这些参数的计算。我得到两个错误

1: nan值和2: TypeError:输入不正确: N=35不能超过M=22

我检查了各个元素，它们似乎可以工作，但不是全局拟合。任何帮助都将不胜感激。我试图在这里使用之前解决的示例中的代码，但显然有些地方不对劲。

##file parser
## coded elsewhere

dataComplete = [[90.47, 91.6, 89.83, 68.58, 56.68, 41.93, 1.09, 17.0, 9.19, 1.84, 0.1, 0.38, -0.37, 0.47, 0.65, -0.01, 0.36, 0.18, 0.29, 0.12, 0.87, -1.26], [108.07, 115.49, 102.48, 104.09, 110.15, 72.64, 49.11, 37.35, 24.94, 14.33, 6.31, 3.07, 9.81, 0.2, 0.74, 0.59, 1.11, 0.44, 0.34, 1.03, 0.06, 1.08], [100.93, 118.82, 89.76, 116.32, 127.06, 97.5, 69.71, 47.17, 44.16, 22.15, 18.6, 9.64, 2.56, 2.02, 0.93, 0.29, 1.11, -0.04, 0.27, -0.23, 1.0, 0.18], [99.87, 110.15, 92.21, 92.59, 96.8, 92.69, 81.39, 70.09, 44.08, 39.48, 7.66, 21.49, 2.96, 5.04, 1.87, 2.49, 1.26, 0.53, 0.65, 0.97, 1.15, 0.52], [102.14, 97.94, 94.85, 94.48, 101.22, 94.81, 111.42, 94.96, 61.2, 33.17, 29.37, 11.33, 23.67, 15.97, 6.52, 2.4, 2.29, 0.27, 0.08, 0.31, 0.09, -0.09]]

dose = [-5.0, -5.30103, -5.60206, -5.90309, -6.20412, -6.50515, -6.80618, -7.10721, -7.40824, -7.70927, -8.0103, -8.31133, -8.61236, -8.91339, -9.21442, -9.51545, -9.81648, -10.11751, -10.41854, -10.71957, -11.0206, -11.32163]



def objective(params,xdata,ydata):
    resid = 0.0*ydata[:]
    for i in range(5):
        logEc50 = 10**(params['logEc50_%i' %(i+1)].value)
        bValue = params['bValue_%i' %(i+1)].value
        pA2 = params['pA2_%i' %(i+1)].value
        schildSlope = params['schildSlope_%i' %(i+1)].value
        top = params['top_%i' %(i+1)].value
        bottom = params['bottom_%i' %(i+1)].value
        hillSlope = params['hillSlope_%i' %(i+1)].value

        ec50 = 10**(logEc50)
        Antag = 1+(bValue/(10**(-1*pA2)))**schildSlope
        LogEC=np.log((ec50*Antag))

        resid[i, :] =bottom + (top-bottom)/(1+10**((LogEC-dose)*hillSlope))
    return resid[i, :] - ydata[i, :]


params = Parameters()
molar = [1.25e-7,2.5e-8,1e-8,5e-9,0]

for i in range(5):
    params.add('logEc50_%i' % (i+1),value=-6.96)
    params.add('bValue_%i' % (i+1),value=molar[i])
    params.add('pA2_%i' % (i+1),value=7.48)
    params.add('schildSlope_%i' % (i+1),value=1)
    params.add('bottom_%i' % (i+1),value=-0.2)
    params.add('top_%i' % (i+1),value=109.496)
    params.add('hillSlope_%i' % (i+1),value=1)


#data
dataComplete = [[]]*5
dataComplete[0]=conc1
dataComplete[1]=conc2
dataComplete[2]=conc3
dataComplete[3]=conc4
dataComplete[4]=conc5


print dataComplete
dataComplete = np.array(dataComplete)
assert(dataComplete.shape) ==(5,22)

#x data
dose = np.array(dose)

#fitting
result = minimize(objective,params,args=(dose,dataComplete))

lmfit.printfuncs.report_fit(result.params)

Answer 1

我猜你是在用10和对数的很多次幂来解决数字问题。例如，with

logEc50_i ~= 7

然后将其升幂为10的幂

logEc50 = 10**(params['logEc50_%i' %(i+1)].value)

然后将其升幂为10的幂：

ec50 = 10**(logEc50)

(这样做两次只是一个错误吗？)

然后你将这个乘以Antag，这是两个数字~1e-7的比率，应该得到接近于1的值，然后取它的对数：

LogEC=np.log((ec50*Antag))

我认为您真的想使用对数的属性，比如

log(a*b) = log(a) + log(b)
log(x**y) = y*log(x)

这是为了防止数值不稳定。

我没有看到哪个参数是要在数据集之间共享的，但是看起来您确实打算适合5个单独的数据集，每个数据集有7个变量--其中一个参数是否应该约束为所有数据集都具有相同的值？如果你总共只有22个观察值，那么无论是35个变量还是31个变量，都无关紧要，仍然没有足够的总观察值来拟合那么多变量。

但看起来五个数据集中的每一个都有22个观察值。我认为您可能会遇到的是，您需要评估每个数据集的残差，然后将结果连接起来，以便所有五个数据集的resdual数组都有110个元素。