ODE的拟合初始条件：使用lmfit最小化和基于似然的方法

Question

我用一阶微分方程来模拟病毒的传播，用似然法来拟合实验数据。但是，每当运行代码时，每次运行代码时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左右

第一次：

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

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

Answer 1

你的模式对我来说有点难学。我想知道为什么

I0 = [likeParams['I0'].value]

然后在negLogLike函数中使用它，而实际上不要在拟合模型中使用I0的可变值。也许应该是这样

 M3 = odeint(model, parameters['I0'].value, ZERO_DAYS_AXIS, args=(parameters,))

或者是其他的？这似乎给出了一个更令人满意的结果。我用method='Nelder'试了一下：

Original I0 = .00372, My I0   =  0.0010018942980410726
Original beta = .871, My beta =  0.7074140684730617
Original psi3 = .083, My psi3 =  0.08807904093119709