NSW 2020-08-04 local cases

"$VERSION"

using Plots, Dates

Plots.default(

    size = (800,600),

    xtickfontsize = 10,

    ytickfontsize = 10,

    legend = :topleft

)

cases = [0,1,1,2,1,1,3,8,11,10,6,6,11,13,15,12,15,15,7,9,11,9,12,17,16,18,

13,11,9,11]

len = length(cases)

println("len is ",len)

day = [ k for k in 0:(len-1)];

@show day;

println("t = $(len-1) for the date $(Date(2020,07,06)+Day(len-1))")

# Use Mathematica to find the partial derivative of "a Exp[b t] + c"

# D[a Exp[b t] + c, a] == Exp(b t)

# D[a Exp[b t] + c, b] == a Exp(b t) t

# D[a Exp[b t] + c, c] == 1

# Initiliase the parameter vector to our guess of [46023.192, 0.0, -46021.7458]

p = [46023.192, 0.0, -46021.7458]

for counter = 1:0

    global p   # use vector p defined outside the for loop

    R = zeros(len,1) # a len*1 column matrix

    X = zeros(len,3) # a len*2 matrix

    for k = 1:len

        X[k,1] = exp(p[2] * day[k])

        X[k,2] = p[1] * exp(p[2] * day[k]) * day[k]

        X[k,3] = 1

        R[k] = cases[k] - (p[1] * exp( p[2] * day[k]) + p[3])

    end

    # use vec() function to turn a 2*1 column matrix into a 1D vector

    delta = vec( inv(X'*X) * X' * R )

    p = p + delta

end

expo_param = copy(p)  # make a copy of p

expo_param_rounded = round.(expo_param,digits=4)

expo_cases = [ expo_param[1] * exp(expo_param[2] * t) + expo_param[3]    for t in day ];

println("exponential model is $(expo_param_rounded[1]) * exp($(expo_param_rounded[2]) * t) + $(expo_param_rounded[3])")

# Use Mathematica to find the partial derivative of "a Exp[b t] + c"

# D[a Exp[b t] , a] == Exp(b t)

# D[a Exp[b t] , b] == a Exp(b t) t

# Initiliase the parameter vector to our guess of [10.6, 0.08]

p = [10.6, 0.08]

for counter = 1:32

    global p   # use vector p defined outside the for loop

    R = zeros(len,1) # a len*1 column matrix

    X = zeros(len,2) # a len*2 matrix

    for k = 1:len

        X[k,1] = exp(p[2] * day[k])

        X[k,2] = p[1] * exp(p[2] * day[k]) * day[k]

        R[k] = cases[k] - p[1] * exp( p[2] * day[k]) 

    end

    # use vec() function to turn a 2*1 column matrix into a 1D vector

    delta = vec( inv(X'*X) * X' * R )

    p = p + delta

end

simple_expo_param = copy(p)  # make a copy of p

simple_expo_param_rounded = round.(simple_expo_param,digits=4)

simple_expo_cases = [ simple_expo_param[1] * exp(simple_expo_param[2] * t)    for t in day ];

println("simple exponential model is $(simple_expo_param_rounded[1]) * exp($(simple_expo_param_rounded[2]) * t)")

X = zeros(len,2) # a len*2 matrix

for k = 1:len

    X[k,1] = 1.0

    X[k,2] = day[k]

end

# Dont forget to reshape cases from a 1D vector into a len*1 column matrix

p = inv(X'*X) * X' * reshape(cases,len,1)

linear_param = copy(p)

linear_param_rounded = round.(linear_param,digits=4)

linear_cases = [ linear_param[1] + linear_param[2] * t    for t in day ];

println("linear model is $(linear_param_rounded[1]) + $(linear_param_rounded[2]) * t")

X = zeros(len,3) # a len*3 matrix

for k = 1:len

    X[k,1] = 1.0

    X[k,2] = day[k]

    X[k,3] = day[k]^2.0

end

# Dont forget to reshape cases from a 1D vector into a len*1 column matrix

p = inv(X'*X) * X' * reshape(cases,len,1)

quadratic_param = copy(p)  # make a copy of p

quadratic_param_rounded = round.(quadratic_param,digits=4)

quadratic_cases = [ quadratic_param[1] + quadratic_param[2]*t + quadratic_param[3]*t^2    for t in day ];

println("quadratic model is $(quadratic_param_rounded[1]) + $(quadratic_param_rounded[2]) * t + $(quadratic_param_rounded[3]) * t^2")

#

# Now Plot the cases with various models

#

plot(day,cases,linewidth=2,

    label="NSW",markershape=:auto)

plot!(day,

title = "NSW Daily Local Cases",

[expo_cases,linear_cases,quadratic_cases,simple_expo_cases],

label=["Expo" "Linear" "Quadratic" "Simple Expo"],

xlabel="Days since 2020-07-06",

ylabel="Num of daily local cases"

)

function getLogTicks(vector,count=1)

    function poststr(expnum)

        n = Int64(expnum) ÷ 3

        if n == 0

            return ""

        elseif n == 1

            return "k"

        else

            return "×\$10^{$(n*3)}\$"

        end

    end

    min = ceil(log10(minimum(vector)))

    max = floor(log10(maximum(vector)))

    major = [ [k*n for k = 1:count] for n in 10 .^ collect(min:max) ]

    major = hcat(major...)

    major = vec(major)

    major = filter(x->x<=maximum(vector),major)

    #majorText = ["\$10^{$(round(Int64,i))}\$" for i=min:max]

    majorText = ["$(Int64(j*10^(i%3)))$(poststr(i))" for i=min:max for j=1:count]

    majorText = majorText[1:length(major)]

    minor = [j*10^i for i=(min-1):(max+1) for j=(count+1):9]

    minor = minor[findall(minimum(vector) .<= minor .<= maximum(vector))]

    ([major; minor], [majorText; fill("", length(minor))])

end

#

# Now Plot the cases with various models

# using LOGSCALE for the Y-Axis

# We can't take log(0) since it is Minus Infinity

# So we assume any value less than one is one

yticks_scale = getLogTicks(

 [ k<1.0 ? 1.0 : k for k in cases ]

 );

plot(day,

    [ k<1.0 ? 1.0 : k for k in cases ],

    linewidth=2,

    yaxis=:log10, yticks=yticks_scale,

    label="NSW",

    markershape=:circle)

plot!(day,

title = "NSW Daily Local Cases",

[ [ k<1.0 ? 1.0 : k for k in expo_cases ],

[ k<1.0 ? 1.0 : k for k in linear_cases ],

[ k<1.0 ? 1.0 : k for k in quadratic_cases ],

[ k<1.0 ? 1.0 : k for k in simple_expo_cases ] ],

label=["Expo" "Linear" "Quadratic" "Simple Expo"],

xlabel="Days since 2020-07-06",

ylabel="Num of daily local cases LOGSCALE"

)

# Calc RSS for various model

expo_residue = sum([ (expo_cases[k] - cases[k])^2 for k in 1:len ])

println("exponential model residue is ",round(expo_residue,digits=1))

simple_expo_residue = sum([ (simple_expo_cases[k] - cases[k])^2 for k in 1:len ])

println("simple exponential model residue is ",round(simple_expo_residue,digits=1))

linear_residue = sum([ (linear_cases[k] - cases[k])^2 for k in 1:len ])

println("linear model residue is ",round(linear_residue,digits=1))

quadratic_residue = sum([ (quadratic_cases[k] - cases[k])^2 for k in 1:len ])

println("quadratic model residue is ",round(quadratic_residue,digits=1))

cumulativecases = zeros(Int64,len)

cumulativecases[1] = cases[1]

for k = 2:len

    cumulativecases[k] = cumulativecases[k-1] + cases[k]

end

p = [13.3, 0.14, -16.9]

for counter = 1:32

    global p   # use vector p defined outside the for loop

    R = zeros(len,1) # a len*1 column matrix

    X = zeros(len,3) # a len*2 matrix

    for k = 1:len

        X[k,1] = exp(p[2] * day[k])

        X[k,2] = p[1] * exp(p[2] * day[k]) * day[k]

        X[k,3] = 1.0

        R[k] = cumulativecases[k] - (p[1] * exp( p[2] * day[k]) + p[3])

    end

    # use vec() function to turn a 2*1 column matrix into a 1D vector

    delta = vec( inv(X'*X) * X' * R )

    p = p + delta

end

expo2_param = copy(p)  # make a copy of p

expo2_param_rounded = round.(expo2_param,digits=4)

expo2_cases = [ expo2_param[1] * exp(expo2_param[2] * t) + expo2_param[3]    for t in day ];

println("exponential2 model is $(expo2_param_rounded[1]) * exp($(expo2_param_rounded[2]) * t) + $(expo2_param_rounded[3])")

plot(day,cumulativecases,

linewidth=2,

title = "NSW Cumulative Local Cases",

markershape=:auto,

label="NSW",

xlabel="Days since 2020-07-06",

ylabel="Cumulative local cases")

plot!(day,

expo2_cases,

label="Expo"

)

yticks_scale = getLogTicks(

 [ k<1.0 ? 1.0 : k for k in cumulativecases ]

 );

plot(day,

[ k<1.0 ? 1.0 : k for k in cumulativecases ],

linewidth=2,

yaxis=:log, yticks=yticks_scale,

title = "NSW Cumulative Local Cases",

markershape=:auto,

label="NSW",

xlabel="Days since 2020-07-06",

ylabel="Cumulative local cases LOGSCALE")

plot!(day,

[ k<1.0 ? 1.0 : k for k in expo2_cases ],

label="Expo"

)

Vic_cases = [0,1,1,0,4,6,2,2,6,11,7,6,12,12,24,15,12,16,19,23,25,40,49,

74,64,70,77,65,108,73,127,191,134,165,288,216,273,176,317,427,216,363,

275,374,484,403,300,357,459,532,384]

Vic_len = length(Vic_cases)

Vic_cumulativecases = zeros(Int64,Vic_len)

Vic_cumulativecases[1] = Vic_cases[1]

for k = 2:Vic_len

    Vic_cumulativecases[k] = Vic_cumulativecases[k-1] + Vic_cases[k]

end

ext=6

xmax=len+ext

ymax=maximum(Vic_cumulativecases[1:(1+xmax)])

color1=palette(:default)[1]

color2=palette(:default)[2]

plot(0:xmax,

Vic_cumulativecases[1:(1+xmax)],

xticks=0:(xmax÷4):xmax,

title = "Current NSW vs early Vic Cumulative Local Cases",

markershape=:square,

markercolor=color2,

linecolor=color2,

label="Vic",

xlabel="Days since Day Zero",

ylabel="Cumulative local cases",

annotation=

  [

    (1/14*xmax,62/90*ymax,text("NSW day zero is 2020 Jul 06",:left)),

    (1/14*xmax,55/90*ymax,text("Vic day zero is 2020 Jun 06",:left))

  ]

)

plot!(day,cumulativecases,

markershape=:circle,

markercolor=color1,

linecolor=color1,

linewidth=2, 

label="NSW"

)

ymax=maximum(Vic_cases[1:(1+xmax)])

plot(0:(xmax),

Vic_cases[1:(1+xmax)],

xticks=0:(xmax÷4):xmax,

title = "Current NSW vs early Vic Daily Local Cases",

markershape=:square,

markercolor=color2,

linecolor=color2,

legend=:topleft,label="Vic",

xlabel="Days since Day Zero",

ylabel="Daily local cases",

annotation=[(1/14*xmax,62/90*ymax,text("NSW day zero is 2020 Jul 06",:left)),

(1/14*(xmax),55/90*ymax,text("Vic day zero is 2020 Jun 06",:left))

]

)

plot!(day,cases,

markershape=:circle,

markercolor=color1,

linecolor=color1,

linewidth=2,

label="NSW"

)

Vic_cases = [0,1,1,0,4,6,2,2,6,11,7,6,12,12,24,15,12,16,19,23,25,40,49,

74,64,70,77,65,108,73,127,191,134,165,288,216,273,176,270,238,317,

295,723,627]

Vic_len = length(Vic_cases)

Vic_cumulativecases = zeros(Int64,Vic_len)

Vic_cumulativecases[1] = Vic_cases[1]

for k = 2:Vic_len

    Vic_cumulativecases[k] = Vic_cumulativecases[k-1] + Vic_cases[k]

end

ext=6

xmax=len+ext

ymax=maximum(Vic_cumulativecases[1:(1+xmax)])

color1=palette(:default)[1]

color2=palette(:default)[2]

plot(0:(xmax),

Vic_cumulativecases[1:(1+xmax)],

xticks=0:(xmax÷4):xmax,

title = "Current NSW vs early Vic Cumulative Local Cases",

markershape=:square,

markercolor=color2,

linecolor=color2,

legend=:topleft,label="Vic",

xlabel="Days since Day Zero",

ylabel="Cumulative local cases",

annotation=[(1/14*xmax,62/90*ymax,text("NSW day zero is 2020 Jul 06",:left)),

(1/14*(xmax),55/90*ymax,text("Vic day zero is 2020 Jun 06",:left))

]

) 

plot!(day,cumulativecases,

markershape=:circle,

markercolor=color1,

linecolor=color1,

linewidth=2,

label="NSW"

) 

ymax=maximum(Vic_cases[1:(1+xmax)])

plot(0:(xmax),

Vic_cases[1:(1+xmax)],

xticks=0:(xmax÷4):xmax,

title = "Current NSW vs early Vic Daily Local Cases",

markershape=:square,

markercolor=color2,

linecolor=color2,

legend=:topleft,label="Vic",

xlabel="Days since Day Zero",

ylabel="Daily local cases",

annotation=[(1/14*xmax,62/90*ymax,text("NSW day zero is 2020 Jul 06",:left)),

(1/14*(xmax),55/90*ymax,text("Vic day zero is 2020 Jun 06",:left))

]

)

plot!(day,cases,

markershape=:circle,

markercolor=color1,

linecolor=color1,

linewidth=2,

label="NSW"

)

yticks_scale = getLogTicks(

 [ k<1.0 ? 1.0 : k for k in Vic_cumulativecases ] , 2

 );

xmax=length(Vic_cumulativecases)-1

ymax=maximum(Vic_cumulativecases)

plot(0:(length(Vic_cumulativecases)-1),

[ k<1.0 ? 1.0 : k for k in Vic_cumulativecases ],

yaxis=:log, yticks=yticks_scale,

title = "Current NSW vs Vic Cumulative Local Cases",

markershape=:square,

markercolor=color2,

linecolor=color2,

legend=:topleft,label="Vic",

xlabel="Days since Day Zero",

ylabel="Cumulative local cases LOGSCALE",

annotation=[(4/14*xmax,10.0^(22/90*log(10,ymax)),text("NSW day zero is 2020 Jul 06",:left)),

(4/14*(xmax),10.0^(15/90*log(10,ymax)),text("Vic day zero is 2020 Jun 06",:left))

]) 

plot!(day,

[ k<1.0 ? 1.0 : k for k in cumulativecases ],

markershape=:circle,

markercolor=color1,

linecolor=color1,

linewidth=2,

label="NSW"

)

yticks_scale = getLogTicks(

 [ k<1.0 ? 1.0 : k for k in Vic_cases ] , 2

 );

xmax=length(Vic_cases)-1

ymax=maximum(Vic_cases)

plot(0:(length(Vic_cases)-1),

[ k<1.0 ? 1.0 : k for k in Vic_cases ],

yaxis=:log, yticks=yticks_scale,

title = "Current NSW vs Vic Daily Local Cases",

markershape=:square,

markercolor=color2,

linecolor=color2,

legend=:topleft,label="Vic",

xlabel="Days since Day Zero",

ylabel="Daily local cases LOGSCALE",

annotation=[(6/14*xmax,10.0^(22/90*log(10,ymax)),text("NSW day zero is 2020 Jul 06",:left)),

(6/14*(xmax),10.0^(15/90*log(10,ymax)),text("Vic day zero is 2020 Jun 06",:left))

])

plot!(day,

[ k<1.0 ? 1.0 : k for k in cases ],

markershape=:circle,

markercolor=color1,

linecolor=color1,

linewidth=2,

label="NSW"

)