# Markov chains¶

In [19]:
p = 0.1

M = [1-p  p
0   1]

Out[19]:
2×2 Array{Float64,2}:
0.9  0.1
0.0  1.0
In [2]:
M = [1-p  p; 0 1]

Out[2]:
2×2 Array{Float64,2}:
0.9  0.1
0.0  1.0
In [4]:
P0 = [1.0  0.0]   # no comma: row vector

Out[4]:
1×2 Array{Float64,2}:
1.0  0.0
In [5]:
P1 = P0 * M

Out[5]:
1×2 Array{Float64,2}:
0.9  0.1
In [6]:
P2 = P1 * M

Out[6]:
1×2 Array{Float64,2}:
0.81  0.19
In [7]:
P3 = P2 * M

Out[7]:
1×2 Array{Float64,2}:
0.729  0.271
In [12]:
function dynamics(P0, M, n)
P = copy(P0)

for i in 1:n
new_P = P * M

# P = copy(new_P)

P, new_P = new_P, P   # swaps "pointers" P and new_P
end

return P
end

Out[12]:
dynamics (generic function with 1 method)
In [11]:
# Alternatives:

# P *= M   # equiv to P = P * M  -- first does P * M, which creates a new vector,
# and then makes P point to that vector -- bad in terms of memory usage

# P .*= M  #

In [13]:
dynamics(P0, M, 1)

Out[13]:
1×2 Array{Float64,2}:
0.9  0.1
In [14]:
dynamics(P0, M, 2)

Out[14]:
1×2 Array{Float64,2}:
0.81  0.19
In [15]:
dynamics(P0, M, 10)

Out[15]:
1×2 Array{Float64,2}:
0.348678  0.651322
In [16]:
dynamics(P0, M, 100)

Out[16]:
1×2 Array{Float64,2}:
2.65614e-5  0.999973
In [17]:
0.9^100

Out[17]:
2.6561398887587544e-5
$$P_0^{(n)} = (1-p)^{n-1}$$

Suppose small chance of re-infection after recovery:

In [20]:
M2 = [1-p  p
0.01   0.99]

Out[20]:
2×2 Array{Float64,2}:
0.9   0.1
0.01  0.99
In [21]:
for n in 1:10
@show dynamics(P0, M2, n)
end

dynamics(P0, M2, n) = [0.9 0.1]
dynamics(P0, M2, n) = [0.811 0.189]
dynamics(P0, M2, n) = [0.73179 0.26821]
dynamics(P0, M2, n) = [0.6612931000000001 0.33870690000000003]
dynamics(P0, M2, n) = [0.598550859 0.40144914100000006]
dynamics(P0, M2, n) = [0.54271026451 0.45728973549000007]
dynamics(P0, M2, n) = [0.49301213541390004 0.5069878645861001]
dynamics(P0, M2, n) = [0.44878080051837105 0.5512191994816291]
dynamics(P0, M2, n) = [0.40941491246135026 0.5905850875386499]
dynamics(P0, M2, n) = [0.37437927209060173 0.6256207279093984]

In [29]:
data2 = [dynamics(P0, M2, n) for n in 1:1000]

Out[29]:
1000-element Array{Array{Float64,2},1}:
[0.9 0.1]
[0.811 0.189]
[0.73179 0.26821]
[0.6612931000000001 0.33870690000000003]
[0.598550859 0.40144914100000006]
[0.54271026451 0.45728973549000007]
[0.49301213541390004 0.5069878645861001]
[0.44878080051837105 0.5512191994816291]
[0.40941491246135026 0.5905850875386499]
[0.37437927209060173 0.6256207279093984]
[0.3431975521606356 0.6568024478393646]
[0.3154458214229657 0.6845541785770345]
[0.2907467810664395 0.7092532189335607]
⋮
[0.09090909090909075 0.9090909090909072]
[0.09090909090909075 0.9090909090909072]
[0.09090909090909075 0.9090909090909072]
[0.09090909090909075 0.9090909090909072]
[0.09090909090909075 0.9090909090909072]
[0.09090909090909075 0.9090909090909072]
[0.09090909090909075 0.9090909090909072]
[0.09090909090909075 0.9090909090909072]
[0.09090909090909075 0.9090909090909072]
[0.09090909090909075 0.9090909090909072]
[0.09090909090909075 0.9090909090909072]
[0.09090909090909075 0.9090909090909072]
In [30]:
data1 = [dynamics(P0, M, n) for n in 1:1000]

