Victorian 2020-06-23 local cases

"$VERSION"

using Plots

cases = [0,1,1,0,4,6,2,2,6,11,7,6,12,12,24,15,12,16]

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

# 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 [ 2.0 , 0.5 ]

p = [2.0, 0.5]

for counter = 1:32

    global p   # use vector p defined outside the for loop

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

    X = zeros(len,2) # a 18x2 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 2x1 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)    for t in day ];

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

X = zeros(len,2) # a 18x2 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 18x1 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 18x2 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 18x1 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"

)

# Calc RSS for various model

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

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

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

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

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

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