Анализ на модела на Лесли

Феноменология

Моделът на Лесли описва и същевременно прогнозира състоянието на една популация. За тази цел популацията се разделя на възрастови групи, всяка от които се характеризира със специфична раждаемост и преживяемост, като последните се разглеждат във вероятностен смисъл.

Означения (описание на състоянието)

m+1 – брой възрастови групи
(Subscript[n, i])(t) – брой индивиди в групата i или векторно: n(t)=(((({{(Subscript[n, 0])(t), (Subscript[n, 1])(t), (Subscript[n, 2])(t), ⋯, (Subscript[n, m-1])(t), (Subscript[n, m])(t)}})))^T)
(Subscript[f, i]) – коефициент на фертилност на групата i
(Subscript[p, i]) – коефициент на преживяемост на групата i

Постановка (прогноза)

След като сме въвели необходимите означения, т.е. сме преценили конкретно по какъв начин моделът ще описва състоянието на популацията, следващата стъпка е съставянето на прогноза. За тази цел трябва да определим как зависи моделът в един следващ времеви момент от състоянието в сегашния.

Естествено е да приемем следните две феноменологични изисквания.

(Subscript[n, 0])(t+1)=(∑_(i=0)^m)(Subscript[f, i])(Subscript[n, i])(t)

(Subscript[n, i])(t+1)=(Subscript[p, i-1])(Subscript[n, i-1])(t),i=1,...,m

Т.е. (1) всяка група допринася с броя и фертилността си към броя новородени в следващия момент и (2) всяка група преминава в следващата със съответна преживяемост.

Съставяне на матрицата

С въвеждането на (1) и (2) моделът сам по себе си е завършен. Обаче тъй като (1) и (2) представляват линейна (дори хомогенна) система уравнения, бихме могли да представим уравненията в матричен вид. Това ще ни даде възможност за далеч по-удобен алгебричен анализ. Всъщност това е естествената (и очакваната) стъпка в изготвянето на един линеен модел.

Представяме системата (1), (2) във вида:

n(t+1)=A n(t)

Тук A=(Subscript[{(Subscript[a, ij])},  m + 1 × m + 1]). Това представяне още се нарича рекурсивна формула на решението.

Сега да заместим.

(({{(Subscript[n, 0])(t+1)}, {(Subscript[n, 1])(t+1)}, {(Subscript[n, 2])(t+1)}, {⋮}, {(Subscript[n, m-1])(t+1)}, {(Subscript[n, m])(t+1)}}))=A(({{(Subscript[n, 0])(t)}, {(Subscript[n, 1])(t)}, {(Subscript[n, 2])(t)}, {⋮}, {(Subscript[n, m-1])(t)}, {(Subscript[n, m])(t)}}))(Overscript[⇔, заместваме])(({{(∑_(i=0)^m)(Subscript[f, i])(Subscript[n, i])(t)}, {(Subscript[p, 0])(Subscript[n, 0])(t)}, {(Subscript[p, 1])(Subscript[n, 1])(t)}, {⋮}, {(Subscript[p, m-2])(Subscript[n, m-2])(t)}, {(Subscript[p, m-1])(Subscript[n, m-1])(t)}}))=A(({{(Subscript[n, 0])(t)}, {(Subscript[n, 1])(t)}, {(Subscript[n, 2])(t)}, {⋮}, {(Subscript[n, m-1])(t)}, {(Subscript[n, m])(t)}}))

Остава да определим елементите на тази матрица. Умножаваме ред по стълб.

(({{(Subscript[f, 0])(Subscript[n, 0])(t)+⋯+(Subscript[f, m])(Subscript[n, m])(t)}, {(Subscript[p, 0])(Subscript[n, 0])(t)}, {(Subscript[p, 1])(Subscript[n, 1])(t)}, {⋮}, {(Subscript[p, m-2])(Subscript[n, m-2])(t)}, {(Subscript[p, m-1])(Subscript[n, m-1])(t)}}))=(({{(Subscript[a, 0,0])(Subscript[n, 0])(t)+(Subscript[a, 0,1])(Subscript[n, 1])(t)+⋯+(Subscript[a, 0,m-1])(Subscript[n, m-1])(t)+(Subscript[a, 0,m])(Subscript[n, m])(t)}, {(Subscript[a, 1,0])(Subscript[n, 0])(t)+(Subscript[a, 1,1])(Subscript[n, 1])(t)+⋯+(Subscript[a, 1,m-1])(Subscript[n, m-1])(t)+(Subscript[a, 1,m])(Subscript[n, m])(t)}, {(Subscript[a, 2,0])(Subscript[n, 0])(t)+(Subscript[a, 2,1])(Subscript[n, 1])(t)+⋯+(Subscript[a, 2,m-1])(Subscript[n, m-1])(t)+(Subscript[a, 2,m])(Subscript[n, m])(t)}, {⋮}, {(Subscript[a, m-1,0])(Subscript[n, 0])(t)+(Subscript[a, m-1,1])(Subscript[n, 1])(t)+⋯+(Subscript[a, m-1,m-1])(Subscript[n, m-1])(t)+(Subscript[a, m-1,m])(Subscript[n, m])(t)}, {(Subscript[a, m,0])(Subscript[n, 0])(t)+(Subscript[a, m,1])(Subscript[n, 1])(t)+⋯+(Subscript[a, m,m-1])(Subscript[n, m-1])(t)+(Subscript[a, m ,m])(Subscript[n, m])(t)}}))