Out[30]:
1000-element Array{Array{Float64,2},1}:
[0.9 0.1]
[0.81 0.19]
[0.7290000000000001 0.271]
[0.6561000000000001 0.34390000000000004]
[0.5904900000000002 0.40951000000000004]
[0.5314410000000002 0.46855900000000006]
[0.47829690000000014 0.5217031000000001]
[0.43046721000000016 0.5695327900000001]
[0.38742048900000015 0.6125795110000002]
[0.34867844010000015 0.6513215599000002]
[0.31381059609000017 0.6861894039100002]
[0.28242953648100017 0.7175704635190002]
[0.25418658283290013 0.7458134171671003]
⋮
[5.569828659391132e-46 1.0]
[5.012845793452018e-46 1.0]
[4.511561214106817e-46 1.0]
[4.060405092696135e-46 1.0]
[3.654364583426521e-46 1.0]
[3.2889281250838693e-46 1.0]
[2.9600353125754822e-46 1.0]
[2.664031781317934e-46 1.0]
[2.397628603186141e-46 1.0]
[2.157865742867527e-46 1.0]
[1.942079168580774e-46 1.0]
[1.7478712517226966e-46 1.0]
In [34]:
data1 = [dynamics(big.(P0), big.(M), n) for n in 1:1000]

Out[34]:
1000-element Array{Array{BigFloat,2},1}:
[0.90000000000000002220446049250313080847263336181640625 0.1000000000000000055511151231257827021181583404541015625]
[0.8100000000000000399680288865056359482888058144019096323303533017413935457540219 0.1900000000000000127675647831893003381312806238275281893325883254353483864385055]
[0.7290000000000000539568389967826091957912766296712997674362064523351291651732287 0.2710000000000000212607709215717481435429990544451144271503457758014564845125599]
[0.6561000000000000647482067961391318336711984918800389746448261819435354487933681 0.3439000000000000307032177460087049524868890987062116998820265147454492010642484]
[0.5904900000000000728417326456565242114419731566738477862317211245261222302138236 0.4095100000000000408201250579054445261384825763896687746093535231340378275185354]
[0.5314410000000000786690712573090471188041558507812798257323740177853629880026913 0.4685590000000000513821762915256407794092748945007577577634818870146132957424792]
[0.4782969000000000826025248201745004937135297269725610041997998920008010514653191 0.5217031000000000621991735894056350129871378448962380305968851077916431851084736]
[0.4304672100000000849625969578937729844736299766426462354471702259111979680073881 0.5695327900000000731145072263572656973913541231853294639337026777693425456771542]
[0.3874204890000000860246294198674462079631533156113804064877557442170927926453398 0.6125795110000000840003399615874056263229387050078437861247453355847159978375624]
[0.3486784401000000860246294198674472691467562798720989536662792732586003248024988 0.6513215599000000947534186390708366676142336105244459737303456858722294799575158]
[0.3138105960900000851643831256687738470270556516915066053765721364513341694115703 0.686189403910000105291435743004599220727652988463310255866203810827737981710096]
[0.2824295364810000836159397961111608085515712665585679347789497881048156686678768 0.7175704635190001155498728013237926919871638486415139138888697407975251277006205]
[0.2541865828329000815255413012083827940214825141245363272683324710651904864375911 0.7458134171671001254792656521119932894235880688528270038585841350507207907265375]
⋮
[5.569828659391124188757232448600998115875120406720079361515822720750214970941043e-46 1.000000000000000277555756156289196735666137413709077605670722188951490677920293]
[5.012845793452011893556549621202792525466626795416519806644992344549027708261368e-46 1.000000000000000277555756156289196735666137413764775892264633433930939012842856]
[4.511561214106810815508431034798220818118629560567011591909004464629444717913867e-46 1.000000000000000277555756156289196735666137413814904350199153555649192918447771]
[4.060405092696129834134370669462537998509359477167786999054726023304548474423564e-46 1.000000000000000277555756156289196735666137413860019962340221666308696797249348]
[3.654364583426516940880038066846011759012171690033884547005518375419375884064593e-46 1.000000000000000277555756156289196735666137413900624013267182966904018115552218]
[3.288928125083865327935228278058167389363600983227166710680775219851709870291419e-46 1.000000000000000277555756156289196735666137413937167659101448138341398346668092]
[2.960035312575478868170580066359433577795468506786371843796054741641770108270761e-46 1.00000000000000027755575615628919673566613741397005694035228679344647249485135]
[2.664031781317931047079509214219866476214087741115930712745515340264835273823332e-46 1.000000000000000277555756156289196735666137413999657293478041583771327974377361]
[2.397628603186138001524946731844619918581113546294769882089849115550226050624717e-46 1.000000000000000277555756156289196735666137414026297611291220895720957777495733]
[2.157865742867524254610501653802225321181669905531464529211891267177208635426397e-46 1.000000000000000277555756156289196735666137414050273897323082277067158484692734]
[1.942079168580771877063696124049864626198472790712206069857473076485666316656767e-46 1.000000000000000277555756156289196735666137414071852554751757520811119617121449]
[1.747871251722694732480146683709954880910050439730411075058778923064513102511265e-46 1.000000000000000277555756156289196735666137414091273346437565240659827082663578]

