Macierze Lesliego i Markowa – K. Leśniak
1
Model Lesliego
Wyodrębniamy w populacji k grup wiekowych. Po każdej
jednostce czasu następują narodziny i zgony oraz starze-
nie (przechodzenie do następnej grupy wiekowej). Chcemy
symulować zmiany stanu liczebności w poszczególnych gru-
pach wiekowych.
Oznaczmy:
0
6 m
i
– liczba potomstwa pojawiającego się co jednost-
kę czasu u osobnika z i-tej grupy wiekowej, i = 1, . . . , k,
[0, 1] 3 s
i
– przeżywalność osobników z i-tej grupy wie-
kowej dożywających przynależności do następnej (i + 1)-ej
grupy wiekowej (jaki procent odobników i-tej grupy dożywa
awansu do (i + 1)-ej grupy wiekowej), i = 1, . . . , k − 1,
x
n,i
– liczba osobników z i-tej grupy wiekowej w n-tej
jednostce czasu, i = 1, . . . , k.
2
Macierze Lesliego i Markowa – K. Leśniak
Struktura wiekowa populacji w n-tej jednostce czasu
x
n
= [x
n,1
, x
n,2
, . . . , x
n,k
]
T
podlega zależnościom:
x
n+1,1
=
k
X
i=1
m
i
· x
n,i
– liczba nowych osobników (wszystkie grupy wiekowe mogą
wydawać potomstwo),
x
n+1,i+1
= s
i
· x
n,i
,
i = 1, 2, . . . , (k − 1),
– liczba osobników, które postarzały się awansując do na-
stępnej grupy wiekowej (pozostałe umarły).
Macierzowo: x
n+1
= M · x
n
, gdzie
M =
m
1
, m
2
, . . . , m
k−1
, m
k
s
1
,
0, . . . ,
0,
0
0,
s
2
, . . . ,
0,
0
. . . . . . . . .
. . .
. . .
0,
0, . . . , s
k−1
,
0
.
Powyższy wzór rekurencyjny prowadzi do wzoru „ jawnego”:
x
n
= M
n
· x
0
.
Macierze Lesliego i Markowa – K. Leśniak
3
Przykład. W pewnej populacji wyróżniają się dwie grupy
wiekowe:
1 – osobniki niedojrzałe, które nie mogą się rozmnażać
(m
1
= 0 o przeżywałności s (s
1
= s), oraz
2 – osobniki dojrzałe o płodności m, które mogą się roz-
mnażać (m
2
= m).
Macierz Lesliego: M =
m
1
m
2
s
1
0
=
0 m
s 0
.
Przez indukcję można zauważyć, że:
M
2n+1
=
0
m
n+1
· s
n
s
n+1
· m
n
0
, M
2n
=
m
n
· s
n
0
0
s
n
· m
n
.
Zatem
x
2n+1
= M
2n+1
· x
0
=
0
m
n+1
· s
n
s
n+1
· m
n
0
·
x
0,1
x
0,2
=
=
m
n+1
· s
n
· x
0,2
s
n+1
· m
n
· x
0,1
= (m · s)
n
· [m · x
0,2
, s · x
0,1
]
T
,
x
2n
= M
2n
· x
0
=
m
n
· s
n
0
0
s
n
· m
n
·
x
0,1
x
0,2
= (m · s)
n
· x
0
.
Zachowanie asymptotyczne:
x
n
=
x
n,1
x
n,2
−→
n→∞
[
0
0
] ,
gdy 0 < ms < 1,
[
∞
∞
] , gdy ms > 1.
4
Macierze Lesliego i Markowa – K. Leśniak
Przykład (cykliczne zmiany struktury wieku).
M =
0 0 4
1
3
0 0
0
3
4
0
. Łatwo dostrzec, że: M
3
=
"
1 0 0
0 1 0
0 0 1
#
.
Zachowanie asymptotyczne:
x
3n
= x
0
= [x
01
, x
02
, x
03
]
T
,
x
3n+1
= x
1
= [x
11
, x
12
, x
13
]
T
=
4 x
03
,
1
3
x
01
,
3
4
x
02
T
,
x
3n+2
= x
2
= [x
21
, x
22
, x
23
]
T
=
3 x
02
,
4
3
x
03
,
1
4
x
01
T
.
Cykl powtarza się co trzy jednostki czasu:
x
0
, x
1
, x
2
, x
0
, x
1
, . . .
Macierze Lesliego i Markowa – K. Leśniak
5
Przykład (stabilna struktura wieku).
M =
2
1
1
0.25
0
0
0
0.(3) 0
.
W tabeli zbieramy dane o procentowym udziale w liczeb-
ności populacji poszczególnych grup wiekowych po trzech
iteracjach.
liczebność
początkowa
udział procentowy
[x
n,i
/
P
3
i=1
x
n,i
, i = 1, 2, 3]
[1, 0, 0]
[49, 772%, 26, 941%, 23, 288%]
[0, 20, 0]
[49, 514%, 27, 185%, 23, 301%]
[0, 0, 8]
[50, 000%, 24.999%, 24, 999%]
[100, 100, 100]
[49, 749%, 26, 936%, 23, 318%]
[250, 10000, 100] [49, 561%, 27, 138%, 23, 300%]
[100000, 5000, 10] [49, 771%, 26, 942%, 23, 288%]
[200, 45000, 6000] [49, 534%, 27, 133%, 23, 334%]
Struktura wiekowa populacji w stosunkowo krótkim czasie
stabilizuje się i zależy od takich parametrów specyficznych
dla populacji jak śmiertelność i rozrodczość.
6
Macierze Lesliego i Markowa – K. Leśniak
Macierze Markowa
M =
p
11
p
12
p
13
. . . p
1s
p
21
p
22
p
23
. . . p
2s
. . . . . . . . . . . . . . .
p
s1
p
s2
p
s3
. . . p
ss
= [p
ij
]
1
6i,j6s
— macierz stochastyczna (= macierz Markowa
= macierz przejścia), gdy
∀
i,j
0
6 p
ij
6 1,
∀
i
s
X
j=1
p
ij
= 1.
Terminologia:
• {1, . . . , s} – zbiór stanów układu;
• p
ij
– prawdopodobieństwo przejścia
ze stanu i w stan j, i, j ∈ {1, . . . , s};
• M = [p
ij
]
i,j∈{1,...,s}
– macierz przejścia
dla jednorodnego łańcucha Markowa;
• p
(n)
ij
– prawdopodobieśtwo z jakim cząstka, która jest w
stanie i znajdzie się w stanie j po n cyklach.
Macierze Lesliego i Markowa – K. Leśniak
7
Przykład. Niech (1) czerwony i (2) niebieski będą stana-
mi jakie przyjmuje pewna cząstka. Zmienia ona swój stan
zgodnie z macierzą przejścia
M =
1
3
2
3
3
4
1
4
.
W szczególności:
p
12
=
2
3
– prawdopodobieństwo zmiany
z czerwonego na niebieski,
p
21
=
3
4
– prawdopodobieństwo zmiany
z niebieskiego na czerwony.
(a) Wyznaczyć prawdopodobieństwo p
(2)
11
znalezienia się
czerwonej cząstki w stanie czerwoności po dwóch cyklach.
1
1
2
1
J
J
J
J
J
J
^
J
J
J
J
J
J
^
p
11
p
11
p
21
p
12
?
p
(2)
11
II cykl
I cykl
p
(2)
11
= p
11
· p
11
+ p
12
· p
21
=
1
3
·
1
3
+
2
3
·
3
4
=
11
18
1
3
= p
11
.
8
Macierze Lesliego i Markowa – K. Leśniak
(b) Wyznaczyć wszystkie prawdopodobieństwa p
(3)
ij
prze-
miany ze stanu i w stan j po trzech cyklach.
j
1
2
i
J
J
J
J
J
J
^
J
J
J
J
J
J
^
p
1j
p
i1
p
2j
p
i2
?
p
(2)
ij
II cykl
I cykl
p
(2)
ij
= p
i1
· p
1j
+ p
i2
· p
2j
=
2
X
k=1
p
ik
· p
kj
,
macierz przejścia po dwóch cyklach:
[p
(2)
ij
]
i,j
= M · M = M
2
.
Macierze Lesliego i Markowa – K. Leśniak
9
j
1
2
i
J
J
J
J
J
J
^
J
J
J
J
J
J
^
p
1j
p
(2)
i1
p
2j
p
(2)
i2
?
p
(3)
ij
III cykl
II cykl
p
(3)
ij
= p
(2)
i1
· p
1j
+ p
(2)
i2
· p
2j
=
2
X
k=1
p
(2)
ik
· p
kj
,
macierz przejścia po trzech cyklach:
[p
(3)
ij
]
i,j
= [p
(2)
ik
]
i,k
· [p
kj
]
k,j
= M
2
· M = M
3
.
10
Macierze Lesliego i Markowa – K. Leśniak
(c) Prawdopodobieństwo pozostawania czerwonym przez
n cykli wynosi
p
11
· p
11
· . . . · p
11
|
{z
}
n
razy
= p
n
11
.
Prawdopodobieństwo pozostawania niezmiennie czerwonym
wynosi
lim
n→∞
p
n
11
= lim
n→∞
1
3
n
= 0.
Co dzieje się z prawdopodobieństwem p
(n)
11
powrotu do
stanu czerwoności po n cyklach?
j
1
2
i
J
J
J
J
J
J
^
J
J
J
J
J
J
^
p
1j
p
(n−1)
i1
p
2j
p
(n−1)
i2
?
p
(n)
ij
cykl n
(n − 1)
cykl
p
(n)
ij
= p
(n−1)
i1
· p
1j
+ p
(n−1)
i2
· p
2j
=
2
X
k=1
p
(n−1)
ik
· p
kj
,
macierz przejścia po n cyklach:
[p
(n)
ij
]
i,j
= [p
(n−1)
ik
]
i,k
· [p
kj
]
k,j
= M
n−1
· M = M
n
.
Macierze Lesliego i Markowa – K. Leśniak
11
Aby znaleźć p
(n)
11
musimy wyznaczyć potęgę M
n
.
Oznaczmy: M
n
=
α
n
β
n
γ
n
δ
n
. Wówczas
α
n+1
β
n+1
γ
n+1
δ
n+1
=
α
n
β
n
γ
n
δ
n
·
α
1
β
1
γ
1
δ
1
,
skąd
α
n
+ β
n
= 1
(1)
α
n+1
=
1
3
α
n
+
3
4
β
n
(2)
(*)
β
n+1
=
2
3
α
n
+
1
4
β
n
α
1
=
1
3
, β
1
=
2
3
(3)
γ
n
+ δ
n
= 1
γ
n+1
=
1
3
γ
n
+
3
4
δ
n
δ
n+1
=
2
3
γ
n
+
1
4
δ
n
γ
1
=
3
4
, δ
1
=
1
4
.
Wystarczy rozwiązać podukład (*). Podstawiając (1) do
(2) dostajemy równanie (L-1):
α
n+1
=
1
3
α
n
+
3
4
(1 − α
n
) = −
5
12
α
n
+
3
4
,
którego rozwiązanie ogólnym jest
α
n
=
−
5
12
n
· α
0
+
1 −
−
5
12
n
1 −
−
5
12
·
3
4
.
12
Macierze Lesliego i Markowa – K. Leśniak
Uwzględniając warunek początkowy (3):
1
3
= α
1
=
−
5
12
α
0
+
3
4
⇔ α
0
= 1
dostajemy
α
n
=
−
5
12
n
+
1 −
−
5
12
n
1 −
−
5
12
·
3
4
.
Zatem
p
(n)
11
= α
n
−→
n→∞
3
4
1 +
5
12
=
9
17
.
Podobnie
p
(n)
12
= β
n
= 1 − α
n
−→
n→∞
8
17
.
Przyjmując za (α
n
, β
n
) odpowiednio ciągi (γ
n
, δ
n
) otrzymu-
jemy rozwiązanie ogólne tego samego kształtu: δ
n
= 1 − γ
n
,
γ
n
=
−
5
12
n
· γ
0
+
1 −
−
5
12
n
1 −
−
5
12
·
3
4
.
Po uwzględnieniu warunków początkowych
3
4
= γ
1
=
−
5
12
n
· γ
0
+
1 −
−
5
12
n
1 −
−
5
12
·
3
4
⇔ γ
0
= 0
Macierze Lesliego i Markowa – K. Leśniak
13
uzyskujemy
γ
n
=
1 −
−
5
12
n
1 −
−
5
12
·
3
4
.
Zatem jak poprzednio
p
(n)
21
= γ
n
−→
n→∞
3
4
1 +
5
12
=
9
17
,
p
(n)
22
= δ
n
−→
n→∞
8
17
.
W ten sposób otrzymaliśmy opis asymptotycznego zacho-
wania się macierzy przejścia po n cyklach:
M
n
−→
n→∞
9
17
8
17
9
17
8
17
.
Mówimy, że macierz M jest ergodyczna, a prawdopodo-
bieństwa graniczne
lim
n→∞
p
(n)
11
= lim
n→∞
p
(n)
21
=
9
17
,
lim
n→∞
p
(n)
12
= lim
n→∞
p
(n)
22
=
8
17
,
nazywane są prawdopodobieństwami ergodycznymi.
14
Macierze Lesliego i Markowa – K. Leśniak
(d) Niech π = [π
1
, π
2
] będzie początkowym rozkładem
prawdopodobienstwa: zanim zadziałał mechanizm zmian ko-
lorów cząstka znajdowała się w stanie (1) z prawdopodo-
bieństwem π
1
, a w stanie (2) z prawdopodobieństwem π
2
.
Wyznaczyć rozkład prawdopodobieństwa stanu cząstki π
(1)
=
[π
(1)
1
, π
(1)
2
] po jednym cyklu.
x
1
2
1
2
1
2
J
J
J
J
J
^
J
J
J
J
J
^
J
J
J
J
J
J
^
p
12
p
11
p
22
p
21
π
1
π
2
stan
po 1 cyklu
stan
początkowy
π
(1)
1
= π
1
· p
11
+ π
2
· p
21
,
π
(1)
2
= π
1
· p
12
+ π
2
· p
22
,
czyli
[π
(1)
1
, π
(1)
2
] = [π
1
, π
2
] ·
p
11
p
12
p
21
p
22
= [π
1
, π
2
] · M.
Macierze Lesliego i Markowa – K. Leśniak
15
(e) Przez indukcję rozkład prawdopodobieństwa π
(n)
=
[π
(n)
1
, π
(n)
2
] stanu cząstki po n cyklach opisuje zależność
[π
(n)
1
, π
(n)
2
] = [π
1
, π
2
] ·
p
11
p
12
p
21
p
22
n
= [π
1
, π
2
] · M
n
.
Uzasadnienie: ćwiczenie.
Niezależnie od rozkładu początkowego π
π
(n)
= π · M
n
−→
n→∞
π
∗
,
gdzie π
∗
= [π
∗
1
, π
∗
2
],
π
∗
j
= lim
n→∞
p
(n)
ij
=
9
17
, j = 1
8
17
, j = 2
niezależne od i = 1, 2 prawdopodobieństwa graniczne zna-
lezione w punkcie (c). Tym samym rozkład graniczny π
∗
wy-
stępuje we wszystkich wierszach macierzy granicznej lim
n→∞
M
n
.
16
Macierze Lesliego i Markowa – K. Leśniak
Przykład.
M =
0.2
0.4
0.4
0.25 0.25 0.5
0.(3) 0.(3) 0.(3)
.
W tabeli zbieramy dane z czterech kroków iteracji wyko-
nanych dla trzech rozkładów początkowych.
n π
(n)
π
(n)
π
(n)
0 [0, 0, 1]
[0.5, 0.25, 0.25]
[.667, .333, 0]
1 [.333, .333, .333] [.246, .346, .408] [.217, .350, .433]
2 [.261, .327, .411] [.272, .321, .407] [.275, .318, .406]
3 [.271, .323, .405] [.271, .325, .406] [.270, .325, .404]
4 [.270, .324, .405] [.271, .324, .406] [.270, .324, .406]
π
∗
= [0.(270), 0.(324), 0.(405)].
Macierze Lesliego i Markowa – K. Leśniak
17
Sam rozkład ergodyczny π
∗
łatwo znaleźć nie iterując ani
macierzy ani rozkładów, gdyż jest on stacjonarny.
Twierdzenie 1
(stacjonarność rozkładu ergodycznego)
Jeżeli łańcuch Markowa M jest ergodyczny tzn.
n-te rozkłady prawdopodobieństwa π
n
= π M
n
zbiegają
do pewnego rozkładu π
∗
niezależnie od rozkładu począt-
kowego π, to
π
∗
= π
∗
· M
jest jedynym rozkładem stacjonarnym.
Dowód. Stacjonarność:
π
(n)
= π · M
n
−→
n→∞
π
∗
,
π
∗
←−
n→∞
π
(n+1)
= π M
n+1
= (π M
n
) M −→
n→∞
π
∗
M.
Jedyność:
ˆ
π = ˆ
π M ⇒
ˆ
π = ˆ
π M = (ˆ
π M ) M = ˆ
π M
2
= . . . = ˆ
π M
n
−→
n→∞
π
∗
⇒ ˆ
π = π
∗
.
18
Macierze Lesliego i Markowa – K. Leśniak
Przykład. M =
1
3
2
3
3
4
1
4
, π = [π
1
, π
2
], π
i
> 0,
P
i
π
i
= 1.
[π
1
, π
2
] = [π
1
, π
2
] M =
1
3
π
1
+
3
4
π
2
,
2
3
π
1
+
1
4
π
2
⇒
π
1
=
9
17
,
π
2
= 1 − π
1
=
8
17
.
Jednak nawet jedyność rozkładu stacjonarnego nie gwa-
rantuje jego ergodyczności.
Przykład (łańcuch okresowy).
M =
0 1
1 0
, π = [π
1
, π
2
], π
i
> 0,
P
i
π
i
= 1.
[π
1
, π
2
] = [π
1
, π
2
] M = [π
2
, π
1
]
⇒
π
1
= π
2
=
1
2
.
Znaleziony rozkład nie jest ergodyczny, bo
M
n
=
1 0
0 1
, n – parzyste,
M,
n – nieparzyste.
Jak się przekonać, że znaleziony rozkład stacjonarny jest
ergodyczny?
Macierze Lesliego i Markowa – K. Leśniak
19
Jako że znajdowanie jawnej postaci n-tej potęgi macie-
rzy jest pracochłonne (nawet z użyciem postaci Jordana),
powstaje pytanie, czy nie można określić ergodyczności ma-
cierzy Markowa bezpośrednio na podstawie jej współczyn-
ników.
Takiego kryterium dostarcza
Twierdzenie 2 (ergodyczne dla łańcuchów)
Niech macierz Markowa M = [p
ij
]
i,j∈{1,...,s}
spełnia ∀
i,j
p
ij
> 0.
Wówczas
∀
j
∃ π
∗
j
> 0 ∀
i
lim
n→∞
p
(n)
ij
= π
∗
j
,
s
P
j=1
π
∗
j
= 1, π
∗
=
π
∗
j
s
j=1
= π
∗
M .
Przypomnijmy, że M
n
=
"
p
(n)
ij
#
i,j∈{1,...,s}
, a więc
M
n
−→
n→∞
[π
∗
, . . . , π
∗
|
{z
}
s
razy
]
T
,
π
(n)
= π M
n
−→
n→∞
π
∗
.
20
Macierze Lesliego i Markowa – K. Leśniak
Dowód. Ustalmy 0 < ε < min
i,j
p
ij
i przyjmijmy oznaczenia:
α
(n)
j
= min
i=1,...,s
p
(n)
ij
[!]
> 0,
β
(n)
j
= max
i=1,...,s
p
(n)
ij
> α
(n)
j
.
[!] Jeśli ∀
ij
p
ij
> 0, to ∀
n
∀
ij
p
(n)
ij
> 0
(ćwiczenie).
Oczywiście
[p
(n+1)
ij
]
i,j
= M
n+1
= M
n
· M = [p
il
]
i,l
· [p
(n)
lj
]
l,j
,
p
(n+1)
ij
=
X
l
p
il
· p
(n)
lj
.
Stąd
p
(n+1)
ij
=
X
l
p
il
− ε · p
(n)
jl
!
· p
(n)
lj
+ ε ·
X
l
p
(n)
jl
· p
(n)
lj
6
6 β
(n)
j
·
X
l
p
il
− ε · p
(n)
jl
!
+ ε · p
(2n)
jj
=
= β
(n)
j
·
X
l
p
il
|
{z
}
=1
− ε ·
X
l
p
(n)
jl
|
{z
}
=1
+ ε · p
(2n)
jj
,
czyli
β
(n+1)
j
= max
i
p
(n+1)
ij
6 β
(n)
j
· (1 − ε) + ε · p
(2n)
jj
.
Analogicznie
α
(n+1)
j
> α
(n)
j
· (1 − ε) + ε · p
(2n)
jj
(*).
Łącznie:
β
(n+1)
j
− α
(n+1)
j
6 (1 − ε) ·
β
(n)
j
− α
(n)
j
!
,
Macierze Lesliego i Markowa – K. Leśniak
21
skąd
β
(n)
j
− α
(n)
j
6 (1 − ε)
n−1
·
β
(1)
j
− α
(1)
j
!
−→
n→∞
0.
Przechodząc w (*) z ε → 0 widzimy, że
α
(n)
j
6 α
(n+1)
j
6 1,
tzn.
α
(n)
j
!
∞
n=1
– rosnący i ograniczony. Istnieją więc granice
π
∗
j
= lim
n→∞
α
(n)
j
,
przy czym π
∗
j
> 0, bo
lim
n→∞
α
(n)
j
> α
(n)
j
> α
(1)
j
= min
i
p
ij
> min
ij
p
ij
> ε > 0.
Ostatecznie
p
(n)
ij
− π
∗
j
6 β
(n)
j
− α
(n)
j
−→
n→∞
0,
a tym samym istnieją lim
n→∞
p
(n)
ij
= π
∗
j
> 0 oraz
X
j
π
∗
j
=
X
j
lim
n→∞
p
(n)
ij
= lim
n→∞
X
j
p
(n)
ij
= 1.
Uwaga: Wystarczy założyć, że dla pewnego n
0
potęga
M
n
0
ma wszystkie wyrazy dodatnie: p
(n
0
)
ij
> 0.