using LinearAlgebra
using Plots
using Interact
using Random: seed!
using Distributions: wsample
using Printf
The WebIO Jupyter extension was not detected. See the WebIO Jupyter integration documentation for more information.
P = [
0.0 0.0 1.0 0.0;
1.0 0.0 0.0 0.0;
0.0 0.6 0.0 1.0;
0.0 0.4 0.0 0.0;
]
4×4 Matrix{Float64}:
0.0 0.0 1.0 0.0
1.0 0.0 0.0 0.0
0.0 0.6 0.0 1.0
0.0 0.4 0.0 0.0
function plot_circle!(x,y;r=10, ic=0)
t = range(0, 2π, 101)
plot!(x .+ r * cos.(t), y .+ r * sin.(t), label="")
plot!(annotate=(x, y+0.8r, "$ic"))
end
plot_circle! (generic function with 1 method)
xc = [-1, 1, -1, 1] * 20
yc = [1, 1, -1, -1] * 20
4-element Vector{Int64}:
20
20
-20
-20
# function for plotting the markov chain diagram
function plot_chain()
plot()
for ic=1:4
plot_circle!(xc[ic], yc[ic], ic=ic)
end
for ii=1:4
if P[ii,ii] != 0
@show "Bug"
#add_loop!()
end
for jj=1:4
if ii != jj && P[ii,jj] != 0
# @show ii, jj, P[ii,jj]
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
# @show rd, frac
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],
annotate=(xm+xo, ym+yo, "$(P[ii,jj])"), arrow=:arrow, label="")
end
end
end
plot!(aspect_ratio=1, axis=:on, title="$totalstep steps")
end
#plot_chain()
plot_chain (generic function with 1 method)
# initial conditions for 100 points
xg = repeat(-4.5:4.5, 1, 10)
yg = xg'
xg = xg[:]
yg = yg[:]
node = 2*ones(Int, length(xg)) # all start in one node
totalstep = 0
plot_chain()
scatter!(xc[node] + xg[:], yc[node] + yg[:], label="")
seed!(0)
function node_update!(node)
global totalstep += 1
for kk=1:length(node)
node[kk] = wsample(1:4, P[:,node[kk]])
end
end
node_update! (generic function with 1 method)
function run_and_plot(;niter=1)
for iter = 1:niter
node_update!(node)
plot_chain()
scatter!(xc[node] + xg[:], yc[node] + yg[:], label="")
end
for ii=1:4
tmp = sum(node .== ii) / length(node)
plot!(annotate=(xc[ii]+10, yc[ii]-10, "$tmp"))
end
plot!()
end
run_and_plot()
Which state has the lowest probability in equilibrium?
@manipulate for niter = ["1", "5", "100"], go=button("go")
run_and_plot(niter=parse(Int, niter))
end
function matprint(V) # function to nicely print the eigenvector matrix
for i=1:size(V,1)
for j=1:size(V,2)
v = V[i,j]
@printf("%5.2f", real(v))
print(imag(v) < 0 ? " -" : " +")
@printf("%5.2fı ", abs(imag(v)))
end
println()
end
end
matprint (generic function with 1 method)
(d, V) = eigen(P)
Eigen{ComplexF64, ComplexF64, Matrix{ComplexF64}, Vector{ComplexF64}}
values:
4-element Vector{ComplexF64}:
-0.5325796049075049 + 0.0im
-0.2337101975462476 - 0.8345303914291096im
-0.2337101975462476 + 0.8345303914291096im
0.9999999999999998 + 0.0im
vectors:
4×4 Matrix{ComplexF64}:
0.383539+0.0im 0.146985+0.524853im 0.146985-0.524853im -0.562544+0.0im
-0.720154+0.0im -0.62892-0.0im -0.62892+0.0im -0.562544+0.0im
-0.204265+0.0im 0.403654-0.245327im 0.403654+0.245327im -0.562544+0.0im
0.54088+0.0im 0.0782812-0.279526im 0.0782812+0.279526im -0.225018+0.0im
[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
matprint(V)
0.38 + 0.00ı 0.15 + 0.52ı 0.15 - 0.52ı -0.56 + 0.00ı -0.72 + 0.00ı -0.63 + 0.00ı -0.63 + 0.00ı -0.56 + 0.00ı -0.20 + 0.00ı 0.40 - 0.25ı 0.40 + 0.25ı -0.56 + 0.00ı 0.54 + 0.00ı 0.08 - 0.28ı 0.08 + 0.28ı -0.23 + 0.00ı
v = Real.(V[:,4])
πss = v / sum(v)
4-element Vector{Float64}:
0.2941176470588236
0.2941176470588233
0.29411764705882365
0.11764705882352948
(Do after covering irreducible matrix, aperiodic matrix, strongly connected, etc.)
The transition matrix P in this example is (choose most specific correct answer):
# P^10