与似然法一起使用的Scipy微分进化

问题描述

在我的工作中，我正在尝试使用基于可能性的方法将模型拟合到数据中。我以前在稍微不同的模型上使用过此代码，但是由于某种原因，在与该模型一起使用时，它一直引发此错误：“ RuntimeError：类似地图的可调用对象必须为f（func，iterable）形式，返回与“可迭代”长度相同的数字序列”。我不确定这是模型还是代码有问题，但是如果您可以帮助我了解此错误消息的含义以及如何解决，我将不胜感激。

import pandas as pd
import numpy as np
from scipy.integrate import odeint
from scipy.optimize import curve_fit
from lmfit import minimize,Parameters,Parameter,report_fit
#==============================================================================
''' Make new data matrix,same as csv except infected cells are one total for convience '''
DataFrame = pd.read_csv('Cell_count_data_Tromas_2014.csv') # Read data from file

TROMAS_DATA = np.empty((DataFrame.shape[0],5),int)
for i in range(DataFrame.shape[0]): # Number of rows
    for j in range(5):              # Number of columns desired
        TROMAS_DATA[i][j] = DataFrame.iloc[i,j]
        if j == 4:
            TROMAS_DATA[i][j] = DataFrame.ix[i,'Venus_only'] + DataFrame.ix[i,'BFP_only'] + DataFrame.ix[i,'Mixed']
'''
Col 0: Days post infection
Col 1: Leaf number
Col 2: Replicate plant number
Col 3: Number of unifected cells
Col 4: Number of total infected cells
'''
#==============================================================================
''' Make axis for negative log likelihood '''
ZERO_DAYS_AXIS = [0,3,5,7,10]
#==============================================================================
''' Parameter list for model '''
likeParams = Parameters()
likeParams.add('I0',value = .00372,min = .0000001,max = 1.0000)
likeParams.add('b',value = .5,min = .0001,max = 20.0000)
likeParams.add('x5',max = 20.0000)
likeParams.add('x6',max = 20.0000)
likeParams.add('x7',max = 20.0000)
likeParams.add('psi3',max = 1.0000)
likeParams.add('psi5',max = 1.0000)
likeParams.add('psi6',max = 1.0000)
likeParams.add('psi7',max = 1.0000)
#==============================================================================
def model(Mk,t,parameters):
    M3 = Mk[0]
    M5 = Mk[1]
    M6 = Mk[2]
    M7 = Mk[3]

    try: # Get parameters
        b = parameters['b'].value
        x3 = parameters['x3'].value
        x5 = parameters['x5'].value
        x6 = parameters['x6'].value
        psi3 = parameters['psi3'].value
        psi5 = parameters['psi5'].value
        psi6 = parameters['psi6'].value
        psi7 = parameters['psi7'].value
    except KeyError:
        b,x3,x5,x6,psi3,psi5,psi6,psi7 = parameters

    if (M3 < psi3):
        S3 = (1 - (M3 / psi3))
    else:
        S3 = 0
    if (M5 < psi5):
        S5 = (1 - (M5 / psi5))
    else:
        S5 = 0
    if (M6 < psi6):
        S6 = (1 - (M6 / psi6))
    else:
        S6 = 0
    if (M7 < psi7):
        S7 = (1 - (M7 / psi7))
    else:
        S7 = 0

    dM3dt = b * M3 * S3 + x3 * S3 * M7
    dM5dt = b * M5 * S5 + x5 * S5 * M7
    dM6dt = b * M6 * S6 + x6 * S6 * M7 
    dM7dt = b * M7 * S7

    return [dM3dt,dM5dt,dM6dt,dM7dt]
#==============================================================================
''' Compute negative log likelihood of Tromas' data given the model,see eq. (3) pg. 11 '''
def negLogLike(parameters):
    # Solve ODE system to get model values; parameters are not yet fitted
    Lk0 = [0,parameters['I0'].value]
    MM = odeint(model,Lk0,ZERO_DAYS_AXIS,args=(parameters,))

    nll = 0
    epsilon = 10**-10
    for t in range(4):          # Iterate through days
        for p in range(5):      # Iterate through replicates
            for k in range(4):  # Iterate through leaves
                Vktp = TROMAS_DATA[20 * t + 4 * p + k][4]          # Number of infected cells
                Aktp = TROMAS_DATA[20 * t + 4 * p + k][3] + Vktp   # Total number of cells observed
                Iktp = MM[t + 1][k]                                # Frequency of cellular infection

                if (Iktp <= 0):
                    Iktp = epsilon
                elif (Iktp >= 1):
                    Iktp = 1 - epsilon

                nll += Vktp * np.log(Iktp) + (Aktp - Vktp) * np.log(1 - Iktp)
    
    return [-nll]
#==============================================================================
''' Miminize negative log likelihood with differential evolution algorithm '''
result_likelihood = minimize(negLogLike,likeParams,method = 'differential_evolution')
report_fit(result_likelihood)

解决方法

我知道了。我没有在likeParams中定义“ x3”。我猜想该错误会传播并吐出为差异演化错误。