Scipy differential evolution used with likelihood method

245 views Asked by At

In my work, I'm trying to fit a model to data using a likelihood based method. I've had this code work before for slightly different models but for some reason it has been throwing this error when used with this model: "RuntimeError: The map-like callable must be of the form f(func, iterable), returning a sequence of numbers the same length as 'iterable'". I'm not sure if it's the model or something wrong with the code, but if you could help me understand what this error message means and how to fix it I would appreciate it.

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)
1

There are 1 answers

0
Innocuous Rift On

I figured it out. I didn't define "x3" in likeParams. I guess that error propagated and got spat out as the differential evolution error.