<jupyter_start><jupyter_text><jupyter_code>pip install tf-utils
import warnings
warnings.filterwarnings("ignore")
from tf_utils import *

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
dtf = pd.read_csv("https://raw.githubusercontent.com/CSSEGISandData/COVID-19/master/csse_covid_19_data/csse_covid_19_time_series/time_series_covid19_confirmed_global.csv", sep=",")
dtf.head()<jupyter_output><empty_output><jupyter_text>Curve fitting for **India**<jupyter_code>## groupby country
dtf = dtf.drop(['Province/State','Lat','Long'], axis=1).groupby("Country/Region").sum().T

## convert index to datetime
dtf.index = pd.to_datetime(dtf.index, infer_datetime_format=True)

## create total cases column
dtf = pd.DataFrame(index=dtf.index, data=dtf["India"].values, columns=["total"])

## create daily changes column
dtf["new"] = dtf["total"] - dtf["total"].shift(1)
dtf["new"] = dtf["new"].fillna(method='bfill')

dtf.tail(10)
fig, ax = plt.subplots(nrows=2, ncols=1, sharex=True, figsize=(13,7))
ax[0].scatter(dtf.index, dtf["total"].values, color="black")
ax[0].set(title="total cases")
ax[1].bar(dtf.index, dtf["new"].values, color="black")
ax[1].set(title="new cases")
plt.show()

'''
Linear function: f(x) = a + b*x
'''
def f(x):
    return 10 + 2500*x

y_linear = f(x=np.arange(len(dtf)))
'''
Exponential function: f(x) = a + b^x
'''
def f(x):
    return 10 + 1.091**x

y_exponential = f(x=np.arange(len(dtf)))
'''
Logistic function: f(x) = a / (1 + e^(-b*(x-c)))
'''
def f(x): 
    return 350000 / (1 + np.exp(-0.5*(x-20)))

y_logistic = f(x=np.arange(len(dtf)))
fig, ax = plt.subplots(figsize=(13,5))
ax.scatter(dtf["total"].index, dtf["total"].values, color="black")
ax.plot(dtf["total"].index, y_linear, label="linear", color="red")
ax.plot(dtf["total"].index, y_exponential, label="exponential", color="green")
ax.plot(dtf["total"].index, y_logistic, label="logistic", color="blue")
ax.legend()
plt.show()<jupyter_output><empty_output><jupyter_text>-> Data starts low and grows exponentially ... looks similar to a exponential function. But I know for fact that   this phenomenon has an upper limit, because the virus can't infect more than the total population of the country.   Therefore the logistic function is more appropriate to model this data.<jupyter_code>'''
Gaussian function: f(x) = a * e^(-0.5 * ((x-μ)/σ)**2)
'''
def f(x):
    return 12000 * np.exp(-0.5 * ((x-145)/21)**2)

y_gaussian = f(x=np.arange(len(dtf)))
fig, ax = plt.subplots(figsize=(13,5))
ax.bar(dtf["new"].index, dtf["new"].values, color="black")
ax.plot(dtf["new"].index, y_gaussian, color="red")
plt.show()
from scipy.optimize import curve_fit
'''
Logistic function: f(x) = capacity / (1 + e^-k*(x - midpoint) )
'''
def logistic_f(X, c, k, m):
    y = c / (1 + np.exp(-k*(X-m)))
    return y
## optimize from scipy
X=np.arange(len(dtf["total"]))
y=dtf["total"].values
logistic_model, cov = curve_fit(logistic_f,
                                X, 
                               y , 
                                maxfev=10000,
                                p0=[np.max(dtf["total"]), 1, 1])
## print the parameters
logistic_model
'''
Gaussian function: f(x) = a * e^(-0.5 * ((x-μ)/σ)**2)
'''
def gaussian_f(X, a, b, c):
    y_prime = a * np.exp(-0.5 * ((X-b)/c)**2)
    return y_prime
## optimize from scipy
X1=np.arange(len(dtf["new"]))
y1=dtf["new"].values
gaussian_model, cov = curve_fit(gaussian_f,
                                X1, 
                                y1, 
                                maxfev=10000,
                                p0=[1,np.mean(dtf["new"]),1])
## print the parameters
gaussian_model
'''
Plot parametric fitting fuction.
'''
def utils_plot_parametric(dtf, zoom=30, figsize=(15,5)):
    ## interval
    dtf["residuals"] = dtf["ts"] - dtf["model"]
    dtf["conf_int_low"] = dtf["forecast"] - 1.96*dtf["residuals"].std()
    dtf["conf_int_up"] = dtf["forecast"] + 1.96*dtf["residuals"].std()
    fig, ax = plt.subplots(nrows=1, ncols=2, figsize=figsize)
    
    ## entire series
    dtf["ts"].plot(marker=".", linestyle='None', ax=ax[0], title="Parametric Fitting", color="black")
    dtf["model"].plot(ax=ax[0], color="green")
    dtf["forecast"].plot(ax=ax[0], grid=True, color="red")
    ax[0].fill_between(x=dtf.index, y1=dtf['conf_int_low'], y2=dtf['conf_int_up'], color='b', alpha=0.3)
   
    ## focus on last
    first_idx = dtf[pd.notnull(dtf["forecast"])].index[0]
    first_loc = dtf.index.tolist().index(first_idx)
    zoom_idx = dtf.index[first_loc-zoom]
    dtf.loc[zoom_idx:]["ts"].plot(marker=".", linestyle='None', ax=ax[1], color="black", 
                                  title="Zoom on the last "+str(zoom)+" observations")
    dtf.loc[zoom_idx:]["model"].plot(ax=ax[1], color="green")
    dtf.loc[zoom_idx:]["forecast"].plot(ax=ax[1], grid=True, color="red")
    ax[1].fill_between(x=dtf.loc[zoom_idx:].index, y1=dtf.loc[zoom_idx:]['conf_int_low'], 
                       y2=dtf.loc[zoom_idx:]['conf_int_up'], color='b', alpha=0.3)
    plt.show()
    
    return dtf[["ts","model","residuals","conf_int_low","forecast","conf_int_up"]]
'''
Forecast unknown future.
:parameter
    :param ts: pandas series
    :param f: function
    :param model: list of optim params
    :param pred_ahead: number of observations to forecast (ex. pred_ahead=30)
    :param freq: None or str - 'B' business day, 'D' daily, 'W' weekly, 'M' monthly, 'A' annual, 'Q' quarterly
    :param zoom: for plotting