Приравняваме съответните елементи отляво и отдясно, т.е. получаваме система. Тогава прирявняваме и коефициентите пред променливите (Subscript[n, i])(t). Тогава получаваме:

(Subscript[f, 0])=(Subscript[a, 0,0]);(Subscript[f, 1])=(Subscript[a, 0,1]);⋯ (Subscript[f, m-1])=(Subscript[a, 0,m-1]);(Subscript[f, m])=(Subscript[a, 0,m])
(Subscript[p, 0])=(Subscript[a, 1,0]);0=(Subscript[a, 1,1]);⋯ 0=(Subscript[a, 1,m-1]);0=(Subscript[a, 1,m])
0=(Subscript[a, 2,0;])(Subscript[p, 1])=(Subscript[a, 2,1]);0=(Subscript[a, 2,2]);⋯ 0=(Subscript[a, 2,m-1]);0=(Subscript[a, 2,m])
⋮
0=(Subscript[a, m-1,0]);⋯ 0=(Subscript[a, m-1,m-3]);(Subscript[p, m-2])=(Subscript[a, m-1,m-2]);0=(Subscript[a, m-1,m-1]);0=(Subscript[a, m-1,m])
0=(Subscript[a, m,0]);⋯ 0=(Subscript[a, m,m-2]);(Subscript[p, m-1])=(Subscript[a, m,m-1]);0=(Subscript[a, m,m])

Т.е. получаваме че A има вида:

(({{(Subscript[f, 0]), (Subscript[f, 1]), (Subscript[f, 2]), ⋯, (Subscript[f, m-2]), (Subscript[f, m-1]), (Subscript[f, m])}, {(Subscript[p, 0]), 0, 0, ⋯, 0, 0, 0}, {0, (Subscript[p, 1]), 0, ⋯, 0, 0, 0}, {0, 0, (Subscript[p, 2]), ⋯, 0, 0, 0}, {⋮, ⋮, ⋮, ⋱, ⋮, ⋮, ⋮}, {0, 0, 0, ⋯, (Subscript[p, m-1]), 0, 0}, {0, 0, 0, ⋯, 0, (Subscript[p, m]), 0}}))

Общо решение

Превид рекурентната зависимост (3) имаме:

n(1)=A n(0)⇔n(2)=A(A n(0))⇔n(3)=A(A(A n(0)))

Т.е. получаваме:

n(t)=(A^t)n(0)

Резултатът (5) се нарича формула за общото решение, а n(0) – начални данни, които в нашия случай представляват началното разпределение на популацията.

Съставяне на матрицата в Mathematica

Да изчистим глобалния контекст от евентуално съществуващи предишни променливи.

Remove["Global`*"]

Нека съставим функция даваща матрицата на Лесли от произволна размерност.

Обърнете внимание, че индексите, които ще използваме са с единица по-големи, тъй като нулевият индекс е запазен за главата на израза, т.е. ще пишем (Subscript[f, 1]) вместо (Subscript[f, 0]) и т.н.

Трябват ни два цикъла, затова и ще полваме два пъти Table. След това проверяваме дали редът (i) не е първият. Ако е, ще нанасяме параметрите на фертилност по колони (j). Ако не е първият ред, трябва да нанесем параметрите на преживяемост точно под главния диагонал, т.е. там където индекса на реда е с единица по-голям от този на колоната (i=j+1), а на останалите места да оставим нули.

ClearAll@leslie
leslie[m_]:=Table[
If[i==1,f[[j]],
If[i==j+1,p[[j]],0]
]
,{j,1,m},{i,1,m}
]

Да пробваме.

