Opdracht 3 (1 oktober 2014)

Opdracht 3 (1 oktober 2014)
Drie onderdelen; stuur de drie bijbehorende bestanden naar [email protected]
Deadline maandag 13 oktober 12:00 uur.
Inleiding
Polynomiale interpolatie wordt ook gebruikt om de afgeleide f 0 (x) numeriek te berekenen
in een getal x zonder de formule van f 0 (·) eerst analytisch af te leiden. De functie f (·)
wordt natuurlijk wel gegeven als een formule of via een algoritme; dwz f (x) kun je voor
elke x direct berekenen. Zie paragraaf 4.3 (tweede deel) voor de theorie.
• Je wilt f 0 (x) berekenen voor een gegeven x. Kies een h > 0 klein, en zet x0 = x − h
en x1 = x + h. Beschouw het 1-ste graads polynoom p1 (t) dat de twee punten
(x0 , f (x0 )) en (x1 , f (x1 )) interpoleert. Bereken dan f 0 (x) als p01 (x). Je krijgt:
f 0 (x) ≈ p01 (x) =
f (x1 ) − f (x0 )
1
=
f (x + h) − f (x − h) .
x1 − x0
2h
(1)
• Je wilt f 0 (x) berekenen voor een gegeven x. Kies een h > 0 klein, en zet
x0 = x − 2h, x1 = x − h, x2 = x + h, x3 = x + 2h.
Beschouw het 3-de graads polynoom p3 (t) dat de vier punten (x0 , f (x0 )), (x1 , f (x1 )),
(x2 , f (x2 )) en (x3 , f (x3 )) interpoleert. Bereken dan f 0 (x) als p03 (x). Je krijgt:
f 0 (x) ≈ p03 (x) =
2
1
(f (x + h) − f (x − h)) −
(f (x + 2h) − f (x − 2h)).
3h
12h
(2)
• Zo kan ook de tweede afgeleide numeriek benaderd worden door (1) tweemaal
toepassen (met h/2):
f 0 (x + h/2) − f (x − h/2)
1 f (x + h) − f (x) f (x) − f (x − h) ≈
−
h
h
h
h
1
= 2 f (x + h) − 2f (x) + f (x − h) .
(3)
h
f 00 (x) ≈
Beschrijving van het probleem
[Zie ook opdracht 2].
Zij X een continue stochastische variabele met kansdichtheidsfunctie
f (x|α, λ) =
λ(λx)α−1 e−λx
,
Γ(α)
x > 0,
waarin vormparameter α > 0 en schaalparameter λ > 0. Dit is de bekende Gamma
verdeling, met verwachting E[X] = α/λ.
1
Opdracht (a)
R∞
Schrijf een Matlabprogramma om de integralen 0 f (x|α, λ) dx numeriek te berekenen
voor de vier combinaties (α = 0.5, λ = 1), (α = 0.5, λ = 3), (α = 2.5, λ = 1), (α =
2.5, λ = 3). Merk op dat er in alle gevallen 1 moet uitkomen (kansdichtheden). De
numerieke berekeningen zullen niet precies 1 geven, want je moet de bovengrens ∞
vervangen door een eindig getal M (groot genoeg), en in sommige gevallen moet je de
ondergrens 0 vervangen door een positief getal δ (klein genoeg). Voor alle gevallen:
gebruik de trapz functie van Matlab (met n trapezia), en de quad of quadl functie (zie
slides of boek Pratap voor hoe je die functies programmeert). Experimenteer eerst zelf
met diverse n, M en δ (indien van toepassing) voordat je die definitief hebt gekozen,
eventueel verschillend voor verschillende integralen. Schrijf je Matlabprogramma in een
scriptfile genaamd opdr3a jeachternaam.m (vervang jeachternaam door je achternaam
met kleine letters).
Opdracht (b)
In het bestand
http://personal.vu.nl/a.a.n.ridder/numprog/src/xdata.txt
staan n = 500 waarnemingen van X bij schaalparameter λ = 1 maar bij een onbekende
vormparameter α. Noteer L(α) voor de likelihood, dus
L(α) =
n
Y
xα−1 e−xi
i
i=1
Γ(α)
Noteer verder `(α) = log L(α) voor de log likelihood, en `0 (α) =
van de log likelihood. Bereken `0 (α) op drie manieren:
d
dα `(α)
voor de afgeleide
(i). Zoals in opdracht 2b, dus met de exacte formule voor `0 (α).
(ii). Als p01 (α) zoals in (1) (dus gebruik makend van de formule voor de log likelihood
`(α)).
(iii). Als p03 (α) zoals in (2).
Doe dit voor h = 0.01 en α ∈ [0.1, 6]. Gevraagd wordt dan om de grafieken van de
foutfuncties `0 (α) − p01 (α) en `0 (α) − p03 (α) te plotten. Schrijf je Matlabprogramma in een
scriptfile genaamd opdr3b jeachternaam.m (vervang jeachternaam door je achternaam
in kleine letters).
Opdracht (c)
In dit onderdeel ga je nog eens het maximum van de log likelihood `(α) bepalen als het
nulpunt van de afgeleide `0 (α). Maar pas nu de Newton methode toe. Iteratieformule:
αk+1 = αk −
2
`0 (αk )
.
`00 (αk )
Berekening van `0 (α) doe je als in (2), en de berekening van `00 (α) als in (3). Begin met
α0 = 0.1. Uitvoer: het maximum α∗ en bijbehorende `0 (α∗ ) in 10 decimalen, en het
aantal iteraties tot dat antwoord. Schrijf je Matlabprogramma in een scriptfile genaamd
opdr3c jeachternaam.m (vervang jeachternaam door je achternaam met kleine letters).
3

Download Report