我正在使用一阶微分方程对病毒传播进行建模，并正在使用似然方法将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,2,7378,102],8978,66],4,8105,116],5,11941,70],[5,12349,214],8958,154,],10661,156],11924,305],9848,417],[7,8854,429],11098,886],12757,470],10233,594],8347,721],[10,8895,586],9554,999],11021,940],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,max = 1.0000)
#====================================================================
''' Initial conditions '''
I0 = [likeParams['I0'].value]
#====================================================================
''' Time axis '''
ZERO_DAYS_AXIS = [0,7,10]
DAYS_AXIS = [3,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左右。

第一次：

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

第二次：

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

第三次：

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

如果您能帮助我学习如何正确地适应初始条件，我将不胜感激。 谢谢！