GitHub

Markov chain

This example illustrates a Markov chain using the Julia language.

In this demo, the elements of the transition matrix $P$ are $p_{ij} = p(X_{k+1} = i | X_{k} = j)$.

This page comes from a single Julia file: markov-chain.jl.

You can access the source code for such Julia documentation using the 'Edit on GitHub' link in the top right. You can view the corresponding notebook in nbviewer here: markov-chain.ipynb, or open it in binder here: markov-chain.ipynb.

Setup

Add the Julia packages used in this demo. Change false to true in the following code block if you are using any of the following packages for the first time.

if false
    import Pkg
    Pkg.add([
        "InteractiveUtils"
        "LaTeXStrings"
        "LinearAlgebra"
        "MIRTjim"
        "Plots"
        "Random"
        "StatsBase"
    ])
end

Tell Julia to use the following packages. Run Pkg.add() in the preceding code block first, if needed.

using LinearAlgebra: eigen
using MIRTjim: prompt
using Plots: annotate!, default, plot, plot!, scatter, scatter!
using Plots: gif, @animate, Plots
using Random: seed!
using StatsBase: wsample

default(
 markerstrokecolor = :auto, label="",
 labelfontsize=8, legendfontsize=8, size = (1,1) .* 600,
)
seed!(0);

The following line is helpful when running this jl-file as a script; this way it will prompt user to hit a key after each image is displayed.

isinteractive() && prompt(:prompt);

Transition matrix

p = 0.4
P = [
 0 0.0 1 0;
 1 0.0 0 0;
 0 1-p 0 1;
 0  p  0 0;
];

Define plot helpers

color = [:red, :green, :blue, :purple];
xc = [-1, 1, -1, 1] * 20 # node centers
yc = [1, 1, -1, -1] * 20

function plot_circle!(x, y, ic=0; r=10)
    t = range(0, 2π, 101)
    plot!(x .+ r * cos.(t), y .+ r * sin.(t), color=color[ic], width=2)
    annotate!(x, y+0.7r, ("$ic", 14, color[ic]))
end;

function for plotting the Markov chain diagram

function plot_chain(steps::Int)

    plot(xtick=nothing, ytick=nothing, axis=:off, aspect_ratio=1,)
    plot_circle!.(xc, yc, 1:4)

    for ii in 1:4
        if P[ii,ii] != 0
            @show "Bug"
            #add_loop!()
        end
        for jj in 1:4
            ((ii == jj) || P[ii,jj] == 0) && continue
            xi = xc[ii]
            xj = xc[jj]
            yi = yc[ii]
            yj = yc[jj]
            xd = xi - xj
            yd = yi - yj
            phi = atan(yd, xd) + π/2
            rd = sqrt(xd^2 + yd^2)
            frac = (rd - 1*10) / rd
            xs = frac * xj + (1-frac) * xi
            xe = frac * xi + (1-frac) * xj
            ys = frac * yj + (1-frac) * yi
            ye = frac * yi + (1-frac) * yj
            xm = (xi + xj) / 2
            ym = (yi + yj) / 2
            xo = 5 * cos(phi)
            yo = 5 * sin(phi)
            plot!([xs, xe], [ys, ye], arrow=:arrow, color=:black, width=2)
            annotate!(xm+xo, ym+yo, ("$(P[ii,jj])", 13))
        end
    end
    title = steps > 0 ? "$steps steps" : "Initial state"
    plot!(; title)
end;

Initial conditions

This block starts the simulation with an initial grid of 100 particles

xg = repeat(-4.5:4.5, 1, 10)
yg = xg'
xg = xg[:]
yg = yg[:]
node = fill(2, length(xg)) # all particles start in node (state) X_0 = 2
plot_chain(0)
scatter!(xc[node] + xg[:], yc[node] + yg[:])
Example block output

Use wsample (weighted sampling) for transitions

function node_update!(node)
    for kk in 1:length(node)
        node[kk] = wsample(1:size(P,1), P[:,node[kk]]) # random process
    end
end;

function run_and_plot(iter::Int)
    node_update!(node)
    plot_chain(iter)
    scatter!(xc[node] + xg[:], yc[node] + yg[:])
    for ii in 1:4
        tmp = sum(node .== ii) / length(node)
        annotate!(xc[ii]+10, yc[ii]-10, ("$tmp", 7, color[ii]))
    end
    plot!()
end;

Simulate

anim1 = @animate for iter in [1:20; 30:10:100]
    run_and_plot(iter)
end
gif(anim1; fps = 4)
Example block output

Clicker question 1

The transition matrix P in this example is (choose most specific correct answer):

  • A. Square
  • B. Nonnegative
  • C. Irreducible
  • D. Primitive
  • E. Positive"

Clicker question 2 (later)

Which state has the lowest probability in equilibrium?

  • A 1
  • B 2
  • C 3
  • D 4
  • E None: they are all equally likely"

Eigenvectors:

