L6maple.mws, Markovin ketjut

ke 2.10.02   Luento 11

Lähteet

Lay Ch 4.9 ss. 288 - 296. (ei prujuja)

Alustukset

>    restart: with(LinearAlgebra):with(linalg):

Warning, the previous binding of the name GramSchmidt has been removed and it now has an assigned value

Warning, the protected names norm and trace have been redefined and unprotected

>    alias(rref=ReducedRowEchelonForm):alias(ref=GaussianElimination):
alias(Id=IdentityMatrix):

 Markovin ketjut

Lay Ch 4.9

Mallinnetaan moninaisia ilmiöitä mm. biologiassa, liiketaloudessa, kemiassa, fysiikassa, insinööritieteissä,...

Malli kuvaa koetta tai mittaussarjaa, joka toistetaan monta kertaa samalla tavalla.

Jokainen koetulos on vektori, jonka alkiot riippuvat vain edellisen tulosvektorin komponenteista.

>    restart: with(LinearAlgebra):

Esim. Kaupugin väestöjakauma

Kaupungin väestö koostukoon 1) kantakaupunkilaisista ja 2) esikaupunkilaisista.

Oletetaan, että vuotuinen muuttoliike tapahtuu seuraavan matriisin määräämällä tavalla:

>    M:=<<.95,.05>|<.03,.97>>;

M := Matrix(%id = 4806268)

>    Kantakaupungista:=M[1..-1,1];

Kantakaupungista := Vector(%id = 4810464)

>    Esikaupungista:=M[1..-1,2];

Esikaupungista := Vector(%id = 4917576)

>    Kantakaupunkiin:=M[1,1..-1];

Kantakaupunkiin := Vector(%id = 4621388)

>    Esikaupunkiin:=M[2,1..-1];

Esikaupunkiin := Vector(%id = 4622148)

>    X[0]:=<.6,.4>;

X[0] := Vector(%id = 4622188)

>    X[1]:=M.X[0];

X[1] := Vector(%id = 5039432)

>    for j from 0 to 10 do X[j+1]:=M.X[j]: end do:

>    seq(X[j],j=0..10);

