Fit Taiwan Influenza Data By Ordinary Least Squares

FitInfluData.py

#
# desc : fet influenza data from Taiwan CDC Open Data
# auth : Jiankai Wang
# date : 2017/2/23
# plat : Python 2.7.12 | Anaconda 4.2.0 (64-bit)
# veri : 0.0.1
#


import sys
import json
import urllib2
import numpy as np
from numpy import pi
from mpl_toolkits.mplot3d import axes3d
import matplotlib.pyplot as plt
from matplotlib import cm
from scipy import optimize

np.random.seed(0)

###############################################
# desc : read the data and parse into json object
################################################
url = 'https://od.cdc.gov.tw/opendataplatform/?s=influlinechart&v=a1'
rawData = ""

try:   
    response = urllib2.urlopen(url)
    for line in response.readlines():
        rawData += line.strip()

except:
    # fetch data is failure
    sys.exit()

jsonData = json.loads(rawData)

###############################################
# desc : prepare data for multiple purpose drawing
###############################################

# 1 ... 52
allWeek = np.arange(1,53,1)
allClins = np.zeros(53)

fig, ax1 = \
    plt.subplots(1, 1, figsize=(8, 8), subplot_kw={'projection': '3d'})
  
crtYear = -1
colorList = {"his" : "blue", "fit" : "green", "ml" : "red"}

X = np.zeros((53 ,4)) # week
Y = np.zeros((53 ,4)) # year
Z = np.zeros((53 ,4)) # ratio

for item in jsonData:
    
    week = (int)(item['week'])
    year = (int)(item['year'])
    
    if week > 52 or year > 2016:
        continue        
    
    ratio = (float)(item['influop'])
    
    X[week][year-2013] = week
    Y[week][year-2013] = year-2013
    Z[week][year-2013] = ratio
    
    allClins[week] += ratio

# plot the last one year        
X = np.delete(X, 0, 0)
Y = np.delete(Y, 0, 0)
Z = np.delete(Z, 0, 0)

# remove the first one week[0]
allClins = allClins[1:]

def ave(x):
    return (float)(x) / (float)(4.0)

allClins = np.array(map(ave, allClins))

# fit data
def eventFun(pa1, pa2, pa3, pa4, pa5, week):
    usual = pa1*np.cos(week*2*pi/float(52))
    peak = pow(pa2 * week * np.e, (-1) * pa3 * week)
    peakW = week / pa4*1.1 - pa5 * week
    return(peak + peakW - usual)

# target function
fitfunc = lambda p, x: eventFun(p[0], p[1], p[2], p[3], p[4], x)
# Distance to the target function
errfunc = lambda p, x, y: fitfunc(p, x) - y
# Initial guess for the parameters
p0 = [3e-2, 2.5e-2, 9e-2, 2e11, -2e-2] 
# get parameters
p1, success = optimize.leastsq(errfunc, p0[:], args=(allWeek, allClins))

def fitLeastSqData(x):
    return fitfunc(p1, x)

fitData = np.array(map(fitLeastSqData, allWeek))

FX = X[:,0:2]
FY = np.zeros((52 ,2))
FY[:,0] = np.full((1,52), -1)
FY[:,1] = np.full((1,52), -2)
FZ = np.zeros((52 ,2))
FZ[:,0] = allClins
FZ[:,1] = fitData

###############################################
# desc : prepare data for drawing a 3d plot
###############################################

# Give the first plot only wireframes of the type y = c
ax1.plot_wireframe(X, Y, Z, rstride=0, cstride=1, color=colorList["his"])
ax1.plot_wireframe(FX, FY, FZ, rstride=0, cstride=1, color=colorList["fit"])

ax1.set_xlabel('X - Week')
ax1.set_ylabel('Y - Year')
ax1.set_zlabel('Z - Ratio')

ax1.set_title("Out-patient Clinicis for Influenza")

plt.tight_layout()
plt.show()