Modeling Logistic Growth

In [5]:
import math
import numpy as np
import pandas as pd
import scipy.optimize as optim
import matplotlib.pyplot as plt

My Logistic Growth Curve

In [6]:
def my_logistic(t):
    return 1000 / (1 + 999 * math.exp(-2 * t))
In [7]:
x = np.linspace(0, 10, 11)
y = [my_logistic(i) for i in x]
plt.plot(x, y)
Out[7]:
[<matplotlib.lines.Line2D at 0xf7156f0>]
In [8]:
pd.DataFrame({'Time':x, 'Infections':y})
Out[8]:
Time Infections
0 0.0 1.000000
1 1.0 7.342147
2 2.0 51.820659
3 3.0 287.664369
4 4.0 748.992325
5 5.0 956.613256
6 6.0 993.899378
7 7.0 999.169992
8 8.0 999.887590
9 9.0 999.984785
10 10.0 999.997941

Comparison Exponential vs Logistic

In [9]:
def exponential(t):
    return 1*2**t

plt.plot(x, y)
plt.plot(x, exponential(x))
plt.title('Logistic Growth vs Exponential Growth')
plt.legend(['Logistic', 'Exponential'])
plt.xlabel('Time')
plt.ylabel('Infections')
plt.show()

Fit Logistic Curve On China Data

In [23]:
# Import the data
data = pd.read_csv('C://Users//jkorstan//Desktop//full_data_logistic.csv', sep=';')
data = data['total_cases']
data = data.reset_index(drop=False)
data.columns = ['Timestep', 'Total Cases']
data.head(10)
Out[23]:
Timestep Total Cases
0 0 278
1 1 310
2 2 574
3 3 835
4 4 1297
5 5 1985
6 6 2761
7 7 4537
8 8 5997
9 9 7736
In [11]:
# Define funcion with the coefficients to estimate
def my_logistic(t, a, b, c):
    return c / (1 + a * np.exp(-b*t))
In [12]:
# Randomly initialize the coefficients
p0 = np.random.exponential(size=3)
p0
Out[12]:
array([0.71591732, 1.6457033 , 1.20268737])
In [13]:
# Set min bound 0 on all coefficients, and set different max bounds for each coefficient
bounds = (0, [100000., 3., 1000000000.])
In [14]:
# Convert pd.Series to np.Array and use Scipy's curve fit to find the best Nonlinear Least Squares coefficients
x = np.array(data['Timestep']) + 1
y = np.array(data['Total Cases'])

(a,b,c),cov = optim.curve_fit(my_logistic, x, y, bounds=bounds, p0=p0)
In [15]:
# Show the coefficients
a,b,c
Out[15]:
(61.150558865123635, 0.19286193341255997, 81802.17465309602)
In [16]:
# Redefine the function with the new a, b and c
def my_logistic(t):
    return c / (1 + a * np.exp(-b*t))

Plot the Real Data vs the Fitted Model

In [24]:
plt.scatter(x, y)
plt.plot(x, my_logistic(x))
plt.title('Logistic Model vs Real Observations of China Coronavirus')
plt.legend([ 'Logistic Model', 'Real data'])
plt.xlabel('Time')
plt.ylabel('Infections')
Out[24]:
Text(0,0.5,'Infections')

Compute the moment of fastest growth

In [32]:
# The time step at which the growth is fastest
t_fastest = np.log(a) / b
t_fastest
Out[32]:
21.327894668265905
In [33]:
# First way to find the y of the fastest growth moment
y_fastest = c / 2
y_fastest
Out[33]:
40901.08732654801
In [34]:
# Second way to find the y of the fastest growth moment
my_logistic(t_fastest)
Out[34]:
40901.08732654801