This example is used for fitting world covid-19 cases number
import numpy as np
import pandas as pd
from datetime import datetime
from lmfit import Minimizer, Parameters, report_fit
import chart_studio.plotly as py
import cufflinks as cf
def Cauchy_cumulative_hazard_fit(x,loc,scale,decaybase):
decayterm=np.power(decaybase,(x-loc))
decayterm[np.where((x-loc)<0)]=1.0
z=(x-loc)/(scale*decayterm)
CHF=-np.log(0.5 - np.arctan(z)/np.pi)
return CHF
# define objective function: returns the array to be minimized
def Cauchy_cumulative_hazard_residual(params, x, data):
"""Model Cauch cumulative hazard function"""
amp = params['amp']
loc = params['loc']
scale = params['scale']
decaybase=params['decaybase']
decayterm=np.power(decaybase,(x-loc))
decayterm[np.where((x-loc)<0)]=1.0
z=(x-loc)/(scale*decayterm)
CHF=-np.log(0.5 - np.arctan(z)/np.pi)
model = amp * CHF
return model - data
data00=pd.read_csv('worldTot.csv')
data00=data00.rename(columns={"Unnamed: 0":"date"})
dataWorld=data00.set_index('date')
start_date=dataWorld.index[0]
end_date=dataWorld.index[-1]
dateData=pd.date_range(start=start_date,end=end_date)
dataWorld['Date']=dateData
dataWorld=dataWorld.set_index('Date')
data=dataWorld['cases']
forecastDays=60
dateForecast= pd.date_range(start=end_date,periods=forecastDays+1)[1:]
dateObsForecast=dateData.append(dateForecast)
# normalize case data
last=data[-1]
data=data/last
# set x values interval = 1day
dataLen=data.count()
x = np.linspace(1, dataLen, dataLen)
# create a set of Parameters
params = Parameters()
params.add('amp', value=1)
params.add('scale', value=1, min=0.1)
params.add('loc', value=100, min=0)
params.add('decaybase',value=1,min=0.8,max=1.05) # shparly increase value>1
# do fit, here with the default leastsq algorithm
minner = Minimizer(Cauchy_cumulative_hazard_residual, params, fcn_args=(x, data))
result = minner.minimize()
# calculate final result
forecastdays=forecastDays
final = data + result.residual
x1=np.linspace(1,dataLen+forecastdays,dataLen+forecastdays)
scale=result.params['scale'].value
amp=result.params['amp'].value
loc=result.params['loc'].value
decaybase=result.params['decaybase'].value
y1=amp*Cauchy_cumulative_hazard_fit(x1,loc,scale,decaybase) # forecast
nn=np.int(np.ceil(loc))
y1[:nn]=final.iloc[:nn] # replace data before
# write error report
report_fit(result)
output = pd.DataFrame({'date' : [],'Forecast':[],'Cases': [],'Fitting':[],'Increase':[]})
output
output['date']=dateObsForecast
output['Forecast']=y1*last
output['Cases'].iloc[:dataLen]=data.values*last
output['Fitting'].iloc[:dataLen]=final.values*last
output['Increase'].iloc[1:]=(y1[1:]-y1[:-1])*last
output=output.set_index('date')
@interact
def plot_ProjectedWorldCOVID19():
fig=output.iplot(asFigure=True,
mode='lines+markers',
size=3,secondary_y = 'Increase',
secondary_y_title='Increase',
xTitle='Date',
yTitle='Cases',
title='Projected COVID-19 Cases over the World',
theme='solar',
filename='COVID19-World',
sharing='public')
fig.show()
url=py.iplot(fig, filename='COVID19-World')
print(url)
Like this:
Like Loading...