Estimating parameters with scipy minimize with unorthodox observed data

Question

I wish to minimize the parameters beta and gamma in this model. However, my observed data isnt in the form of a time series. The values I want to estimate are for when two certain trajectories have equilibrium values. Namely, when equilibrium values for I (prevalence) and J_diff (incidence) reach 0.4 and 0.3 respectively. My code is as follows:

def peak_infections(x):

    # Total population, N.
    N = 1
    # Initial number of infected and recovered individuals, I0 and R0.
    I0, R0 = 0.001, 0
    # Everyone else, S0, is susceptible to infection initially.
    beta = x[0]
    gamma = x[1]
    U0 = N - I0 - R0
    J0 = I0
    Lf0, Ls0 = 0, 0
    # Contact rate, beta, and mean recovery rate, gamma, (in 1/years).
    beta, gamma = 15, 2/5
    mu, muTB, sigma, rho = 1/80, 1/6, 1/6, 0.03
    u, v, w = 0.083, 0.88, 0.0006

    # A grid of time points 
    times = np.arange(0,20,2.5)

    def deriv(y, times, N, beta, gamma, mu, muTB, sigma, rho, u, v, w):
        U, Lf, Ls, I, R, cInc = y
        b = (mu * (U + Lf + Ls + R)) + (muTB * I)
        lamda = beta * I
        clamda = 0.2 * lamda
        dU = b - ((lamda + mu) * U)
        dLf = (lamda*U) + ((clamda)*(Ls + R)) - ((u + v + mu) * Lf)
        dLs = (u * Lf) - ((w + clamda + mu) * Ls)
        dI = w*Ls + v*Lf - ((gamma + muTB + sigma) * I) + (rho * R)
        dR = ((gamma + sigma) * I) - ((rho + clamda + mu) * R)
        cI = w*Ls + v*Lf + (rho * R)
        return dU, dLf, dLs, dI, dR, cI

    # Initial conditions are S0, I0, R0
    # Integrate the SIR equations over the time grid, t.
    solve = odeint(deriv, (U0, Lf0, Ls0, I0, R0, J0), times, args=(N, beta, gamma, mu, muTB, sigma, rho, u, v, w))
    U, Lf, Ls, I, R, cInc = solve.T

    return I

def residual(x):

    # Total population,  N.
    StartingPop = 1
    prev= 0.4/StartingPop
    return np.sum((peak_infections(x) - prev) ** 2)


x0 = [12, 0.4] #estimates for beta and gamma starting point
res = minimize(residual, x0, method="Nelder-Mead", options={'fatol':1e-04}).x
print(res)

However, where I attempt the minimizing as res, it simply returns the initial estimates in x0 that I gave it. How do I correct this code to include in the residual function, that this must be optimised for when I and J_diff reach their equilibrium states for 0.4 and 0.3?

Answer 1

You are overwriting your input arguments to the function 'peak_infections'. beta and gamma are being assigned the values of x[0] and x[1], respectively. But a few lines later, they are being reassigned as 15 and 2/5. No matter what you pass to the function, the result is the same. Just delete the line where you assign those values to 15 and 2/5 and you will get a result.