{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Markov chains"
]
},
{
"cell_type": "code",
"execution_count": 19,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"2×2 Array{Float64,2}:\n",
" 0.9 0.1\n",
" 0.0 1.0"
]
},
"execution_count": 19,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"p = 0.1\n",
"\n",
"M = [1-p p\n",
" 0 1]"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"2×2 Array{Float64,2}:\n",
" 0.9 0.1\n",
" 0.0 1.0"
]
},
"execution_count": 2,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"M = [1-p p; 0 1]"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"1×2 Array{Float64,2}:\n",
" 1.0 0.0"
]
},
"execution_count": 4,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"P0 = [1.0 0.0] # no comma: row vector"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"1×2 Array{Float64,2}:\n",
" 0.9 0.1"
]
},
"execution_count": 5,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"P1 = P0 * M"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"1×2 Array{Float64,2}:\n",
" 0.81 0.19"
]
},
"execution_count": 6,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"P2 = P1 * M"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"1×2 Array{Float64,2}:\n",
" 0.729 0.271"
]
},
"execution_count": 7,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"P3 = P2 * M"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"dynamics (generic function with 1 method)"
]
},
"execution_count": 12,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"function dynamics(P0, M, n)\n",
" P = copy(P0)\n",
" \n",
" for i in 1:n\n",
" new_P = P * M\n",
" \n",
" # P = copy(new_P)\n",
" \n",
" P, new_P = new_P, P # swaps \"pointers\" P and new_P\n",
" end\n",
" \n",
" return P\n",
"end"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [],
"source": [
"# Alternatives:\n",
"\n",
"# P *= M # equiv to P = P * M -- first does P * M, which creates a new vector, \n",
" # and then makes P point to that vector -- bad in terms of memory usage\n",
"\n",
"# P .*= M # "
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"1×2 Array{Float64,2}:\n",
" 0.9 0.1"
]
},
"execution_count": 13,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"dynamics(P0, M, 1)"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"1×2 Array{Float64,2}:\n",
" 0.81 0.19"
]
},
"execution_count": 14,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"dynamics(P0, M, 2)"
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"1×2 Array{Float64,2}:\n",
" 0.348678 0.651322"
]
},
"execution_count": 15,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"dynamics(P0, M, 10)"
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"1×2 Array{Float64,2}:\n",
" 2.65614e-5 0.999973"
]
},
"execution_count": 16,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"dynamics(P0, M, 100)"
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"2.6561398887587544e-5"
]
},
"execution_count": 17,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"0.9^100 "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$P_0^{(n)} = (1-p)^{n-1}$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Suppose small chance of re-infection after recovery:"
]
},
{
"cell_type": "code",
"execution_count": 20,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"2×2 Array{Float64,2}:\n",
" 0.9 0.1\n",
" 0.01 0.99"
]
},
"execution_count": 20,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"M2 = [1-p p\n",
" 0.01 0.99]"
]
},
{
"cell_type": "code",
"execution_count": 21,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"dynamics(P0, M2, n) = [0.9 0.1]\n",
"dynamics(P0, M2, n) = [0.811 0.189]\n",
"dynamics(P0, M2, n) = [0.73179 0.26821]\n",
"dynamics(P0, M2, n) = [0.6612931000000001 0.33870690000000003]\n",
"dynamics(P0, M2, n) = [0.598550859 0.40144914100000006]\n",
"dynamics(P0, M2, n) = [0.54271026451 0.45728973549000007]\n",
"dynamics(P0, M2, n) = [0.49301213541390004 0.5069878645861001]\n",
"dynamics(P0, M2, n) = [0.44878080051837105 0.5512191994816291]\n",
"dynamics(P0, M2, n) = [0.40941491246135026 0.5905850875386499]\n",
"dynamics(P0, M2, n) = [0.37437927209060173 0.6256207279093984]\n"
]
}
],
"source": [
"for n in 1:10\n",
" @show dynamics(P0, M2, n)\n",
"end"
]
},
{
"cell_type": "code",
"execution_count": 29,
"metadata": {
"scrolled": false
},
"outputs": [
{
"data": {
"text/plain": [
"1000-element Array{Array{Float64,2},1}:\n",
" [0.9 0.1]\n",
" [0.811 0.189]\n",
" [0.73179 0.26821]\n",
" [0.6612931000000001 0.33870690000000003]\n",
" [0.598550859 0.40144914100000006]\n",
" [0.54271026451 0.45728973549000007]\n",
" [0.49301213541390004 0.5069878645861001]\n",
" [0.44878080051837105 0.5512191994816291]\n",
" [0.40941491246135026 0.5905850875386499]\n",
" [0.37437927209060173 0.6256207279093984]\n",
" [0.3431975521606356 0.6568024478393646]\n",
" [0.3154458214229657 0.6845541785770345]\n",
" [0.2907467810664395 0.7092532189335607]\n",
" ⋮\n",
" [0.09090909090909075 0.9090909090909072]\n",
" [0.09090909090909075 0.9090909090909072]\n",
" [0.09090909090909075 0.9090909090909072]\n",
" [0.09090909090909075 0.9090909090909072]\n",
" [0.09090909090909075 0.9090909090909072]\n",
" [0.09090909090909075 0.9090909090909072]\n",
" [0.09090909090909075 0.9090909090909072]\n",
" [0.09090909090909075 0.9090909090909072]\n",
" [0.09090909090909075 0.9090909090909072]\n",
" [0.09090909090909075 0.9090909090909072]\n",
" [0.09090909090909075 0.9090909090909072]\n",
" [0.09090909090909075 0.9090909090909072]"
]
},
"execution_count": 29,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"data2 = [dynamics(P0, M2, n) for n in 1:1000]\n"
]
},
{
"cell_type": "code",
"execution_count": 30,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"1000-element Array{Array{Float64,2},1}:\n",
" [0.9 0.1]\n",
" [0.81 0.19]\n",
" [0.7290000000000001 0.271]\n",
" [0.6561000000000001 0.34390000000000004]\n",
" [0.5904900000000002 0.40951000000000004]\n",
" [0.5314410000000002 0.46855900000000006]\n",
" [0.47829690000000014 0.5217031000000001]\n",
" [0.43046721000000016 0.5695327900000001]\n",
" [0.38742048900000015 0.6125795110000002]\n",
" [0.34867844010000015 0.6513215599000002]\n",
" [0.31381059609000017 0.6861894039100002]\n",
" [0.28242953648100017 0.7175704635190002]\n",
" [0.25418658283290013 0.7458134171671003]\n",
" ⋮\n",
" [5.569828659391132e-46 1.0]\n",
" [5.012845793452018e-46 1.0]\n",
" [4.511561214106817e-46 1.0]\n",
" [4.060405092696135e-46 1.0]\n",
" [3.654364583426521e-46 1.0]\n",
" [3.2889281250838693e-46 1.0]\n",
" [2.9600353125754822e-46 1.0]\n",
" [2.664031781317934e-46 1.0]\n",
" [2.397628603186141e-46 1.0]\n",
" [2.157865742867527e-46 1.0]\n",
" [1.942079168580774e-46 1.0]\n",
" [1.7478712517226966e-46 1.0]"
]
},
"execution_count": 30,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"data1 = [dynamics(P0, M, n) for n in 1:1000]\n"
]
},
{
"cell_type": "code",
"execution_count": 34,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"1000-element Array{Array{BigFloat,2},1}:\n",
" [0.90000000000000002220446049250313080847263336181640625 0.1000000000000000055511151231257827021181583404541015625]\n",
" [0.8100000000000000399680288865056359482888058144019096323303533017413935457540219 0.1900000000000000127675647831893003381312806238275281893325883254353483864385055]\n",
" [0.7290000000000000539568389967826091957912766296712997674362064523351291651732287 0.2710000000000000212607709215717481435429990544451144271503457758014564845125599]\n",
" [0.6561000000000000647482067961391318336711984918800389746448261819435354487933681 0.3439000000000000307032177460087049524868890987062116998820265147454492010642484]\n",
" [0.5904900000000000728417326456565242114419731566738477862317211245261222302138236 0.4095100000000000408201250579054445261384825763896687746093535231340378275185354]\n",
" [0.5314410000000000786690712573090471188041558507812798257323740177853629880026913 0.4685590000000000513821762915256407794092748945007577577634818870146132957424792]\n",
" [0.4782969000000000826025248201745004937135297269725610041997998920008010514653191 0.5217031000000000621991735894056350129871378448962380305968851077916431851084736]\n",
" [0.4304672100000000849625969578937729844736299766426462354471702259111979680073881 0.5695327900000000731145072263572656973913541231853294639337026777693425456771542]\n",
" [0.3874204890000000860246294198674462079631533156113804064877557442170927926453398 0.6125795110000000840003399615874056263229387050078437861247453355847159978375624]\n",
" [0.3486784401000000860246294198674472691467562798720989536662792732586003248024988 0.6513215599000000947534186390708366676142336105244459737303456858722294799575158]\n",
" [0.3138105960900000851643831256687738470270556516915066053765721364513341694115703 0.686189403910000105291435743004599220727652988463310255866203810827737981710096]\n",
" [0.2824295364810000836159397961111608085515712665585679347789497881048156686678768 0.7175704635190001155498728013237926919871638486415139138888697407975251277006205]\n",
" [0.2541865828329000815255413012083827940214825141245363272683324710651904864375911 0.7458134171671001254792656521119932894235880688528270038585841350507207907265375]\n",
" ⋮\n",
" [5.569828659391124188757232448600998115875120406720079361515822720750214970941043e-46 1.000000000000000277555756156289196735666137413709077605670722188951490677920293]\n",
" [5.012845793452011893556549621202792525466626795416519806644992344549027708261368e-46 1.000000000000000277555756156289196735666137413764775892264633433930939012842856]\n",
" [4.511561214106810815508431034798220818118629560567011591909004464629444717913867e-46 1.000000000000000277555756156289196735666137413814904350199153555649192918447771]\n",
" [4.060405092696129834134370669462537998509359477167786999054726023304548474423564e-46 1.000000000000000277555756156289196735666137413860019962340221666308696797249348]\n",
" [3.654364583426516940880038066846011759012171690033884547005518375419375884064593e-46 1.000000000000000277555756156289196735666137413900624013267182966904018115552218]\n",
" [3.288928125083865327935228278058167389363600983227166710680775219851709870291419e-46 1.000000000000000277555756156289196735666137413937167659101448138341398346668092]\n",
" [2.960035312575478868170580066359433577795468506786371843796054741641770108270761e-46 1.00000000000000027755575615628919673566613741397005694035228679344647249485135]\n",
" [2.664031781317931047079509214219866476214087741115930712745515340264835273823332e-46 1.000000000000000277555756156289196735666137413999657293478041583771327974377361]\n",
" [2.397628603186138001524946731844619918581113546294769882089849115550226050624717e-46 1.000000000000000277555756156289196735666137414026297611291220895720957777495733]\n",
" [2.157865742867524254610501653802225321181669905531464529211891267177208635426397e-46 1.000000000000000277555756156289196735666137414050273897323082277067158484692734]\n",
" [1.942079168580771877063696124049864626198472790712206069857473076485666316656767e-46 1.000000000000000277555756156289196735666137414071852554751757520811119617121449]\n",
" [1.747871251722694732480146683709954880910050439730411075058778923064513102511265e-46 1.000000000000000277555756156289196735666137414091273346437565240659827082663578]"
]
},
"execution_count": 34,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"data1 = [dynamics(big.(P0), big.(M), n) for n in 1:1000]\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**Stationary state**"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Unchanging probability vector -- **eigenvector** $v$ of the transition matrix:\n",
"\n",
"$$ M v = \\lambda v$$"
]
},
{
"cell_type": "code",
"execution_count": 35,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"1×2 Array{Float64,2}:\n",
" 0.0909091 0.909091"
]
},
"execution_count": 35,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"v = data2[end]"
]
},
{
"cell_type": "code",
"execution_count": 38,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"true"
]
},
"execution_count": 38,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"v * M2 == v"
]
},
{
"cell_type": "code",
"execution_count": 40,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"true"
]
},
"execution_count": 40,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"v * M2 == 1.0 .* v # eigenvector with eigenvalue 1"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"`dynamics` multiplies repeatedly by the matrix"
]
},
{
"cell_type": "code",
"execution_count": 41,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"-2984622845537545263"
]
},
"execution_count": 41,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"3^100"
]
},
{
"cell_type": "code",
"execution_count": 42,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"1×2 Array{Float64,2}:\n",
" 0.031698 0.968302"
]
},
"execution_count": 42,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"v * M^10"
]
},
{
"cell_type": "code",
"execution_count": 44,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"1×2 Array{Float64,2}:\n",
" 0.0909091 0.909091"
]
},
"execution_count": 44,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"v * M2^100"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Continuous limit"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Back to the recovery model. We have talking about the number of people that recover each *day*"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$A \\rightarrow B$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Chemical reactions: **Continuous time**: $t \\in \\mathbb{R}$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Probability of recovering at time $n$ is $p_n := p (1-p)^{n-1}$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**Cumulative** Probability of having recovered at some time <= $n$:"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$s_N := \\sum_{n=1}^N p_n$"
]
},
{
"cell_type": "code",
"execution_count": 46,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"20-element Array{Float64,1}:\n",
" 0.1\n",
" 0.09000000000000001\n",
" 0.08100000000000002\n",
" 0.0729\n",
" 0.06561\n",
" 0.05904900000000001\n",
" 0.05314410000000001\n",
" 0.04782969000000001\n",
" 0.04304672100000001\n",
" 0.03874204890000001\n",
" 0.03486784401000001\n",
" 0.031381059609000006\n",
" 0.028242953648100012\n",
" 0.02541865828329001\n",
" 0.02287679245496101\n",
" 0.02058911320946491\n",
" 0.018530201888518418\n",
" 0.016677181699666577\n",
" 0.015009463529699918\n",
" 0.013508517176729929"
]
},
"execution_count": 46,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"p = 0.1\n",
"\n",
"probs = [p * (1 - p)^(n-1) for n in 1:20]"
]
},
{
"cell_type": "code",
"execution_count": 47,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"20-element Array{Float64,1}:\n",
" 0.1\n",
" 0.19\n",
" 0.271\n",
" 0.3439000000000001\n",
" 0.40951000000000004\n",
" 0.46855900000000006\n",
" 0.5217031000000001\n",
" 0.5695327900000001\n",
" 0.6125795110000002\n",
" 0.6513215599000002\n",
" 0.6861894039100002\n",
" 0.7175704635190002\n",
" 0.7458134171671003\n",
" 0.7712320754503903\n",
" 0.7941088679053513\n",
" 0.8146979811148162\n",
" 0.8332281830033346\n",
" 0.8499053647030012\n",
" 0.864914828232701\n",
" 0.878423345409431"
]
},
"execution_count": 47,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"s = cumsum(probs)"
]
},
{
"cell_type": "code",
"execution_count": 48,
"metadata": {},
"outputs": [],
"source": [
"using Plots"
]
},
{
"cell_type": "code",
"execution_count": 50,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 50,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"plot(s, m=:o)"
]
},
{
"cell_type": "code",
"execution_count": 58,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"41-element Array{Float64,1}:\n",
" 0.0\n",
" 0.0\n",
" 0.1\n",
" 0.1\n",
" 0.19\n",
" 0.19\n",
" 0.271\n",
" 0.271\n",
" 0.3439000000000001\n",
" 0.3439000000000001\n",
" 0.40951000000000004\n",
" 0.40951000000000004\n",
" 0.46855900000000006\n",
" ⋮\n",
" 0.7712320754503903\n",
" 0.7941088679053513\n",
" 0.7941088679053513\n",
" 0.8146979811148162\n",
" 0.8146979811148162\n",
" 0.8332281830033346\n",
" 0.8332281830033346\n",
" 0.8499053647030012\n",
" 0.8499053647030012\n",
" 0.864914828232701\n",
" 0.864914828232701\n",
" 0.878423345409431"
]
},
"execution_count": 58,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"ys = [0; reduce(vcat, [ [s[n], s[n]] for n in 1:20 ])]\n",
"\n",
"pop!(ys)\n",
"pushfirst!(ys, 0)"
]
},
{
"cell_type": "code",
"execution_count": 59,
"metadata": {},
"outputs": [],
"source": [
"xs = [0; reduce(vcat, [ [n, n] for n in 1:20 ])];"
]
},
{
"cell_type": "code",
"execution_count": 68,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 68,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"plot(xs, ys, leg=false)\n",
"scatter!(s)\n",
"\n",
"ylabel!(\"cumulative probability of recovery up to time n\")\n",
"xlabel!(\"number of days, n\")"
]
},
{
"cell_type": "code",
"execution_count": 70,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 70,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"plot(s, m=:o)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Think of our discrete-time process in steps, not of days, but of a time $\\delta$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Instead of \"cumulative probability that have decayed by day $n$\", want to talk about\n",
"\"cumulative probability that have decayed by time $n \\delta$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"What is the probability $p(\\delta)$ of recovering in time $\\delta$"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "Julia 1.4.1",
"language": "julia",
"name": "julia-1.4"
},
"language_info": {
"file_extension": ".jl",
"mimetype": "application/julia",
"name": "julia",
"version": "1.4.1"
}
},
"nbformat": 4,
"nbformat_minor": 4
}