Opdracht 2 (17 september 2014)

Opdracht 2 (17 september 2014)
Drie onderdelen; stuur de drie bijbehorende bestanden naar [email protected] Deadline maandag 29 september 12:00 uur.
Vooraf: Enige Matlab Tips
1. Een tekstbestand met data inlezen en opslaan in een vector of matrix. Het bestand
http://personal.vu.nl/a.a.n.ridder/numprog/src/proefdata.txt
bevat 30 data, weergegeven in 10 rijen en 3 kolommen. Om met deze data te kunnen
rekenen in een matlabprogramma, moet je eerst het bestand downloaden naar je eigen
drive (H:). NB: bewaar als zogeheten platte tekst. Veronderstel nu dat je met Matlab
in een current directory bezig bent waarin ook het bestand proefdata.txt staat. Het
volgende programma leest het bestand en maakt een 10 × 3 matrix aan waarin de
data komen.
function A = datainlezen
fi = fopen(’proefdata.txt’,’r’);
A = fscanf(fi,’%f’,[3,10])’;
fclose(fi);
end
Toelichting:
• Met fopen wijs je het betreffende bestand aan een variabele toe, hier fi, en geef
je aan of je het bestand wil lezen of schrijven.
• Door fscanf ga je inderdaad het bestand lezen (scannen);
• ’r’ betekent lezen (read ).
• ’%f’ betekent dat de data als reele getallen (floats) worden gelezen (dus niet als
gehele getallen).
• De data worden rijgewijs ingelezen en kolomgewijs opgeslagen. Dus door A =
fscanf(fi,’%f’); krijg je ´e´en kolomvector A met achtereenvolgens (onder elkaar) −4.00, −7.92, 5.23, 4.14, 9.56, . . .. En door A = fscanf(fi,’%f’,[3,10]);
wordt A een 3 × 10 matrix waarin de i-de rij de i-de kolom van de data is. Dus
die moet je transponeren om het goede resultaat te krijgen (zie de apostrofe aan
het eind van fscanf.
• Weet je wel dat de data in drie kolommen staan, maar je weet niet hoe lang die
kolommen zijn, kies dan een voldoende grote lengte, bv (in dit geval) 100: A =
fscanf(fi,’%f’,[3,100])’; het resultaat is toch een 10 × 3 matrix.
2. De lijst met invoerargumenten van een functie kan ook vectoren of matrices bevatten.
Bijvoorbeeld, de matrix A van tip 1 wordt gebruikt in het volgende (deel)programma:
1
function main
A = datainlezen;
x = linspace(1,10,101);
y = funf(x,A);
plot(x,y);
end
function y = funf(x,A)
h1 = prod(prod(abs(A)));
h2 = max(max(A));
y = log(h1)./x + x.^(sqrt(h2));
end
Toelichting:
• De functie datainlezen staat natuurlijk wel in het volledige programma.
Q Q
• prod(prod(abs(A))) berekent i j |aij | en max(max(abs(A))) berekent maxi,j aij .
• Een enkele aanroep van prod (prod(A)) geeft een rijvector van de kolomsgewijze
producten. Idem voor max.
• Een alternatief programma is het product en het maximum in main te berekenen
en in de argumenten van funf aan te roepen. Dan hoeft niet meer de matrix A
in de argumentenlijst:
function main
A = data_inlezen;
p = prod(prod(abs(A)));
m = max(max(A));
x = linspace(1,10,101);
y = funf(x,p,m);
plot(x,y);
end
function y = funf(x,a,b)
y = log(a)./x + x.^(sqrt(b));
end
P
3. Als x een vector is van lengte n, en je wilt ni=1 log(xi ); in Matlab sum(log(x));
4. In Matlab probeer je zoveel gevectoriseerde berekeningen te doen. Bijvoorbeeld bij
plotten:
function grafiekplot(a,b)
x = linspace(a,b,501);
y = fun(x);
2
plot(x,y);
end
function y = fun(x) % x mag vector zijn
y = ...;
end
In sommige gevallen is dat deels niet mogelijk. Bijvoorbeeld de functie
n
X
f (x) =
τix .
i=1
Hierin is τ = (τ1 , . . . , τn ) een rij gegeven getallen.
function grafiekplot(a,b,tau)
x = linspace(a,b,501);
y = zeros(1,501);
for i=1:501
y(i) = fun(x(i),tau);
end
plot(x,y);
end
function y = fun(x,tau) % x is scalar; tau is vector
y = sum(tau.^x);
end
Toelichting: τ is een vector; elk element wordt tot de macht x (getal) berekend
(vandaar de punt), en dan worden alle τix gesommeerd. Je kunt niet tau.^x doen
met zowel τ als x een vector.
5. Matlab kent de volgende functies:
gamma(x) voor Γ(x); zoals bekend:
∞
Z
Γ(x) =
tx−1 e−t dt.
0
gammaln(x) voor log Γ(x);
d
psi(x) voor dx
log Γ(x) =
d
dx
Γ(x) /Γ(x), ook wel de digamma functie genoemd.
Beschrijving van het probleem
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] = α/λ.
3
Opdracht (a)
Schrijf een Matlabprogramma om in ´e´en figuur de grafieken van de functie f te plotten
op [0.2, 6] voor de combinaties (α = 0.5, λ = 1), (α = 0.5, λ = 3), (α = 2.5, λ = 1),
(α = 2.5, λ = 3). Geef elke grafiek een andere kleur, en geef met een legend aan welke
kleur bij welke combinatie hoort. Schrijf je Matlabprogramma in een scriptfile genaamd
opdr2a jeachternaam.m (vervang jeachternaam door je achternaam in kleinde letters).
Opdracht (b)
In het bestand
http://personal.vu.nl/a.a.n.ridder/numprog/src/xdata.txt
staan 500 waarnemingen van X bij schaalparametert λ = 1 maar bij een onbekende
vormparameter α. De bedoeling is om α te schatten met behulp van de maximum likelihood methode. Noteer L(α) voor de likelihood, `(α) = log L(α) voor de log likelihood,
d
`(α) voor de afgeleide van de log likelihood. In dit onderdeel wordt geen `0 (α) = dα
vraagd een Matlabprogramma te schrijven om de grafieken van `(α) en `0 (α) te plotten.
Werk daartoe eerst de uitdrukkingen van deze twee functies op papier uit voordat je hun
code programmeert. Het plotinterval heeft α als x-as en wel tussen α = 0.5 en α = T .
Experimenteer zelf met T totdat `(α) duidelijk een unimodaal gedrag vertoont. Zet er
ook de ”x-as”bij (bv als de grafiek van de functie f (·) ≡ 0) zodat je kunt zien in de
figuur dat `(α) een maximum heeft waar `0 (α) = 0. Schrijf je Matlabprogramma in een
scriptfile genaamd opdr3b jeachternaam.m (vervang jeachternaam door je achternaam
in kleinde letters). In dat programma heeft T natuurlijk een door jou gevonden waarde.
Onderdeel (c)
In dit onderdeel bepaal je het maximum α∗ van `(α) in het interval [0, T ] op twee
manieren:
(i). Pas de Matlabfunctie fminbnd toe op `(α), zie bv het voorbeeld tijdens het eerste
college.
(ii). Pas bisectie toe op `0 (α). De code van bisectie moet je zelf programmeren en toepassen op `0 (α). Zie bijvoorbeeld de code van de false position methode toegepast
tijdens het college:
http://personal.vu.nl/a.a.n.ridder/numprog/src/falseposvb.m.txt
Rapporteer in de output van beide methoden: de gevonden optimale parameter α∗ in 6
decimalen en de waarde van `(α∗ ) in 6 decimalen. En van methode (ii) (bisectie) druk
ook het aantal iteraties en de lengte van het laatste inklemmingsinterval af. Tenslotte,
vergelijk de twee resultaten ook door α = E[X] te schatten via het gemiddelde van de
data. Schrijf je Matlabprogramma in een scriptfile genaamd opdr3c jeachternaam.m
(vervang jeachternaam door je achternaam).
4

Download Report