CoCalc Public FilesMarkovketten.sagewsOpen with one click!
Authors: Johannes Hatzl, Walter Hochfellner, Lisa Hochfellner
Views : 58
Description: Markov-Ketten-Demo
Compute Environment: Ubuntu 18.04 (Deprecated)
p0=vector([.75, .25])
M=matrix(QQ,[[0.8,0.2],[0.1,0.9]])
g = DiGraph( {'A':{'A':M[0,0].n(digits = 1),'B':M[0,1].n(digits = 1)}, 'B':{'A':M[1,0].n(digits = 1), 'B':M[1,1].n(digits = 1)}}) g.graphplot(prog = "dot", vertex_size=1500, edge_labels=True ,edge_color='lightgray',loop_size=.15, dist=.12, pos={'A': [0, 0], 'B': [ 2, 0]}).show()
# Marktanteile nach 1 Jahr p1 = p0 * M p1
(0.625000000000000, 0.375000000000000)
# Marktanteile nach 2 Jahren p2 = p1 * M p2 = p0 * M^2 p2
(0.537500000000000, 0.462500000000000)
# Marktanteile nach 3 Jahren p3 = p2 * M p3 = p0 * M^3 p3
(0.476250000000000, 0.523750000000000)
p0 * M^4
(0.433375000000000, 0.566625000000000)
p0 * M^5
(0.403362500000000, 0.596637500000000)
p0 * M^6
(0.382353750000000, 0.617646250000000)
p0 * M^10
(0.345103135375000, 0.654896864625000)
p0 * M^20
(0.333665801109573, 0.666334198890427)
p0 * M^50
(0.333333340826938, 0.666666659173062)
var('a,b') eq =[ 0.8 * a + 0.1 * b == a , 0.2 * a + 0.9 * b == b , a + b == 1 ] solve(eq,a,b)
(a, b) [[a == (1/3), b == (2/3)]]
# Alternative Formel zur Berechnung des Gleichgewichtszustandes v = M.eigenvectors_left()[0][1][0] p = v/sum(v) p
(1/3, 2/3)
# Kontrolle, ob Gleichgewicht vorliegt p.numerical_approx()
(0.333333333333333, 0.666666666666667)
M = matrix(QQ, [[.87, 0 , .1 , 0 , .03, 0 ], [.85, 0 , 0 , .11, 0 , .04 ], [0 , .82, 0 , 0 , .13, .05 ], [0 , 0 , .82, 0 , 0 , 0.18], [0 , 0 , 0 , .8 , 0 , .2 ], [0 , 0 , 0 , 0 , .78, .22 ]])
# Kontrolle, ob Übergangsmatrix korrekt ist (alle Zeilensummen = 1) sum(M.transpose())
(1, 1, 1, 1, 1, 1)
G=DiGraph(M.numerical_approx(digits=2),format='weighted_adjacency_matrix') G.plot(layout = "circular", prog = "dot",edge_labels=True,edge_color='gray',loop_size=.15, dist=.12)
p0 = vector([.63, .05, .06, .06, .08, .12])
# Verteilung nach 1 Jahr p1 = p0 * M p1.numerical_approx(digits=3)
(0.591, 0.0492, 0.112, 0.0695, 0.120, 0.0582)
# Verteilung nach 2 Jahren p2 = p0 * M^2 p2.numerical_approx(digits=3)
(0.556, 0.0920, 0.116, 0.102, 0.0777, 0.0570)
# Verteilung nach 5 Jahren p5 = p0 * M^5 p5.numerical_approx(digits=3)
(0.592, 0.0947, 0.116, 0.0753, 0.0721, 0.0502)
# Gleichgewichtszustand v = M.eigenvectors_left()[0][1][0] p = v/sum(v) p.numerical_approx(digits=3)
(0.614, 0.0939, 0.115, 0.0648, 0.0681, 0.0446)
# Kontrolle p * M.numerical_approx(digits=3)
(0.614, 0.0939, 0.115, 0.0648, 0.0681, 0.0446)
# Bei Erstanmeldung wird man in Prämienklasse~3 eingestuft. p0 = vector([0,0,0,1,0,0])
# Wie groß ist die Wahrscheinlichkeit, nach 5~Jahren ebenfalls in Klasse~3 eingestuft zu sein? p5 = p0 * M^5 p5.numerical_approx(digits=3)
(0.433, 0.229, 0.0910, 0.0724, 0.115, 0.0592)
# Mit welcher Wahrscheinlichkeit befindet man sich nach 5~Jahren in Stufe~0? p5[0].numerical_approx(digits=3)
0.433
# Mit welcher Wahrscheinlichkeit befindet man sich nach 25 Jahren in Stufe~0? p25 = p0 * M^25 p25[0].numerical_approx(digits=3)
0.614
def step(pos): r=zero_vector(QQ,pos.length()) for i in range(pos.length()): for j in range(pos[i]): X=GeneralDiscreteDistribution(M[i]) x=X.get_random_element() r[x]+=1 return r def color_dict(pos): d={} s=max(pos) for i in range(pos.length()): d[Color('lightgreen').blend(Color('pink'), fraction=pos[i]/s)]=[] for i in range(pos.length()): d[Color('lightgreen').blend(Color('pink'), fraction=pos[i]/s)]+=[i] return d def plot_pos(p): d = color_dict(p) G=DiGraph(M.numerical_approx(digits=2),format='weighted_adjacency_matrix') return G.plot(layout = "circular", prog = "dot",edge_labels=True,edge_color='lightgray',loop_size=.15, dist=.12, vertex_colors=d)
# Wähle Startposition (rot markiert) q = vector([0,0,0,1,0,0]) plot_pos(q)
# Simuliere einzelnen 'Schritt' (1 Jahr) q=step(q) plot_pos(q)
def simulate_steps(q,n): plots=[] z=[q] for i in range(n): q=step(q) plots.append(plot_pos(q)) z[0]=q print('Distribution after '+str(n)+' steps: ' +str(z[0])) print( z[0]/sum(z[0]).n(digits=2)) a=animate(plots) return a.apng(savefile='markov.png')
# Simuliere 30 Schritte q = vector([0,0,0,1,0,0]) simulate_steps(q,30)
Distribution after 30 steps: (0, 0, 0, 0, 0, 1) (0.00, 0.00, 0.00, 0.00, 0.00, 1.0)
q = vector([0,0,0,0,300,700]) plot_pos(q)
# Simuliere einzelnen 'Schritt' (1 Jahr) q=step(q) print(q) plot_pos(q)
(0, 0, 0, 242, 519, 239)
# Simuliere 10 Schritte q = vector([0,0,0,0,100,900]) simulate_steps(q,10)
Distribution after 10 steps: (548, 101, 125, 89, 80, 57) (0.55, 0.10, 0.12, 0.089, 0.080, 0.057)
# Theoretische Gleichgewichtsverteilung p.numerical_approx(digits=3)
(0.614, 0.0939, 0.115, 0.0648, 0.0681, 0.0446)