leslie[7]//MatrixForm

Работи. Тъй като Mathematica не знае какво е f и p, а по подразбиране счита символите за скалари (атоми по-точно), се оплаква от получения резултат.

Нека заглушим съобщенията partd сQuiet.

ClearAll@leslie
leslie[m_]:=Quiet[
Table[
If[i==1,
f[[j]],
If[i==j+1,
p[[j]],
0]
]
,{j,1,m},{i,1,m}
]
,Part::partd
]

Нека вземем два примерни f1 и p1.

f1=RandomReal[{0.1,0.9},7];
p1=RandomReal[{0.1,0.9},6];

Достатъчно е да заместим f и p в получения резултат за матрицата.

leslie[7]/.{f->f1,p->p1}//MatrixForm

Наистина става. Това е удобен начин да видим вида на матрицата за произволна размерност без да задаваме стойности на параметрите. Да ползваме MatrixPlot, за да покаже визуално Mathematica как изглежда матрицата.

MatrixPlot@%

Такъв един начин на получаване на матрицата, макар директен, е непрактичен, неефикасен и крие възможности за грешка. По-добре е параметрите да се предават като променливи. Да редактираме дефиницията за leslie. Единствено променяме списъка с параметри и заменяме m с Length[f].

leslie[f_,p_]:=Table[
If[i==1,f[[j]],
If[i==j+1,p[[j]],0
]
]
,{j,1,Length[f]},{i,1,Length[f]}
]

Да пробваме.

leslie[f1,p1]//MatrixForm

При това така построена програмата работи и в случаите, когато p има броя на елементи на f (а не с един по-малко).

leslie[f1,p1~Join~{0}]//MatrixForm

Използваме едно име, но функциите имат различни сигнатури и затова Mathematica може да ги разграничава. Въпреки че едно такова обяснение е удобно за един конвенционален програмист, истината е, че става дума за „пасване“ на шаблони и изрази.

Да видим сега какво стои зад символа leslie.

?leslie

Собствени стойности

Нека вземем за пример матрицата на Лесли от 3-та размерност.

leslie[3]//MatrixForm

Да помолим Mathematica да ни даде собствените стойности.

Eigenvalues@%

В първия момент резултатът може да произведе странни реакции у читателите, но всъщност Mathematica просто връща несметнати стойностите на първия, втория и третия корен на характеристичния полином, т.е. заявява, че такива не може да намери. В общия случай не можем да намерим корени на полином от произволна степен. Има формули за намиране на корените на полином от 3-та, 4-та степен и 5-та степен (вж. [1]), но прилагайки ги към разглеждания модел, не получаваме корените в удобна за съставяне на изводи форма. Които иска може да пробва. Освен това Solve също може да бъде полезен.

Какво ни остава?
Първо – анализ на характеристичния полином (който ще забележим, че има специфичен вид).
Второ – дали няма да можем да намерим собствените вектори?
Трето –  не бива да забравяме и стойността на експерименталните резултати.

Характеристичен полином

Нека чрез CharacteristicPolynomial да видим формата на характеристичните полиноми за няколко размерности.

{
{cp[1]},
{cp[2]},
{cp[3]},
{cp[4]}
}/.{cp[k_]->CharacteristicPolynomial[leslie[k],λ]}//Grid//Quiet

Горното записване не е най-коретният възможен начин, а по-скоро цели да покаже абстракността на Mathematica. По-коректно е да ползваме Table:

Table[{CharacteristicPolynomial[leslie[k],λ]},{k,4}]//Grid//Quiet

Забелязва се, че полиномът има следният вид:

(((-λ))^n)+(((-1))^(n-1))(∑_(i=0)^(m-1))(f[k+1](∏_(j=1)^i)p[j])

Според теоремата на Перон-Фробениус (вж. [2])за всяка реална неразложима квадратна матрица със неотрицателни стойности съществува положителен корен (Subscript[λ, 1]), такъв че:

|(Subscript[λ, k])|≤(Subscript[λ, 1]), където k≠1

Т.е. полиномът има еднократен положителен доминиращ корен. Нещо повече:

|(Subscript[λ, k])|<(Subscript[λ, 1]),където k≠1, ако (Subscript[f, i])≠0∧(Subscript[f, i+1])≠0 за някое i=1,...,m-1

Т.е. ако има две поредни възрастови групи, които са фертилни – почти винаги изпълнено за една реална популация, то (Subscript[λ, 1]) е стриктно доминиращ (т.е. неравенството е строго).

По въпроса дали матрицата е неразложима ще говорим по-долу.

