Each individual infects $c$ new people each day.
$I_{n+1} = I_n + (c I_n) = I_n(1 + c)$
So
$I_{n+1} = \lambda I_n$, where $\lambda := 1 + c$ -- growth rate
and $I_0 = 1$
I0 = 1
1
I_0 = 1
1
I₀ = 1 # \_0 <TAB>
1
c = 0.01 # average no. of people that each individual infects on each day
0.01
λ = 1 + c # \lambda <TAB>
1.01
λ = 1 + c;
;
suppresses output
I_1 = λ * I_0
1.01
I_2 = λ * I_1
1.0201
I_3 = λ * I_2
1.030301
T = 10 # final time
I = zeros(T)
10-element Array{Float64,1}: 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
I
10-element Array{Float64,1}: 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
I[1]
0.0
I[1:2]
2-element Array{Float64,1}: 0.0 0.0
I = zeros(Int64, T)
10-element Array{Int64,1}: 0 0 0 0 0 0 0 0 0 0
I = zeros(T)
10-element Array{Float64,1}: 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
I[1] = I_0 # Could use OffsetArrays.jl -- enables arbitrary indexing
1
for n in 1:T-1
I[n+1] = λ * I[n]
@show n, I[n]
end
# for loops do not return anything so running a for loop does not output
(n, I[n]) = (1, 1.0) (n, I[n]) = (2, 1.01) (n, I[n]) = (3, 1.0201) (n, I[n]) = (4, 1.030301) (n, I[n]) = (5, 1.04060401) (n, I[n]) = (6, 1.0510100501) (n, I[n]) = (7, 1.061520150601) (n, I[n]) = (8, 1.0721353521070098) (n, I[n]) = (9, 1.08285670562808)
I
10-element Array{Float64,1}: 1.0 1.01 1.0201 1.030301 1.04060401 1.0510100501 1.061520150601 1.0721353521070098 1.08285670562808 1.0936852726843609
using Plots
┌ Info: Precompiling Plots [91a5bcdd-55d7-5caf-9e0b-520d859cae80] └ @ Base loading.jl:1260
plot(I, m=:o, label="I[n]", legend=:topleft)
T = 20
I = zeros(I)
I[1] = I_0
for n in 1:T-1
I[n+1] = λ * I[n]
@show n, I[n]
end
plot(I)
(n, I[n]) = (1, 1.0) (n, I[n]) = (2, 1.01) (n, I[n]) = (3, 1.0201) (n, I[n]) = (4, 1.030301) (n, I[n]) = (5, 1.04060401) (n, I[n]) = (6, 1.0510100501) (n, I[n]) = (7, 1.061520150601) (n, I[n]) = (8, 1.0721353521070098) (n, I[n]) = (9, 1.08285670562808)
BoundsError: attempt to access 10-element Array{Float64,1} at index [11] Stacktrace: [1] setindex!(::Array{Float64,1}, ::Float64, ::Int64) at ./array.jl:825 [2] top-level scope at ./In[30]:6
function run_infection(T=20) # default
I = zeros(T)
I[1] = I_0
for n in 1:T-1
I[n+1] = λ * I[n]
end
return I
end
run_infection (generic function with 2 methods)
methods(run_infection)
run_infection(10)
10-element Array{Float64,1}: 1.0 1.01 1.0201 1.030301 1.04060401 1.0510100501 1.061520150601 1.0721353521070098 1.08285670562808 1.0936852726843609
I_result = run_infection(10);
I_result
10-element Array{Float64,1}: 1.0 1.01 1.0201 1.030301 1.04060401 1.0510100501 1.061520150601 1.0721353521070098 1.08285670562808 1.0936852726843609
Always separate data generation from plotting!
I_result = run_infection(20)
plot(I_result, m=:o)
using Interact
Unable to load WebIO. Please make sure WebIO works for your Jupyter client. For troubleshooting, please see the WebIO/IJulia documentation.
end_T = 1000
@manipulate for T in slider(1:end_T, value=1)
I_result = run_infection(T)
plot(I_result, m=:o)
xlims!(0, end_T)
ylims!(0, 10)
end
plot(run_infection(1000))
plot(run_infection(1000), yscale=:log10)
I_result = run_infection(1000)
plot(log10.(I_result)) # log of each element in I_result
ylabel!("log(I_n)")
$\log(I_n) = a n + b$
Take exponentials of both sides:
$I_n = \exp(an + b) = C e^{a n}$
Solve $I_{n+1} = \lambda I_n$:
$I_n = \lambda^n I_0$
$\log(\lambda^n) = n \log(\lambda)$
Exp growth is unrealistic: Assumes that there are always more people to infect -- wrong (finite population).
$I_{n+1} = I_n + c I_n$
Each individual will be in contact with a fraction $\alpha$ of the population. At each contact there will a probability $p$ that you infect each person.
Number of contacts = $\alpha N$
So $c = p \alpha N$
Original model:
$I_{n+1} = I_n + (p \alpha N) I_n$ -- good approximation when almost everybody still susceptible
New model:
Can only infect uninfected people!
$I_{n+1} = I_n + (p \alpha S_n) I_n$
$I_{n+1} = I_n + [p \alpha (N - I_n)] I_n = f(I_n)$
$I_{n+1} = I_n + \beta(I_n, S_n) I_n = f(I_n)$
p = 0.02
α = 0.01
N = 1000
β(I, S) = p * α * S
function run_infection(T=20) # default
I = zeros(T)
I[1] = I_0
for n in 1:T-1
I[n+1] = I[n] + β(I[n], N - I[n]) * I[n]
end
return I
end
run_infection (generic function with 2 methods)
I = run_infection(20)
20-element Array{Float64,1}: 1.0 1.1998 1.439472095992 1.7269520992073721 2.071746046338255 2.4852368292698026 2.9810489147042514 3.57548136711873 4.288020827141153 5.141947568046585 6.165049156697394 7.390457421815974 8.857625133998434 10.613458656195233 12.713621286504926 15.224018310562577 18.222467825971023 21.80054972443152 26.065606875660343 31.14284507843312
plot(I, m=:o)
plot(I, m=:o, yscale=:log10)
I = run_infection(100)
plot(I, m=:o, yscale=:log10)
p = 0.01
α = 0.1
N = 100
I = run_infection(100)
plot(I, m=:o)
Sigmoid shape -- S-shaped
plot(I, m=:o, yscale=:log10)
"Logistic-type growth"
$I_{n+1} = \lambda I_n$
Idea: Instead of modelling population globally, model each individual
"Patch model": Local patches where population is well mixed
Network model: Links between nodes
Exponential growth:
$I_{n+1} = \lambda_n I_n$ -- growth rate changes in time
rand()
0.2767437763441829
rand()
0.7089008345855983
rand()
0.3789568897206148
www.random.org
randn()
-0.33523291258515997
p = 0.02
α = 0.01
N = 1000
function run_infection(c_average=1.1, T=20) # default
I = zeros(T)
I[1] = I_0
for n in 1:T-1
c = c_average + 0.1*randn()
I[n+1] = I[n] + c * I[n]
end
return I
end
run_infection (generic function with 3 methods)
c_average = 1.1
cs = [c_average + 0.1*randn() for i in 1:100]
100-element Array{Float64,1}: 1.1794683814069813 1.0510617699472489 1.0799811965453343 1.3160413008055867 1.0655215426073297 1.014913328568307 1.2142141455660913 0.9952511944053644 1.1118228511102424 1.2274911767386665 1.097194489799497 1.156157408743076 1.0138054185809025 ⋮ 1.252209922391299 1.0305288124096037 1.1816089517166688 0.9848866930013738 1.175520029029598 1.1695318311722314 1.1237453859981914 1.1925550201637143 0.9331362892813453 1.0848564135882912 1.184898469695791 1.0756350756110384
scatter(cs)