(d, V) = eigen(P)
round.(V; digits=3)
4×4 Matrix{ComplexF64}:
  0.384+0.0im   0.147+0.525im   0.147-0.525im  -0.563+0.0im
  -0.72+0.0im  -0.629-0.0im    -0.629+0.0im    -0.563+0.0im
 -0.204+0.0im   0.404-0.245im   0.404+0.245im  -0.563+0.0im
  0.541+0.0im   0.078-0.28im    0.078+0.28im   -0.225+0.0im

Eigenvalues:

[d abs.(d)] # exactly one λ=1 and only one |λ| = 1
4×2 Matrix{ComplexF64}:
 -0.53258+0.0im       0.53258+0.0im
 -0.23371-0.83453im  0.866638+0.0im
 -0.23371+0.83453im  0.866638+0.0im
      1.0+0.0im           1.0+0.0im

Plot eigenvalues in complex plane

scatter(real(d), imag(d), color=:blue,
 xaxis = ("Re(λ)", -1:1),
 yaxis = ("Im(λ)", -1:1),
 framestyle = :origin,
 size = (400,400),
)
tmp = range(0, 2π, 301)
plot!(cos.(tmp), sin.(tmp), color=:black)
Example block output
prompt()

Steady-state distribution

v = real(V[:,4])
πss = v / sum(v) # normalize
[πss; "check:"; 1 / (p + 4 - 1); p / (p + 4 - 1)]
7-element Vector{Any}:
 0.2941176470588236
 0.2941176470588233
 0.29411764705882365
 0.11764705882352948
  "check:"
 0.2941176470588235
 0.11764705882352941

For insight: $4^2 - 2 ⋅ 4 + 2 = N^2 - 2N + 2$ in Ch8

P^10
4×4 Matrix{Float64}:
 0.432  0.288   0.216  0.16
 0.216  0.432   0.16   0.48
 0.288  0.1936  0.432  0.216
 0.064  0.0864  0.192  0.144

Approximate limiting distribution

P^200
4×4 Matrix{Float64}:
 0.294118  0.294118  0.294118  0.294118
 0.294118  0.294118  0.294118  0.294118
 0.294118  0.294118  0.294118  0.294118
 0.117647  0.117647  0.117647  0.117647

Reproducibility

This page was generated with the following version of Julia:

using InteractiveUtils: versioninfo
io = IOBuffer(); versioninfo(io); split(String(take!(io)), '\n')
12-element Vector{SubString{String}}:
 "Julia Version 1.12.1"
 "Commit ba1e628ee49 (2025-10-17 13:02 UTC)"
 "Build Info:"
 "  Official https://julialang.org release"
 "Platform Info:"
 "  OS: Linux (x86_64-linux-gnu)"
 "  CPU: 4 × AMD EPYC 7763 64-Core Processor"
 "  WORD_SIZE: 64"
 "  LLVM: libLLVM-18.1.7 (ORCJIT, znver3)"
 "  GC: Built with stock GC"
 "Threads: 1 default, 1 interactive, 1 GC (on 4 virtual cores)"
 ""

And with the following package versions

import Pkg; Pkg.status()
Status `~/work/book-la-demo/book-la-demo/docs/Project.toml`
  [6e4b80f9] BenchmarkTools v1.6.3
  [aaaa29a8] Clustering v0.15.8
  [35d6a980] ColorSchemes v3.31.0
  [3da002f7] ColorTypes v0.12.1
  [c3611d14] ColorVectorSpace v0.11.0
  [717857b8] DSP v0.8.4
  [72c85766] Demos v0.1.0 `~/work/book-la-demo/book-la-demo`
  [e30172f5] Documenter v1.16.0
  [4f61f5a4] FFTViews v0.3.2
  [7a1cc6ca] FFTW v1.10.0
  [587475ba] Flux v0.16.5
  [a09fc81d] ImageCore v0.10.5
  [71a99df6] ImagePhantoms v0.8.1
  [b964fa9f] LaTeXStrings v1.4.0
  [7031d0ef] LazyGrids v1.1.0
  [599c1a8e] LinearMapsAA v0.12.0
  [98b081ad] Literate v2.20.1
  [7035ae7a] MIRT v0.18.2
  [170b2178] MIRTjim v0.26.0
  [eb30cadb] MLDatasets v0.7.18
⌅ [efe261a4] NFFT v0.13.7
  [6ef6ca0d] NMF v1.0.3
  [15e1cf62] NPZ v0.4.3
  [0b1bfda6] OneHotArrays v0.2.10
  [429524aa] Optim v1.13.2
  [91a5bcdd] Plots v1.41.1
  [f27b6e38] Polynomials v4.1.0
  [2913bbd2] StatsBase v0.34.8
  [d6d074c3] VideoIO v1.4.0
  [b77e0a4c] InteractiveUtils v1.11.0
  [37e2e46d] LinearAlgebra v1.12.0
  [44cfe95a] Pkg v1.12.0
  [9a3f8284] Random v1.11.0
Info Packages marked with ⌅ have new versions available but compatibility constraints restrict them from upgrading. To see why use `status --outdated`

This page was generated using Literate.jl.