Дали можем да оценим в какви граници са възможните стойности на този корен. Отиваме на адрес http://functions.wolfram.com, търсим информация за root и забелязваме, че има на разположение неравенство, оценяващо граничните стойности на дължината на корените (т.е. на модула, а не самите корени) на полином, използвайки неговите коефициенти. Всъщност това е следствие от теоремата на Руше (вж. [3]).

|(Subscript[z, k])|≤1+max(|((Subscript[a, 0])/(Subscript[a, n]))|,|((Subscript[a, 1])/(Subscript[a, n]))|,⋯,|((Subscript[a, n-1])/(Subscript[a, n]))|)

Където (Subscript[a, i]) е коефециентът пред i-тата степен на z.

Т.е. получаваме, че за водещата собствена стойност (Subscript[λ, 1]) (както отбелязахме (Subscript[λ, 1]) е неотрицателна) на матрицата имаме:

0≤(Subscript[λ, 1])≤2

Собствени вектори

Въпреки, че не знаем точната стойност на собствените стойности, намирането на собствените вектори не е невъзможно. Обаче това тук няма да правим.

Асимптотично решение

Взимайки предвид важното свойство (8), бихме могли да запишем формулата за общото решение (5) и рекурентния му вид (3) в асимптотична форма.

Предвид че собствените стойности за матрицата на Лесли са с кратност единица (може да се покаже), то за матрицата на Лесли A имаме познатото разлагане:

A=X D (X^(-1)),

където D е диагонална матрица. Имаме още:

(A^t)=(({{(Subscript[x, 1]), (Subscript[x, 2]), ⋯, (Subscript[x, m])}}))(({{(Subscript[λ, 1]^t), 0 , ⋯, 0 }, {0, (Subscript[λ, 2]^t),  , ⋮}, {⋮,  , ⋱, 0 }, {0, …, 0, (Subscript[λ, m]^t)}}))(({{(Subscript[(OverBar[x]), 1])}, {(Subscript[(OverBar[x]), 2])}, {⋮}, {(Subscript[(OverBar[x]), m])}})),

където (Subscript[x, 1]),...,(Subscript[x, m]) са собствените вектори на A във вид на вектор-стълбове, съставляващи матрицата X; (Subscript[(OverBar[x]), 1]),...,(Subscript[(OverBar[x]), m]) са вектор-редове, съставляващи (X^(-1)), λ са собствените стойности на A. Тогава може да се провери и докаже, че за (11) още е в сила:

(A^t)=(∑_(k=1)^m)(Subscript[λ, k]^t)(Subscript[x, k])(Subscript[(OverBar[x]), k])

Предвид свойство (8) има смисъл изнасянето на (Subscript[λ, 1]) пред скоби:

(A^t)=(∑_(k=1)^m)(Subscript[λ, k]^t)(Subscript[x, k])(Subscript[(OverBar[x]), k])=(Subscript[λ, 1]^t)(Subscript[x, 1])(Subscript[(OverBar[x]), 1])+(Subscript[λ, 1]^t)(∑_(k=2)^m)((Subscript[λ, k]^t)/(Subscript[λ, 1]^t))(Subscript[x, k])(Subscript[(OverBar[x]), k])

В общия случай (Subscript[λ, k]) могат да бъдат комплексни. Тогава по формулата на Моавър (не забравяме, че (Subscript[λ, 1]) е реално) имаме:

(A^t)=(Subscript[λ, 1]^t)(Subscript[x, 1])(Subscript[(OverBar[x]), 1])+(Subscript[λ, 1]^t)(∑_(k=2)^m)(((((|(Subscript[λ, k])|)/(Subscript[λ, 1]))))^t)(cos(t ϕ)+i sin(t ϕ))(Subscript[x, k])(Subscript[(OverBar[x]), k])

Но от (8) очевидно следва:

((|(Subscript[λ, k])|)/(Subscript[λ, 1]))<1

Тогава след граничен преход при t -> ∞ и предвид факта, че cos(t ϕ)+i sin(t ϕ) е ограничено (в комплексен смисъл) имаме:

(A^t)(Underscript[⟶, t->∞])(Subscript[λ, 1]^t)(Subscript[x, 1])(Subscript[(OverBar[x]), 1])

Вече намерили граничната стойност на матрицата A, можем да изведем такива стойности и за векторите на състоянието n.

n(t)=(A^t)n(0)(Underscript[⟶, t->∞])(Subscript[λ, 1]^t)(Subscript[x, 1])(Subscript[(OverBar[x]), 1])n(0)=(Subscript[λ, 1])(Subscript[λ, 1]^(t-1))(Subscript[x, 1])(Subscript[(OverBar[x]), 1])n(0)(Underscript[←, t->∞])n(t-1)

Т.е. за рекурентната връзка (3) получихме асимптотичното решение:

n(t)(Underscript[⟶, t->∞])λ n(t-1)

Тогава (Subscript[λ, 1]) може да се разглежда като асимптотична скорост на модела (в случая – на възпроизводство на популацията), при това:

при (Subscript[λ, 1])<1 популацията ще загине след достатъчно дълго време

при (Subscript[λ, 1])>1 популацията ще се пренасели (броят ще расте неограничено)

при (Subscript[λ, 1])=1 броят на индивидите в популацията ще е в равновесие.

От (17) получихме асимптотичния вид на общото решение

n(t)(Underscript[⟶, t->∞])(Subscript[λ, 1]^t)(Subscript[x, 1])(Subscript[c, 1]), където (Subscript[c, 1])=(Subscript[(OverBar[x]), 1])n(0), т.е. (Subscript[c, 1])е число

Оттук и изводите:

за достатъчно големи t асимптотичното решение (19) може да се разглежда като приближено решение

за достатъчно големи t векторът на състоянието n ще бъде точно правопропорционален на (Subscript[x, 1]). Това означава, че пропорциите между отделните групи се запазват във времето, т.е.:

(((Subscript[n, i])((Subscript[t, 1])))/((Subscript[n, k])((Subscript[t, 1]))))=(((Subscript[n, i])((Subscript[t, 2])))/((Subscript[n, k])((Subscript[t, 2]))))

Този факт се нарича устойчиво възрастово разпределение.

Life-Cycle диаграма

Да направим и LifeCycle диаграма по дадени f и p. За тази цел се използва GraphPlot. Като аргумент се подава списък с дъгите от графа във вида

{...,{3->4,"дъга 4"},{3->1,"дъга 5"},...}

Нека напишем функция. Използваме Module, за да може част от променливите да останат локални. Да забележим, че имаме два основни вида дъги:

Дъги, отговарящи на коефициентите на преживяемост, които свързват всяка група със следващата.

Дъги, отговарящи на коефициентите на фертилност, които свързват всички групи с първата.

Също така съм задал някои опции за GraphPlot – големина на изображението, насоченост, големина на стрелките, надписване на дъгите и метод на чертане.

ClearAll@leslieDiagram
leslieDiagram[f_,p_]:=Module[
{
edgeList={},
n=Length[f] (* брой възрастови групи *)
},
(* въвеждаме връзките на преживяемост между групите *)
(* U1 може да ползваме и стрингове за да надпишем стойностите *)
edgeList=edgeList~Join~Table[{i->i+1,p[[i]]},{i,1,n-1}];
(* въвеждаме връзките на фертилност на групите *)
(* U2 може да ползваме и стрингове за да надпишем стойностите *)
edgeList=edgeList~Join~Table[{i->1,f[[i]]},{i,1,n}];
(* Чертаем ориентиран граф чрез GraphPlot, използвайки информацията в edgeList *)
GraphPlot[edgeList,ImageSize->400,DirectedEdges->{True,"ArrowheadsSize"->0.03},VertexLabeling->True,Method->"RadialDrawing"]
]

Да пробваме.

leslieDiagram[f1,p1]

Сега нека си припомним, че остана да докажем неразложимост на матрицата. Една квадратна n×n матрица е разложима, ако чрез елементарни преобразувания можем да я трансформираме в блочна горнотриъгълна форма: (({{(Subscript[A, 1]), (Subscript[A, 12])}, {0, (Subscript[A, 2])}})), където (Subscript[A, 1]),(Subscript[A, 2]) са квадратни подматрици (вж. [4]). Друг удобен критерий е следният: Една матрицата е разложима, тогава и само тогава, когато е матрица на съседство на граф, който не е силно свързан (вж. [5]). Оттук следва, че една матрица е неразложима, ако има пътища (в двете посоки) между всички точки на съответният й граф. Разглеждайки последната диаграма лесно бихме се убедили в този факт.

Източници

http://functions.wolfram.com/ElementaryFunctions/Root/03/01/

http://en.wikipedia.org/wiki/Perron-Frobenius_theorem

http://en.wikipedia.org/wiki/Properties_of_polynomial_roots

http://homepage.newschool.edu/het/essays/math/matrix.htm

http://mathworld.wolfram.com/ReducibleMatrix.html