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 mi – liczba potomstwa pojawiającego się co jednost- kę czasu u osobnika z i-tej grupy wiekowej, i = 1, . . . , k,
[0, 1] 3 si – 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,
xn,i – liczba osobników z i-tej grupy wiekowej w n-tej jednostce czasu, i = 1, . . . , k.
Struktura wiekowa populacji w n-tej jednostce czasu xn = [xn,1, xn,2, . . . , xn,k]T
podlega zależnościom:
xn+1,1 = Xk
i=1mi · xn,i
– liczba nowych osobników (wszystkie grupy wiekowe mogą wydawać potomstwo),
xn+1,i+1 = si · xn,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: xn+1 = M · xn, gdzie
M =
m1, m2, . . . , mk−1, mk s1, 0, . . . , 0, 0
0, s2, . . . , 0, 0 . . . . . . . . . .
0, 0, . . . , sk−1, 0
.
Powyższy wzór rekurencyjny prowadzi do wzoru „ jawnego”:
xn = Mn · x0.
Przykład. W pewnej populacji wyróżniają się dwie grupy wiekowe:
1 – osobniki niedojrzałe, które nie mogą się rozmnażać (m1 = 0 o przeżywałności s (s1 = s), oraz
2 – osobniki dojrzałe o płodności m, które mogą się roz- mnażać (m2 = m).
Macierz Lesliego: M =
m1 m2 s1 0
=
0 m s 0
. Przez indukcję można zauważyć, że:
M2n+1 =
0 mn+1 · sn sn+1 · mn 0
, M2n =
mn · sn 0 0 sn · mn
. Zatem
x2n+1 = M2n+1 · x0 =
0 mn+1 · sn sn+1 · mn 0
·
x0,1 x0,2
=
=
mn+1 · sn · x0,2 sn+1 · mn · x0,1
= (m · s)n · [m · x0,2, s · x0,1]T, x2n = M2n· x0 =
mn · sn 0 0 sn · mn
·
x0,1 x0,2
= (m · s)n· x0. Zachowanie asymptotyczne:
xn =
xn,1 xn,2
n→∞−→
[00] , gdy 0 < ms < 1, [∞∞] , gdy ms > 1.
Przykład (cykliczne zmiany struktury wieku).
M =
0 0 4
1
3 0 0 0 34 0
. Łatwo dostrzec, że: M3 =
"1 0 0 0 1 0 0 0 1
#
. Zachowanie asymptotyczne:
x3n = x0 = [x01, x02, x03]T, x3n+1 = x1 = [x11, x12, x13]T =
4 x03, 1
3 x01, 3 4 x02
T
, x3n+2 = x2 = [x21, x22, x23]T =
3 x02, 4
3 x03, 1 4 x01
T
. Cykl powtarza się co trzy jednostki czasu:
x0, x1, x2, x0, x1, . . .
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
[xn,i/ P3i=1 xn,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ść.
Macierze Markowa
M =
p11 p12 p13 . . . p1s p21 p22 p23 . . . p2s . . . . ps1 ps2 ps3 . . . pss
= [pij]16i,j6s
— macierz stochastyczna (= macierz Markowa
= macierz przejścia), gdy
∀i,j 0 6 pij 6 1, ∀i Xs
j=1pij = 1.
Terminologia:
• {1, . . . , s} – zbiór stanów układu;
• pij – prawdopodobieństwo przejścia ze stanu i w stan j, i, j ∈ {1, . . . , s};
• M = [pij]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.
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:
p12 = 23 – prawdopodobieństwo zmiany z czerwonego na niebieski,
p21 = 34 – 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
JJ^
J J
J J
JJ^
p11 p11
p21 p12
?
p(2)11
II cykl I cykl
p(2)11 = p11 · p11 + p12 · p21 = 1 3 · 1
3 + 2 3 · 3
4 = 11
18 1
3 = p11.
(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
JJ^
J J
J J
JJ^
p1j pi1
p2j pi2
?
p(2)ij
II cykl I cykl
p(2)ij = pi1 · p1j + pi2 · p2j = X2
k=1pik · pkj, macierz przejścia po dwóch cyklach:
[p(2)ij ]i,j = M · M = M2.
j
1 2
i
J J
J J
JJ^
J J
J J
JJ^
p1j p(2)i1
p2j p(2)i2
?
p(3)ij
III cykl II cykl
p(3)ij = p(2)i1 · p1j + p(2)i2 · p2j = X2
k=1p(2)ik · pkj, macierz przejścia po trzech cyklach:
[p(3)ij ]i,j = [p(2)ik ]i,k · [pkj]k,j = M2 · M = M3.
(c) Prawdopodobieństwo pozostawania czerwonym przez n cykli wynosi
p11 · p11 · . . . · p11
| {z }
n razy
= p11n.
Prawdopodobieństwo pozostawania niezmiennie czerwonym wynosi
n→∞lim p11n = limn→∞
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
JJ^
J J
J J
JJ^
p1j p(n−1)i1
p2j p(n−1)i2
?
p(n)ij
cykl n (n − 1)cykl
p(n)ij = p(n−1)i1 · p1j + p(n−1)i2 · p2j = X2
k=1p(n−1)ik · pkj, macierz przejścia po n cyklach:
[p(n)ij ]i,j = [p(n−1)ik ]i,k · [pkj]k,j = Mn−1 · M = Mn.
Aby znaleźć p(n)11 musimy wyznaczyć potęgę Mn. Oznaczmy: Mn =
α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 = 13 αn + 34 βn (2)
(*) βn+1 = 23 αn + 14 βn
α1 = 13, β1 = 23 (3) γn + δn = 1
γn+1 = 13 γn + 34 δn δn+1 = 23 γn + 14 δn γ1 = 34, δ1 = 14.
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.
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
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:
Mn n→∞−→
9 17
8 17 9 17
8 17
.
Mówimy, że macierz M jest ergodyczna, a prawdopodo- bieństwa graniczne
n→∞lim p(n)11 = limn→∞p(n)21 = 179 ,
n→∞lim p(n)12 = limn→∞p(n)22 = 178 ,
nazywane są prawdopodobieństwami ergodycznymi.
(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), π2(1)] po jednym cyklu.
x
1 2 1 2
1 2
J J
J J J
^
J J
J J J
^
J J
J J
JJ^
p12
p11 p21 p22
π1 π2
stan
po 1 cyklu stan
początkowy
π1(1) = π1 · p11 + π2 · p21, π2(1) = π1 · p12 + π2 · p22, czyli
[π1(1), π2(1)] = [π1, π2] ·
p11 p12 p21 p22
= [π1, π2] · M.
(e) Przez indukcję rozkład prawdopodobieństwa π(n) = [π1(n), π2(n)] stanu cząstki po n cyklach opisuje zależność
[π1(n), π2(n)] = [π1, π2] ·
p11 p12 p21 p22
n
= [π1, π2] · Mn. Uzasadnienie: ćwiczenie.
Niezależnie od rozkładu początkowego π π(n) = π · Mn n→∞−→ π∗, gdzie π∗ = [π1∗, π2∗],
πj∗ = limn→∞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 limn→∞Mn.
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)].
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 = π Mn zbiegają do pewnego rozkładu π∗ niezależnie od rozkładu począt- kowego π, to
π∗ = π∗ · M
jest jedynym rozkładem stacjonarnym.
Dowód. Stacjonarność:
π(n) = π · Mn n→∞−→ π∗,
π∗ n→∞←− π(n+1) = π Mn+1 = (π Mn) M −→n→∞ π∗M.
Jedyność:
π = ˆˆ π M ⇒
π = ˆˆ π M = (ˆπ M ) M = ˆπ M2 = . . . = ˆπ Mn n→∞−→ π∗
⇒ ˆπ = π∗.
Przykład. M =
1 3
2 3 3 4
1 4
, π = [π1, π2], πi > 0, Pi π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, Pi πi = 1.
[π1, π2] = [π1, π2] M = [π2, π1]
⇒ π1 = π2 = 1 2. Znaleziony rozkład nie jest ergodyczny, bo
Mn =
1 0 0 1
, n – parzyste, M, n – nieparzyste.
Jak się przekonać, że znaleziony rozkład stacjonarny jest ergodyczny?
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 = [pij]i,j∈{1,...,s}
spełnia ∀i,j pij > 0.
Wówczas
∀j ∃ πj∗ > 0 ∀i n→∞lim p(n)ij = πj∗,
Ps
j=1πj∗ = 1, π∗ = πj∗s
j=1 = π∗ M . Przypomnijmy, że Mn =
"
p(n)ij
#
i,j∈{1,...,s}, a więc Mn n→∞−→ [π∗, . . . , π∗
| {z }
s razy
]T, π(n) = π Mn n→∞−→ π∗.
Dowód. Ustalmy 0 < ε < min
i,j pij i przyjmijmy oznaczenia:
α(n)j = min
i=1,...,sp(n)ij > 0,[!]
βj(n) = max
i=1,...,sp(n)ij > α(n)j .
[!] Jeśli ∀ij pij > 0, to ∀n ∀ij p(n)ij > 0 (ćwiczenie).
Oczywiście
[p(n+1)ij ]i,j = Mn+1 = Mn · M = [pil]i,l · [p(n)lj ]l,j, p(n+1)ij = X
l pil · p(n)lj . Stąd
p(n+1)ij = X
l pil − ε · p(n)jl ! · p(n)lj + ε · X
l p(n)jl · p(n)lj 6 6 βj(n) · X
l
pil − ε · p(n)jl ! + ε · p(2n)jj =
= βj(n) ·
X
l pil
| {z }
=1
− ε · X
l p(n)jl
| {z }
=1
+ ε · p(2n)jj , czyli
βj(n+1) = max
i p(n+1)ij 6 βj(n) · (1 − ε) + ε · p(2n)jj . Analogicznie
α(n+1)j > α(n)j · (1 − ε) + ε · p(2n)jj (*).
Łącznie:
βj(n+1) − αj(n+1) 6 (1 − ε) · βj(n) − α(n)j !,
skąd
βj(n) − α(n)j 6 (1 − ε)n−1 · βj(1) − α(1)j ! n→∞−→ 0.
Przechodząc w (*) z ε → 0 widzimy, że α(n)j 6 α(n+1)j 6 1, tzn. αj(n)
!∞
n=1 – rosnący i ograniczony. Istnieją więc granice πj∗ = limn→∞α(n)j ,
przy czym πj∗ > 0, bo
n→∞lim α(n)j > α(n)j > α(1)j = min
i pij > minij pij > ε > 0.
Ostatecznie
p(n)ij − πj∗
6 βj(n) − α(n)j n→∞−→ 0, a tym samym istnieją lim
n→∞p(n)ij = πj∗ > 0 oraz
X
j πj∗ = X
j
n→∞lim p(n)ij = limn→∞ X
j p(n)ij = 1.
Uwaga: Wystarczy założyć, że dla pewnego n0 potęga Mn0 ma wszystkie wyrazy dodatnie: p(nij0) > 0.
