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]:
5 10 15 20 0.2 0.4 0.6 0.8 y1
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]:
0 5 10 15 20 0.0 0.2 0.4 0.6 0.8 number of days, n cumulative probability of recovery up to time n
In [70]:
plot(s, m=:o)
Out[70]:
5 10 15 20 0.2 0.4 0.6 0.8 y1

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 [ ]: