0
$\begingroup$

I'm trying to fit a simple logistic growth model to dummy data using Python's Scipy package. The code is shown below, along with the output that I get. The correct output is shown below it. I'm not quite sure what's going wrong here.

import scipy.optimize as optim
from scipy.integrate import odeint
import numpy as np
import pandas as pd

N0 = 0.37
parsic = [.25, 12.9]

df_yeast = pd.DataFrame({'cd': [9.6, 18.3, 29., 47.2, 71.1, 119.1, 174.6, 257.3, 350.7, 441., 513.3, 559.7, 594.8, 629.4, 640.8, 651.1, 655.9, 659.6], 'td': np.arange(18)})

def logistic_de(t, N, r, K):
    return r*N*(1 - N/K)

def logistic_solution(t, r, K):
    return odeint(logistic_de, N0, t, (r, K), tfirst=True).ravel()

params, _ = optim.curve_fit(logistic_solution, df_yeast['td'], df_yeast['cd'], p0=parsic)

N1 = odeint(logistic_de, N0, np.linspace(0, 20, 10000), (params[0], params[1]), tfirst=True)

plt.plot(np.linspace(0, 20, 10000), N1)
plt.scatter(df_yeast['td'], df_yeast['cd'])
plt.ylabel('num yeast')
plt.xlabel('time')

My output: enter image description here

Correct output: enter image description here

$\endgroup$
2
  • $\begingroup$ Add the initial value N0 to the parameters that get optimized. Compare the first data point 9.6 to your initial value of 0.37. $\endgroup$ Sep 23 '21 at 6:30
  • $\begingroup$ That works quite well and results in the parameters r=0.5476140280399281, K=662.6552616132678, N0=9.10156146739931. $\endgroup$ Sep 23 '21 at 6:50
2
$\begingroup$

To repeat the answers you got at your cross-post https://stackoverflow.com/questions/69292456/logistic-growth-curve-using-scipy-is-not-quite-right

The logistic curve is rather rigid in its symmetry. Forcing the initial value to be some prescribed constant effectively removes a degree of freedom in the family of curves that are available to the fitting process.

Thus include N0 in the set of parameters, do not forget to unpack it for the computation for the plot, and you will get a fitted solution that looks like your second graph with parameters r=0.5476140280399281, K=662.6552616132678, N0=9.10156146739931

Changes in code were

parsic = [.25, 12.9, N0]
...
def logistic_solution(t, r, K, N0):
    return odeint(logistic_de, N0, t, (r, K), tfirst=True, atol=1e-6, rtol=1e-9).ravel()
...
t = np.linspace(0, 20, 1000)
N1 = logistic_solution(t, *params)
plt.plot(t, N1)
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.