Victorian 2020-07-29 local cases

"$VERSION"

using Plots, Dates

Plots.default(

    size = (800,600),

    xtickfontsize = 10,

    ytickfontsize = 10,

    legend = :topleft

)

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,

427,216,363,275,374,484,403,300,357,459,532,384,295]

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,06,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 [8.3, 0.09, -13.9]

p = [8.3, 0.09, -13.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

        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="Vic",markershape=:auto)

plot!(day,

title = "Victorian Daily Local Cases",

[expo_cases,linear_cases,quadratic_cases,simple_expo_cases],

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

xlabel="Days since 2020-06-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="Vic",

    markershape=:circle)

plot!(day,

title = "Victorian 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-06-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 =  [97.0, 0.0881, -241.0]

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 = "Victorian Cumulative Local Cases",

markershape=:auto,

label="Vic",

xlabel="Days since 2020-06-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 = "Victorian Cumulative Local Cases",

markershape=:auto,

label="Vic",

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

ylabel="Cumulative local cases LOGSCALE")

plot!(day,

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

label="Expo"

)

import Statistics:mean

function movingaverage(X::Vector{T},numofele::Int) where T <: Real

    BackDelta = div(numofele,2)

    ForwardDelta = isodd(numofele) ? div(numofele,2) : div(numofele,2) - 1

    len = length(X)

    Y = similar(X)

    for n = 1:len

        lo = max(1,n - BackDelta)

        hi = min(len,n + ForwardDelta)

        Y[n] = mean(X[lo:hi])

    end

    return Y

end

function geometricmovingaverage(X::Vector{T},numofele::Int) where T <: Real

    XX = similar(X)

    for k in eachindex(X)

        XX[k] = log(X[k])

    end

    BackDelta = div(numofele,2)

    ForwardDelta = isodd(numofele) ? div(numofele,2) : div(numofele,2) - 1

    len = length(XX)

    Y = similar(XX)

    for n = 1:len

        lo = max(1,n - BackDelta)

        hi = min(len,n + ForwardDelta)

        Y[n] = exp(  mean(XX[lo:hi])  )

    end

    return Y

end

fcumcases = Float64.(cumulativecases)

growth = similar(fcumcases); growth[1]=1.0

for k = 2:length(fcumcases)

    growth[k]=1.0 + (fcumcases[k]-fcumcases[k-1])/fcumcases[k-1]

end

avggrowth = geometricmovingaverage(growth,7)

doubling = similar(fcumcases); doubling[1]=1.0

for k = 2:length(avggrowth)

    doubling[k]=log(2.0)/log(avggrowth[k])

    if doubling[k]===Inf || doubling[k] <= 1.0

        doubling[k]=1.0

    end

end

rounded_doubling = round.(doubling,digits=1);

rounded_growth = round.(growth,digits=2);

# Perform a Plot of doubling time

plot(day,doubling,

legend=:none,

linewidth=2,

title="Doubling Time for Vic Cumulative Cases",

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

ylabel="Doubling Time in Days",

)

plot(day,avggrowth,

legend=:none,

linewidth=2,

ylims=(0.9,1.6),

title="Average growth for Vic Cumulative Cases",

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

ylabel="Average daily growth",

)