Victorian 2020-07-12 local cases

"$VERSION"

using Plots

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]

len = length(cases)

println("len is ",len)

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

@show day;

# 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 [6.5, 0.3, -8.0]

p = [6.5, 0.3, -8.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

        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])")

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,

    legend=:topleft,label="Vic",markershape=:auto)

plot!(day,

title = "Victorian Daily Local Cases",

[expo_cases,linear_cases,quadratic_cases],

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

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

ylabel="Num of daily local cases"

)

function getLogTicks(vec)

    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(vec)))

    max = floor(log10(maximum(vec)))

    major = 10 .^ collect(min:max)

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

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

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

    ([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 ],

    yaxis=:log10, yticks=yticks_scale,

    legend=:topleft,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 ] ],

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

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))

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 = [17.0, 0.1, -22.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,

title = "Victorian Cumulative Local Cases",

markershape=:auto,

legend=:topleft,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 ],

yaxis=:log, yticks=yticks_scale,

title = "Victorian Cumulative Local Cases",

markershape=:auto,

legend=:topleft,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"

)