Подгонка начального условия ODE с использованием lmfit Minimum и подхода, основанного на правдоподобии

Я использую дифференциальные уравнения первого порядка для моделирования распространения вирусов и использую вероятностный подход, чтобы подогнать ODE к экспериментальным данным. Однако всякий раз, когда код запускается, подогнанное начальное значение ODE изменяется каждый раз при запуске кода, в то время как все другие параметры в модели каждый раз сходятся к одному и тому же значению. Мой код показан ниже.

import numpy as np
from scipy.integrate import odeint
from lmfit import minimize, Parameters, Parameter, report_fit
import matplotlib.pyplot as plt
#====================================================================
''' Data on which the model will be fitted '''
TROMAS_DATA = [[3,3,1,7219,33],[3,3,2,7378,102],[3,3,3,8978,66],[3,3,4,8105,116],[3,3,5,11941,70],
               [5,3,1,12349,214],[5,3,2,8958,154,],[5,3,3,10661,156],[5,3,4,11924,305],[5,3,5,9848,417],
               [7,3,1,8854,429],[7,3,2,11098,886],[7,3,3,12757,470],[7,3,4,10233,594],[7,3,5,8347,721],
               [10,3,1,8347,721],[10,3,2,8895,586],[10,3,3,9554,999],[10,3,4,11021,940],[10,3,5,9416,917]]
#====================================================================
''' Parameter list '''
likeParams = Parameters()
likeParams.add('I0', value = .00372, min = .0000001, max = 1.0000)
likeParams.add('b', value = .871, min = .0001, max = 1.0000)
likeParams.add('psi3', value = .083, min = .0001, max = 1.0000)
#====================================================================
''' Initial conditions '''
I0 = [likeParams['I0'].value]
#====================================================================
''' Time axis '''
ZERO_DAYS_AXIS = [0, 3, 5, 7, 10]
DAYS_AXIS = [3, 5, 7, 10]
t = np.linspace(0, 10, 10000)
#====================================================================
''' Define the model '''
def model(M3, t, parameters):

    try: # Get parameters
        b = parameters['b'].value
        psi3 = parameters['psi3'].value
    except KeyError:
        b, psi3 = parameters

    if (M3 < psi3):
        S3 = (1 - (M3 / psi3))
    else:
        S3 = 0

    dM3dt = b * M3 * S3

    return dM3dt
#====================================================================
''' Compute negative log likelihood of Tromas' data given the model '''
def negLogLike(parameters):
    # Solve ODE system to get model values; parameters are not yet fitted
    M3 = odeint(model, I0, 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
            Vktp = TROMAS_DATA[5 * t + p][4]          # Number of infected cells
            Aktp = TROMAS_DATA[5 * t + p][3] + Vktp   # Total number of cells observed
            Iktp = M3[t + 1]                          # 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')
#====================================================================
''' Compare previous and fitted values '''
print("Original I0 = .00372, My I0   = ", result_likelihood.params['I0'].value)
print("Original beta = .871, My beta = ", result_likelihood.params['b'].value)
print("Original psi3 = .083, My psi3 = ", result_likelihood.params['psi3'].value)

Чтобы помочь проиллюстрировать проблему, результат трех прогонов кода показан ниже. Я также знаю, что «правильное» значение для I0 должно быть примерно 0,003.

1-й раз:

Original I0 = .00372, My I0   =  0.9711066426171252
Original beta = .871, My beta =  0.44374676021094367
Original psi3 = .083, My psi3 =  0.11330842894454747

2-й раз:

Original I0 = .00372, My I0   =  0.09183577426999114
Original beta = .871, My beta =  0.4437477735310379
Original psi3 = .083, My psi3 =  0.11330747932246728

3-й раз:

Original I0 = .00372, My I0   =  0.47857244584676956
Original beta = .871, My beta =  0.44374801498547956
Original psi3 = .083, My psi3 =  0.11330824660617532

Если бы вы могли помочь мне научиться правильно соответствовать начальному состоянию, я был бы очень признателен. Спасибо!