Vector(%id = 4622188), Vector(%id = 4773216), Vector(%id = 5072648), Vector(%id = 5105720), Vector(%id = 5220768), Vector(%id = 5194488), Vector(%id = 12118240), Vector(%id = 5118700), Vector(%id = 121...
Vector(%id = 4622188), Vector(%id = 4773216), Vector(%id = 5072648), Vector(%id = 5105720), Vector(%id = 5220768), Vector(%id = 5194488), Vector(%id = 12118240), Vector(%id = 5118700), Vector(%id = 121...
Vector(%id = 4622188), Vector(%id = 4773216), Vector(%id = 5072648), Vector(%id = 5105720), Vector(%id = 5220768), Vector(%id = 5194488), Vector(%id = 12118240), Vector(%id = 5118700), Vector(%id = 121...

>    for j from 10 to 20 do X[j+1]:=M.X[j]: end do:

>    seq(X[j],j=11..20);

Vector(%id = 12298988), Vector(%id = 12224764), Vector(%id = 12255000), Vector(%id = 12414584), Vector(%id = 12404556), Vector(%id = 12495472), Vector(%id = 12457968), Vector(%id = 12574280), Vector(%i...
Vector(%id = 12298988), Vector(%id = 12224764), Vector(%id = 12255000), Vector(%id = 12414584), Vector(%id = 12404556), Vector(%id = 12495472), Vector(%id = 12457968), Vector(%id = 12574280), Vector(%i...

>    for j from 20 to 100 do X[j+1]:=M.X[j]: end do:

>    seq(X[j],j=[30,40,50,60,70,80,90,100]);

Vector(%id = 12978324), Vector(%id = 12684372), Vector(%id = 13754744), Vector(%id = 13928984), Vector(%id = 14456036), Vector(%id = 14683196), Vector(%id = 14443172), Vector(%id = 5194536)
Vector(%id = 12978324), Vector(%id = 12684372), Vector(%id = 13754744), Vector(%id = 13928984), Vector(%id = 14456036), Vector(%id = 14683196), Vector(%id = 14443172), Vector(%id = 5194536)

>   

Lähestyykö tasapainotilaa . Vaikuttiko alkupisteen valinta  asiaan? Kokeillaan aivan toista.

>    X[0]:=<0.8,0.2>;

X[0] := Vector(%id = 12186056)

>    for j from 0 to 100 do X[j+1]:=M.X[j]: end do:

>    seq(X[j],j=[50,60,70,80,90,100]);

Vector(%id = 13916308), Vector(%id = 14398756), Vector(%id = 14798396), Vector(%id = 15130116), Vector(%id = 15369748), Vector(%id = 15612984)
Vector(%id = 13916308), Vector(%id = 14398756), Vector(%id = 14798396), Vector(%id = 15130116), Vector(%id = 15369748), Vector(%id = 15612984)

Eipä näyttäisi juurikaan vaikuttavan .

>    q:=<.375,.625>;q:=convert(q,rational); M:=convert(M,rational);

q := Vector(%id = 15761252)

q := Vector(%id = 13359852)

M := Matrix(%id = 15793600)

>    M.q;

Vector(%id = 4331192)

q on kiintopiste (tasapainotila, "steady state").

Tästä prosessi ei etene mihinkään.

Esim. Vaalit

>    restart: with(LinearAlgebra):alias(Id=IdentityMatrix):

>    <<Demokraateilta>|<Republikaaneilta>|<Liberaaleilta>>;
P:=<<.7,.2,.1>|<.1,.8,.1>|<.3,.3,.4>>;

Matrix(%id = 12150168)

P := Matrix(%id = 12207548)

>    X[0]:=<.55,.40,.05>;

X[0] := Vector(%id = 12311952)

>    for j from 0 to 10 do X[j+1]:=P.X[j]: end do:

>    seq(X[j],j=0..10);

Vector(%id = 12311952), Vector(%id = 12267328), Vector(%id = 12482284), Vector(%id = 12569972), Vector(%id = 12536060), Vector(%id = 12627320), Vector(%id = 12593612), Vector(%id = 12704304), Vector(%i...
Vector(%id = 12311952), Vector(%id = 12267328), Vector(%id = 12482284), Vector(%id = 12569972), Vector(%id = 12536060), Vector(%id = 12627320), Vector(%id = 12593612), Vector(%id = 12704304), Vector(%i...
Vector(%id = 12311952), Vector(%id = 12267328), Vector(%id = 12482284), Vector(%id = 12569972), Vector(%id = 12536060), Vector(%id = 12627320), Vector(%id = 12593612), Vector(%id = 12704304), Vector(%i...

>   

>   

>   

Taspainotila "Steady state" (kiintopiste)

>    restart: with(LinearAlgebra):with(linalg):

Warning, the previous binding of the name GramSchmidt has been removed and it now has an assigned value

Warning, the protected names norm and trace have been redefined and unprotected

>    alias(rref=ReducedRowEchelonForm):alias(ref=GaussianElimination):
alias(Id=IdentityMatrix):

>   

Kaupunkimatriisi

>    M:=<<.95,.05>|<.03,.97>>; M:=convert(M,rational);

M := Matrix(%id = 14664252)

M := Matrix(%id = 14614904)

               M x = x  <==> (M - I) x = 0

Tehtävänä on siis etsiä vektori x#0 siten, että (M - I) x = 0,

ts. (HY):lle ei-triv. ratk.

>    M1:=M-Id(2);

M1 := Matrix(%id = 14828704)

>    rref(M1);

Matrix(%id = 14887272)

>    x2:=t: x1:=3/5*t;

x1 := 3/5*t

>    T:=solve(x1+x2=1,t);

T := 5/8

>    x:=subs(t=T,<x1,x2>);

x := Vector(%id = 14987140)

>    convert(%,float);

Vector(%id = 14755212)

>   

>   

Normaalisti lasketaan kaikki liukuluvuilla.

Huom!  Tämähän oli puhdaslinjainen ominaisvektorilasku . Taspainovektorihan  tarkoittaa ominaisarvoa 1 vastaavaa ominaisvetoria.

Näyttäisi siis siltä, että stokastisilla matriiseilla on ominaisarvo 1, ainakin tämän esimerkin valossa.

Asia saa vahvistuksen harj. 4:n AV-tehtävässä.