Stationary state

Unchanging probability vector -- eigenvector $v$ of the transition matrix:

$$M v = \lambda v$$
In [35]:
v = data2[end]

Out[35]:
1×2 Array{Float64,2}:
0.0909091  0.909091
In [38]:
v * M2 == v

Out[38]:
true
In [40]:
v * M2 == 1.0 .*  v  # eigenvector with eigenvalue 1

Out[40]:
true

dynamics multiplies repeatedly by the matrix

In [41]:
3^100

Out[41]:
-2984622845537545263
In [42]:
v * M^10

Out[42]:
1×2 Array{Float64,2}:
0.031698  0.968302
In [44]:
v * M2^100

Out[44]:
1×2 Array{Float64,2}:
0.0909091  0.909091

## Continuous limit¶

Back to the recovery model. We have talking about the number of people that recover each day

$A \rightarrow B$

Chemical reactions: Continuous time: $t \in \mathbb{R}$

Probability of recovering at time $n$ is $p_n := p (1-p)^{n-1}$

Cumulative Probability of having recovered at some time <= $n$:

$s_N := \sum_{n=1}^N p_n$

In [46]:
p = 0.1

probs = [p * (1 - p)^(n-1) for n in 1:20]

Out[46]:
20-element Array{Float64,1}:
0.1
0.09000000000000001
0.08100000000000002
0.0729
0.06561
0.05904900000000001
0.05314410000000001
0.04782969000000001
0.04304672100000001
0.03874204890000001
0.03486784401000001
0.031381059609000006
0.028242953648100012
0.02541865828329001
0.02287679245496101
0.02058911320946491
0.018530201888518418
0.016677181699666577
0.015009463529699918
0.013508517176729929
In [47]:
s = cumsum(probs)

Out[47]:
20-element Array{Float64,1}:
0.1
0.19
0.271
0.3439000000000001
0.40951000000000004
0.46855900000000006
0.5217031000000001
0.5695327900000001
0.6125795110000002
0.6513215599000002
0.6861894039100002
0.7175704635190002
0.7458134171671003
0.7712320754503903
0.7941088679053513
0.8146979811148162
0.8332281830033346
0.8499053647030012
0.864914828232701
0.878423345409431
In [48]:
using Plots

In [50]:
plot(s, m=:o)

Out[50]:
In [58]:
ys = [0; reduce(vcat, [ [s[n], s[n]] for n in 1:20 ])]

pop!(ys)
pushfirst!(ys, 0)

Out[58]:
41-element Array{Float64,1}:
0.0
0.0
0.1
0.1
0.19
0.19
0.271
0.271
0.3439000000000001
0.3439000000000001
0.40951000000000004
0.40951000000000004
0.46855900000000006
⋮
0.7712320754503903
0.7941088679053513
0.7941088679053513
0.8146979811148162
0.8146979811148162
0.8332281830033346
0.8332281830033346
0.8499053647030012
0.8499053647030012
0.864914828232701
0.864914828232701
0.878423345409431
In [59]:
xs = [0; reduce(vcat, [ [n, n] for n in 1:20 ])];

In [68]:
plot(xs, ys, leg=false)
scatter!(s)

ylabel!("cumulative probability of recovery up to time n")
xlabel!("number of days, n")

Out[68]:
In [70]:
plot(s, m=:o)

Out[70]:

Think of our discrete-time process in steps, not of days, but of a time $\delta$

Instead of "cumulative probability that have decayed by day $n$", want to talk about "cumulative probability that have decayed by time $n \delta$

What is the probability $p(\delta)$ of recovering in time $\delta$

In [ ]: