faculteit Wiskunde en Natuurwetenschappen
Modelleren van de zang van de geelgors
Gebruikmakend van Markovketens
Bacheloronderzoek Wiskunde
Juli 2013 Student: F.H. Dijkstra Eerste Begeleider: dr. E.C. Wit Tweede Begeleider: dr. A.C.D. van Enter
Bron van het plaatje op de voorkant: http://www.duikonder.nl/wp-content /uploads/2011/06/Geelgors 20090512 3320.jpg, (06-2013)
Samenvatting
De zang van de geelgors kan worden gezien als een opeenvolging van in- dividuele noten. Op deze manier kan de zang van de geelgors beschreven worden door een Markovproces. Met behulp van een Markovproces is het mogelijk om een voorspelling te doen wat de geelgors gaat fluiten. In deze scriptie wordt de zang van de geelgors met verschillende modellen gesimu- leerd door middel van Markovketens. Vervolgens worden de verschillende modellen getoetst en met elkaar vergeleken.
Inhoudsopgave
Inleiding 6
1 Markovketens 8
1.1 Markovprocessen . . . 9
1.2 Discrete Markovketens . . . 10
1.3 Classificatie van toestanden . . . 16
2 Goodness of fit 19 2.1 Maximum-likelihood-schatter . . . 19
2.2 Likelihood ratio . . . 22
2.3 Statistiche toetsen . . . 23
3 Het model 24 3.1 De zang van de geelgors . . . 24
3.1.1 Hoe vogels geluid maken . . . 24
3.1.2 Het typerende geluid van de geelgors . . . 25
3.1.3 De vijfde symfonie van Beethoven . . . 28
3.2 Data . . . 29
3.2.1 Verwerken van de data . . . 29
3.2.2 Het maken van modellen op basis van de data . . . 31
3.3 Uitkomsten . . . 33
3.3.1 Statistiche toetsen . . . 33
3.4 Verbetering van het model . . . 35
3.4.1 Ingebedde Markovketen . . . 35
3.4.2 Markovketens met variabele lengte . . . 38
Conclusie 40
Bibliografie 41
A Matlab codes 43
B R codes 44
B.1 Nulde orde Markovketens . . . 44
B.2 Codes voor eerste orde Markovketens . . . 45
B.3 Tweede orde Markovketens . . . 46
B.4 Derde orde Markovketens . . . 50
B.5 Histogrammen voor de verschillende zangtypes . . . 73
B.6 Ingebedde Markovketens . . . 82
B.7 Markovketen met variabele lengte . . . 122
Inleiding
“Bird-watching is either the most scientific of sports or the most sporting of sciences” (E.M. Nicholson, 1931)[14]
Wetenschappers worden al jaren ge¨ıntrigeerd door vogels. Hierdoor is er in het verleden al veel onderzoek gedaan naar vogels, maar ook tegenwoor- dig ontdekken wetenschappers nog steeds nieuwe onderzoeksgebieden met betrekking tot vogels. Zoals de distributie van kolonies keizerspingu¨ıns op Antarctica, het collectieve gedrag van spreeuwen in een spreeuwenzwerm en de biodiversiteit van vogels.
Omdat ik al van jongs af aan ge¨ınteresseerd ben in vogels leek het mij mooi dit onderwerp bij mijn scriptie te betrekken. Deze scriptie gaat over de zang van de geelgors. De zang van dit vogeltje is zo typerend dat men niet een groot vogelaar hoeft te wezen om haar geluid te herkennen. De zang is in zo’n grote mate typerend dat ook Beethoven hierdoor ge¨ınspireerd zou zijn voor de vijfde symfonie van Beethoven1.
De zang van de geelgors kan worden gezien als een opeenvolging van in- dividuele noten. Op deze manier kan de zang van de geelgors beschreven worden door een Markovproces. Met behulp van zo’n Markovproces is het mogelijk om een voorspelling te doen over wat de geelgors gaat fluiten.
Het doel van deze scriptie is om de zang van de geelgors met verschillende modellen te simuleren en te toetsen welk model we het beste kunnen kiezen om de zang te simuleren. Ook willen we dat het typische geluid van de geelgors terug te vinden is in de gesimuleerde liedjes. De wiskundige theorie die achter het maken van een model ligt is een Markovketen, deze theorie wordt in hoofdstuk 1 uitgelegd.
1We komen hier in hoofdstuk 3 op terug.
In het tweede hoofdstuk wordt uitgelegd hoe je een likelihood van een Markovketen kan berekenen en hoe je daarmee een model kan toetsen.
In het derde hoofdstuk wordt eerst enige achtergrond informatie gegeven omtrent de zang van de geelgors. Vervolgens wordt besproken hoe de data verkregen is en hoe hieruit de zang van de geelgors gesimuleerd kan worden.
Tot slot van dit hoofdstuk worden de verschillende modellen getest en be- sproken welk model het beste de zang van een geelgors kan simuleren.
In de conclusie wordt uiteindelijk geformuleerd dat het beste model een tweede orde ingebedde Markovketen van variabele lengte blijkt te zijn. Als met dit model liedjes worden gesimuleerd is het typische geluid van de geel- gors terug te vinden.
Aan het eind zijn acht appendices toegevoegd. De eerste appendix bevat een Matlab code waarin de limietverdeling van een Markovketen wordt be- rekend. De tweede tot en met achtste appendix bevatten codes waarmee de zang van een geelgors gesimuleerd en getest kan worden. Dit is gedaan door verschillende Markovordes te gebruiken. In de scriptie wordt aangegeven wanneer voor een bepaalde afbeelding of matrix gebruik gemaakt is van een appendix.
Hoofdstuk 1
Markovketens
In dit hoofdstuk zullen we de Markovketens bespreken. Er zullen definities, theorie¨en en voorbeelden worden behandeld. Om zo de kennis te verkrijgen die nodig is om het vervolg van de scriptie te begrijpen.
Figuur 1.1: Andrei Andreevich Markov (1856-1922) [15].
Andrei Andreevich Markov was een Russisch wiskundige. Hij is bekend door zijn werk in de getaltheorie, analyse en probabilieteitstheorie. Naar hem zijn de Markov-ongelijkheid, de Markovprocessen en de Markovketens vernoemd. Hij heeft grote bijdragen geleverd voor het nieuwe onderzoeks- gebied Markovketens, waar we in dit hoofdstuk verder op in zullen gaan [2].
1.1 Markovprocessen
Het is vaak mogelijk om het gedrag van een systeem wiskundig of fysisch te beschrijven. Een systeem kan dan gemodelleerd worden door de verschil- lende toestanden van het systeem te beschrijven en de overgangen hiertussen.
Bij een Markovproces is de volgende toestand van het systeem alleen afhan- kelijk van de huidige toestand van het model en niet van de geschiedenis.
Dit wordt ook wel de Markoveigenschap genoemd.
Voorbeeld 1.1
Beschouw de vogelzang van een geelgors1 die verschillende noten fluit. De verschillende noten representeren de verschillende toestanden van het sys- teem. Laten we aannemen dat welke noot de geelgors als volgende fluit alleen afhankelijk is van de noot die de geelgors op dat moment fluit2. Als we een model van dit voorbeeld willen maken dan zijn we ge¨ınteresseerd wat de kans is naar een volgende noot met als gegeven de net gefloten noot. Als we de kansen van en naar elke noot weten, dan kunnen we een voorspelling doen welke noten de geelgors gaat fluiten.
Een Markovproces is een stochastisch proces met stochastische variabelen.
Laat { X(t), t ∈ T } zo’n stochastische variabele zijn. Waarbij de toestands- ruimte S = {1, ..., i, ..., j, ...} de verschillende toestanden aangeeft die X(t) kan aannemen. De toestandsruimte S kan zowel continu zijn als discreet. Als een toestandsruimte discreet is dan wordt een Markovproces een Markovke- ten genoemd. Voor het verdere verloop zijn we alleen ge¨ınteresseerd in Markovketens en zullen we Markovprocessen achterwege laten. Waarbij de stochastische variabele verandert in: { Xn, n = 0, 1, 2... } en de toestands- ruimte in: S = {x1, ..., xi, ..., xm} met m verschillende toestanden die Xn kan aannemen.
1Dit is een zangvogel die onder andere in Nederland voorkomt.
2Dit is de Markov eigenschap.
1.2 Discrete Markovketens
Een Markovketen { Xn, n = 0, 1, 2... } is een stochastisch proces die voldoet aan de Markoveigenschap [20]
P (Xn+1= xn+1|Xn= xn, Xn−1= xn−1, ..., X1 = x1) = P (Xn+1= xn+1|Xn= xn) ∀n, ∀xn.
Als dit het geval is dan is Xneen eerste orde Markovketen. Voor een tweede orde Markovketen geldt
P (Xn+1= xn+1|Xn= xn, Xn−1= xn−1, ..., X1 = x1) = P (Xn+1= xn+1|Xn= xn, Xn−1= xn−1) ∀n, ∀xn.
Hierbij is de volgende toestand van het systeem alleen afhankelijk van de huidige toestand van het model en de vorige toestand van het model. Op deze manier kunnen nog hogere orde Markovketens gemaakt3. Bij hogere orde Markovketens wordt er dus verder terug gekeken in de tijd.
In het vervolg zullen we de notatie iets versimpelen. We zullen de toestanden xi veranderen in enkele letters zoals i en j. Als we nu ´e´en stap in de tijd doen, van n naar n + 1 kunnen we de transitiekans berekenen met
ˆ
pij = P (Xn+1= j|Xn= i).
Dit zijn de (´e´enstaps) transitiekansen van de Markovketen. Om de k-staps transitiekansen te berekenen gebruiken we
ˆ
p(k)ij = P (Xn+k = j|Xn= i).
De matrix ˆP die geconstrueerd wordt door ˆpij in rij i en kolom j te zetten (∀
i,j) wordt de (´e´enstaps) transitiematrix genoemd. De elementen van deze matrix moeten voldoen aan de volgende eigenschappen:
• ˆP (i, j) ≥ 0, ∀i, j.
• P
alle j
P (i, j) = 1, ∀i.ˆ
Een matrix die aan deze eigenschappen voldoet is een stochastische matrix.
3Een nulde orde Markovketen is onafhankelijk van vorige toestanden.
Definitie 1.1 (Homogene Markovketen) Een Markovketen is homogeen wanneer geldt [18]
P (Xn+1= j|Xn= i) = P (X1 = j|X0 = i)∀n, ∀i, j.
Met andere woorden, als de transitiekansen onafhankelijk zijn van de tijd, en deze kansen niet zullen veranderen in de Markovketen. In het vervolg zullen we alleen naar homogene Markovketens kijken.
De kans bestaat dat de geelgors (teruggaand op voorbeeld 1.1) meerdere ma- len dezelfde noot herhaald. De Markovketen zal dan een aantal tijdperiodes in dezelfde toestand bevinden. Het aantal tijdperiodes dat een Markovketen achtereenvolgend in dezelfde toestand blijft wordt de verblijftijd genoemd [6]. De kans van het verlaten van toestand i is
X
i6=j
ˆ
pij = 1 − ˆpii.
De kans om in dezelfde toestand te blijven is ˆpii en de kans om de toestand te verlaten is 1 − ˆpii. Hieruit volgt dat de verblijftijd Vi met h tijdstappen een geometrische verdeling heeft en gelijk is aan
P (Vi = h) = (1 − ˆpii)ˆph−1ii .
De verwachtingswaarde en de variantie van deze geometrische verdeling zijn [5]
E[Vi] = 1 1 − ˆpii
; V ar[Vi] = pˆii (1 − ˆpii)2.
Waarmee we zouden kunnen uitrekenen hoe vaak de geelgors een bepaalde noot achter elkaar gaat zingen.
We kunnen nu een Markovketen op een andere manier simuleren. We ke- ken steeds naar de verandering van elke tijdstap, waardoor het mogelijk was om in dezelfde toestand te blijven. We kunnen nu ook kijken naar de Markovketen waarbij we alleen kijken naar de tijdstappen waar de toestan- den daadwerkelijk veranderen. Hierdoor zou de verblijftijd dus altijd gelijk zijn aan 1. We kunnen deze Markovketen maken door de transitiekansen te veranderen in
ˆ
qij = pˆij
(1 − ˆpii) voor j 6= i en ˆqii= 0.
Deze nieuwe Markovketen wordt de ingebedde Markovketen genoemd [20].
Stelling 1.1 (Chapman Kolmogorov vergelijkingen)
Neem een Markovketen Xm met transitiekans ˆpij. Dan geldt ∀i,j [20]
ˆ
p(k)ij = X
alle z
ˆ
p(v)iz pˆ(k−v)zj voor 0 < v < k.
Bewijs
ˆ
p(k)ij = P (Xk= j|X0 = i)
= X
alle z
P (Xk= j, Xv = z|X0 = i) voor 0 < v < k
= X
alle z
P (Xk= j|Xv = z, X0= i)P (Xv = z|X0 = i) voor 0 < v < k
= X
alle z
P (Xk= j|Xv = z)P (Xv = z|X0= i) voor 0 < v < k
= X
alle z
ˆ
p(v)iz pˆ(k−v)zj voor 0 < v < k.
Als we deze Chapman Kolmogorov vergelijking opschrijven in matrix notatie wordt dit
Pˆ(k)= ˆP(v)Pˆ(k−v).
Omdat bij een homogene Markovketen de transitiematrix niet veranderd wordt dit
Pˆk= ˆPvPˆk−v.
Met de Chapman Kolmogorov vergelijking kunnen we sneller een k-stap transitiekansen berekenen. Dit gaan we duidelijk maken in een voorbeeld:
Voorbeeld 1.2
Stel dat een geelgors drie verschillende noten fluit. Een a, b en een c met transitiematrix
Pˆ1=
0.3 0.1 0.6 0.05 0.45 0.5 0.1 0.2 0.7
.
Bij een (´e´enstaps) transitiematrix staat in de ide rij en de jde kolom voor het aantal transities van toestand i naar toestand j gedeeld door het aantal transities van de ide toestand.
We kunnen nu de kans uitrekenen dat de geelgors over 2 noten een b fluit, gegeven dat de geelgors nu een a fluit. We krijgen dan
ˆ
paapˆab+ ˆpabpˆbb+ ˆpacpˆcb= 0.3 × 0.1 + 0.1 × 0.45 + 0.6 × 0.2 = 0.1950.
a
a
a
a
a
b b
b
b c
c c
c
0.3
0.3
0.1
0.1
0.05
0.45 0.5 0.6
0.6
0.1 0.25
0.7 0.2
Figuur 1.2: Grafische weergave4, hoe je kan berekenen om van noot a naar noot b te gaan met ´e´en tussennoot.
We kunnen dit ook op een andere manier zien, namelijk door eerst de twee- staps transitiematrix uit te rekenen. We krijgen dan
Pˆ2=
0.1550 0.1950 0.6500 0.0875 0.3075 0.6050 0.1100 0.2400 0.6500
.
Als we vervolgens kijken naar (i,j)=(1,2) zien we dat deze gelijk is aan 0.1950 precies als we al eerder hadden uitgerekend.
Laat
πj(n) = P (Xt= j)
de kans zijn dat de Markovketen na n tijdstappen zich in toestand j bevindt.
Vorig voorbeeld zouden we dan kunnen herschrijven als
π(2) = π(0) ˆP2= (1, 0, 0) ˆP2 = (0.1550, 0.1950, 0.6500), met πb(2)=0.1950. In het algemeen zou dan gelden
π(n) = π(n − 1) ˆP = π(0) ˆPn.
Op deze manier kunnen we snel π berekenen na n tijdstappen. Als we bij dit voorbeeld hogere staps transitiematrixen berekenen dan zien we dat deze matrix convergeert naar een matrix waar de rijen identiek zijn. We hebben
n→∞lim Pˆn=
0.1092 0.2521 0.6387 0.1092 0.2521 0.6387 0.1092 0.2521 0.6387
.
De transitiekansen hebben zich dan naar een evenwicht ontwikkeld.
Definitie 1.2 (Stationaire distributie)
Een Stationaire distributie van een Markovketen is een verdeling π die in- variant blijft bij vermenigvuldiging van de transitiematrix ˆP .
Als
πj ∈ <, 0 ≤ πj ≤ 1 en X
alle j
πj = 1, 1,
dan is π een stationaire distributie dan en slechts dan als π ˆP =π[9].
4De tekeningen, zoals deze, zijn gemaakt in Inkscape.
Definitie 1.3 (Limietverdeling) De kansverdeling l waarvoor geldt
l = lim
n→∞π(n) = lim
n→∞π(0) ˆPn= π(0) lim
n→∞
Pˆn,
wordt de limietverdeling van de Markovketen genoemd. Wanneer de toe- standen van de Markovketen ergodisch5 zijn dan bestaat de limietverdeling altijd en is uniek [20].
Definitie 1.4 (Evenwichtsverdeling)
Een limietverdeling l waarvoor geldt dat deze convergeert onafhankelijk van π0, waarvan de elementen van de vector strikt positief zijn en bij elkaar opgeteld ´e´en zijn, wordt een evenwichtsverdeling genoemd [20].
5In de volgende paragraaf zal dit behandeld worden.
1.3 Classificatie van toestanden
In deze paragraaf zullen we de verschillende eigenschappen van individuele toestanden bespreken.
Zo is er een verschil tussen overgangstoestanden, recurrente toestanden, en absorberende toestanden. Overgangstoestanden zijn toestanden waarbij er een kans bestaat dat de Markovketen nooit meer terugkeert naar zo’n toe- stand. Bij recurrente toestanden is de kans dat de Markovketen nooit meer naar zo’n toestand komt nul, waardoor een Markovketen altijd naar zo’n toestand kan terugkeren. Als een Markovketen in een absorberende toe- stand komt, zal de Markovketen nooit meer uit deze toestand komen. We zullen dit verduidelijken met een plaatje:
1 2
3
4 5 6
8 7 9
Figuur 1.3: Verschillende toestanden.
• Toestand 1,2,3,4 en 5 zijn overgangstoestanden. Uiteindelijk zal de Markovketen deze toestanden verlaten en hier niet meer terugkomen.
• Toestand 7,8 en 9 zijn recurrente toestanden. Eenmaal bij ´e´en van deze toestanden aangekomen, zal de Markovketen telkens naar deze toestand terugkeren.
• Toestand 6 is een absorberende toestand. Als de Markovketen bij deze toestand is aangekomen, zal de Markovketen altijd in deze toestand blijven.
Verder zijn de drie toestanden 7, 8 en 9 periodiek met een periode drie. Als een toestand i wordt verlaten en het is alleen mogelijk om terug te keren naar toestand i na een veelvoud van p transities met p > 1, dan wordt toestand i periodiek met een periode p genoemd. Een toestand waarbij p=1 wordt aperiodiek genoemd. Een toestand die zowel aperiodiek als positief recur- rent6 is, wordt ergodisch genoemd. Als alle toestanden in een Markovketen ergodisch zijn, dan is de Markovketen ergodisch.
Definitie 1.5 (Terugkeertijd)
Het minimale aantal transitieovergangen in een Markovketen van toestand i naar dezelfde toestand i wordt de terugkeertijd genoemd. Dit noteren we als fii
Definitie 1.6 (Positief recurrent)
Een toestand i is positief recurrent als i recurrent is en als geldt [13]
E[fii] < ∞.
Recurrente toestanden waarbij de terugkeertijd oneindig is worden nul-recurrent genoemd. Voor een eindige Markovketen geldt:
• Geen enkele toestand is nul-recurrent
• Ten minste ´e´en toestand is positief recurrent
Omdat wij alleen eindige Markovketens zullen behandelen zullen in het ver- volg deze twee eigenschappen bij de Markovketens van toepassing zijn.
6Zie definitie 1.6.
Theorie 1.1
Laat j een toestand zijn van een eindige discrete Markovketen, dan geldt [20]:
• Als j een nul-recurrent of een overgangstoestand is en i een willekeurige toestand van de Markovketen dan
n→∞lim pˆ(n)ij = 0.
• Als j positief recurrent en aperiodiek is dan
n→∞lim pˆ(n)jj > 0.
Definitie 1.7 (Bereikbaar)
Een toestand j wordt bereikbaar vanuit toestand i genoemd als ˆpnij > 0 voor een zeker aantal tijdstappen n [12].
De toestanden i en j communiceren als j bereikbaar is vanuit i en andersom.
Definitie 1.8 (Irreducibel)
De Markovketen is irreducibel als alle toestanden met elkaar communiceren [13].
Hoofdstuk 2
Goodness of fit
Om uiteindelijk te onderzoeken hoe goed ons model is, zullen we enige the- orie behandelen over maximum-likelihood-schatter voor Markovketens.
2.1 Maximum-likelihood-schatter
Neem ˆP de geschatte transitiematrix op basis van de data, L het aantal geobserveerde toestanden van de Markovketen, y de geobserveerde data met y = (y0, y1, · · · , yL) en Y de gemodelleerde data op basis van de geobser- veerde data met Y = (Y0, Y1, · · · , YL). De log-likelihood van een Markovke- ten is dan
lY =y( ˆP(1)) = log[P (Y = y| ˆP )]
= log[PPˆ(Y0 = y0)PPˆ(Y1= y1|Y0 = y0) · · · PPˆ(YL= yL|Y0· · · YL−1)]
=1 log[PPˆ(Y0 = y0)
L−1
Y
i=0
PPˆ(Yi+1= yi+1|Yi = yi)]
= log[PPˆ(Y0 = y0)] +
L−1
X
i=0
log[PPˆ(Yi+1= yi+1|Yi = yi)]
= C +
L−1
X
i=0
log[ ˆPYiYi+1]
= C +X
ij
nijlog[ˆpij].
1Vanwege de Markoveigenschap
Waarbij nij het aantal keer is dat de transitie (i, j) voorkomt. Verder geldt C = log[PPˆ(Y0 = y0)], maar we nemen Y0 als gegeven, hierdoor is PPˆ(Y0= y0) = 1 waardoor C wegvalt.
We kunnen nu de maximum-likelihood-schatter afleiden door middel van de Langrange-multiplicator [19]. Stel dat de Markovketen m verschillende parameters heeft. Dus dat ˆP een matrix is met m2 elementen. Dan hebben we m vergelijkingen waarvoor geldt
X
j
ˆ
pij = 1 ∀i. (2.1)
Hierdoor hebben we m Langrange-multiplicators: λ1, λ2, · · · λm. Hiermee verkrijgen we de Langrangefunctie
Λ(ˆpij, λi) = lY =y( ˆP(1)) +
m
X
i=1
λi(X
j
ˆ pij− 1).
Bij het nemen van de afgeleiden naar λi krijgen we de vergelijkingen (2.1) terug. Bij het nemen van de afgeleide naar ˆpij krijgen we
0 = nij ˆ pij
+ λi
λi= −nij
ˆ pij
ˆ
pij = −nij
λi . Vanwege de vergelijkingen (2.1) hebben we
m
X
j=1
nij
λi = −1
m
X
j=1
nij = −λ.
Waardoor de maximum-likelihood-schatter gelijk is aan ˆ
pij = nij
m
P
j=1
nij .
Dit alles is verkregen onder de voorwaarden dat de Markovketen eerste orde is. De afleiding voor hogere ordes gaat soortgelijks en krijg je hetzelfde antwoord, behalve dat ˆpij ditmaal een m×mu matrix is, waarbij u de orde van de Markovketen is.
Met behulp van deze maximum-likelihood-schatter kunnen we nieuwe lied- jes van de geelgors simuleren. De naam maximum-likelihood-schatter laat al zien dat het om een geschatte maximum-likelihood gaat. De matrix ˆP is een geschatte transitiematrix die kan veranderen naarmate er meer data bij zou komen. Het is daarom niet mogelijk om de ’ware’ transitiematrix2 P te berekenen, omdat deze altijd weer verandert als er meer transities bij- komen. Wel wordt de geschatte transitiematrix preciezer naarmate er meer data gebruikt wordt.
Bij een geschatte transitiematrix kan daarom rekening worden gehouden met zekere foutmarge . De maximum-likelihood-schatter wordt dan
ˆ
pij = nij m
P
j=1
nij
± .
Waarbij nog steeds moet gelden:
X
j
ˆ
pij = 1 ∀i en ˆpij ≥ 0, ∀i, j.
De keuze van kan verschillen, maar een logische keus zou zijn om rekening te houden met de hoeveelheid geobserveerde data. Stel we hebben te maken met een liedje van de geelgors, met L opeenvolgende noten. Voor ˆP geldt dan
Pˆij ∼ Bin(L, Pij) L
Voor een Binomiale verdeling geldt E[ ˆPij] = Pij en V ( ˆPij) = Pij(1 − Pij)
L .
Door deze variantie is een logische keuze om < 1
2L te kiezen. In het vervolg zullen we deze foutmarge niet meenemen in de berekeningen, maar de transitiematrix die we gaan gebruiken is dus een geschatte matrix en niet de matrix met de precieze kansen om van de ene noot naar de andere noot te gaan.
2Met ’ware’ transitiematrix wordt de matrix bedoeld waarin de precieze kansen staan om van de ene toestand in de andere toestand te komen. Zo kan in een geobserveerde data een bepaalde overgang niet voorkomen terwijl de kans op deze overgang niet nul is.
2.2 Likelihood ratio
Laat L( ˆPu, y) de likelihood zijn van de transitiematrix ˆP met orde u. Laat verder L( ˆPu?, y) de likelihood zijn van de transitiematrix ˆP met orde u? waarvoor geldt u > u?. Dan geeft de likelihood ratio [1]
λ = L( ˆPu, y) L( ˆPu?, y)
een manier om te bepalen wat de goodness of fit voor het model is. Waarbij de goodness of fit van een model beschrijft hoe goed het model past bij de geobserveerde waarden. De log-likelihood ratio logλ = l( ˆPu, y) − l( ˆPu?, y) is in de praktijk makkelijker toe te passen. Als je de log-likelihood ratio met twee vermenigvuldigt krijg je de log-likelihood statistiek. Deze log-likelihood statistiek is asymptotisch chi-kwadraat verdeeld, in formulevorm wordt dit
D = 2[l( ˆPu, y) − l( ˆPu?, y)] ∼ χ2dfu−df
u?
Hierbij is dfu het aantal parameters van ˆPu en is te berekenen met3 dfu= mu(m − 1)
waarbij m het aantal toestanden is van de Markovketen met orde u.
3Iets soortgelijks geldt voor u?
2.3 Statistiche toetsen
Met behulp van de log-likelihood kunnen we testen welke modellen het liedje goed beschrijven en welk model we het beste kunnen kiezen [1]. Stel we heb- ben twee modellen die niet even groot zijn, in dit geval door de verschillende orde tussen de twee modellen. Stel verder dat het kleinere model gebaseerd is op een Markovketen van orde u? en het grotere model gebaseerd is op een Markovketen van orde u. Er geldt dan u > u?. Dan kunnen we de test opstellen
H0 : Het kleinere model met orde u? moeten we kiezen.
H1 : Het grotere model met orde u moeten we kiezen.
We kunnen H0tegen H1testen door de log-likelihood statistiek te gebruiken.
Als beide modellen de data goed beschrijven moet er gekeken worden of de extra parameters in het grotere model de verbetering ten opzichte van het kleinere model waard zijn. De nulhypothese H0 wordt verworpen als D zich in het kritische gebied bevindt. D bevindt zich in het kritische gebied als er minder dan 100×α % van de χ2dfu−df
u? verdeling meer is dan D4.
4Van nu af aan zullen we α = 5 kiezen
Hoofdstuk 3
Het model
De zang van de geelgors is een zeer typerend geluid. Hierdoor is het niet nodig om een heel scala vogelgeluiden te weten om de zang van de geelgors te herkennen. Indien we de voorgaande theorie¨en toepassen op de zang van een geelgors kunnen we de zang van de geelgors simuleren en testen hoe goed deze gemodelleerde zang is. Dit zullen we in dit hoofdstuk behandelen. Ook zullen we kijken of het typerende geluid van de geelgors is terug te vinden in de gemodelleerde zang.
3.1 De zang van de geelgors
3.1.1 Hoe vogels geluid maken Vogels maken niet op dezelfde manier geluid als zoogdieren doen [3]. Zoog- dieren maken geluid met hun stem- banden in het strottenhoofd, die zich bovenin de luchtpijp bevindt.
Figuur 3.1: Versimpelde weergave van de syrinx [16].
Vogels daarentegen gebruiken de sy- rinx. Dit orgaan bevindt zich pre- cies op de plek waar de luchtpijp wordt opgesplitst naar de primaire bronchi¨en, oftewel aan het eind van de luchtpijp. Stembanden bestaan uit spieren en kraakbeentjes die de bijbehorende stemplooien laten tril-
len. Deze trilling wordt overgebracht op de luchtstroom die door de luchtpijp stroomt bij het uitademen. Bij de syrinx van de vogel wordt een resonerende ruimte die in verbinding staat met elastische membranen gebruikt om de luchtstroom te laten trillen.
Er wordt bij vogelgeluiden onderscheid gemaakt tussen vogelzang en vogel- roep. De roep is korter dan de zang met een maximum van drie elementen.
De zang komt alleen voor bij zangvogels die hierbij spieren gebruiken die de niet-zangvogels niet hebben. Zangvogels gebruiken ook de roep, denk bijvoorbeeld aan de bedelroep van een jonge koolmees, de alarmroep van de merel of de vluchtroep van de vink. Bij het modelleren van het geluid van de geelgors zullen we ons alleen richt op de zang. Dit omdat in de roep wei- nig verschillende elementen voorkomen omdat deze maar uit drie elementen kan bestaan. Hierdoor is dit statistisch minder interessant en is de vogel hierin minder uniek dan in zijn zang. Zangvogels kunnen dus een ander geluid maken dan niet-zangvogels, maar ook elke zangvogelsoort maakt een ander geluid. Dit verschil in zang komt voornamelijk door het verschillende leerproces en het verschil in grootte van het gedeelte in de hersenen die de zang regelt.
3.1.2 Het typerende geluid van de geelgors
Ter introductie, een aantal quotes over de zang van de geelgors uit verschil- lende vogelboekjes:
“Roep een onzuiver stuuff. Heeft ook diverse meer gesmoorde, korte, tikkende roepjes, pt...pt,pt,pittiLITT...; ook een scherp tsit en een fijn, lang- gerekt tsiih. Zang welbekende strofe van 5 − 8 snel herhaalde korte noten met verschillend einde, si-si-si-si-si-si-SUUUUUU : voorlaatste noot vaak hoger en laatste lager, sre-sre-sre-sre-sre-sreSIIII-suuuu; klank soms meer herinnerend aan Krekelzanger, dre-dre-dre-...” [7]
“Zang kenmerkend, beginnend met zes snelle, hoge, vaak insektachtige, sjirpende noten, gevolgd door een langgerekte, hese, weemoedige eindtoon tsietsie-tsie-tsie-tsie-tsie-tzuuuuuuh. Roept indien opgeschrikt een kort tsik, soms gevolgd door een droog, rollend prullulu... ” [11]
“De zang van de geelgors klinkt van een uitkijkpost of op telegraaf- draden. Het is een eentonige serie klanken gevold door een vaak tweetonige uithaal. De vluchtroep is een metaalachtig tsik.” [10]
Uit bovenstaande stukjes is al duidelijk te merken dat het moeilijk is om vogelgeluiden op schrift te zetten. Zo hoort de ene ornitholoog een roep als een tzuuuh en de ander als suuu. Ondanks dat het moeilijk is om de vogelgeluiden in fonetische spelling op te schrijven kunnen we nog wel wat meer zeggen over de zang van de geelgors.
Bij de zang van de geelgors worden drie verschillende zangtypes onderschei- den [17]. Waarbij het verschil tussen de zangtypes afhangt van het einde van de zangtypes. Het meest eenvoudige stuk bestaat uit een serie korte een- tonige noten met op het eind een langgerekte noot waarbij geen modulatie in frequentie optreedt. Deze laatste noot wordt door sommige vogelboeken beschreven als een suuuu noot. We noemen dit zangtype voor het verdere verloop van de scriptie zangtype A (figuur 3.2A).
In het tweede zangtype is er voor de suuuu noot een hoge korte noot. Bij deze noot is er ook zo goed als geen modulatie in frequentie. Deze hoge noot wordt fonetisch soms als si opgeschreven. We noemen dit zangtype voor het verdere verloop van de scriptie zangtype B (figuur 3.2B ).
Het derde zangtype is hetzelfde als het eerste zangtype, behalve dat er nu nog een hoge noot achter komt. Deze noot verandert in frequentie. Deze noot wordt fonetisch geschreven als siii. We noemen dit zangtype voor het verdere verloop van de scriptie zangtype C (figuur 3.2C ).
Als een geelgors enige tijd aan het fluiten is, kunnen alle drie de zangtypes erin voorkomen, ook is het mogelijk dat een geelgors het suuuu ,si-suuuu of het suuuu-siii gedeelte weglaat waardoor alleen de korte eentonige noten gezongen worden. We noemen dit zangtype voor het verdere verloop van de scriptie zangtype D (figuur 3.2D ).
Figuur 3.2: Sonograms die de verschillende zangtypes van de zang van de geelgors laten zien1.
1Verkregen met SAP.
3.1.3 De vijfde symfonie van Beethoven
Een leuk weetje is de inspiratie die Beethoven misschien door de geelgors heeft gekregen. Onder andere Carl Czerny, componist en leerling van Beet- hoven, beweerde dat de eerste vier noten van de vijfde symfonie van Beetho- ven door Beethoven ge¨ımiteerd zou zijn van de geelgors. Zo fluit de geelgors vaak vijf of meer korte noten gevolgd door een lange lage noot. Vergelijk- baar met de dun dun dun dunnnn van Beethoven.
Anton Felix Schindler, een biograaf van Beethoven, vertelde dat de ’Scene at the brook’ ook ontleend zou zijn aan de geelgors. Dit is niet gemakkelijk te horen omdat alleen de stijgende toonhoogte van de zang van de geelgors terug te vinden is:
“On a bright, sunny dat in April, 1823, Beethoven took Schindler for a long ramble through the scenes in which he had composed his Fifth and Sixth symphonies. Schindler writes(1, pp. 153-154): “After we had loo- ked at the bath-house and its adjacent garden at Heiligenstadt and he had given expression to many agreeable recollections touching his creations, we continued our walk towards the Kahlenberg in the direction past Grinzing.
Passing through the pleasant meadow-valley between Heiligenstadt and the latter village, which is traversed by a gently murmuring brook which hurries down from a near-by mountain and is bordered with high elms, Beethoven repeatedly stopped and let his glances roam, full of happiness, over the glo- rious landscape. Then seating himself on the turf and leaning against an elm, Beethoven asked me if there were any yellowhammers2 to be heard3 in the trees around us. But all was still. He then said: ’Here I composed the
“Scene by the Brook” and the yellowhammers up there, the quils, nighting- ales and cuckoos round about, composed with me.’ To my question why he had not also put the yellowhammers into the scene4, he drew out his sketchbook and wrote:
’That’s the composer up there,’ he remarked, ’hasn’t she a more important role to play than the others? They are meant only for a joke.’ And really the
2Een yellohammer is de engelse benaming van een geelgors.
3Beethoven was rond 1819 volledig doof geworden.
4Beethoven had bij de zesde symfonie geschreven dat de nachtegaal bij de fluit hoorde, de kwartel bij de hobo en de koekoek bij de klarinet.
entrance of this figure in G major gives the tone-picture a new charm. Spea- king now of the whole work and its parts, Beethoven said that the melody of this variation from the species of the yellowhammers was pretty plainly imi- tated in the scale written down in Andante rhythm and the same pitch.” [21]
In een iets ouder boek staat het nog duidelijker beschreven:
“But the note of the yellowhammer, both in England and in Austria, is not an arpeggio- cannot in any way be twisted into one, or represented by one. It is a quick succession of the same note, ending with a longer one, sometimes rising above the preceding note, but more frequently falling. In fact, Schindler himself tells us that it was the origin of the mighty theme which opened the C minor Symphony!5” [8]
3.2 Data
3.2.1 Verwerken van de data
De verschillende vogelgeluiden zijn verkregen van de website: http://www.xeno- canto.org/. Vervolgens is dit omgezet in bruikbare data door middel van de programma’s Sound Analysis Pro(SAP) en Melodyne Singletrack (MS). De eigenschappen van beide programma’s zal worden uitgelegd met behulp van de data van een grijze boswinterkoning. Dit is een vogeltje uit Zuid-Amerika die telkens een reeks van noten herhaald. Met behulp van SAP kunnen we het lied van deze vogel omzetten in een sonogram.
Figuur 3.3: Een sonogram van vijftien noten van een grijze boswinterkoning.
5De hele titel van de vijfde symfonie is:’The Symphony No. 5 in C minor ’
In dit geval fluit de vogel steeds vijf verschillende noten om vervolgens in precies dezelfde volgorde die vijf noten te herhalen. In de sonogram hier- boven (figuur 3.3) is duidelijk te zien dat het steeds om dezelfde vijf noten gaat. Met het programma MS kunnen we bepalen wat de daadwerkelijke muzieknoot van een gefloten noot is (figuur 3.4).
Figuur 3.4: De bijbehorende muzieknoten van de gefloten noten van een grijze boswinterkoning.
Ondanks dat figuur 3.3 en figuur 3.4 hetzelfde gefloten gedeelte laten zien, zien we bij het programma SAP vijf verschillende noten waar we bij het pro- gramma MS acht verschillende noten zien. Dit komt omdat het programma MS fijner differentieert.
Figuur 3.5: Het vergelijken van de programma’s SAP en MS.
In figuur 3.5 worden de zesde en elfde noot van eerder genoemd sonogram
(figuur 3.3) met elkaar vergeleken. Zo wordt de zesde noot van SAP in MS omgezet in een Bb8 en een A7. Hierbij staat de hoofdletter voor de muziek- noot en geeft de kleine letter b aan dat het gaat om een noot die met een halve toon is verlaagd, dus een Bb wordt een Bes. Verder geeft het getal aan, uit hoeveelste octaaf de noot komt.
Het hele liedje van de grijze boswinterkoning duurt in totaal 122 secon- den. Met het programma SAP worden in dit liedje vijf verschillende noten gebruikt terwijl MS zevenentwintig verschillende noten noteert. Hierdoor zien we dat het uitmaakt welk programma we gebruiken voor het maken van een model. Zo kunnen we het liedje van de grijze boswinterkoning op basis van de verkregen data van SAP al met een eerste orde Markovketen exact voorspellen. Dan zijn het namelijk vijf verschillende noten die elkaar afwisselend opvolgen volgens het patroon in figuur 3.3. Als we het liedje om- zetten met MS kunnen we het liedje pas met een achtste orde Markovketen exact voorspellen. In het vervolg zullen we SAP gebruiken om het liedje om te zetten in noten. Dit omdat het model hiervan niet te groot wordt en dat SAP beter het aantal verschillende noten kan laten zien. In figuur 3.5 zien we namelijk dat de vogel volgens SAP twee dezelfde noten fluit, maar doordat twee gefloten nooit precies hetzelfde zijn ziet MS hier zelfs vier verschillende noten in. Door de sonogram kunnen we echter zien dat het ongeveer om
´e´en en dezelfde noot gaat. Daarom zullen we een gefloten liedje eerst met SAP omzetten in een m aantal verschillende noten, vervolgens kunnen we deze noten met MS omzetten in muzieknoten. Hierbij wordt gekeken welke muzieknoot het vaakst in MS voorkomt bij een bepaalde noot in SAP.
3.2.2 Het maken van modellen op basis van de data
Nu kunnen we de verschillende liedjes van de geelgors analyseren. Het lang- ste liedje die te vinden is op http://www.xeno-canto.org/ duurt 408 secon- den. Bij dit liedje fluit de geelgors 453 noten waarvan er 9 verschillend zijn6. Deze 9 verschillende noten verschillen in frequentie en tijdsduur zoals in fi- guur 3.6 te zien is. Dit liedje gebruiken we in het verdere verloop van dit verslag, omdat dit liedje de meeste data bevat7.
6Volgens SAP.
7De link voor dit liedje is http://www.xeno-canto.org/123750.
Figuur 3.6: Negen verschillende noten van een gefloten liedje van de geel- gors8.
De eerste drie noten van het liedje nemen we als gegeven, zodat we een nulde, eerste, tweede en derde orde transitiematrix kunnen maken. De eerste orde transitiematrix van dit liedje is9
Pˆ1=
0.538 0.231 0 0.077 0.154 0 0 0 0
0.500 0 0 0 0.500 0 0 0 0
0 0 0.316 0.684 0 0 0 0 0
0 0 0 0.838 0.010 0 0 0.152 0
0 0.024 0.012 0 0.602 0.337 0 0.024 0
0 0 0 0 0 0.692 0.209 0.099 0
0 0 0 0 0 0 0.740 0.260 0
0.023 0.023 0.045 0.045 0.364 0 0 0 0.500
0.045 0 0.455 0 0.500 0 0 0 0
10.
Elke toestand in Markovketen is zowel positief recurrent als aperiodiek waar- door de Markovketen ergodisch is. Volgens theorie 1.1 zou moeten gelden dat lim
n→∞pˆ(n)jj > 0. De limietverdeling voor deze Markovketen is l=
0.023 0.012 0.043 0.222 0.185 0.202 0.162 0.100 0.050
. Dit is onafhankelijk van π0 waardoor bovenstaande limietverdeling de even- wichtsverdeling is. Dus lim
n→∞pˆ(n)jj > 0 is inderdaad waar en verder is de Markovketen irreducibel.
8De twee grafieken zijn verkregen met SAP met behulp van clustering.
9De getallen in de matrix zijn in drie decimalen nauwkeurig.
10Zie bijlage: A.
3.3 Uitkomsten
3.3.1 Statistiche toetsen
Nu we de nulde orde, eerste orde, tweede orde en derde orde Markovke- ten hebben gemaakt, kunnen we op basis van deze Markovketens een nieuw liedje van de geelgors simuleren. Voor deze verschillende ordes hebben we elk twintigduizend nieuwe liedjes gesimuleerd om zo het ´e´en en ander erover te kunnen zeggen. Zo hebben we gecontroleerd of de likelihood van deze twintigduizend liedjes overeenkomt met bijbehorende Markovketen.
Figuur 3.7: Vergelijken van de likelihood van twintigduizend gesimuleerde liedjes met de likelihood van het originele liedje11.
In figuur 3.7 wordt de likelihood van twintigduizend gesimuleerde liedjes vergeleken met de likelihood van het originele liedje. Hierbij zijn de twin- tigduizend liedjes gesimuleerd met een nulde orde Markovketen, en is de likelihood bij zowel de twintigduizend gesimuleerde liedjes als het originele liedje met behulp van een nulde orde Markovketen berekend. Uit dit plaatje valt af te lezen dat de likelihood van het originele liedje het gemiddelde is van deze twintigduizend gesimuleerde liedjes. Voor de hogere ordes is hetzelfde resultaat af te lezen12. Dit is precies wat wiskundig te verwachten valt. Het
11Zie bijlage B.1.
12Voor de bijbehorende code zie bijlage B.2, B.3 en B.4.
is ook noodzakelijk dat deze twintigduizend gesimuleerde liedjes gemiddeld dezelfde eigenschappen hebben als het model, dit omdat we verderop achter meer eigenschappen van de verschillende modellen willen komen door mid- del van twintigduizend liedjes te simuleren.
Zonder dat we al die twintigduizend liedjes simuleren kunnen we ook al wat zeggen over de verschillende modellen. We kunnen immers de verschil- lende likelihood uitrekenen, en daarmee kunnen we met de log-likelihood statistiek de verschillende modellen met elkaar vergelijken. Hieronder staat de tabel met de verschillende Markovorders ingevuld13.
Orde kleine model Orde grote model df χdf ;0.95 D Nulhypothese
0 1 64 83.675 1059.999 verwerpen
0 2 640 699.963 1144.966 verwerpen
0 3 5824 6002.652 1208.972 behouden
1 2 576 632.942 84.967 behouden
1 3 5760 5937.674 148.973 behouden
2 3 5184 5352.614 64.006 behouden
Hierbij geldt steeds dat de nulhypothese is dat we het model moeten kiezen met de kleinere orde. Uit deze tabel valt af te lezen dat een eerste orde Markovketen het beste model is. Daarna komt een tweede orde Markovketen gevolgd door een nulde orde Markovketen. Het model dat we echt niet moeten kiezen zou een derde orde Markovketen zijn. Dit komt omdat dit model erg veel parameters nodig heeft.
13Deze en volgende getallen zijn in drie decimalen nauwkeurig.
3.4 Verbetering van het model
3.4.1 Ingebedde Markovketen
Als we nogmaals kijken naar de eerste orde transitiematrix van paragraaf 3.2.2 dan zien we dat de kansen op de hoofddiagonaal behoorlijk hoog zijn.
Met andere woorden, de kans dat er meerdere keren dezelfde noot herhaald wordt is vrij groot. We kunnen daarom misschien een verbeterd model ma- ken met behulp van een ingebedde Markovketen. We maken dan eerst een ingebedde Markovketen om vervolgens handmatig de herhalende noten toe te voegen. Er worden dan E[Vi]−1 = pˆii
1 − ˆpii
herhalende noten toegevoegd14. In paragraaf 3.1.2 zijn vier verschillende zangtypes van de geelgors ge¨ıntroduceerd.
In het liedje van de geelgors dat we het meest hebben onderzocht, dit omdat dit liedje het langst duurde en dus de meeste data bevatte, hebben we geteld hoe vaak de verschillende zangtypes voorkomen.
In de vogelzang van deze geelgors kwam zangtype
• A 23 keer voor.
• B 0 keer voor.
• C 22 keer voor.
• D 3 keer voor.
Deze aantallen kunnen we vervolgens vergelijken met het aantal keer dat deze zangtypes voorkomen in twintigduizend gesimuleerde liedjes. Als we dit doen voor zangtype C krijgen we figuur 3.8.
14Een beter model is om het aantal herhalende noten met dezelfde verwachtings- waarde Poisson verdeeld te laten wezen. Maar omdat we voor het verdere verloop niet ge¨ınteresseerd zijn naar het aantal herhalende noten zullen we dit niet toepassen.
Figuur 3.8: Twintigduizend keer een eerste orde ’normale’ Markovketen gesimuleerd.
Als we vervolgens het model iets aanpassen naar een ingebedde Markovketen krijgen we figuur 3.915.
Figuur 3.9: Twintigduizend keer een eerste orde ingebedde Markovketen gesimuleerd.
15Bij beide figuren geldt dat de rode lijn het aantal keren is dat zangtype C in de originele vogelzang voorkwam(22 keer). Voor de code zie bijlage B.2, B.5 en B.6.
Als we de twee figuren met elkaar vergelijken zien we dat het model van fi- guur 3.9 het liedje veel beter weergeeft dan het model van figuur 3.8. Ook de andere zangtypes en andere Markovordes gaven soortgelijke figuren16. Hier- uit concluderen we dat een ingebedde Markovketen een verbetering geeft aan het model. Dit zullen we nogmaals laten zien door onderstaande tabel.
Waarin voor elke Markovorde en voor elk zangtype het gemiddeld aantal keren weergegeven wordt dat een bepaald zangtype voorkomt bij een ’nor- male’ Markovketen en bij een ingebedde Markoveketen17.
Markovorde Zangtype µn σn µi σi
0 A 0.183 0.423 0 0
0 C 0.0102 0.101 0 0
0 D 1.862 1.712 0 0
1 A 24.0467 3.374 22.469 3.842
1 C 11.999 2.917 22.193 3.164
1 D 2.995 1.734 2.99525 1.701
2 A 25.745 3.240 22.154 3.990
2 C 16.952 3.122 22.462 3.678
2 D 3.028 1.753 2.946 2.281
3 A 28.114 3.423 22.329 4.203
3 C 14.439 3.219 22.327 3.987
3 D 3.0564 1.810 2.976 2.352
In deze tabel zijn µn en σn de gemiddelde en de standaardafwijking van een ’normale’ Markovketen. Verder zijn µi en σi de gemiddelde en de stan- daardafwijking van een ingebedde Markovketen.
We zien in deze tabel dat als je van een nulde orde Markovketen naar een eer- ste Markovketen gaat, dit een enorme verbetering is. Verder zien we dat een nog hogere orde Markovketen niet een significante verbetering geeft. Deze twee punten hadden we ook al gezien bij de log-likelihood statistiek. Wat we in deze tabel verder nog zien is de enorme verbetering die een ingebedde Markovketen geeft.
16Zie voor de codes de bijlage B
17Voor elke rij zijn dus twintigduizend liedjes gesimuleerd met een ’normale’ Markovke- ten en met een ingebedde Markovketen.
3.4.2 Markovketens met variabele lengte
We kunnen ook kijken of Markovketens met variabele lengte een verbete- ring aan het model geven. Een Markovketen met variabele lengte is een Markovketen waarbij voor elke transitieovergang wordt gekeken of het voor- delig is om voor die overgang een hogere orde te gebruiken. Hierdoor zou het kunnen dat niet elke transitieovergang dezelfde Markovorde heeft.
Voorbeeld 3.1
Stel dat een vogel drie verschillende noten fluit, namelijk noot a,b en c. Stel verder dat deze vogel deze drie noten altijd in de volgende volgorde fluit
abcbabcbabcb...
Dan is de eerste orde transitiematrix
Pˆ1=
0 1 0
0.5 0 0.5
0 1 0
.
Hierdoor kan een eerste orde Markovketen niet een liedje simuleren die altijd volledig overeenkomt met het originele liedje. Een tweede orde Markovketen kan dit echter wel, maar dan hebben we een 9×3 transitiematrix nodig. Met behulp van een Markovketen met variabele lengte kunnen we een kleinere transitiematrix maken, waarmee we nog steeds liedjes kunnen simuleren die volledig overeenkomen met het originele liedje. Deze transitiematrix is dan
Pˆ2=
a b c
a 0 1 0
ab 0 0 1 cb 1 0 0
c 0 1 0
.
We hebben voor het liedje van de geelgors gekeken of een Markovketen met variabele lengte een significante verbetering geeft. We hebben dit bekeken ten opzichte van een ’normale’ Markovketen en een ingebedde Markovke- ten. Bij de ’normale’ Markovketen blijkt dat er geen significante verbetering plaats vindt. Bij de ingebedde Markovketen daarentegen blijkt dat er nog wel degelijk een significante verbetering mogelijk is. Als we de toets doen met
H0 : W e moeten een eerste orde ingebedde M arkovketen kiezen.
H1 : W e moeten een tweede orde ingebedde M arkovketen met variabele lengte kiezen.
Dan krijgen we p-waarde=0.026818 wat minder is dan de grenswaarde van 5% die we eerder in dit verslag als grens hadden gekozen. Daarom kun- nen we de nulhypothese verwerpen en aannemen dat een tweede orde inge- bedde Markovketen met variabele lengte een beter model is. Deze ingebedde Markovketen met variabele lengte gebruikt voor elke noot een eerste orde Markovketen behalve voor ´e´en noot waar die een tweede orde Markovketen gebruikt.
18Zie bijlage B.7.
Conclusie
Het doel van deze scriptie is om de zang van de geelgors met verschillende modellen te simuleren en te toetsen welk model we het beste kunnen kiezen om de zang te simuleren. Zo’n model is verkregen door gebruik te maken van de wiskundige theorie van Markovketens. Om vervolgens de verschillende modellen te toetsen is gekeken naar de goodness of fit van desbetreffende Markovketen. Verder stellen we nog een extra conditie aan het model, na- melijk dat het typische geluid van de geelgors terug te vinden moet zijn in de gesimuleerde liedjes.
Eerst hebben we gekeken naar ’normale’ Markovketens waar alleen de orde verschillend is. Hieruit blijkt dat een eerste orde Markovketen het beste model is. Vervolgens hebben we gekeken of een ingebedde Markovketen een verbetering geeft aan het model. Dit blijkt inderdaad zo te zijn. Door een
’normale’ Markovketen te veranderen in een ingebedde Markovketen wordt het typische geluid van de geelgors goed gesimuleerd. Goed gesimuleerd omdat de verschillende zangtypes van de geelgors in de juiste hoeveelheid voorkomen. De verschillende zangtypes komen dan dus gemiddeld in de ge- simuleerde liedjes van de geelgors even vaak voor als in het originele liedje.
Een Markovketen met variabele lengte zorgt voor een nog verdere verbete- ring. Het beste model blijkt daarom een tweede orde ingebedde Markovketen van variabele lengte te zijn.
Bibliografie
[1] Barnett A.G. en Dobson A.J., An Introduction to Generalized Linear Models, p.79-266,(2008)
[2] Basharin G.P.,Langville A.N. en Naumov V.A., The Life and Work of A.A. Markov, Linear Algebra and its Applications, p.1-12,(2004) http://langvillea.people.cofc.edu/MarkovReprint.pdf ?referrer=webcluster, (06-2013)
[3] Casteleyn C., Michiels A., Simoens P. en Van den Broeck W., Geluids- productie bij vogels, Vlaams Diergeneeskundig Tijdschrift, 2011,80, p.
15-23 (2011)
[4] CRAN (1997 ff.),The Comprehensive R Archive Network, http://cran.R- project.org/ (07-2013)
[5] Dekking F.M., Kraaikamp C., Lopuha¨a H.P en Meester L.E., A Mondern Introduction to Probability and Statistics,(2010)
[6] Gerardo R. en Sericola B., Sojourn Times in finite Markov Processes, (1989)
[7] Grant P.J. en Svensson L.,Vogelgids van Europa, p.372 , (2005) [8] Grove G., Beethoven and his Nine Symphonies, p.211, (1896) [9] Haigh J., Introduction to Markov Chains - the finite case, (2010) [10] Hayman P., Zakgids Vogels, p.44, (1999)
[11] Johnsson L, Complete gids Vogels van Nederland, p.526, (1998)
[12] Kallenberg L.C.M., Lecture notes: Besliskunde 2 (Universiteit Leiden) (2012) http://www.math.leidenuniv.nl/ spieksma/colleges/besliskunde/
Besliskunde2-12.pdf (06-2013)
[13] L´evˆeque O., Lecture notes on Markov Chains, (National University of Ireland) (2011) http://www.hamilton.ie/ollie/Downloads/Mar1.pdf (06- 2013)
[14] Nicholson E.M., The Art of Bird-Watching, (1931)
[15] Plaatje van Markov: http://estadisticamigable.blogspot.nl/2010/08/andrei- marcov-1856-1922-en-el-157.html, (06-2013)
[16] Plaatje van Syrinx: http://members.chello.nl/p.dinjens/PDFbestanden/zang.pdf (06-2013)
[17] Panov E.N., Roubtsov A.S. en Mozikov D.G., Hybridization between Yellowhammer and Pine Bunting in Russia, Dutch birding 25: 17-31, (2003)
[18] Savat B., Monte Carlo Markov Chain methoden, (2010)
[19] Shalizi C., Lecture notes on Chaos, Complexity, and Inference: Maxi- mum Likelihood Estimation for Markov Chains (Carnegie Mellon
University) (2009) http://www.stat.cmu.edu/ cshalizi/462/lectures/06/markov- mle.pdf (06-2013)
[20] Stewart W.J., Probability, Markov Chains, Queues, and Simulation, p.193-382, (2009)
[21] Thayer A.W., Thayer’ life of Beethoven, p.437, (1991)
Bijlage A
Matlab codes
In deze bijlage staat de gebruikte Matlab code.
G=zeros(9)
G(1,1)=7 G(1,2)=3 G(1,4)=1 G(1,5)=2 G(2,1)=3 G(2,5)=3 G(3,3)=6 G(3,4)=13 G(4,4)=83 G(4,5)=1 G(4,8)=15 G(5,2)=2 G(5,3)=1 G(5,5)=50 G(5,6)=28 G(5,8)=2 G(6,6)=63 G(6,7)=19 G(6,8)=9 G(7,7)=54 G(7,8)=19 G(8,1)=1 G(8,2)=1 G(8,3)=2 G(8,4)=2 G(8,5)=16 G(8,9)=22 G(9,1)=1 G(9,3)=10 G(9,5)=11
K=[G(1,:)/13;G(2,:)/6;G(3,:)/19;G(4,:)/99;G(5,:)/83;G(6,:)/91;
G(7,:)/73;G(8,:)/44;G(9,:)/22]
x=zeros(1,9) x(1,1)=1 K100=K^100 x*K100
Bijlage B
R codes
In deze bijlage staan de gebruikte R codes1.
B.1 Nulde orde Markovketens
y=rep(NA,20000) for (i in 1:20000)
{ y[i] = loglik ; i=i+1 require("seqinr")
alp <- c("1","2","3.1","3.2","4.0","4.1","4.2","5","6")
zeroOrder <- c(13/450,6/450,19/450,99/450,83/450,91/450,75/450,44/450,22/450) names(zeroOrder) <- alp
zeroOrderSeq <- sample(alp,450,rep=T,prob=zeroOrder)
zeroOrderFreq <- count(zeroOrderSeq,1,alphabet=alp,freq=TRUE) z=rep(NA,9)
for(i in 1:9){z[i]=(zeroOrderFreq[i]*450)*log(zeroOrder[i]);i=i+1 } loglik=sum(z)}
k=-2*(y)
hist(k,breaks=100,col="lightblue",xlab="Log-likelihood ",ylab="Frequentie", main="Tweeduizend log-likelihoods van een nulde orde Markovketen")
1Voor de codes heb ik gebruik gemaakt van de site: http://tata-box- blog.blogspot.nl/2012/04/introduction-to-markov-chains-and.html
B.2 Codes voor eerste orde Markovketens
y=rep(NA,20000) for (i in 1:20000)
{ y[i] = loglik ; i=i+1
alp <- c("1","2","3.1","3.2","4.0","4.1","4.2","5","6") require("seqinr")
firstOrderMat <- matrix(NA,nr=9,nc=9)
colnames(firstOrderMat) <- rownames(firstOrderMat) <- alp firstOrderMat[1,] <- c(7/13,3/13,0,1/13,2/13,0,0,0,0) firstOrderMat[2,] <- c(3/6,0,0,0,3/6,0,0,0,0)
firstOrderMat[3,] <- c(0,0,6/19,13/19,0,0,0,0,0) firstOrderMat[4,] <- c(0,0,0,83/99,1/99,0,0,15/99,0)
firstOrderMat[5,] <- c(0,2/83,1/83,0,50/83,28/83,0,2/83,0) firstOrderMat[6,] <- c(0,0,0,0,0,63/91,19/91,9/91,0)
firstOrderMat[7,] <- c(0,0,0,0,0,0,54/73,19/73,0)
firstOrderMat[8,] <- c(1/44,1/44,2/44,2/44,16/44,0,0,0,22/44) firstOrderMat[9,] <- c(1/22,0,10/22,0,11/22,0,0,0,0)
inProb <- c(1,0,0,0,0,0,0,0,0) ; names(inProb) <- alp generateFirstOrderSeq <- function
(lengthSeq,alphabet, initialProb,bfirstOrderMatrix) {outputSeq <- rep(NA,lengthSeq)
outputSeq[1] <- sample(alphabet,1,prob=initialProb) for(i in 2:length(outputSeq)){
prevNuc <- outputSeq[i-1]
currentProb <- firstOrderMatrix[prevNuc,]
outputSeq[i] <- sample(alp,1,prob=currentProb)}
return(outputSeq)}
firstOrderSeq <- generateFirstOrderSeq(451,alp,inProb,firstOrderMat)}
B.3 Tweede orde Markovketens
y=rep(NA,20000) for (i in 1:20000)
{alp <- c(1,2,3,4,5,6,7,8,9)
alpa<-c(1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25, 26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42,43,44,45,46,47,48,49,50, 51,52,53,54,55,56,57,58,59,60,61,62,63,64,65,66,67,68,69,70,71,72,73,74,75, 76,77,78,79,80,81)
require("seqinr")
secOrderMat <- matrix(NA,nr=81,nc=9) colnames(secOrderMat) <- alp
rownames(secOrderMat)<-alpa
secOrderMat[1,] <- c(4/7,2/7,0,0,1/7,0,0,0,0) secOrderMat[2,] <- c(1,0,0,0,0,0,0,0,0)
secOrderMat[3,] <- c(0,0,0,0,0,0,0,0,0) secOrderMat[4,] <- c(0,0,0,1,0,0,0,0,0) secOrderMat[5,] <- c(0,0,0,0,1,0,0,0,0) secOrderMat[6,] <- c(0,0,0,0,0,0,0,0,0) secOrderMat[7,] <- c(0,0,0,0,0,0,0,0,0) secOrderMat[8,] <- c(0,0,0,0,0,0,0,0,0) secOrderMat[9,] <- c(0,0,0,0,0,0,0,0,0)
secOrderMat[10,] <- c(1/3,1/3,0,0,1/3,0,0,0,0) secOrderMat[11,] <- c(0,0,0,0,0,0,0,0,0)
secOrderMat[12,] <- c(0,0,0,0,0,0,0,0,0) secOrderMat[13,] <- c(0,0,0,0,0,0,0,0,0) secOrderMat[14,] <- c(0,0,0,0,1,0,0,0,0) secOrderMat[15,] <- c(0,0,0,0,0,0,0,0,0) secOrderMat[16,] <- c(0,0,0,0,0,0,0,0,0) secOrderMat[17,] <- c(0,0,0,0,0,0,0,0,0) secOrderMat[18,] <- c(0,0,0,0,0,0,0,0,0) secOrderMat[19,] <- c(0,0,0,0,0,0,0,0,0) secOrderMat[20,] <- c(0,0,0,0,0,0,0,0,0) secOrderMat[21,] <- c(0,0,0,1,0,0,0,0,0) secOrderMat[22,] <- c(0,0,0,1,0,0,0,0,0) secOrderMat[23,] <- c(0,0,0,0,0,0,0,0,0) secOrderMat[24,] <- c(0,0,0,0,0,0,0,0,0)
secOrderMat[25,] <- c(0,0,0,0,0,0,0,0,0) secOrderMat[26,] <- c(0,0,0,0,0,0,0,0,0) secOrderMat[27,] <- c(0,0,0,0,0,0,0,0,0) secOrderMat[28,] <- c(0,0,0,0,0,0,0,0,0) secOrderMat[29,] <- c(0,0,0,0,0,0,0,0,0) secOrderMat[30,] <- c(0,0,0,0,0,0,0,0,0)
secOrderMat[31,] <- c(0,0,0,67/83,1/83,0,0,15/83,0) secOrderMat[32,] <- c(0,0,0,0,1,0,0,0,0)
secOrderMat[33,] <- c(0,0,0,0,0,0,0,0,0) secOrderMat[34,] <- c(0,0,0,0,0,0,0,0,0)
secOrderMat[35,] <- c(0,0,1/14,2/14,1/14,0,0,0,10/14) secOrderMat[36,] <- c(0,0,0,0,0,0,0,0,0)
secOrderMat[37,] <- c(0,0,0,0,0,0,0,0,0) secOrderMat[38,] <- c(0,0,0,0,1,0,0,0,0) secOrderMat[39,] <- c(0,0,0,1,0,0,0,0,0) secOrderMat[40,] <- c(0,0,0,0,0,0,0,0,0)
secOrderMat[41,] <- c(0,2/50,0,0,28/50,18/50,0,2/50,0) secOrderMat[42,] <- c(0,0,0,0,0,23/28,4/28,1/28,0) secOrderMat[43,] <- c(0,0,0,0,0,0,0,0,0)
secOrderMat[44,] <- c(0,1/2,0,0,1/2,0,0,0,0) secOrderMat[45,] <- c(0,0,0,0,0,0,0,0,0) secOrderMat[46,] <- c(0,0,0,0,0,0,0,0,0) secOrderMat[47,] <- c(0,0,0,0,0,0,0,0,0) secOrderMat[48,] <- c(0,0,0,0,0,0,0,0,0) secOrderMat[49,] <- c(0,0,0,0,0,0,0,0,0) secOrderMat[50,] <- c(0,0,0,0,0,0,0,0,0)
secOrderMat[51,] <- c(0,0,0,0,0,40/63,15/83,8/63,0) secOrderMat[52,] <- c(0,0,0,0,0,0,17/19,2/19,0) secOrderMat[53,] <- c(0,0,0,0,3/9,0,0,0,6/9) secOrderMat[54,] <- c(0,0,0,0,0,0,0,0,0) secOrderMat[55,] <- c(0,0,0,0,0,0,0,0,0) secOrderMat[56,] <- c(0,0,0,0,0,0,0,0,0) secOrderMat[57,] <- c(0,0,0,0,0,0,0,0,0) secOrderMat[58,] <- c(0,0,0,0,0,0,0,0,0) secOrderMat[59,] <- c(0,0,0,0,0,0,0,0,0) secOrderMat[60,] <- c(0,0,0,0,0,0,0,0,0)
secOrderMat[61,] <- c(0,0,0,0,0,0,37/54,17/54,0) secOrderMat[62,] <- c(1/19,0,1/19,0,1/11,0,0,0,6/19) secOrderMat[63,] <- c(0,0,0,0,0,0,0,0,0)
secOrderMat[64,] <- c(0,0,0,1,0,0,0,0,0) secOrderMat[65,] <- c(0,0,0,0,1,0,0,0,0) secOrderMat[66,] <- c(0,0,1,0,0,0,0,0,0) secOrderMat[67,] <- c(0,0,0,1,0,0,0,0,0)
secOrderMat[68,] <- c(0,0,0,0,9/16,7/16,0,0,0) secOrderMat[69,] <- c(0,0,0,0,0,0,0,0,0)
secOrderMat[70,] <- c(0,0,0,0,0,0,0,0,0) secOrderMat[71,] <- c(0,0,0,0,0,0,0,0,0)
secOrderMat[72,] <- c(1/22,0,10/22,0,11/22,0,0,0,0) secOrderMat[73,] <- c(1,0,0,0,0,0,0,0,0)
secOrderMat[74,] <- c(0,0,0,0,0,0,0,0,0)
secOrderMat[75,] <- c(0,0,4/10,6/10,0,0,0,0,0) secOrderMat[76,] <- c(0,0,0,0,0,0,0,0,0)
secOrderMat[77,] <- c(0,0,1/11,0,7/11,3/11,0,0,0) secOrderMat[78,] <- c(0,0,0,0,0,0,0,0,0)
secOrderMat[79,] <- c(0,0,0,0,0,0,0,0,0) secOrderMat[80,] <- c(0,0,0,0,0,0,0,0,0) secOrderMat[81,] <- c(0,0,0,0,0,0,0,0,0)
inProb <- c(0,1,0,0,0,0,0,0,0); names(inProb) <- alp secProb<- c(1,0,0,0,0,0,0,0,0)
generateSecOrderSeq <- function
(lengthSeq,alphabet, initialProb,secOrderMatrix) {outputSeq <- rep(NA,lengthSeq)
outputSeq[1] <- sample(alp,1,prob=inProb) outputSeq[2]<-sample(alp,1,pro=secProb)
for(i in 3:length(outputSeq)){
prevNuc <- outputSeq[i-1]
prevprevNuc <- outputSeq[i-2]
currentProb <- secOrderMat[(prevprevNuc-1)*9+prevNuc,]
outputSeq[i] <- sample(alp,1,prob=currentProb)}
return(outputSeq)}
secOrderSeq <- generateSecOrderSeq(452,alp,inProb,secOrderMat) secOrderSeq1 <-n2s(secOrderSeq, levels =
c("1", "2", "3.1","3.2","4.0","4.1","4.2","5" , "6"), base4 = FALSE) y[i]=loglik; i=i+1 }
B.4 Derde orde Markovketens
y=rep(NA,20000) for (i in 1:20000)
{ alp <- c(1,2,3,4,5,6,7,8,9)
alpa<-c(1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25, 26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42,43,44,45,46,47,48,49,50, 51,52,53,54,55,56,57,58,59,60,61,62,63,64,65,66,67,68,69,70,71,72,73,74,75, 76,77,78,79,80,81,82,83,84,85,86,87,88,89,90,91,92,93,94,95,96,97,98,99,100, 101,102,103,104,105,106,107,108,109,110,111,112,113,114,115,116,117,118,119, 120,121,122,123,124,125,126,127,128,129,130,131,132,133,134,135,136,137,138, 139,140,141,142,143,144,145,146,147,148,149,150,151,152,153,154,155,156,157, 158,159,160,161,162,163,164,165,166,167,168,169,170,171,172,173,174,175,176, 177,178,179,180,181,182,183,184,185,186,187,188,189,190,191,192,193,194,195, 196,197,198,199,200,201,202,203,204,205,206,207,208,209,210,211,212,213,214, 215,216,217,218,219,220,221,222,223,224,225,226,227,228,229,230,231,232,233, 234,235,236,237,238,239,240,241,242,243,244,245,246,247,248,249,250,251,252, 253,254,255,256,257,258,259,260,261,262,263,264,265,266,267,268,269,270,271, 272,273,274,275,276,277,278,279,280,281,282,283,284,285,286,287,288,289,290, 291,292,293,294,295,296,297,298,299,300,301,302,303,304,305,306,307,308,309, 310,311,312,313,314,315,316,317,318,319,320,321,322,323,324,325,326,327,328, 329,330,331,332,333,334,335,336,337,338,339,340,341,342,343,344,345,346,347, 348,349,350,351,352,353,354,355,356,357,358,359,360,361,362,363,364,365,366, 367,368,369,370,371,372,373,374,375,376,377,378,379,380,381,382,383,384,385, 386,387,388,389,390,391,392,393,394,395,396,397,398,399,400,401,402,403,404, 405,406,407,408,409,410,411,412,413,414,415,416,417,418,419,420,421,422,423, 424,425,426,427,428,429,430,431,432,433,434,435,436,437,438,439,440,441,442, 443,444,445,446,447,448,449,450,451,452,453,454,455,456,457,458,459,460,461, 462,463,464,465,466,467,468,469,470,471,472,473,474,475,476,477,478,479,480, 481,482,483,484,485,486,487,488,489,490,491,492,493,494,495,496,497,498,499, 500,501,502,503,504,505,506,507,508,509,510,511,512,513,514,515,516,517,518, 519,520,521,522,523,524,525,526,527,528,529,530,531,532,533,534,535,536,537, 538,539,540,541,542,543,544,545,546,547,548,549,550,551,552,553,554,555,556, 557,558,559,560,561,562,563,564,565,566,567,568,569,570,571,572,573,574,575, 576,577,578,579,580,581,582,583,584,585,586,587,588,589,590,591,592,593,594, 595,596,597,598,599,600,601,602,603,604,605,606,607,608,609,610,611,612,613, 614,615,616,617,618,619,620,621,622,623,624,625,626,627,628,629,630,631,632,
