Victorian 2020-06-30 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]

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 [0.5, 0.3, 4.0]

p = [0.5, 0.3, 4.0]

for counter = 1:32

    global p   # use vector p defined outside the for loop

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

    X = zeros(len,3) # a lenx2 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 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) + 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 lenx2 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 20x1 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 lenx3 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 lenx1 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: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 20x1 column matrix

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

)