1. No category

Toni Fadjukoff
CAUCHYN JAKAUMAN ESTIMOINTI
NUMEERISIN MENETELMIN
Kandidaatintyö
Tarkastaja: Simo Ali-Löytty
II
TIIVISTELMÄ
TAMPEREEN TEKNILLINEN YLIOPISTO
Teknisluonnontieteellinen koulutusohjelma
TEKIJÄN NIMI: Toni Fadjukoff
Kandidaatintyö, 18 sivua
Joulukuu 2011
Pääaine: Matematiikka
Tarkastaja: Simo Ali-Löytty
Avainsanat: Cauchyn jakauma, Newtonin menetelmä
Cauchyn jakauma on hyvin paksuhäntäinen todennäköisyysjakauma. Jakaumaa voidaan käyttää joustavissa menetelmissä, sillä se sallii satunnaisotoksissa isotkin poikkeamat ilman että sen parametrit muuttuvat. Cauchyn jakauma on tunnettu myös
teoreettisista ominaisuuksistaan, sillä jakaumalle ei pystytä laskemaan odotusarvoa.
Tässä työssä esitetään useita Cauchyn jakauman ominaisuuksia. Jakauman parametrien merkitys selvitetään. Cauchyn kertymäfunktio ja affiinimuunnos lasketaan auki.
Näytetään myös, miksi jakaumalla ei ole odotusarvoa. Lisäksi johdetaan Cauchyn
jakaumalle suurimman uskottavuuden estimaatti. Suurimman uskottavuuden estimaatit lasketaan käyttäen muunneltua Newtonin algoritmia. Newtonin algoritmin
tarvitsevat derivaatat lasketaan sekä analyyttisesti että numeerisesti.
Estimaattorin tarkkuutta tutkitaan ajamalla lukuisat määrät simulaatioita. Numeerisen ratkaisijan nopeutta ja tarkkuutta verrataan MATLAB-ohjelmistoon sisältyviin valmisratkaisijoihin. Testauksessa havaitaan itse kirjoitettu algoritmi nopeammaksi kuin yleiskäyttöiset optimointirutiinit. Lisäksi selvitetään mahdollisia ongelmia optimointialgoritmin toteuttamisessa.
III
SISÄLLYS
1. Johdanto . . . . . . . . . . . . . . . . . . . . . . . .
2. Teoria . . . . . . . . . . . . . . . . . . . . . . . . . .
2.1 Cauchyn jakauma . . . . . . . . . . . . . . . .
2.2 Suurimman uskottavuuden estimointi . . . . .
2.2.1 Analogia Bayesiläisiin menetelmiin . . . . .
2.2.2 Log-uskottavuus . . . . . . . . . . . . . . .
2.3 Estimaatin arvon laskeminen . . . . . . . . . .
2.3.1 Newtonin menetelmä . . . . . . . . . . . .
2.3.2 Newtonin menetelmä numeerisesti . . . . .
2.3.3 Numeerinen osittaisderivointi . . . . . . . .
2.3.4 Newtonin alkuarvaus . . . . . . . . . . . .
3. Testaus . . . . . . . . . . . . . . . . . . . . . . . . .
3.1 Satunnaislukujen generointi . . . . . . . . . . .
3.2 Vertailtavat algoritmit . . . . . . . . . . . . . .
3.3 Newtonin menetelmän toteuttaminen . . . . . .
3.3.1 Ongelmia algoritmin kanssa . . . . . . . .
3.4 Suurimman uskottavuuden estimaatin tarkkuus
3.5 Numeeristen menetelmien nopeus . . . . . . . .
4. Yhteenveto . . . . . . . . . . . . . . . . . . . . . . .
Lähteet . . . . . . . . . . . . . . . . . . . . . . . . . .
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
1
2
2
4
5
5
6
7
8
9
10
11
11
12
13
14
16
17
18
19
IV
TERMIT JA SYMBOLIT
f (x; θ, σ)
F (x; θ, σ)
Cauchy(θ, σ)
Normal(θ, v 2 )
L(x1 , . . . , xn ; θ, σ)
l(x1 , . . . , xn ; θ, σ)
Lneg (x1 , . . . , xn ; θ, σ)
lneg (x1 , . . . , xn ; θ, σ)
Mk (x)
�k
i=1 xi
qr
limx→x0
R
Z+
hT
∇f (y)
H(y)
p(θ, σ|x)
p(x|θ, σ)
O(h2 )
Cauchyn tiheysfunktio, jonka sijaintiparametri on θ ja skaalausparametri σ
Cauchyn kertymäfunktio, jonka sijaintiparametri on θ ja skaalausparametri σ
Cauchyn jakauma, jonka sijaintiparametri on θ ja skaalausparametri σ
normaalijakauma, jonka odotusarvo on θ ja varianssi v 2
uskottavuusfunktio otoksesta x1 , . . . , xn
uskottavuusfunktion logaritmi
−L(x1 , . . . , xn ; θ, σ)
−l(x1 , . . . , xn ; θ, σ)
todennäköisyysjakauman k:s momentti
tulo alkioista {xi |1 ≤ i ≤ k}
otoskvantiili, pienin luku jolla r:s osa otoksesta on pienempi
raja-arvo, kun x lähestyy x0 :a
reaalilukujen joukko
positiivisten kokonaislukujen joukko
vektorin h transpoosi
funktion f gradientti pisteessä y
Hessen matriisi kohdassa y
θ:n ja σ:n tiheysfunktio kun otos on x
x:n tiheysfunktio, kun parametrit ovat θ ja σ
h:n polynomi, joissa kaikkien termien aste vähintään 2
1
1.
JOHDANTO
Mallintamisessa käytetään usein apuna todennäköisyysjakaumia. Monissa mittauksissa virhettä arvioidaan jotakin jakaumaa noudattavana satunnaismuuttujana. Suosittu ja paljon tutkittu normaalijakauma sopii hyvin monenlaiseen dataan, mutta
poikkeuksellisen paljon poikkeavia mittauksia, ulkolaisia, ei normaalijakaumaa noudattavilla mittauksilla pitäisi olla lainkaan. Tosimaailmassa mittauksissa kuitenkin
usein havaitaan ulkolaisia, eikä edes kovin harvoin. Tällaisiin malleihin siis sopii
todennäköisyysjakauma, joka sallii poikkeamat mittauksissa.
Menetelmissä, joita kutsutaan myös robusteiksi menetelmiksi, käytetään paksuhäntäisiä jakaumia mallintamaan data, jossa esiintyy ulkolaisia. Cauchyn jakauma
on eräs todennäköisyysjakauma, jolla on erityisen paksu häntä. Se on sopiva tällaisen
datan mallintamiseen. Lisäksi Cauchyn jakauma tunnetaan myös erityispiirteestään:
jakaumalla ei ole odotusarvoa eikä muitakaan momentteja.
Tässä kandidaatintyössä tutkitaan Cauchyn jakaumaa ja sen parametrien estimointia numeerisin menetelmin. Teoria-luvussa esitellään Cauchyn jakauman ominaisuuksia. Cauchyn parametreille esitetään suurimman uskottavuuden estimaatti,
jolle esitetään myös Bayesiläinen analogia. Suurimman uskottavuuden estimaatin
arvon laskeminen on optimointitehtävä, joka voidaan laskea useilla numeerisilla menetelmillä, joista Newtonin menetelmä käydään läpi hieman tarkemmin. Newtonin
menetelmän laskentaan esitetään kaksi tapaa. Menetelmässä käytetyille derivaatoille
esitetään analyyttinen ratkaisu sekä numeerinen erotusosamäärään perustuva ratkaisu.
Testaus-luvussa tutkitaan käytännössä eri menetelmiä laskea suurimman uskottavuuden estimaatti. Teoria-luvussa esitettyjen kahden menetelmän lisäksi vertailua tehdään myös tietokoneohjelmisto MATLABiin sisältyviin optimointirutiineihin.
Menetelmien laskenta-aikoja ja tarkkuutta verrataan keskenään. Lisäksi pohditaan
hyviä ja huonoja puolia eri menetelmissä. Viimeisessä luvussa tehdään yhteenveto
teoriasta ja testauksesta.
2
2.
TEORIA
2.1
Cauchyn jakauma
Yksiulotteinen Cauchyn jakauma on todennäköisyysjakauma, joka noudattaa koko
reaalilukuavaruudessa tiheysfunktiota
�
�
�2 �−1
1
x−θ
f (x; θ, σ) =
1+
, σ ∈ (0, ∞), θ ∈ R.
πσ
σ
(2.1)
Huomataan, että (2.1) todella on tiheysfunktio, sillä se saa aina positiivisia arvoja,
�
ja lisäksi f (x) d x = limx→∞ F (x) = 1 (2.5).
Parametriä θ kutsutaan sijaintiparametriksi. Tiheysfunktion huippu on kohdassa
x = θ. Tiheysfunktion maksimiarvo on
max f (x) = f (θ) =
1
πσ
(2.2)
Jakauman parametriä σ kutsutaan yleisesti skaalausparametriksi. Skaalaa voidaan
havainnollistaa, kun lasketaan f (θ + σ):
�
� σ �2 �−1
1
1
f (θ + σ) =
1+
=
πσ
σ
2πσ
(2.3)
Mitä suurempi σ, sitä leveämpi huippu tiheysfunktiolla on.
Cauchy-jakauman affinimuunnos on seuraava. Oletetaan että x ∼ Cauchy(θ, σ).
Nyt y = ax + b, a > 0 noudattaa tiheysfunktiota [9]
1
f (y) = f
a
�
x−b
a
�
�
�
�2 �−1
1
x − b − θa
=
1+
,
πσa
σa
(2.4)
jolloin y ∼ Cauchy(θa + b, σa).
Cauchyn jakauma on erikoistapaus Studentin t-jakaumasta. Se vastaa Studentin
t-jakaumaa vapausasteilla 1. Cauchyn jakauman muoto on suurinpiirtein samanlainen kuin normaalijakauman, mutta sen hännillä on enemmän massaa. Jakaumien tiheysfunktiot ovat verrattavissa kuvassa 2.1. Cauchy-jakauman paksuhäntäisyys verrattuna normaalijakaumaan näkyy selvästi logaritmisella asteikolla kuvassa 2.2.
2. Teoria
3
Cauchy(θ, σ)
f (θ; θ, σ)
1
f (θ; θ, σ)
2
Normal(θ, v 2 )
θ−σ θ θ+σ
θ + 2v 2 ln 1/2 θ − 2v 2 ln 1/2
Kuva 2.1: Cauchyn jakauma ja normaalijakauma
θ + 2v 2 ln 1/2
θ+σ
x
ln(f (θ; θ, σ))
ln( 12 f (θ; θ, σ))
Cauchy(θ, σ)
ln(f (x))
Normal(θ, v 2 )
Kuva 2.2: Paksuhäntäisyys semilog-asteikolla
Cauchyn jakaumalle kertymäfunktio saadaan integroimalla tiheysfunktiota:
�
�
�2 �−1
1
x−θ
F (x; θ, σ) = lim
1+
dx
y→−∞ y πσ
σ
�
�
�
�
��
y−θ
1
x−θ
=
σ arctan
− lim arctan
y→−∞
πσ
σ
σ
�
�
1
x−θ
1
=
arctan
+
π
σ
2
�
x
(2.5)
Kertymäfunktion kuvaaja on piirretty kuvaan 2.3. Myös kertymäfunktiosta on havaittavissa jakauman paksuhäntäisyys.
Todennäköisyysjakauman k:s momentti on määritelty
�
Mk (x) = xk f (x) d x, k ∈ Z+ .
(2.6)
Cauchyn jakauman momentit eivät ole äärellisiä. Tämä voidaan osoittaa tarkastele-
2. Teoria
4
1
π
1
arctan(x) +
1
2
Kuva 2.3: Kertymäfunktio
malla momenttien itseistä suppenemista:
�
σ
1
|Mk (x)| =
|x|k
dx
2
π σ + (x − θ)2
�
σ ∞
xk
dx
>
π max(1,θ) σ 2 + (x − θ)2
�
σ ∞
x
≥
dx
2
π max(1,θ) σ + (x − θ)2
�
σ ∞
x� + θ
=
d x�
π max(1−θ,0) σ 2 + x�2
�� ∞
�
� ∞
σ
x�
θ
�
�
=
dx +
dx
2
�2
2
�2
π
max(1−θ,0) σ + x
max(1−θ,0) σ + x
�
�
� x ��
σ
1
θ
2
2
=
lim
log(x + σ ) + arctan
π x→∞ 2
σ
σ
�
��
1
θ
max(1 − θ, 0)
2
2
− log(max(1 − θ, 0) + σ ) − arctan
2
σ
σ
Koska log(x2 ) → ∞, kun x → ∞, niin integraali ei suppene. Tällöin momentti ei
myöskään suppene ja sitä ei ole määritelty. Cauchyn jakaumalle ei siis ole määritelty
yhtään äärellistä momenttia. Ensimmäinen momentti on odotusarvo, ja Cauchyn
jakaumalle ei ole siis myöskään odotusarvoa olemassa.
2.2
Suurimman uskottavuuden estimointi
Cauchyn jakauman estimoinnissa puhutaan ongelmasta, jossa yritetään määrätä parametrit θ ja σ realisoituneiden satunnaismuuttujien xi , i = 1, . . . , n
xi ∼ Cauchy(θ, σ),
(2.7)
perusteella. Jos muuttujat xi ovat riippumattomia, niin uskottavuusfunktio on
L(x1 , . . . , xn ; θ, σ) =
n
�
i=1
f (xi ; θ, σ).
(2.8)
2. Teoria
5
Uskottavuusfunktion arvo lasketaan otoksen x1 , . . . , xn ja parametrien θ ja σ suhteen. Uskottavuusfunktion arvo on itse asiassa todennäköisyys, jolla Cauchyn jakaumaa parametrein θ, σ noudattavat satunnaismuuttujat saavat arvot x1 , . . . , xn . Suurimman uskottavuuden menetelmässä etsitään uskottavuusfunktion maksimiarvoa
jakauman parametrien suhteen. Tätä arvoa kutsutaan suurimman uskottavuuden
estimaatiksi.
2.2.1
Analogia Bayesiläisiin menetelmiin
Suurimman uskottavuuden estimoinnille voidaan löytää vastine Bayesiläisestä maailmasta. Bayesin kaava kirjoitetaan seuraavassa muodossa
p(θ, σ|x1:n ) =
p(x1:n |θ, σ)p(θ, σ)
.
p(x1:n )
(2.9)
Nimittäjä p(x1:n ) on itseasiassa aina vakio estimoitavien parametrien suhteen. Tällöin sen voidaan ajatella olevan vain normalisointivakio tehtävän kannalta. Tekijää
p(θ, σ) kutsutaan prioriksi. Jos otetaan vakiopriori, siten että p(θ, σ) ∝ 1 todennäköisyyden määritelyalueella σ > 0, θ ∈ R, ja muuttujat xi riippumattomiksi, niin
saadaan johdettua suurimman uskottavuuden kaava.
p(θ, σ|x1:n ) ∝ p(x1:n |θ, σ)
= p(x1 |θ, σ)p(x2 |θ, σ) . . . p(xn |θ, σ)
n
�
=
p(xi |θ, σ)
(2.10)
i=1
Näin posteriorijakaumaksi saadaan klassista uskottavuusfunktiota vastaava jakauma.
2.2.2
Log-uskottavuus
Uskottavuusfunktion arvon laskemiseksi käytetään usein log-muunnosta. Koska logaritmifunktio on aidosti kasvava, sen maksimiarvo on samassa pisteessä kuin alkuperäisen funktion maksimiarvo. Log-uskottavuusfunktiossa tiheysfunktioiden tulon
sijasta tutkitaan tulon logaritmia, joka muuttuu logaritmien summaksi:
l(x1 , . . . , xn ; θ, σ) = ln(
n
�
i=1
f (xi ; θ, σ)) =
n
�
ln(f (xi ; θ, σ))
(2.11)
i=1
Log-muunnosta käytetään, koska log-uskottavuuden laskeminen on monin tavoin
mielekkäämpää kuin tavallisen uskottavuusfunktion. Suurimman uskottavuuden estimaattori on log-uskottavuusfunktion maksimikohdassa.
2. Teoria
6
Cauchyn jakaumalle log-uskottavuusfunktio (2.11) on
l(x1 , . . . , xn ; θ, σ) =
n
�
i=1
− ln(π) + ln(σ) − ln(σ 2 + (xi − θ)2 ).
(2.12)
Osittaisderivoimalla log-uskottavuusfunktio (2.12) molempien parametrien suhteen saadaan
n
�
∂
xi − θ
l(x1 , . . . , xn ; θ, σ) = 2
(2.13)
2 + (x − θ)2
∂θ
σ
i
i=1
ja
n
∂
n �
2σ
l(x1 , . . . , xn ; θ, σ) = −
.
2
2
∂σ
σ
σ
+
(x
−
θ)
i
i=1
(2.14)
Log-uskottavuuden osittaisderivaatan nollakohdalle ei ole analyyttistä ratkaisua
olemassa enää, kun otoksen kooksi kasvaa viisi, sillä silloin yhtälön aste kasvaa viideksi tai suuremmaksi. Analyyttinen kaava on olemassa vain neljännen tai pienemmän asteen yhtälöille. [4]
Nollakohdat voi ratkaista numeerisesti myös suuremmille otoksille. Molempien
parametrien yhtäaikainen suurimman uskottavuuden estimaatti on todistetusti olemassa ja yksikäsitteinen määrittelyalueellaan [2].
2.3
Estimaatin arvon laskeminen
Numeeriseen ratkaisuun on useampia tapoja. Yksi yksinkertaisimmista tavoista on
Newtonin iterointi. Newtonin iteroinnilla voidaan etsiä funktion nollakohtia. Menetelmällä voidaan etsiä funktion ääriarvoa, kun etsitäänkin funktion derivaatan
nollakohtia.
Ennen menetelmän esittelyä on syytä esitellä muutama käsite. Funktion saa lokaalin maksimiarvon pisteessä y, jos
∇f (y) = 0
(2.15)
hT H(y)h < 0 kaikille h �= 0
Vektorin h osoittama suunta on kasvusuunta, jos
0 < hT ∇f (y).
Funktion lokaali maksimipiste löytyy iteroimalla kasvusuuntaan päin. [7]
(2.16)
2. Teoria
2.3.1
7
Newtonin menetelmä
Merkitään y = (θ, σ)T , missä θ on jakauman sijaintiparametri ja σ skaalausparametri. Halutaan löytää piste y + h, missä gradientin arvo on 0. Tällaisessa pisteessä
funktio voi saavuttaa ääriarvonsa. Newtonin menetelmän johtaminen käy Taylorin
sarjan avulla.
∇f (y + h) = ∇f (y) + H(y)h + O(h2 )
≈ ∇f (y) + H(y)h
(2.17)
Ratkaistaan yhtälö (2.17) h:n suhteen. Oletuksella H(y) ei ole singulaarinen saadaan
h = −[H(y)]−1 ∇f (y).
(2.18)
Newtonin algoritmissa lasketaan iteratiivisesti uutta estimaattia kaavalla
yk+1 = yk − [H(y)]−1 ∇f (y).
(2.19)
Jos algoritmille annettava alkuarvo y0 on tarpeeksi lähellä funktion f (y) nollakohtaa,
Newtonin iterointi suppenee neliöllisesti [3].
Vain jos Hessen matriisi on negatiivisesti definiitti pisteessä y, Newtonin menetelmä suppenee kohti lokaalia maksimia. Negatiivisesti definiitille matriisille pätee,
että hT H(y)h < 0, h �= 0, jolloin
−hT ∇f (y) = hT H(y)h < 0.
(2.20)
Silloin myös hT ∇f (y) > 0 ja h:n suuntaan funktion arvo kasvaa. Tällaisessa tapauksessa Newtonin algoritmia voidaan käyttää ratkaisun etsimiseen.
Käyttämällä linjahaku-menetelmää Newtonin menetelmän lisäksi voidaan hakua
tehostaa lisäämällä apukerroin β ∈ (0, 1] yhtälöön (2.19):
yk+1 = yk − β[H(y)]−1 ∇f (y).
(2.21)
Joillakin arvoilla iteraatioaskel (2.18) saattaa muodostua liian pitkäksi ja iterointi
karkaa kauas. Tällöin viivahakumenetelmän periaatteiden mukaisesti voidaan iteraatioaskelta yrittää myös muilla suunnan h pisteillä.
Newtonin menetelmällä optimoitava funktio on meidän tapauksessamme loguskottavuusfunktio (2.12). Funktion gradientti saadaan yhdistämällä yhtälöt (2.13)
2. Teoria
8
ja (2.14).
∇f (y) =
�
2
n
σ
�
�n
−
xi −θ
i=1 σ 2 +(xi −θ)2
�n
2σ
i=1 σ 2 +(xi −θ)2
(2.22)
Kaksiulotteisessa tapauksessa Hessen matriisi voidaan kirjoittaa
H(y) =
�
∂ 2 f (x;y)
dθ 2
∂ 2 f (x;y)
∂σ∂θ
∂ 2 f (x;y)
dθ∂σ
∂ 2 f (x;y)
∂σ 2
�
.
(2.23)
Hessen matriisin alkioina olevat osittaisderivaatat ovat Cauchyn jakauman loguskottavuusfunktiolle
n
�
∂ 2 f (x; θ, σ)
2(−σ 2 + (yi − θ)2 )
=
∂θ2
(σ 2 + (xi − θ)2 )2
i=1
n
�
∂ 2 f (x; θ, σ)
−4σ(xi − θ)
=
∂θ∂σ
(σ 2 + (xi − θ)2 )2
i=1
n
�
∂ 2 f (x; θ, σ)
−4σ(xi − θ)
=
∂σ∂θ
(σ 2 + (xi − θ)2 )2
i=1
(2.24)
(2.25)
(2.26)
n
∂ 2 f (x; θ, σ)
−n � 4σ 2 − 2(σ 2 + (xi − θ)2 )
=
+
∂σ 2
σ2
(σ 2 + (xi − θ)2 )2
i=1
2.3.2
(2.27)
Newtonin menetelmä numeerisesti
Newtonin algoritmin analyyttinen versio on käypä menetelmä, mutta sen käyttö vaatii derivaattojen laskemista. Analyyttisten derivaattojen ratkaisu saattaa olla työlästä ja se on tehtävä jokaiselle funktiolle erikseen. Newtonin algoritmia voi kuitenkin
käyttää myös numeerisesti laskettujen derivaattojen avulla. Derivaatan määritelmä
on
f (x + h) − f (x)
.
h→0
h
f � (x) = lim
(2.28)
Derivaattaa voidaan arvioida numeerisesti laskemalla suoraan erotusosamäärän määritelmästä pienellä h:n arvolla, joskin parempi approksimaatio saadaan laskemalla
molemman puoleiset derivaatat yhteen
f � (x) ≈
f (x + h) − f (x − h)
,
2h
(2.29)
missä h on pieni luku [3].
Numeerista derivaatan laskemista voidaan soveltaa myös useaan otteeseen. Toisen
2. Teoria
9
asteen derivaatta saadaan kätevästi esimerkiksi:
f � (x) − f � (x − h)
f (x) ≈
h
�
�
1 f (x + h) − f (x) f (x) − f (x − h)
=
−
h
h
h
f (x + h) − 2f (x) + f (x − h)
=
h2
��
2.3.3
Numeerinen osittaisderivointi
Newtonin algoritmin ratkaisuun tarvitaan toisen asteen osittaisderivaattoja. Voimme ajatella funktiomme kolmen muuttujan funktiona.
f (x, θ, σ) = f (x|θ, σ)
Osittaisderivoinnissa kohdellaan tiettyä funktion parametria derivoitavana muuttujana. Esimerkiksi θ:n suhteen osittaisderivaatta voidaan laskea:
∂
f (x, θ + h, σ) − f (x, θ, σ)
f (x, θ, σ) ≈
∂θ
h
Näin ollen osittaisderivaatta molempien muuttujien suhteen on
�
�
∂2
∂ f (x, θ + h, σ) − f (x, θ, σ)
f (x, θ, σ) ≈
∂θ∂σ
∂σ
h
f (x, θ + h, σ + h) − f (x, θ + h, σ) − f (x, θ, σ + h) + f (x, θ, σ)
=
h2
Toisia osittaisderivaattoja voidaan siis arvioida numeerisesti seuraavilla lausekkeilla:
∂ 2 f (x, θ, σ)
∂θ2
2
∂ f (x, θ, σ)
∂θ∂σ
∂ 2 f (x, θ, σ)
∂σ∂θ
∂ 2 f (x, θ, σ)
∂σ
f (x, θ + h, σ) − 2f (x, θ, σ) + f (x, θ − h, σ)
(2.30)
h2
f (x, θ + h, σ + h) − f (x, θ + h, σ) − f (x, θ, σ + h) + f (x, θ, σ)
=
(2.31)
h2
f (x, θ + h, σ + h) − f (x, θ + h, σ) − f (x, θ, σ + h) + f (x, θ, σ)
=
(2.32)
h2
f (x, θ, σ + h) − 2f (x, θ, σ) + f (x, θ, σ − h)
=
(2.33)
h2
=
Kun gradientti (2.22) lasketaan osittaisderivaattojen (2.29) avulla ja yhtälöt (2.30)–
(2.33) sijoitetaan Hessen matriisin (2.23) alkioiksi, voidaan Newtonin algoritmia ratkaista ilman analyyttisten derivaattojen tuntemista.
2. Teoria
2.3.4
10
Newtonin alkuarvaus
Newtonin menetelmän nopea suppeneminen riippuu menetelmälle annetuista alkuarvoista. Hyvä alkuarvo θ-parametrille on mediaani. Noin puolet näytteistä saa
pienemmän arvon kuin θ, ja noin puolet saa suuremman arvon. Parametrille σ löydetään sopiva alkuarvo tutkimalla kertymäfunktiota:
1
F (θ + σ) = arctan
π
�
θ+σ−σ
σ
�
+
1
1
1
3
= arctan(1) + =
2
π
2
4
(2.34)
Vastaavasti F (θ − σ) = 14 . Tällöin saadaan hyvä alkuestimaatti lausekkeesta
σ
ˆ=
q3/4 − q1/4
.
2
(2.35)
11
3.
TESTAUS
Edellisessä kappaleessa esitellyt funktiot on ohjelmoitu MATLABilla. Testaamme
estimaattien laskenta-aikoja sekä suurimman uskottavuuden estimaatin tarkkuutta
erikokoisilla otoksilla Cauchyn jakaumasta. Lisäksi mukana on muutama mustan
laatikon algoritmi.
Tässä luvussa käyttämäni ja toteuttamani ohjelmat ovat saataville osoitteesta:
http://www.students.tut.fi/~lambergt/kandi/
3.1
Satunnaislukujen generointi
Ensimmäiseksi on tarpeellista luoda Cauchyn jakaumasta satunnaislukuja testausdataa varten. Hyvä tapa luoda satunnaislukuja jatkuvasta jakaumasta on käänteiskertymämenetelmä. Menetelmässä luodaan tasan jakautunut satunnaisluku u välille [0, 1). Sen jälkeen saadaan laskettua jakauman mukainen satunnaisluku kaavalla
x = F −1 (u). [10]
Kertymäfunktion (2.5) käänteisfunktio on
F
−1
� �
��
1
(u) = σ tan π u −
+ θ, u ∈ [0, 1).
2
(3.1)
Cauchy-jakautunutta dataa on esitelty kuvassa 3.1.
40
35
30
25
20
15
10
5
0
−50
−40
−30
−20
−10
0
Kuva 3.1: Cauchy(0,1)-jakautuneen datan histogrammi, 60 alkiota
10
3. Testaus
3.2
12
Vertailtavat algoritmit
Sekä analyyttistä että numeerista derivaattaa käyttävät ratkaisut on ohjelmoitu
käyttäen MATLAB-ohjelmointiympäristöä. Lisäksi vertailuun on otettu mukaan kaksi valmisrutiinia, joilla voi löytää suurimman uskottavuuden estimaatin.
MATLAB-ympäristön optimointi-työkalulaatikkoon sisältyvä fminunc etsii annetulle funktiolle minimin. Koska haluamme etsiä funktion maksimia, käytämme
funktioita
Lneg (x; θ, σ) = −L(x; θ, σ)
lneg (x; θ, σ) = −l(x; θ, σ).
(3.2)
(3.3)
Funktio fminunc käyttää sisäisesti eri algoritmeja. Rutiinia testattiin kolmella eri tavalla. Ilman analyyttistä gradienttia ja Hessen matriisia rutiini etsii ääriarvoa BFGS
(Broyden–Fletcher–Goldfarb–Shanno) kvasi-Newton -menetelmällä. Menetelmässä
Hessen matriisia päivitetään joka aika-askeleella sen täydellisen uudelleenlaskemisen sijaan kuten Teoria-luvussa esitellyssä menetelmässä [5]. Kun funktiolle annetaan parametrina analyyttisen gradientin laskeva funktio, käytettävä algoritmi on
sisäisesti-heijastuva Newton-menetelmä, jota esitellään Colemanin paperissa [1]. Tämä menetelmä oli käytettävissä analyyttisen Hessen matriisin kanssa tai ilman. [8]
Toinen valmisalgoritmi on myös optimointi-työkalulaatikkoon sisältyvä rutiini
fminsearch. Se etsii myös funktion minimin, jolloin optimoitava funktio on myös
(3.2) tai (3.3). Funktio fminsearch käyttää algoritminaan Nelder-Mead simplexmenetelmää [6]. Menetelmässä ei käytetä funktion derivaattoja ollenkaan.
3. Testaus
3.3
13
Newtonin menetelmän toteuttaminen
Ohjelmoituna Newtonin algoritmi näyttää algoritmilta 1. Ainoa ero toteuttamani
analyyttisen Newtonin iteraation ja täysin numeerisen iteroinnin välillä on rivin 3
toteutus. Ensimmäisessä menetelmässä käytettiin analyyttista gradienttia ja Hessen
matriisia ja toisella menetelmällä approksimaatiot laskettiin erotusosamäärien avulla.
Algoritmi 1 Newtonin algoritmi
1: (θ, σ)T ← (q1/2 , (q3/4 − q1/4 )/2)T
2: while β > εβ do
3:
(θn , σn )T ← (θ, σ)T − β[H(θ, σ)]−1 ∇f (θ, σ)
4:
if f (θn , σn ) > f (θ, σ) then
5:
(θ, σ)T ← (θn , σn )T
6:
else
7:
β ← β/2
8:
end if
9: end while
Kuvassa 3.2 nähdään Newtonin menetelmän suppeneminen tasossa. Data sisältää 20 pistettä ja se generoitiin jakaumasta Cauchy(30,20). Suurimman uskottavuuden estimaatti asettui pisteeseen (24.51, 17.85). Useimmissa tapauksissa muutama
iteraatio riittää Newtonin menetelmällä suppenemiseen.
Newton iteration, N=20
30
25
σ
20
15
10
5
0
5
10
15
20
25
θ
30
35
40
45
50
Kuva 3.2: Newtonin menetelmän suppeneminen θ, σ-tasossa
3. Testaus
3.3.1
14
Ongelmia algoritmin kanssa
Simulaatioissa havaittiin, että algoritmi 1 ei aina supennut samaan pisteeseen kuin
MATLABin referenssimenetelmät. Tarkempi tarkastelu osoitti, että välillä yhtälön
(2.18) antama suunta ei ollutkaan kasvusuunta.
Tilannetta tutkittiin piirtämällä Newtonin suppenemissuunta monessa pisteessä
tasolle. Kuvassa 3.3 nähdään, mihin suuntaan algoritmi lähtee suppenemaan kussakin pisteessä. Vihreän taustavärin alueella Hessen matriisi on negatiivisesti definiitti. Taulukossa 3.1 olevalla otoksella ei saada käyttämäämme alkuarvausta sopivalle
alueelle. Algoritmin alkupisteessä Hessen matriisi ei siis ole negatiivisesti definiitti.
0.6
Alkuarvaus
Ratkaisu
0.5
σ
0.4
0.3
0.2
0.1
0
4.2
4.25
4.3
4.35
4.4
4.45
θ
4.5
4.55
4.6
Kuva 3.3: Newtonin menetelmän iterointisuunta
[4.5427, -33.8699, 4.4489, 5.9635, 4.3210, 4.4862, 4.5230, 5.9778]
Taulukko 3.1: Ei-suoraan ratkeavan tapauksen otos
4.65
4.7
3. Testaus
15
Ratkaisuksi muodostui muodostaa hybridi-menetelmä, jossa tarkastellaan jokaisella iteraatiolla Hessen matriisin definiittisyyttä. Jos matriisi ei ole negatiivisesti
definiitti, käytetään gradienttimenetelmää. Jos matriisi on negatiivisesti definiitti,
iteroidaan Newtonin menetelmällä. Hybridi-menetelmä on esitetty esimerkiksi lähteessä [7].
Algoritmi 2 Hybridi-algoritmi: Newton-, gradientti- ja viivahakumenetelmät
1: (θ, σ) ← (q1/2 , (q3/4 − q1/4 )/2)
2: while |∇f (θ, σ)| < tol ∧ Iter < MaxIter do
3:
l ← f (θ, σ)
4:
if H on negatiivisesti definiitti then
5:
δ ← −[H(θ, σ)]−1 ∇f (θ, σ)
6:
B ← B1 (δ)
7:
else
8:
δ ← −∇f (θ, σ)
9:
B ← B2 (δ)
10:
end if
11:
(θ0 , σ0 ) ← (θ, σ)
12:
for β ∈ B do
13:
(θ� , σ � ) ← (θ, σ) + βδ
14:
if f (θ� , σ � ) > l then
15:
(θ, σ) ← (θ� , σ � )
16:
l ← f (θ� , σ � )
17:
end if
18:
end for
19:
if (θ, σ) = (θ0 , σ0 ) then
20:
virhe(’Ei löydetä parempaa siirtymää.’)
21:
end if
22:
if |(θ, σ) − (θ0 , σ0 )| < dtol then
23:
palauta(’Siirtymä on tarpeeksi pieni.’)
24:
end if
25:
Iter ← Iter + 1
26: end while
27: if Iter < MaxIter then
28:
palauta(’Maksimi löytynyt.’)
29: else
30:
virhe(’Algoritmi ei supennut.’)
31: end if
Algoritmi 2 vastaa toimivaa toteutusta. Tärkein parannus edelliseen versioon on
siis Hessen matriisin definiittisyyden tutkiminen. Riveillä 12-18 tarkastellaan parasta siirtymää nykyisestä pisteestä. Kertoimen β saamat arvot riippuvat sekä käytettävästä iteraatiomenetelmästä että siirtymän δ suuruusluokasta.
Joukot B1 (δ) ja B2 (δ) ovat tarkemmin määritelty ohjelmakoodissa. Gradienttimenetelmää käytettäessä pitää tutkimaan suurempia siirtymiä (B1 ) kuin Newtonin
3. Testaus
16
meneltemällä (B2 ). Tarvittavat siirtymät testattiin ajamalla ohjelmaa useasti.
Simulaatioissa uusi algoritmi pystyi ratkaisemaan yli sata tuhatta estimaattia
samalla tarkkuudella kuin MATLABin valmisoptimointialgoritmit ilman yhtäkään
virhettä.
3.4
Suurimman uskottavuuden estimaatin tarkkuus
Taulukossa 3.2 esitellään suurimman uskottavuuden estimaatin tarkkuuksia erikokoisille otoksille. Jokaista otoskokoa luotiin 10 eri otosta. Näistä laskettiin suurimman uskottavuuden estimaatti Newtonin menetelmällä. Jokaisesta 10 otoksen otoksesta lasketaan keskiarvo, keskivirhe ja 95%-kvantiili virheelle.
n Keskiarvo
10
24.45
20
18.44
30
22.24
40
18.95
50
21.32
100
20.89
500
20.34
1000
20.05
10000
19.91
θ
Keskivirhe
10.37
11.32
7.14
6.35
3.16
2.35
1.71
1.21
0.31
95%-virhe
28.72
27.63
14.98
12.19
5.96
3.65
3.84
2.87
0.62
Keskiarvo
27.23
33.24
29.94
26.63
29.23
30.47
29.54
30.34
30.10
σ
Keskivirhe
7.69
7.95
5.14
5.21
1.93
2.62
1.65
1.29
0.31
95%-virhe
18.54
19.97
9.86
12.12
4.94
6.42
3.48
2.29
0.57
Taulukko 3.2: Jakaumasta Cauchy(20,30) generoitujen otosten tunnuslukuja
Taulukosta 3.2 havaitaan MLE-estimaatin hyvä tarkkuus. Kun jakauman koko
kasvaa jo muutaman kymmenen otoksen kokoiseksi, antaa suurimman uskottavuuden menetelmä melko tarkkoja arvoja jakauman parametreille.
3. Testaus
3.5
17
Numeeristen menetelmien nopeus
Eri menetelmistä testattiin myös niiden nopeutta. Jokainen menetelmä ratkaisee periaatteessa saman suurimman uskottavuuden ongelman, kuitenkin niiden käyttämä
menetelmä on hieman erilainen.
Taulukossa (3.3) on listattu eri menetelmien laskenta-aikoja suurten otosten kanssa. Nopein ratkaisija on itse kirjoitettu algoritmi. Tämä selittyy silläkin, että se on
erikoistunut siihen mitä tekee. Valmiit optimointialgoritmit tekevät monia tarkistuksia parametreihinsä rutiinin alussa sekä algoritmin laskemisen aikana. Ne ovat
paljon yleiskäyttöisempiä ja sopivat monien tehtävien ratkaisuun.
Menetelmä
Newton (an.)
Newton (num.)
fminsearch
fminunc 1 1
fminunc 2 2
fminunc 3 3
Aika (s)
5
8
12
23
41
200
Taulukko 3.3: Eri menetelmien laskenta-aikoja, 1400 peräkkäistä testiä, otoskoko 7–20
Suurimman uskottavuuden estimaatin löytäminen on nopea toimenpide. Newtonin menetelmän suppeneminen on nopeaa, ja gradienttimenetelmää ei useimmille
otoksille edes tarvitse käyttää.
Omaa toteutusta algoritmista voisi nopeuttaa vielä vaihtamalla ongelma minimointitehtäväksi käyttäen funktiota (3.3) optimoitavana funktiona. Tällöin Newtonin menetelmää voi käyttää optimointiin kun Hessen matriisi on positiivisesti definiitti. Cholesky-hajotelman avulla voi nopeasti tarkistaa positiivisen definiittisyyden
ja saada tarvittavia välituloksia käänteismatriisin muodostamiseksi. [7]
1
Ilman analyyttistä gradienttiä, BGFS
Analyyttinen gradientti, sisäisesti-heijastuva Newton
3
Analyyttinen gradientti ja analyyttinen Hessen matriisi
2
18
4.
YHTEENVETO
Tässä kandidaatintyössä kävimme läpi Cauchyn jakauman ominaisuuksia. Selvitimme jakauman parametrien merkityksen sekä jakauman affiinimuunnoksen. Laskimme jakauman kertymäfunktion ja sen käänteisfunktion. Osoitimme, ettei Cauchyn
jakaumalla ole odotusarvoa eikä muitakaan momentteja.
Esitimme suurimman uskottavuuden menetelmän sekä suurimman uskottavuuden estimaatin ratkaisun Cauchyn jakaumalle. Newtonin menetelmä esiteltiin ja
sille annettiin kaksi eri toteutustapaa.
Näitä toteutustapoja verrattiin keskenään useaan valmiiseen optimointirutiiniin
MATLAB-ohjelmiston kirjastosta. Huomattiin, että itse kirjoitetut erikoistuneet ohjelmanpätkät olivat kaikkein nopeimpia ratkaisussa. Valmisohjelmien käyttö on silti
helpompaa, sillä niille riittää antaa ongelma, itse ratkaisu on jo valmiiksi tehty. Lisäksi kunnollisen optimointirutiinin kirjoittaminen vaatii paljon työtä. Jo valmiilta
näyttävä ratkaisija usein epäonnistuu joillain arvoilla.
Varsinainen suurimman uskottavuuden estimaatti on hyvin tarkka arvio Cauchyn jakauman parametreista useimmilla otoksilla. Cauchyn jakauma sietää hyvin
ulkolaisia. Vain alle kymmenen alkion otoksissa saattaa olla selvä ulkolainen, mutta
estimaatin arvo on silti hyvin lähellä satunnaisesti luodun otoksen paratmetreja.
19
LÄHTEET
[1] Thomas F. Coleman and Yuying Li. On the convergence of interior-reflective
Newton methods for nonlinear minimization subject to bounds. Math. Program., 67:189–224, marraskuu 1994.
[2] J. B. Copas. On the unimodality of the likelihood for the Cauchy distribution:
Some comments. Biometrika, 1975.
[3] L. Elden, L. Wittmeyer-Koch, and H. B. Nielsen. Introduction to numerical
computation - analysis and MATLAB illustrations, kesäkuu 2004.
[4] T. S. Ferguson. Maximum likelihood estimates of the parameters of the Cauchy
distribution for samples of size 3 and 4. Jour. Amer. Statist. Assoc., 1978.
[5] Osmo Kaleva. Numeerinen analyysi 1. http://butler.cc.tut.fi/~kaleva/
Numa1.pdf, kesäkuu 2001. Haettu 18.11.2011.
[6] Jeffrey C. Lagarias, James A. Reeds, Margaret H. Wright, and Paul E. Wright.
Convergence properties of the Nelder–Mead simplex method in low dimensions.
SIAM J. on Optimization, 9:112–147, toukokuu 1998.
[7] K. Madsen, H. B. Nielsen, and O. Tingleff. Methods for non-linear least squares
problems (2nd ed.), 2004.
[8] Mathworks. Find minimum of unconstrained multivariable function. http:
//www.mathworks.se/help/toolbox/optim/ug/fminunc.html, 2011. Haettu
18.11.2011.
[9] Ilkka Mellin.
Todennäköisyyslaskenta.
https://noppa.tkk.fi/noppa/
kurssi/mat-1.2600/materiaali/todennakoisyyslaskenta__osa_2_
_satunnaismuuttujat_ja_todennakoisyysjakaumat.pdf,
elokuu
2008.
Haettu 18.11.2011.
[10] Keijo Ruohonen. Tilastomatematiikka. http://math.tut.fi/~ruohonen/TM.
pdf, 2011. Haettu 18.11.2011.