1. No category

HENRI NURMINEN
KUULUVUUSALUEEN ESTIMOINTI
METROPOLISIN JA HASTINGSIN ALGORITMILLA
Kandidaatintyö
Tarkastaja: professori Robert Piché
Palautettu 21.6.2011
II
TIIVISTELMÄ
TAMPEREEN TEKNILLINEN YLIOPISTO
Teknis-luonnontieteellinen koulutusohjelma
NURMINEN, HENRI: Kuuluvuusalueen estimointi Metropolisin ja Hastingsin
algoritmilla
Kandidaatintyö, 26 sivua, 0 liitesivua
Kesäkuu 2011
Pääaine: Matematiikka
Tarkastaja: professori Robert Piché
Avainsanat: kuuluvuusalue, Markovin ketjut, Monte Carlo -menetelmät, Metropolisin ja Hastingsin algoritmi
Tässä työssä pyritään estimoimaan langattoman verkon tukiaseman kuuluvuusaluetta käyttämällä Metropolisin ja Hastingsin algoritmia. Mittausdatana käytetään joukkoa pisteitä,
joissa kuuluvuusalue on kuultu. Yksi mittauspiste sisältää pisteen paikan kaksiulotteisessa koordinaatistossa sekä pisteessä kuultujen kuuluvuusalueitten tunnistenumerot.
Työssä esitellään Markovin ketjujen sekä Metropolisin ja Hastingsin algoritmin keskeisin
teoria. Tärkein tulos on suurten lukujen laki Metropolisin ja Hastingsin algoritmille.
Kuuluvuusaluetta mallinnetaan satunnaisjakaumana, jonka tiheysfunktio kuvaa tukiaseman havaitsemisen todennäköisyyttä. Satunnaismuuttujan jakaumaksi oletetaan kahden
tasajakauman mikstuuri, koska osa mittauksista saattaa olla virheellisiä eli niin sanottuja
ulkolaisia. Virheettömien mittausten jakauma on ellipsin muotoinen ja ulkolaisten ympyrän muotoinen.
Tässä työssä myös mittausjakauman parametreja pidetään satunnaismuuttujina, joten työ
pohjautuu bayesläiseen näkemykseen tilastotieteestä. Metropolisin ja Hastingsin algoritmilla estimoidaan parametrien jakaumien odotusarvoja sekä variansseja. Bayesläisen
päättelyn avulla esitetään myös algoritmi olemassa olevan kuuluvuusalueen päivittämiseen, kun on saatu uusia havaintoja.
Algoritmia testattiin simuloidulla datalla eräissä kriittisissa tilanteissa. Näitä olivat algoritmin sietokyky virheellisille havainnoille, algoritmin toiminta tapauksessa, jossa suuri
osa havainnoista sijaitsee pienessä kuuluvuusalueen osassa, sekä kuuluvuusalueen aikainvariantti päivittäminen uusilla havainnoilla. Algoritmi todettiin laskennallisesti raskaaksi
sekä joissakin tapauksissa jonkin verran epäluotettavaksi. Toisaalta algoritmi on joustava, ja siinä on paljon muokattavia parametrejä, joilla algoritmia voidaan virittää. Sopivilla
parametreilla algoritmi toimii esimerkiksi testatuissa kriittisissä tilanteissa erittäin hyvin.
III
ALKUSANAT
Tämä kandidaatintyö on tehty Tampereen teknillisen teknillisen yliopiston Matematiikan
laitokselle ja sen henkilökohtaisen paikannuksen algoritmien tutkimusryhmälle.
Kiitän työn ohjaajaa ja tarkastajaa professori Robert Pichéä ja TkT Simo Ali-Löyttyä
tutkimukseen ja raporttiin liittyvistä neuvoista ja kommenteista. Lisäksi haluan kiittää
muita paikannusryhmän jäseniä tähän työhön vaikuttaneista ajatuksista ja kommenteista.
Tampere, 21. kesäkuuta 2011
Henri Nurminen
IV
SISÄLLYS
1. Johdanto . . . . . . . . . . . . . . . . . . . . . . . . . .
2. Teoreettinen tausta . . . . . . . . . . . . . . . . . . . . .
2.1 Markovin ketjujen teoriaa . . . . . . . . . . . . . .
2.1.1 Markovin ketjun määritelmä . . . . . . . . . .
2.1.2 Invariantti jakauma . . . . . . . . . . . . . . .
2.1.3 Suurten lukujen laki ja keskeinen raja-arvolause
2.2 Metropolisin ja Hastingsin algoritmi . . . . . . . . .
2.2.1 Erilaisia ehdokasjakaumia . . . . . . . . . . . .
2.2.2 MH-algoritmin testaaminen . . . . . . . . . . .
3. Kuuluvuusalueen estimointi . . . . . . . . . . . . . . . .
3.1 Kuuluvuusalueen malli . . . . . . . . . . . . . . . .
3.2 Alkuarvot ja ehdokasjakaumat . . . . . . . . . . . .
3.3 Uskottavuusfunktio . . . . . . . . . . . . . . . . . .
3.4 Posteriorijakauma . . . . . . . . . . . . . . . . . .
3.5 Keskipisteen ja kovarianssimatriisin laskeminen . .
3.6 Posteriorin aikainvariantti päivittäminen . . . . . . .
4. Algoritmin testausta . . . . . . . . . . . . . . . . . . . .
4.1 Yleistä . . . . . . . . . . . . . . . . . . . . . . . .
4.2 Virheellisen havainnon tunnistus . . . . . . . . . . .
4.3 Havaintopisteen moninkertaistaminen . . . . . . . .
4.4 Posteriorin aikainvariantti päivittäminen . . . . . . .
5. Johtopäätökset . . . . . . . . . . . . . . . . . . . . . . .
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
1
4
4
4
5
6
8
11
11
12
12
12
15
16
17
18
20
20
20
22
24
25
V
MERKINNÄT JA SYMBOLIT
A∪B
A∩B
A\B
A⊂B
min A
max A
f :A→B
x
X
X=Y
T
x∈A
{a1 ,a2 , . . . ,an }
N
N0
R
Rn
∞
E
(ai )∞
i=0
(ai )∞
i=0 ⊂ A
(Xi )∞
i=0
(Xi )∞
i=0 ⊂ E
limn→∞ Xn
!n
i=0 Xi
¯
Xn
∀
∃
⇒
⇔
joukkojen A ja B leikkaus
joukkojen A ja B yhdiste
joukkojen A ja B erotus
joukko A on joukon B osajoukko
joukon A pienin alkio
joukon A suurin alkio
funktio, jonka määrittelyjoukko on A ja maalijoukko B
skalaari tai vekrori
satunnaismuuttuja tai satunnaisvektori
satunnaismuuttujalla X on nollamittaista joukkoa lukuun
ottamatta sama tiheysfunktio kuin satunnaismuuttujalla Y
matriisin tai vektorin transpoosi
x on joukon A alkio
alkioitten a1 ,a2 , . . . ,an muodostama joukko
luonnollisten lukujen joukko
N ∪{0}
reaalilukujen joukko
n-ulotteisten reaalivektorien joukko
ääretön
tila-avaruus
jono, jonka alkiot ovat a0 ,a1 , . . . .
jonon alkiot kuuluvat joukkoon A
stokastinen prosessi
stokastinen prosessi, jonka alkioitten
arvoalue on E
stokastisen prosessin (Xi )∞
i=0 raja-arvo
satunnaismuuttujien X0 ,X1 , . . . ,Xn summa
satunnaismuuttujien X0 ,X1 , . . . ,Xn keskiarvo
kaikkikvanttori
olemassaolokvanttori
implikaatio
ekvivalenssi
VI
P(X ∈ A)
P(X ∈ A | Y ∈ B)
E(X)
µπ
Σπ
X ∼π
Tas(a,b)
N(µ,σ 2 )
w ← N(µ,σ 2 )
eig(Σ)
ln
logit
-xH
S(f )
S(X)
"
f (x) dx
A
K(x,y)
π
q(x | y)
r(x,x# )
α(x,x# )
∝
p
$
y1 : n
XF 1 ,XF 2
XF 1,1 ,XF 1,2
S
G
todennäköisyys, että satunnaismuuttuja X kuuluu joukkoon A
todennäköisyys, että satunnaismuuttuja X kuuluu joukkoon A
ehdolla satunnaismuuttuja Y kuuluu joukkoon B
satunnaismuuttujan X odotusarvo
jakauman π odotuarvo
jakauman π kovarianssimatriisi
satunnaismuuttuja X noudattaa tiheysfunktion π
määrittelemää jakaumaa
tasajakauma välillä (a,b)
yksiulotteinen normaalijakauma odotusarvolla µ ja varianssilla σ 2
normaalijakaumasta generoitu luku w
matriisin Σ ominaisarvojen joukko
luonnollinen logaritmi
logit-funktio
euklidinen normi
Heavisiden funktio
funktion f kantaja
satunnaismuuttujan X kantaja
funktion f integraali skalaari- tai vektorimuuttujan x suhteen
joukon A yli
Markovin ketjun siirtymätodennäköisyys pisteestä x pisteeseen y
Markovin ketjun invariantti jakauma, Metropolisin ja Hastingsin
algoritmin kohdejakauma, posteriorijakauma
Metropolisin ja Hastingsin algoritmin ehdokasketjun
siirtymätodennäköisyys pisteestä y pisteeseen x
MH-suhde Metropolisin ja Hastingsin algoritmissa,
nykyinen arvo x, ehdokasarvo x#
ehdokasarvon hyväksymisen todennäköisyys Metropolisin
ja Hastingsin algoritmissa
suoraan verrannollisuus estimoitavien muuttujien suhteen
priorijakauma
uskottavuusfunktio
mittaukset numeroilla 1,2, . . . ,n
ellipsin polttopisteitä kuvaavat satunnaismuuttujat
ellipsin polttopisteen ens. ja toinen komponentti
ellipsin isoakselin pituutta kuvaava satunnaismuuttuja
mittauksen oikeellisuustodennäköisyyttä kuvaava
satunnaismuuttuja
1
1.
JOHDANTO
Nykyisin on käytössä lukuisia liikkuvissa päätteissä toimivia sovelluksia, joissa vaaditaan päätteen maantieteellisen paikan tuntemista. Liikkuvan päätteen paikannus suoritetaan yleensä mittaamalla satelliitin radiolähetyksiä. Esimerkiksi GPS (Global Positioning
System) on tällainen paikannusjärjestelmä. Jos tällaista järjestelmää ei ole käytettävissä
tai liikkuva pääte sijaitsee katvealueella, esimerkiksi sisätiloissa tai tiheään rakennetuilla taajama-alueilla, tarvitaan muita tapoja paikantaa pääte. Muuta kuin satelliittipohjaista
paikannusinformaatiota onkin nykyisin useimmiten saatavilla.
Yksi tapa paikantaa liikkuva pääte on kuuluvuusaluepaikannus, johon tämä työ liittyy. Kuuluvuusalue on maantieteellinen alue, jossa päätelaite kykenee havaitsemaan kuuluvuusalueeseen liittyvän tukiaseman lähettämän langattoman radiosignaalin. Käytännössä kuuluvuusaluepaikannuksessa voidaan käyttää esimerkiksi matkapuhelinverkkojen, tvtai radioverkkojen tai WLAN-verkkojen tukiasemia. Esimerkiksi matkapuhelimet ovat
usein yhteydessä puhelinverkon lisäksi WLAN-verkkoihin.
Tässä työssä ei oleteta kuuluvuusalueita tunnetuiksi, mutta niiden sijaintia ja kokoa
voidaan tietyssä mielessä havainnoida. Tukiasemien yksilöimiseen käytetään ID-numeroita (identity, yksikäsitteinen tunnus). Tietyllä kuuluvuusalueella sijaitsevan liikkuvan
päätteen vastaanottama radiosignaali sisältää kuuluvuusalueen ID:n, ja jos käytössä on tällöin jokin paikannusmenetelmä, pääte voi raportoida sijaintinsa ja havaitsemansa kuuluvuusalueen ID:n tietokantaa ylläpitävälle järjestelmälle. Näin kerätyistä sijainti–ID-vektoreista käytetään tässä työssä nimitystä sormenjälki (fingerprint), ja käytössä on siis suuren määrän sormenjälkiä sisältävä tietokanta.
Kuuluvuusalueitten sijainteja, kokoja ja muotoja ei siis voida tietää tarkasti, ja lisäksi
ne saattavat vaihdella ajan kuluessa jonkin verran esimerkiksi teknisistä taikka luonnon
tai ympäröivän maaston muutosten aiheuttamista syistä. Näin ollen kuuluvuusaluetta on
mielekästä mallintaa todennäköisyysjakaumalla, josta sormenjäljet ovat näytteitä. Kaikkien kuuluvuusalueitten sijainnit, koot ja muodot sisältävää tietokantaa kutsutaan radiokartaksi (radiomap). Tämän kandidaatintyön tarkoituksena on estimoida radiokartan kuuluvuusalueitten parametrit käyttäen havaintoina sormenjälkitietokantaa. Lisäksi halutaan
pystyä päivittämään olemassa olevaa radiokarttaa uusien havaintojen avulla tuntematta
koko sormenjälkitietokantaa.
Käytännössä kuuluvuusalueen jakauman mallinnusta rajoittavat muun muassa liikkuvan päätteen muistin määrä ja laskentakapasiteetti, joten kuuluvuusaluejakauman mal-
1. Johdanto
2
lin on oltava parametrimäärältään ja matemaattisilta ominaisuuksiltaan mahdollisimman
yksinkertainen. Ellipsi on määriteltävissä viidellä parametrilla, ja sillä voidaan approksimoida varsin yleisesti erityisesti kahden akselin suhteen symmetrisiä konvekseja alueita.
Kuuluvuusalueitten muodot eivät kaikissa tapauksissa ole tällaisia, mutta ne ovat silti varsin toimivia yleistyksiä. Ellipsin muotoista aluetta parhaiten mallintavia jakaumia ovat
esimerkiksi ellipsin muotoinen tasajakauma ja moniulotteinen normaali- tai t-jakauma.
Kuuluvuusaluehavaintojen lisäksi sormenjälkitietokanta saattaa sisältää virheellisiä havaintoja. Tämä voi näkyä muista havainnoista huomattavasti poikkeavien havaintojen
määränä, joka on suurempi, kuin sen pitäisi oletetun mallin mukaan olla. Muista havainnoista huomattavasti pokkeavat havainnot ovat usein vaikuttavia, eli ne saattavat häiritä
jakauman estimointia oleellisesti. Siksi voidaan ajatella, että virheelliset havainnot tulevat omasta erillisestä todennäköisyysjakaumastaan ja että havainnolla on tietty todennäköisyys olla virheellinen. Tässä työssä kuuluvuusaluetta mallinnetaan kaksiulotteisen
ellipsin alueelle tasajakautuneen jakauman sekä tietyn virhejakauman mikstuurina.
Koska tietyn kuuluvuusalue-ellipsin parametrien estimaateissa on epävarmuutta ja koska ne ovat lisäksi muuttuvia, ongelmaa on hyvin mielekästä lähestyä bayesläisen päättelyn
kautta. Bayesläisessä tilastotieteessä mittausten lisäksi koko taustalla oleva ilmiö nähdään
satunnaisprosessina ja tarkoituksena on estimoida kiinteitten parametriarvojen sijaan jakaumia. Tässä tapauksessa siis myös kuuluvuusalue-ellipsin sijainti-, koko- ja muotoparametreja tarkastellaan satunnaismuuttujina, joitten jakaumia pyritään estimoimaan.
Bayesläisessä päättelyssä estimoitavalle satunnaismuuttujalle on muodostettava priorijakauma. Priorijakauma sisältää kaiken satunnaismuuttujasta ennen estimointiprosessia saadun tiedon. Saaduista havainnoista lasketaan uskottavuusfunktion (likelihood) arvo, jonka tarkoitus on painottaa priorijakaumaa havaintojen mukaisesti. Kertomalla priori uskottavuudella saadaan posteriorijakauma, joka on bayesläinen ratkaisu estimoitavan
satunnaismuuttujan jakaumalle ehdolla saadut havainnot. Tässä työssä osoitetaan lisäksi,
että saataessa uusia havaintoja vanhaa posterioria voidaan käyttää priorina estimoitaessa uutta posterioria ja uskottavuutta laskettaessa tarvitsee tällöin huomioida vain uudet
mittaukset.
Estimointimenetelmänä tässä työssä käytetään Metropolisin ja Hastingsin algoritmia,
johon tästä eteenpäin viitataan lyhenteellä MH-algoritmi. Se on numeerinen, satunnaislukujen generointiin perustuva näytteistysalgoritmi, jolla siis näytteistetään jakaumaa eli
generoidaan jakauman mukaisesti jakautunut joukko lukuja tai vektoreita. Algoritmilla
voidaan periaatteessa estimoida mitä tahansa jakaumaa, kunhan sen tiheysfunktion arvot
voidaan laskea. Tiheysfunktio saa lisäksi olla normalisoimaton, eli sen ei tarvitse integroitua tai summautua ykköseksi, mikä on merkittävä etu, sillä bayesläisessä päättelyssä jakaumien analyyttisistä muodoista saattaa tulla hyvin monimutkaisia.
MH-algoritmi kuuluu Markovin ketjuun perustuviin Monte Carlo -menetelmiin (MCMCmenetelmät, Markov chain Monte Carlo). Markovin ketju on tietyt ehdot täyttävä satun-
1. Johdanto
3
naismuuttujien jono, ja MH-algoritmin käytännön toteutus perustuu tällaisen jonon pseudosatunnaislukujen avulla generoituun realisoitumaan. Pseudosatunnaisluvut ovat lukuja,
jotka generoidaan täysin deterministisillä algoritmeilla, joitten on osoitettu tuottavan määrättyä jakaumaa hyvin tarkasti noudattavia lukujonoja. Markovin ketjuja ja MCMC-menetelmiä on tutkittu paljon koneellisen laskentakapasiteetin nopean kasvun myötä, mutta tässä työssä esitellään vain MH-algoritmin kannalta keskeisin osa Markovin ketjujen
teoriaa. MH-algoritmin tuottama näyte ei ole korreloimaton, joten perinteisestä suurten
lukujen lauseesta johdetaan Markovin ketjuille oma versio.
4
2.
TEOREETTINEN TAUSTA
2.1
Markovin ketjujen teoriaa
2.1.1
Markovin ketjun määritelmä
Olkoon joukko E ⊂ Rn Borel-joukko, jota kutsumme tila-avaruudeksi.
Määritelmä 2.1.1. Tila-avaruuden E stokastinen prosessi on jono, jonka alkiot ovat satunnaismuuttujia, joiden arvojoukko on tila-avaruus E.
Merkitään stokastista prosessia (Xi )∞
i=0 . Tällöin satunnaismuuttujaa X0 kutsukaan stokastisen prosessin alkutilaksi ja satunnaismuuttujaa Xi stokastisen prosessin i:nneksi tilaksi. Jos tilaa Xi kutsutaan prosessin nykyiseksi tilaksi, niin tila Xi−1 on prosessin edellinen tila ja tila Xi+1 on prosessin seuraava tila.
Markovin ketju on stokastisen prosessin erikoistapaus, jossa ketjun nykyisen tilan jakauma riippuu aikaisemmista tiloista vain edellisen tilan kautta.
Määritelmä 2.1.2. Markovin ketju on stokastinen prosessi (Xi )∞
i=0 , jolle pätee
∀i ∈ N0 : (Xi+1 | X0 : i ) = (Xi+1 | Xi ),
(2.1)
missä satunnaismuuttujien yhtäsuuruus tarkoittaa, että niitten tiheysfunktiot ovat samat
nollamittaista joukkoa lukuun ottamatta.
Markovin ketju on aikahomogeeninen, jos ehdollisen satunnaismuuttujan Xi+1 | Xi jakauma on riippumaton indeksistä i.
Määritellään nyt jatkon kannalta tärkeä kantajan käsite. Satunnaismuuttujalle kantaja
on joukko, jossa on satunnaismuuttujan kaikki ”todennäköisyysmassa”.
Määritelmä 2.1.3. Olkoon funktion f : E → [0,∞) kantaja joukko
S(f ) = {x ∈ E|f (x) > 0}.
(2.2)
Kutsutaan lisäksi satunnaismuuttujan kantajaksi sen tiheysfunktion kantajaa.
Aikahomogeenisen Markovin ketjun määrittelevät alkutila X0 ja funktio K : E ×E →
2. Teoreettinen tausta
5
[0,∞), jolle pätee kaikilla indekseillä i
#
∀x ∈ S(Xi ) : ∀A ⊂ E : P(Xi+1 ∈ A | Xi = x) =
∀x ∈ E \ S(Xi ) : ∀y ∈ E : K(x,y) = 0
"
A
K(x,y) dy
.
(2.3)
Tällöin pätee
∀x ∈ S(Xi ) :
$
E
K(x,y) dy = P(Xi+1 ∈ E | Xi = x) = 1,
joten K(x,·) on tiheysfunktio. Funktiota K kutsutaan Markovin ketjun siirtymätodennäköisyydeksi eli ytimeksi.
2.1.2
Invariantti jakauma
Markovin ketjulla saattaa olla invariantti jakauma. Jos ketjun alkutila tällöin noudattaa
invarianttia jakaumaa, niin ketju säilyttää sen. Formalisoidaan tämä seuraavaksi määritelmäksi.
Määritelmä 2.1.4. Olkoon K aikahomogeenisen Markovin ketjun siirtymätodennäköisyys. Tällöin tiheysfunktio π : E → R on Markovin ketjun invariantti jakauma, jos pätee
$
∀y ∈ E :
π(x)K(x,y) dx = π(y).
(2.4)
E
Kääntyvä ketju säilyttää jakauman myös toiseen suuntaan.
Määritelmä 2.1.5. Aikahomogeeninen Markovin ketju on kääntyvä, jos on olemassa jakauma λ : E → R siten, että
∀x ∈ E : ∀y ∈ E : λ(x)K(x,y) = λ(y)K(y,x).
(2.5)
Tällöin sanotaan, että ketju on kääntyvä jakauman λ suhteen.
Todistetaan nyt MCMC-algoritmien kannalta keskeinen tulos.
Lause 1. Jakauman λ suhteen kääntyvällä Markovin ketjulla on invariantti jakauma, ja
sen invariantti jakauma on λ.
Todistus.
∀y ∈ E :
$
$
λ(x)K(x,y) dx =
λ(y)K(y,x) dx
E
$
= λ(y) K(y,x) dx = λ(y),
E
E
sillä K(y,·) on tiheysfunktio.
2. Teoreettinen tausta
2.1.3
6
Suurten lukujen laki ja keskeinen raja-arvolause
Suurten lukujen laki ja keskeinen raja-arvolause takaavat, että tietyin ehdoin invariantin
jakauman omaavan Markovin ketjun avulla voidaan estimoida kyseisen invariantin jakauman tunnuslukuja. Koska Markovin ketjujen satunnaismuuttujat eivät yleisesti ole riippumattomia, niille on johdettava omat versionsa näistä lauseista.
Selvästikin estimoivan Markovin ketjun tilojen on voitava saada arvoja mistä tahansa
invariantin jakauman kantajan osajoukosta. Tällaista ketjua sanotaan pelkistymättömäksi.
Määritelmä 2.1.6. Markovin ketju (Xi )∞
i=0 , jolla on invariantti jakauma π, on pelkistymätön jakauman π suhteen, jos
∀x0 ∈ E : ∀A ⊂ E : ∃m ∈ N :
$
A
π(t) dt > 0 ⇒ P(Xm ∈ A | X0 = x0 ) > 0. (2.6)
Suurten lukujen lakia varten ketjulle määritellään lisäksi seuraava ominaisuus.
Määritelmä 2.1.7. Markovin ketju (Xi )∞
i=0 , jolla on invariantti jakauma π, on Harrispalautuva jakauman π suhteen, jos sille pätee kaikilla alkutiloilla x0 ∈ E seuraava yhtälö:
$
∀B ⊂ E : ∃n ∈ N :
π(x) dx > 0 ⇒ P(Xn ∈ B | X0 = x0 ) = 1.
(2.7)
B
Pelkistyvyys luonnollisesti seuraa Harris-palautuvuudesta, mutta tavallisesti suurten
lukujen laissa esitetään molemmat ehdot, koska laki pätee tietyin erityisehdoin (ks. [4,
s.32–33]) myös ilman Harris-palautuvuutta. Kutsutaan funktiota f : E → R π-integroituvaksi, jos
$
E
f (x)π(x) dx < ∞,
missä π on tiheysfunktio. Nyt esitetään todistamatta tämän luvun päätulos.
Lause 2 (Suurten lukujen laki Markovin ketjuille). Olkoon (Xi )∞
i=0 aikahomogeeninen
Markovin ketju. Oletetaan, että seuraavat ehdot patevät:
1. Ketjulla on invariantti jakauma π.
2. Ketju on pelkistymätön jakauman π suhteen.
3. Ketju on Harris-palautuva jakauman π suhteen.
Tällöin jokaisella π-integroituvalla funktiolla f : E → R pätee
%
n
1 &
f (Xi ) =
∀x0 ∈ E : P lim
n→∞ n + 1
i=0
$
E
f (x)π(x) dx | X0 = x0
'
= 1,
(2.8)
eli funktion arvoista muodostetun ketjun aritmeettinen keskiarvo suppenee melkein varmasti kohti invarianttia jakaumaa π.
2. Teoreettinen tausta
7
Todistus. Sivuutetaan. Katso artikkeli [4, s.32–33].
p
Seuraus 3. Olkoon (Xi )∞
i=0 ⊂ R lauseen 2 ehdot toteuttava Markovin ketju. Merkitään
n
¯n =
X
1 &
Xi .
n + 1 i=0
Tällöin jos ketjun invariantilla jakaumalla π on odotusarvo µπ , niin
P
(
)
¯
lim Xn = µπ = 1.
(2.9)
n→∞
Todistus. Olkoon funktio fk : Rp → R muuttujavektorin k:nnen komponentin palauttava
funktio, ja olkoon X ∼ π satunnaismuuttuja. Tällöin lauseen 2.8 nojalla pätee
%
n
1 &
∀k ∈ {1,2, . . . ,p} : P lim
fk (Xi ) = E(fk (X))
n→∞ n + 1
i=0
'
= 1.
'
=1
Tästä päätellään, että pätee
%
n
1 &
P ∀k ∈ {1,2, . . . ,p} : lim
fk (Xi ) = E(fk (X))
n→∞ n + 1
i=0
eli


 

f
(X
)
E(f
(X))
1
i
1
n
1 &

 

..
..
P  lim

=
 = 1.
.
.
n→∞ n + 1
i=0
fp (Xi )
E(fp (X))
Tämä on ekvivalentti väite yhtälön (2.9) kanssa.
p
Seuraus 4. Olkoon (Xi )∞
i=0 ⊂ R lauseen 2 ehdot toteuttava Markovin ketju. Merkitään
n
Sn2
1 &
¯nX
¯ nT .
=
Xi XiT − X
n + 1 i=0
Tällöin jos ketjun invariantilla jakaumalla π on kovarianssimatriisi Σπ , niin
P
(
lim
n→∞
Sn2
)
= Σπ = 1.
Todistus. Olkoon funktio fij : Rp×p → R muuttujamatriisin i:nnen rivin ja j:nnen sarakkeen alkion palauttava funktio, ja olkoon satunnaismuuttuja X ∼ π. Tällöin lauseen 2.8
nojalla pätee
%
n
1 &
fij (Xk XkT ) = E(fij (XX T ))
∀i,j ∈ {1,2, . . . ,p} : P lim
n→∞ n + 1
k=0
'
= 1,
2. Teoreettinen tausta
8
joten edelleen pätee
%
n
1 &
fij (Xk XkT ) = E(fij (XX T ))
P ∀i,j ∈ {1,2, . . . ,p} : lim
n→∞ n + 1
k=0
'
= 1.
Tästä päätellään, että väite
%
n
1 &
P lim
Xk XkT = E(XX T )
n→∞ n + 1
k=0
'
= 1.
on tosi. Näin ollen edellisen seurauksen nojalla pätee
%
'
n
1 &
¯ n = µπ
1 = P lim
Xi XiT = E(XX T ), lim X
n→∞ n + 1
n→∞
i=0
%
7
'
6
n
1 &
¯nX
¯ nT = E(XX T ) − µπ µT
Xi XiT − X
= P lim
π
n→∞ n + 1
i=0
(
)
= P lim Sn2 = Σπ .
n→∞
Käytännössä ei voida konstruoida äärettömän pitkää Markovin ketjua, mutta esimerkiksi artikkelissa [3, Luku 4] todistetusta Markovin ketjujen keskeisestä raja-arvolauseesta saadaan arvio estimaatin tarkkuudelle. Tarkkojen virhe-estimaattien laskemisen kannalta arvio ei ole yhtä käyttökelpoinen kuin riippumattomien satunnaismuuttujien tapauksessa, mutta se osoittaa, ettei tila-avaruuden dimensio vaikuta ketjun suppenemisnopeuteen
vaan ketjun korreloituneisuus.
2.2
Metropolisin ja Hastingsin algoritmi
MCMC-algoritmit eli Markovin ketjuun perustuvat Monte Carlo -algoritmit ovat satunnaisalgoritmeja, joilla simuloidaan tunnettua jakaumaa muista jakaumista generoitujen
satunnaislukujen avulla muodostetulla Markovin ketjulla. Tunnetuin MCMC-algoritmi on
Metropolisin ja Hastingsin algoritmi.
Kutsutaan kohdejakaumaksi sitä jakaumaa, jota MH-algoritmilla halutaan simuloida.
Olkoon kohdejakauman tiheysfunktio π : E → R. Olkoon ehdokasketju sellainen Markovin ketju (Qi )∞
i=1 ⊂ E, että siitä voidaan generoida pseudosatunnaisia realisoitumia.
Ehdokasketjun siirtymätodennäköisyyttä pisteestä y pisteeseen x merkitään q(x | y) ja jakaumaa q(· | y) kutsutaan ehdokasjakaumaksi.
Tällöin MH-algoritmi Markovin ketjun generoimiseen on seuraava: Olkoon meillä ketjun nykyinen tila x ∈ S(π). Generoidaan ehdokasjakaumasta q(· | x) luku x# , jota kutsu-
2. Teoreettinen tausta
9
taan ehdokasarvoksi. Kutsutaan nyt MH-suhteeksi lukua
r(x,x# ) =
π(x# )q(x | x# )
aπ(x# )q(x | x# )
=
,
π(x)q(x# | x)
aπ(x)q(x# | x)
(2.10)
missä a on mielivaltainen nollasta poikkeava vakio. Koska a on mielivaltainen, tiheysfunktiot π ja q voivat olla normalisoimattomia. Generoidaan sitten luku u standarditasajakaumasta Tas(0,1). Jos pätee r(x,x# ) > u, niin sanotaan, että ehdokasarvo hyväksytään,
ja asetetaan ketjun seuraavan tilan arvoksi ehdokasarvo. Muutoin ehdokasarvo hylätään
ja asetetaan ketjun seuraavan tilan arvoksi nykyinen tila. Havaitaan, että MH-algoritmissa
ehdokasarvon hyväksymisen todennäköisyys on
#
#
α(x,x ) = P(x hyväksytään | x) = min
8
9
π(x# )q(x | x# )
,1 .
π(x)q(x# | x)
(2.11)
MH-algoritmi tuottaa selvästi Markovin ketjun, sillä seuraava tila riippuu ketjun aikaisemmista tiloista vain nykyisen tilan kautta. Ehdokasjakauman siirtymätodennäköisyys ei
muutu algoritmin suorituksen aikana, joten Markovin ketju on aikahomogeeninen. Todistetaan seuraavaksi, että kohdejakauma π on MH-algoritmin tuottaman Markovin ketjun
invariantti jakauma.
Apulause 5. MH-algoritmin tuottama Markovin ketju on kääntyvä kohdejakauman π suhteen.
Todistus. Kun A ja B ovat tila-avaruuden E mielivaltaisia osajoukkoja, ketjun siirtymistodennäköisyys näitten joukkojen välillä on
$ $
K(A,B) =
K(x,y) dx dy
B
A
= P(Xi+1 ∈ B | Xi ∈ A)
= P(Xi+1 ∈ B, hyväksytään | Xi ∈ A) + P(Xi+1 ∈ B, hylätään | Xi ∈ A)
= Khyv (A,B) + Khyl (A,B).
Jos ehdokasarvoa ei hyväksytä, ketjun seuraavan tilan arvo on sama kuin nykyinen, joten
selvästi yhtälö (2.5) pätee. Tarkastellaan nyt tapausta, jossa ehdokasarvo hyväksytään.
Khyv (A,B) = P(Xi+1 ∈ B, hyväksytään | Xi ∈ A)
= P(Xi# ∈ B, hyväksytään | Xi ∈ A)
$ $
=
q(y | x)α(x,y) dx dy,
B
A
missä Xi# on i:nnettä tilaa vastaava ehdokasarvo. Tästä päätellään, että siirtymistodennäköisyystiheys hyväksyttäessä on Khyv (x,y) = q(y | x)α(x,y). Yhtälön (2.5) vasen puoli
2. Teoreettinen tausta
10
voidaan siis kirjoittaa
π(x)Khyv (x,y) = π(x)q(y | x)α(x,y)
9
8
π(y)q(x | y)
,1
= π(x)q(y | x) min
π(x)q(y | x)
= min {π(y)q(x | y), π(x)q(y | x)} .
Vastaavasti voidaan osoittaa, että pätee
π(y)Khyv (y,x) = min {π(x)q(y | x), π(y)q(x | y)}
= min {π(y)q(x | y), π(x)q(y | x)} ,
joten
π(x)Khyv (x,y) ≡ π(y)Khyv (y,x),
eli MH-algoritmin tuottama Markovin ketju on kääntyvä kohdejakauman π suhteen.
MH-algoritmin määritelmän perusteella on selvää, että MH-algoritmiin liittyvä Markovin ketju on pelkistymätön, jos ehdokasjakauman kantajajoukko sisältää kohdejakauman kantajajoukon eli ehdokasjakaumasta voidaan generoida kaikki tila-avaruuden arvot.
Muotoillaan tämä lauseeksi.
Apulause 6. MH-algoritmin tuottama Markovin ketju on pelkistymätön invariantin jakaumansa π suhteen, jos
∀x ∈ E : S(π) ⊂ S(q(· | x)),
(2.12)
missä q on ehdokasketjun siirtymätodennäköisyys.
Yleisesti Harris-palautuvuus on pelkistymättömyyttä vahvempi ehto, mutta seuraavan
lauseen mukaan MH-algoritmin tapauksessa pelkistymättömyys on riittävä ehto Harrispalautuvuudelle.
Apulause 7. Jos MH-algoritmin tuottama Markovin ketju on pelkistymätön invariantin
jakaumansa π suhteen, niin se on Harris-palautuva jakauman π suhteen.
Todistus. Sivuutetaan. Katso esimerkiksi artikkeli [5, Korollaari 2].
Näin olleen voidaan esittää suurten lukujen laki MH-algoritmille.
Lause 8. MH-algoritmin tuottama Markovin ketju toteuttaa lauseen 2 ehdot, eli sille pätee suurten lukujen laki (2.8).
2. Teoreettinen tausta
2.2.1
11
Erilaisia ehdokasjakaumia
Usein ehdokasjakauma on symmetrinen ehdokas- ja nykyisen arvon suhteen eli
q(x|x# ) = q(x# |x),
(2.13)
π(x# )
.
π(x)
(2.14)
jolloin MH-suhde on
r(x,x# ) =
Jos MH-algoritmin tuottaman Markovin ketjun seuraava tila on riippumaton nykyisestä tilasta ja siis kaikista aikaisemmista tiloista, MH-algoritmia kutsutaan riippumattomaksi. Satunnaiskävelyyn perustuvassa MH-algoritmissa taas generoidaan Markovin ketjun
tilasta ja ajanhetkestä riippumattomasta jakaumasta näyte w ja asetetaan ehdokasarvoksi
x# = x + w.
Jos satunnaisluvun w generoiva jakauma on symmetrinen origon suhteen, on ehdokasjakauma symmetrinen.
2.2.2
MH-algoritmin testaaminen
MH-algoritmin testaamiseksi simuloidusta näytteestä on syytä piirtää näytehistoria, jossa on kuvaajat kustakin komponentista Markovin ketjun ajanhetken funktiona. Esimerkki
näytehistoriasta on kuvassa 4.1. Jos kuvaajissa on suuria hyppäyksiä ja ehdokasarvojen
hyväksymisaste on alhainen, algoritmissa on paljon turhaa toistoa. Jos taas hyppäykset
ovat hyvin pieniä ja näytteessä on paljon autokorrelaatiota eli kuvaaja näyttää aaltoilevan, näytteitten lukumäärän on oltava hyvin suuri, jotta simulointi todella vastaisi oikeata
jakaumaa. Varsinkin monihuippuinen jakauma saattaa tällöin tulla simuloiduksi puutteellisesti, kun algoritmi jää pitkäksi aikaa kiinni yhteen tiheysfunktion huippuun eli moodiin. Satunnaiskävelyssä nämä ongelmat voivat liittyä esimerkiksi liian suuriin tai pieniin
ehdokasjakauman variansseihin.
Teoriassa algoritmi simuloi kohdejakaumaa riittävän suurella toistomäärällä riippumatta Markovin ketjun alkuarvosta. Käytännössä ketjun alkutila on jokin määrätty tai
laskettu piste eikä siis kohdejakauman mukaisesti jakautunut satunnaismuuttuja, minkä
vuoksi ketjun alussa saattaa olla muusta ketjusta selvästi poikkeavia arvoja, joista käytetään nimitystä burn-in. Näytehistorioita tutkimalla tulee selvittää, milloin ketju alkaa
simuloida kohdejakaumaa, ja sitten poistaa burn-in lopullisesta näytteestä. [2, Kappale
13.2]
12
3.
KUULUVUUSALUEEN ESTIMOINTI
3.1
Kuuluvuusalueen malli
Kuuluvuusalueet sijaitsevat kaksiulotteisessa koordinaatistossa. Mallinnetaan kuuluvuusaluetta jatkuvalla jakaumalla, jonka tiheysfunktion arvo kuvaa havainnon saamisen todennäköisyyttä kustakin pisteestä. Koska kuuluvuusalueita ei tunneta tarkasti, jakauman parametrit ovat satunnaismuuttujia.
Kuuluvuusalue on mikstuurijakauma f : R2 → R+ , joka on muotoa
f = G · foikeelliset + (1 − G) · fvirheelliset .
(3.1)
Jakauma foikeelliset on ellipsin muotoiselle alueelle jakautunut tasajakauma. Sitä parametrisoivat satunnaismuuttujat ovat polttopisteet XF 1 = [XF 1,1 ; XF 1,2 ] ja XF 2 = [XF 2,1 ; XF 2,2 ]
sekä isoakselin pituus S. Jakauma fvirheelliset on ympyrän muotoinen tasajakauma. Sitä parametrisoivat luvut ovat keskipiste aC sekä säde a, joita ei yksinkertaisuuden vuoksi estimoida vaan niille lasketaan jollakin deterministisellä menetelmällä sellaiset arvot, että
kaikki havaintopisteet ovat ympyrän sisäpuolella. Näin määriteltyä ympyrää kutsutaan
virheympyräksi. Lisäksi satunnaismuuttuja G kuvaa yksittäisen havainnon todennäköisyyttä olla oikeellinen. Olkoon X ellipsin muotoparametreista muodostettu satunnaisvektori


XF 1,1


X =  XF 1,2  .
S
Kuuluvuusaluejakaumasta saadaan joukko vektoreita, joita kutsutaan havainnoiksi tai sormenjäljiksi. Näitten perusteella pyritään estimoimaan MH-algoritmilla parametrien XF 1 ,
XF 2 , S ja G jakaumia. Tila-avaruus on siis E = R4 ×[0,∞) × [0,1]. Parametrien aC ja a
on oltava sellaiset, että kaikki havainnot ovat virheympyrän sisäpuolella.
3.2
Alkuarvot ja ehdokasjakaumat
Asetetaan ellipsin molemmiksi polttopisteiksi mediaanihavainnon piste, jotta mahdollisten virheellisten havaintojen vaikutus alkuarvoihin minimoituisi. Isoakselin pituuden S
alkuarvoksi s0 asetetaan kaksi kertaa havaintojen keskimääräinen etäisyys polttopisteistä.
Kuuluvuusalueen kovarianssimatriisista voi olla käytössä prioriestimaattina matriisi
3. Kuuluvuusalueen estimointi
13
ˆ − , jonka voidaan ajatella kuvaavan oletetun priorijakauman moodia. Tässä työssä prioria
Σ
käytetään vain ketjun alkuarvojen generoimiseen.
Kirjassa [1, Kappale 4.2] osoitetaan, että kaksidimensioisessa tapauksessa yhtälö
(x − µ)T Σ−1 (x − µ) = 1,
(3.2)
missä Σ on symmetrinen positiivisesti definiitti matriisi, määrittelee µ-keskisen ellipsin.
Tällöin matriisilla Σ on similaarimuunnos
6
76
7
:
; a2 0
dT
Σ= d n
,
(3.3)
0 b2
nT
missä d on ellipsin isoakselin suuntainen yksikkövektori, n vektoria d vastaan kohtisuora
yksikkövektori, a isopuoliakselin pituus ja b pikkupuoliakselin pituus. Pituuksien neliöt a2
ja b2 ovat siis matriisin Σ ominaisarvot ja vektorit d ja n niitä vastaavat ominaisvektorit,
jotka ovat ortonormaalit. Yhtälön (3.2) mielessä siis jokaista ellipsiä vastaa matriisi Σ,
jota kutsutaan usein ellipsin kovarianssimatriisiksi.
Olkoon meillä ellipsi ja sen sisäpuolelle tasan jakautunut jakauma. Käytetään tarkasteltavan ellipsin parametreille äskeisiä merkintöjä. Yhtälön (3.2) mielessä yksikköympyrän
kovarianssimatriisi on yksikkömatriisi I, ja tarkastelemamme ellipsi saadaan yksikköympyrästä lineaarimuunnoksella
Tx =
:
d n
;
6
a 0
0 b
7
x,
sillä
6
76
7
6
7
:
; 1/a2
T
a
0
d
0
(T x)T Σ−1 T x = xT
d n
0 b
nT
0
1/b2
6
7
6
7
; a 0
dT :
·
x
d
n
nT
0 b
6
76
76
7
2
a
0
1/a
0
a
0
= xT
x
0 b
0
1/b2
0 b
= xT I −1 x.
Näin ollen lineaarimuunnos T muuntaa yksikköympyrän alueelle tasajakautuneen satunnaismuuttujan tarkastellun ellipsin alueelle tasajakautuneeksi satunnaismuuttujaksi. Tulkitaan nyt T äskeistä lineaarimuunnosta vastaavaksi matriisiksi. Yksinkertaisella laskutoimituksella voidaan osoittaa, että yksikköympyrän alueelle tasajakautuneen satunnaismuuttujan varianssi on 1/4 · I, joten koska ellipsimuunnnos T on lineaarinen, kaksiulot-
3. Kuuluvuusalueen estimointi
14
teisen ellipsin muotoisen tasajakauman kovarianssimatriisi on
1
ΣTas = T · IT T
4
6
76
7
:
; a2 0
dT
1
=
,
d n
4
0 b2
nT
(3.4)
Yhtälön (3.4) nojalla isoakselin pituuden alkuarvoksi aseteˆ − perusteella
taan prioriestimaatin Σ
s0 = 2
<

ˆ − )),
max(eig(4Σ



 
ˆ − ) on priorin kovarianssimatriisin ominaisarvojen
missä eig(Σ
joukko. Kuvasta 3.1 nähdään, että polttopisteille saadaan vastaa- Kuva 3.1: Ellipsin poltvasti alkuarvot yhtälöstä
topisteitten laskeminen
xF 1,F 2,0 = µ ± yd
<
ˆ − )) − min(eig(4Σ
ˆ − )) d,
= µ ± max(eig(4Σ
ˆ − suurinta ominaisarvoa vastaava normalisoitu ominaisvektori.
missä d on matriisin Σ
Ehdokasarvot simuloidaan riippumattomista normaalijakaumista, joitten odotusarvoina käytetään parametrien viimeksi hyväksyttyjä arvoja ja joitten keskihajonnat ovat vakioita. Parametreilla S ja G on normaalijakauman arvoalueesta poikkeavat arvoalueet,
joten stokastisina muuttujina käytetään muuttujia
ln(S − -x#F 1 − x#F 2 -) ja logit(G) = ln
G
,
1−G
missä x#F 1 ja x#F 2 ovat polttopisteille generoidut ehdokasarvot. Tällöin siis ehdokasarvon
s# generoiva jakauma riippuu nykyisen tilan lisäksi aiemmin generoiduista ehdokasarvoista x#F 1 ja x#F 2 . Tämä ratkaisu on tehokkaampi nykyisiä polttopisteitä käyttävään verrattuna, koska aina hylättäviä arvoja, joissa isoakseli on lyhyempi kuin polttopisteitten välimatka, ei generoida. Kyseessä ei ole puhdas satunnaiskävely, sillä parametrit generoidaan
lausekkeilla
2
s# = (s − -x#F 1 − x#F 2 -)ew1 + -x#F 1 − x#F 2 - , w1 ← N(0,σln
s ),
1
2
g# =
1−g −w2 , w2 ← N(0,σlogit(g) ).
1+ g e
Koska satunnaisluvut wi ovat origon suhteen symmetrisesti jakautuneet, nykyisten ja ehdokasarvojen paikkoja yllä olevissa yhtälöissä voidaan vaihtaa, joten molemmat ehdokasjakaumat ovat symmetrisiä yhtälön (2.13) mielessä.
3. Kuuluvuusalueen estimointi
15
Algoritmin toiminnan kannalta ehdokasjakaumien valinta on hyvin kriittinen. Ehdokaskeskihajonnat on kokeilemalla asetettu sellaisiksi, että algoritmi toimisi useimmissa
tapauksissa mahdollisimman hyvin. Jotta ehdokaskeskihajontojen arvioinnissa voidaan
huomioida sovitettavien ellipsien mahdolliset suuruusluokkaerot, käytetyt keskihajonnat
on asetettu riippumaan havaintopisteitten otoskeskihajonnasta. Sovelluksissa sen laskeminen ei välttämättä ole rajoittuneen laskentakapasiteetin vuoksi mahdollista, joten tässäkin voitaisiin mahdollisesti hyödyntää prioritietoa. Ehdokaskeskihajonnoille on asetettu
minimi, jotta algoritmi kävisi läpi koko tila-avaruuden myös hyvin pienillä datamäärillä.
Minimi on pyritty asettamaan siten, että algoritmi toimisi muutamalla lähekkäin sijaitsevalla havaintopisteellä. Koska ei-degeneroituneen normaalijakauman arvoalue on koko
avaruus, apulauseen 6 ehto on voimassa.
Kaikki komponentit päivitetään kerralla. Erilaisia komponenteittaisen päivityksen mahdollisuuksia on tutkittu, mutta ne lisäävät merkittävästi ohjelman ajankulutusta, eikä niitten kuitenkaan todettu oleellisesti parantavan radiokartan tarkkuutta.
3.3
Uskottavuusfunktio
Alkeisgeometriasta tiedetään ellipsin pinta-ala isoakselin pituuden s sekä polttopisteitten
xF 1 ja xF 2 funktiona.
A = π$isopuoliakseli $pikkupuoliakseli
=
?2
>
-xF 1 − xF 2 s ( s )2
−
=π
2
2
2
<
s s2 − -xF 1 − xF 2 -2
=π
.
4
Lisäksi tiedetään, että piste y = [yx ,yy ]T on ellipsin sisäpuolella, jos ja vain jos
-y − xF 1 - + -y − xF 2 - < s.
Näin ollen yksittäisen kuuluvuusalueen oikeellisten havaintojen uskottavuusfunktio on
foikeelliset (y | x) =
H(s − -y − xF 1 - − -y − xF 2 -)
√
,
s s2 −&xF 1 −xF 2 &2
π
4
(3.5)
missä H tarkoittaa Heavisiden funktiota. Funktion muuttuja x on ellipsiparametreista
koostuva vektori. Virheellisten jakauma on ympyrä, jonka sisällä kaikki havainnot ovat.
Koska virheympyrä sovitetaan siten, että kaikki havinnot ovat sen alueella, virheellisten
3. Kuuluvuusalueen estimointi
16
havaintojen uskottavuusfunktio on
fvirheelliset (y | x) =
1
,
πa2
(3.6)
Tässä mallissa siis oletetaan, että myös kuuluvuusalueella voi kaavan (3.7) mukaisesti
olla virheellisiä havaintoja.
Kun siis oletetaan tunnetut havaintopisteet toisistaan riippumattomiksi ehdolla ellipsiparametrit, havaintojoukon y1 : n uskottavuusfunktioksi saadaan
$(y1 : n |x,g)
n
@
=
(gfoikeelliset (yi | x) + (1 − g)fvirheelliset (yi | x))
(3.7)
i=1
=
n
@


− xF 1 - − -yi − xF 2 -)
1
g H(s − -yi √
+ (1 − g) 2 
2
2
πa
s s −&xF 1 −xF 2 &
i=1
π
4
virheympyrän alueella, muualla uskottavuus on nolla. Uskottavuusfunktion arvot ovat
käytännössä hyvin pieniä, joten kertolasku johtaa liukuluvuilla laskettaessa yleensä alivuotovirheeseen. Alivuodon välttämiseksi käytetään tämän funktion logaritmeja. Tämä
suurentaa lukujen itseisarvoja ja muuttaa kertolaskun yhteenlaskuksi.
Virheympyrän parametrien määräämiseksi lasketaan pienin koko havaintopistejoukon
sisältävä sivuiltaan koordinaattiakselien suuntainen suorakulmio, ja käytetään virheympyrän säteenä a tämän suorakulmion lävistäjän puolikasta ja keskipisteenä aC tämän puoliväliä. Asetetaan säteelle lisäksi minimi. Hyvin pienellä virheympyrän säteellä uskottavuusfunktion arvot ovat hyvin suuria, jolloin virheympyrän vaikutus koko kuuluvuusaluejakaumaan on suhteettoman suuri ja oikeellisten havaintojen jakauman siis pieni. Tällä
menetelmällä valittu virheympyrä ei myöskään kasva havaintopisteitten hajontaan nähden
suhteettoman suureksi, jolloin virheellisten jakauman tiheys olisi lähellä nollaa. Tällöin
uskottavusfunktio (3.7) ei siis juurikaan painottaisi virheellisiä havaintoja ja ulkolaiset
voisivat häiritä estimointia.
3.4
Posteriorijakauma
Jos tietystä kuuluvuusalueesta ei ole olemassa aiemmin laskettua radiokarttaa, kuuluvuusalueen sijainnista ei ole prioritietoa. Halutaan mahdollisimman pieni oikeelliset havainnot sisältävä ellipsi. Valitaan siis priorijakaumaksi sijaintiparametreille tasajakauma
ja muotoparametrin logaritmille tasajakauma, joka johtaa muotoparametrin käänteisluvulla kertomiseen. Valitaan siis ellipsiparametrien prioriksi
1
p(x,g) = p(xF 1 ,xF 2 ,s,g) = .
s
(3.8)
3. Kuuluvuusalueen estimointi
17
Priori ei normalisoidu, mutta tällä ei ole merkitystä, koska uskottavuusfunktion kantaja on
rajoitettu joukko. Testattiin myös prioria, joka pyrki pienentämään polttopisteitten välistä
etäisyyttä mikä johti kuitenkin helposti ympyränmuotoiseen kuuluvuusalueeseen. Informatiivisemman priorin käyttö on ilmeisesti perusteltua vain, jos radiokartasta on olemassa
aiemmin laskettua tietoa.
Estimoitava posteriorijakauma on Bayesin säännön nojalla
π(x,g | y1 : n ) ∝ $(y1 : n | x,g)p(x,g),
(3.9)
joten MH-suhde on tässä tapauksessa
r(x# ,x) =
3.5
s
exp (ln $(y1# : n |x# ,g # ) − ln $(y1 : n |x,g)) .
s#
(3.10)
Keskipisteen ja kovarianssimatriisin laskeminen
Muotoillaan analyyttinen lauseke keskipisteen eli kuuluvuusaluejakauman odotusarvon ja
kovarianssimatriisin estimaattoreille olettaen, että tunnetaan estimaattorit ellipsiparametreille. Olkoon nyt Z kuuluvuusaluejakauman mukaisesti jakautunut satunnaismuuttuja.
Sen ellipsiparametrien suhteen ehdollisen jakauman odotusarvo eli keskipisteen estimaattori on
µ = E(Z | X,G)
$
= zf (z | X,G) dz
$
= z(Gfoikeelliset (z | XF 1 ,XF 2 ,S) + (1 − G)fvirheelliset (z)) dz
$
$
= G zfoikeelliset (z | XF 1 ,XF 2 ,S) dz + (1 − G) zfvirheelliset (z) dz
(3.11)
= Gµoikeelliset + (1 − G)µvirheelliset
XF 1 + XF 2
=G·
+ (1 − G) · aC .
2
Ennustejakauman kovarianssimatriisille voidaan johtaa seuraava lauseke.
Σ = var(Z)
= E(ZZ T ) − µµT
$
= zz T (Gfoikeelliset (z) + (1 − G)fvirheelliset (z)) dz − µµT
$
$
T
= G zz foikeelliset (z) dz + (1 − G) zz T fvirheelliset (z) dz − µµT
= G(Σoikeelliset + µoikeelliset µT
oikeelliset )
T
+ (1 − G)(Σvirheelliset + µvirheelliset µT
virheelliset ) − µµ .
(3.12)
3. Kuuluvuusalueen estimointi
18
Kaavan (3.4) nojalla ellipsin alueelle tasajakautuneen satunnaismuuttujan kovarianssimatriisi on
6
 S
7
0
:
;
T
2 A
A
1
,
(3.13)
Σoikeelliset =
A B 
B S C2 ( &XF 1 −XF 2 & )2 
T
4
B
0
−
2
2
missä
XF 1 − XF 2
, ja
-XF 1 − XF 2 6
7
A2
B =
.
−A1
A =
Matemaattisesti perustelluin estimaatti kuuluvuusalueelle saataisiin laskemalla odotusarvo ja kovarianssimatriisi jokaisella MH-algoritmin kierroksella ja ottamalla lopuksi
keskiarvot, mutta laskenta-ajan säästämiseksi estimaatit lasketaan MH-algoritmin simuloimien parametrinäytteitten keskiarvoista.
3.6
Posteriorin aikainvariantti päivittäminen
Olemassa olevaa radiokarttaa halutaan päivittää uusilla havaintopisteillä yn+1:m tuntematta edellisiä havaintopisteitä y1:n ja tallettamalla päivitysten välillä mahdollisimman vähän lukuja. Oletetaan uusien havaintojen tulevan samasta jakaumasta kuin aikaisempienkin havaintopisteitten (aikainvarianttius). Tällöin seuraavan lauseen nojalla olemassa olevaa posteriorijakaumaa voidaan käyttää priorijakaumana uutta posterioria estimoitaessa.
Symboli ∝ tarkoittaa suoraan verrannollisuutta estimoitavien parametrien suhteen.
Lause 9.
π2 (x,g | y1 : m ) ∝ $(yn+1 : m | x,g)π1 (x,g | y1 : n ).
(3.14)
Todistus. Merkitään tässä kirjaimella p aina niitten satunnaismuuttujien tiheysfunktiota, joihon muuttujasymboleilla kulloinkin viitataan. Käytetään ehdollisen tiheysfunktion
määritelmää:
p(x,g | y1 : m ) = p(x,g | y1 : n ,yn+1 : m )
p(x,g,y1 : n ,yn+1 : m )
p(y1 : n ,yn+1 : m )
p(yn+1 : m | x,g,y1 : n )p(x,g | y1 : n )p(y1 : n )
=
.
p(yn+1 : m | y1 : n )p(y1 : n )
=
3. Kuuluvuusalueen estimointi
19
Mittaukset oletetaan toisistaan riippumattomiksi ehdolla jakauman parametrit, joten pätee
p(x,g | y1 : m ) =
p(yn+1 : m | x,g)p(x,g | y1 : n )
p(yn+1 : m | y1 : n )
∝ p(yn+1 : m | x,g)p(x,g | y1 : n ).
Päivitetty posteriori voidaan siis laskea aikaisemmin esiteltyä vastaavalla MH-algoritmilla, jossa priorina on edellinen posteriori π1 (x,g | y1 : n ). Käytetään edelleen uskottavuusfunktiona $(yn+1 : m | x,g) kaavan (3.7) mukaista tasajakautuneen ellipsin ja tasajakautuneen virheympyrän mikstuuria. Virheympyrän parametrit lasketaan jokaisella päivityskerralla uudesta pistejoukosta.
Laskennan yksinkertaistamiseksi ja radiokarttaan tallennettavien parametrien määrän
pienentämiseksi oletetaan, että parametrien XF 1,1 , XF 1,2 , XF 2,1 , XF 2,2 , S ja G yhteisjakauma on normaali. Tällöin prioritiheyksien laskemiseen voidaan käyttää yleisimpiin
tilasto-ohjelmistoihin kuuluvaa moniulotteisen normaalijakauman tiheysfunktiota. Jotta
päivittäminen onnistuisi, priorista on odotusarvona käytettävän vanhan radiokarttaellipsin lisäksi tunnettava radiokarttaparametrien kovarianssimatriisi. Seurauksen 4 nojalla parametrien varianssia voidaan estimoida MH-algoritmin näytehistorioista lasketulla otoskovarianssimatriisilla. Sen laskeminen on yksinkertaista, mutta symmetrisessä 6 × 6matriisissa on 21 vapausastetta, mikä on varsin paljon. Siksi oletetaan parametrit toisistaan riippumattomiksi. Tämä vähentää vapausasteitten määrän kuuteen mutta ei testien
perusteella näytä vaikuttavan tuloksiin merkittävästi. Yksi mahdollisuus voisi olla vakiokovarianssin käyttäminen; päivitysjakauman muotoparametrin karkeahkokaan approksimointi ei välttämättä heikennä tulosta oleellisesti.
20
4.
ALGORITMIN TESTAUSTA
Tässä kandidaatintyössä esiteltävän algoritmin toteuttavat MATLAB-funktiot ovat saatavilla osoitteessa http://math.tut.fi/posgroup/nurminen2011a_koodit.html samoin kuin tässä kappaleessa esiteltyjen testien toteutukset.
4.1
Yleistä
Monimutkainen kuusiulotteinen jakauma vaatii MCMC-algoritmilta varsin paljon. Teoriassa ongelman dimensio ei vaikuta ketjun suppenemiseen, mutta käytännössä jakauman
monimutkaisuus vaikuttaa huomattavasti ketjun estimointiominaisuuksiin.
Vaikka ehdokasjakauma asetettiin riippuvaksi saadun näytteen koosta, optimaalisen
ehdokasjakauman varianssin kaava vaikutti riippuvan tehtävän mittasuhteista. Ehdokasjakaumien varianssien oli oltava suhteellisen pieniä, jotta hyväksymisprosentin pieneneminen ei tekisi algoritmista liian tehotonta. Tällöin esiintyi kuitenkin yksittäistapauksia,
joissa tulosellipsi ei silmämääräisesti tarkastellen selvästikään ollut optimaalinen. Näissä
tapauksissa voitiin siis tulkita, että parametrijakaumassa oli useita huippuja, ja algoritmi
oli juuttunut ”väärään” huippuun.
Pienet ehdokasvarianssit aiheuttivat monissa tapauksissa varsin voimakkaan burn-inilmiön; usein näytehistoriakuvaajan tasaantuminen vaati jopa noin tuhat iteraatiokierrosta.
Ongelmaa pyrittiin korjaamaan kehittämällä ketjun alkuarvojen generointia. Algoritmin
toimintakyvyn havaittiinin olevan varsin herkkä ketjun alkuarvojen suhteen.
Joissain tapauksissa hyväksymistodennäköisyydet olivat varsin pienistä ehdokasjakauman variansseista huolimatta hyvin pieniä. Testien mukaan suhteellisen luotettava raja,
jonka jälkeen ketjun keskiarvona saatu ellipsi näytti useimmissa tapauksissa pysyvän varsin muuttumattomana, oli noin 2500 MH-iteraatiota. Laskennallisesti algoritmi on siis
hyvin raskas.
Kuvasarjassa 4.1 ovat kuvasarjan 4.2 vasemmapuoleisen tilanteen näytehistoriat, joista
ei ole poistettu burn-in-osuuksia. Tässä tapauksessa keskimääräinen hyväksymistodennäköisyys näyttää olevan varsin alhainen, mutta testien perusteella ei voida sanoa, että näin
olisi yleensä.
4.2
Virheellisen havainnon tunnistus
Simuloitiin joukko pisteitä ellipsin muotoisesta tasajakaumasta sekä muutama ulkolainen.
Kuvassa 4.2 vertaillaan MH-algoritmin laskemia ellipsejä ilman näitä ulkolaisia sekä ul-
4. Algoritmin testausta
21
x
x
F1,1
100
400
0
200
−100
0
−200
0
500
1000
1500
2000
2500
−200
0
500
xF2,1
F1,2
1000
1500
2000
2500
1500
2000
2500
1500
2000
2500
xF2,2
200
0
100
−200
0
−100
0
500
1000
1500
2000
2500
−400
0
500
1000
s
g
760
1
740
0.95
720
0.9
700
0
500
1000
1500
2000
2500
0.85
0
500
1000
Kuva 4.1: Esimerkkejä näytehistoriasta.
kolaiset mukaan ottaen. Ellipsin muotoisesta tasajakaumasta generoituun 100 havainnon
pisteistöön lisättiin kolme viiden ulkolaishavainnon ryhmää.
gˆ = 0,952
Ilman ulkolaisia
Ulkolaiset mukana
gˆ = 0,914
Ilman ulkolaisia
Ulkolaiset mukana
gˆ = 0,906
Ilman ulkolaisia
Ulkolaiset mukana
Kuva 4.2: Lisättiin virheellisiä havaintoja viiden havainnon ryhmissä.
Testien perusteella pienet määrät virheellisiä havaintoja eivät yleensä vaikuta ellipsin
odotusarvoon merkittävästi. Näytehistorioitten perusteella virheelliset havainnot vaikuttivat lähinnä oikeellisuustodennäköisyyttä kuvaavaan parametriin G. Tämän parametrin
ollessa pieni kuuluvuusalueen kovarianssin laskennassa painotetaan tavallista enemmän
virheympyrää, mikä suurentaa kovarianssia.
4. Algoritmin testausta
4.3
22
Havaintopisteen moninkertaistaminen
Käytettiin testiaineistona havaintopistejoukkoa, joka sisälsi aluksi kuvasarjassa 4.3 näkyvät 100 pistettä ja johon sovitettiin MH-algoritmilla kuvaajaan piirretty ellipsi. Moninkertaistettiin kuvasarjassa punaisella merkittyä pistettä ja sovitettiin MH-algoritmilla uusi
ellipsi.
Kuvasarjassa 4.3 moninkertaistettu piste on ellipsin keskellä, kun taas kuvasarjassa
4.4 moninkertaistettu piste on ellipsin reunalla. Testien perusteella keskellä oleva moninkertaistettu piste ei varsin huomattavallakaan moninkertaistuksella vaikuta olennaisesti
laskennan tulokseen. Reunalla olevan pisteen moninkertaistaminen vaikuttaa selvästi herkemmin, mutta tällöinkin vaikutus näyttäisi olevan merkittävä vasta alkuperäistä pistejoukkoa suuremmalla moninkertaistuksella.
Todellisuudessa havaintosarjat saattavat olla simuloitua tilannetta paljon epäsäännöllisempiä, ja testeissä havaittiinkin myös tilanteita, joissa moninkertaistettu piste häiritsi
ellipsin estimaattia havaintosarjan kokoon nähden pienemmilläkin moninkertaistuksilla.
50 pistettä lisätty
200 pistettä lisätty
Ilman monink.
Monink. mukana
Monink. piste
500 pistettä lisätty
Ilman monink.
Monink. mukana
Monink. piste
1000 pistettä lisätty
Ilman monink.
Monink. mukana
Monink. piste
Kuva 4.3: Punaisella merkittyä pistettä on moninkertaistettu.
Ilman monink.
Monink. mukana
Monink. piste
4. Algoritmin testausta
23
50 pistettä lisätty
200 pistettä lisätty
Ilman monink.
Monink. mukana
Monink. piste
500 pistettä lisätty
Ilman monink.
Monink. mukana
Monink. piste
1000 pistettä lisätty
Ilman monink.
Monink. mukana
Monink. piste
Kuva 4.4: Punaisella merkittyä pistettä on moninkertaistettu.
Ilman monink.
Monink. mukana
Monink. piste
4. Algoritmin testausta
4.4
24
Posteriorin aikainvariantti päivittäminen
Posteriorin aikainvarianttia päivitystä testattiin neljällä eri varianssimallilla. Vähiten vapausasteita, 12, on mallissa, jossa kaikki ellipsiparametrit oletetaan riippumattomiksi toisistaan ja saadaan siis kuusi erillistä normaalijakaumaa. 14 vapaustasteen mallissa kunkin polttopisteen karteesisille koordinaateille lasketaan korrelaatiot, ja 18 vapausasteen
mallissa kaikki polttopistekoordinaatit voivat korreloida keskenään. Täysin korreloivassa
mallissa on 21 vapausastetta.
Kuvassa 4.5 on radiokartan muodostaminen neljässä vaiheessa korreloimattomalla mallilla (kovarianssimatriisi diagonaalinen). Havaitaan, että ellipsi lähestyy ilman päivitystä
laskettua ellipsiä. Testit osoittivat, että korreloimaton päivitys on likimain yhtä tehokas
kuin monimutkaisemmatkin mallit.
Koko pistejoukko
1 ensimm ist joukkoa
Koko pistejoukko
2 ensimm ist joukkoa
Koko pistejoukko
3 ensimm ist joukkoa
Koko pistejoukko
4 ensimm ist joukkoa
Kuva 4.5: Radiokartan päivitys korreloimattomalla mallilla (vasemmalta oikealle, ylhäältä alas).
4×50 simuloitua havaintopistettä.
Testien perusteella oletus jakauman pysymisestä vakiona on hyvin kriittinen. Jos siis
joissain päivitysvaiheissa pisteitä tulee vain tietystä osasta kuuluvuusaluetta, algoritmi on
epäluotettava. Ellipsin isoakselin lyheneminen kasvattaa uskottavuusfunktion (3.7) arvoa
nopeammin, kuin se pienentää normaalijakauman tiheysfunktion arvoa.
25
5.
JOHTOPÄÄTÖKSET
Tässä kandidaatintyössä esitettiin menetelmä ellipsinmuotoisen kuuluvuusalueen sijaintija muotoparametrien estimoimiseen sormenjälkimittausten perusteella. Menetelmä perustui Metropolisin ja Hastingsin otanta-algoritmiin, joka on Markovin ketjuun perustuva satunnaisalgoritmi (MCMC-algoritmi). Sovelluksessa mittauksia kerättäisiin kuuluvuusalueella sijaitsevien paikannustoiminnolla varustettujen liikkuvien päätteitten avulla. Mittaus
sisältää liikkuvan päätteen paikan sekä sen havaitsemat kuuluvuusalueet.
Kuuluvuusaluetta mallinnettiin ellipsinmuotoisen tasajakauman ja virheellisien havaintojen lähteeksi oletetun tasajakauman mikstuurina. Mittaukset oletettiin siis tästä jakaumasta generoduiksi näytteiksi.
Tässä työssä algoritmin testaamiseen käytettiin ainoastaan simuloituja mittauksia. Yleisesti ottaen algoritmin puutteena oli huomattava laskennallinen raskaus. Ehdokasjakauman varianssit oli asetettava suhteellisen pieniksi, mikä aiheutti monissa tapauksissa varsin voimakkaan burn-in-ilmiön. Joissakin tapauksissa tästä seurasi myös ketjun suppenemisen sellaisia arvoja kohti, jotka eivät silmämääräisen tarkastelun perusteella olleet
optimaalisia ratkaisuja. Algoritmin toiminta vaikutti olevan varsin herkkä sekä ehdokasjakauman varianssien että Markovin ketjulle laskettujen alkuarvojen suhteen. Ainakin alkuarvojen laskennassa onkin syytä hyödyntää prioritietoa, jos sitä on saatavilla.
Erityisesti testattiin algoritmin toimintaa joissakin kriittisiksi arvioiduissa tapauksissa.
Testien perusteella pienet määrät ulkolaishavaintoja eivät useimmissa tapauksissa häiritse
estimaattia olennaisesti vaan näissä tilanteissa virheellisten havaintojen jakauma painottuu. Algoritmi toimii ainakin simulaatioissa varsin hyvin myös tapauksissa, joissa pienelle osalle kuuluvuusaluetta on keskittynyt suuri osa havainnoista, eikä sen laskema ellipsi
siis keskity kovin nopeasti tuon osan ympärille. Algoritmin etuna onkin joustavuus vaikeissa tilanteissa. Suurten lukujen lain nojalla algoritmin avulla on ainakin periaatteessa
mahdollista estimoida optimaaliset parametrit mielivaltaisen tarkasti lisäämällä laskennan määrää. Lisäksi ketjun alkuarvojen ja ehdokasjakaumien varianssien laskentakaavoja
voidaan periaatteessa optimoida tilanteen mukaan.
Testattiin myös kuuluvuusalueen päivittämistä uusilla havainnoilla siten, että välillä
tallennettiin vain kuuluvuusalueen parametrit. Algoritmi toimii tässä tapauksessa hyvin
vain, jos uudet havainnot noudattavat likimain samaa jakaumaa kuin priorin muodostamiseen käytetyt havainnotkin. Jos havaintojen jakaumat poikkeavat, algoritmi on epäluotettava.
26
KIRJALLISUUTTA
[1] Richard A Johnson and Dean W Wichern. Applied Multivariate Statistical Analysis.
Prentice-Hall, 3rd edition, 1982.
[2] John F. Monahan. Numerical Methods of Statistics. Cambridge University Press,
2001.
[3] Esa Nummelin. MC’s for MCMC’ists. International Statistical Review / Revue Internationale de Statistique, 70(2):215–240, 2002.
[4] Gareth O. Roberts and Jeffrey S. Rosenthal. General state space Markov chains and
MCMC algorithms. Probability Surveys, 1:20–71, 2004.
[5] Luke Tierney. Markov Chains for Exploring Posterior Distributions. The Annals of
Statistics, 22(4):1701–1728, 1994.