python - 与似然法一起使用的 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', value = .5, min = .0001, max = 20.0000)
likeParams.add('x6', value = .5, min = .0001, max = 20.0000)
likeParams.add('x7', value = .5, min = .0001, max = 20.0000)
likeParams.add('psi3', value = .5, min = .0001, max = 1.0000)
likeParams.add('psi5', value = .5, min = .0001, max = 1.0000)
likeParams.add('psi6', value = .5, min = .0001, max = 1.0000)
likeParams.add('psi7', value = .5, min = .0001, 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, 0, 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)

标签： pythonscipycurve-fittinglmfitlog-likelihood

python - 与似然法一起使用的 Scipy 差分进化

问题描述

解决方案

推荐阅读