'''
def forecast_curve(ts, f, model, pred_ahead=None, freq="D", zoom=30, figsize=(15,5)):
    ## fit
    X = np.arange(len(ts))
    fitted = f(X, model[0], model[1], model[2])
    dtf = ts.to_frame(name="ts")
    dtf["model"] = fitted
    
    ## index
    index = pd.date_range(start='2020-06-17',periods=pred_ahead,freq=freq)
    index = index[1:]
    ## forecast
    Xnew = np.arange(len(ts)+1, len(ts)+1+len(index))
    preds = f(Xnew, model[0], model[1], model[2])
    dtf = dtf.append(pd.DataFrame(data=preds, index=index, columns=["forecast"]))
    
    ## plot
    utils_plot_parametric(dtf, zoom=zoom)
    return dtf
preds = forecast_curve(dtf["total"], logistic_f, logistic_model, 
                       pred_ahead=30, freq="D", zoom=7)
## Forecast
preds = forecast_curve(dtf["new"], gaussian_f, gaussian_model, pred_ahead=30,freq="D", zoom=7, figsize=(15,5))<jupyter_output><empty_output><jupyter_text>For some reason we are not getting good fit with gaussian function, I guess it's just the matter of parameter scaling. Trying with lmfit, if we can get the results.<jupyter_code>!pip install lmfit
from lmfit import Model
import math
def gaussian(x, amp, cen, wid):
    """1-d gaussian: gaussian(x, amp, cen, wid)"""
    return (amp / (math.sqrt(2*math.pi) * wid)) * np.exp(-(x-cen)**2 / (2*wid**2))


gmodel = Model(gaussian)
result = gmodel.fit(y1, x=X1, amp=1, cen=np.mean(dtf["new"]), wid=1)

print(result.fit_report())

plt.plot(X1, y1, 'bo')
plt.plot(X1, result.init_fit, 'k--', label='initial fit')
plt.plot(X1, result.best_fit, 'r-', label='best fit')
plt.legend(loc='best')
plt.show()
# <end examples/doc_model_gaussian.py><jupyter_output>[[Model]]
    Model(gaussian)
[[Fit Statistics]]
    # fitting method   = leastsq
    # function evals   = 4
    # data points      = 149
    # variables        = 3
    chi-square         = 2.9719e+09
    reduced chi-square = 20355151.4
    Akaike info crit   = 2510.46723
    Bayesian info crit = 2519.47907
##  Warning: uncertainties could not be estimated:
    amp:  at initial value
    cen:  at initial value
    wid:  at initial value
[[Variables]]
    amp:  1.00000000 (init = 1)
    cen:  2553.90604 (init = 2553.906)
    wid:  1.00000000 (init = 1)
{'fraction_non_alphanumeric': 0.1835830618892508, 'fraction_numerical': 0.032052117263843646, 'mean_word_length': 3.8065122103944895, 'pattern_counts': {'":': 0, '<': 17, '<?xml version=': 0, '>': 18, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 2, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '4', 'dir': 'jupyter-structured-clean-dedup', 'id': '247727', 'license': 'MIT', 'n_tokens_mistral': 3054, 'n_tokens_neox': 2867, 'n_words': 595, 'path': 'Time_Series_Forecasting_with_Parametric_Curve_Fitting.ipynb', 'repo_name': 'deepacefic/Covid19-Dashboard', 'size': '7675'}