Scriptie - TU Delft Institutional Repository

Technische Universiteit Delft
Faculteit Elektrotechniek, Wiskunde en Informatica
Delft Institute of Applied Mathematics
Een statistische analyse van Wicksells
lichaampjes probleem
(Engelse titel: A statistical analysis of Wicksell’s
corpuscle problem)
Verslag ten behoeve van het
Delft Institute of Applied Mathematics
als onderdeel ter verkrijging
van de graad van
BACHELOR OF SCIENCE
in
TECHNISCHE WISKUNDE
door
Ramona Esmeralda Boes
Delft, Nederland
Augustus 2014
c 2014 door Ramona Esmeralda Boes. Alle rechten voorbehouden.
Copyright BSc verslag TECHNISCHE WISKUNDE
Een statistische analyse van Wicksells lichaampjes probleem
(Engelse titel: A statistical analysis of Wicksell’s corpuscle problem)
RAMONA ESMERALDA BOES
Technische Universiteit Delft
Begeleider
Prof.dr.ir. G. Jongbloed
Overige commissieleden
Dr.ing. D. Jeltsema
Dr.ir. M. Keijzer
Augustus, 2014
Delft
Samenvatting
Een statistische analyse van Wicksells lichaampjes probleem
In deze bachelorscriptie wordt het lichaampjes probleem van de Zweedse wiskundige Sven Dag
Wicksell geanalyseerd. In verschillende organen in het menselijk lichaam bevinden zich bolvormige lichaampjes, Wicksell bestudeerde de bolletjes in je milt. Dit deed hij aan de hand van
data verkregen door post mortem doorsnijdingen van verschillende milten. Op deze doorsnedes
namen anatomen de cirkelvormige profielen van de bolletjes waar en noteerden hiervan de diameters. Het doel van Wicksell was het vinden van een schatting voor de verdelingsfunctie van
de diameters van de bollen, aan de hand van de diameters van de geobserveerde cirkelprofielen.
In deze scriptie wordt beschreven hoe Wicksell zijn schatter algebra¨ısch afleidde en hoe hij een
kant en klare methode ontwikkelde die ook door niet-wiskundigen toegepast kon worden. Vervolgens wordt deze methode toegepast op de beschikbare data om zo tot de door Wicksell gevonden
verdelingsfunctie te komen.
Ook worden er moderne schatters afgeleid en toegepast en wordt er een zogenoemde histogramschatter geconstrueerd om ook in moderne notatie een verdelingsfunctie van de bollen aan de
hand van de echte data te vinden. Vervolgens worden de resultaten van beide schattingsmethoden vergeleken en kan worden geconcludeerd dat Wicksell in zijn tijd een heel aardige schatting
had gemaakt.
Ten slotte wordt nog de situatie bekeken waarin de lichaampjes niet netjes bolvormig, maar
ellipso¨ıden zijn. Hierbij wordt eerst afgeleid welke relaties van de bollen ook voor ellipso¨ıde
lichamen geldig zijn en vervolgens wordt algebra¨ısch aangetoond dat de doorsnede van een ellipso¨ıde een ellips is
door Ramona Esmeralda Boes
Inhoudsopgave
Samenvatting
1
1 Inleiding
3
2 Bolvormige lichamen
2.1 Het probleem . . . . . . . . . . . . . .
2.2 Wicksells oplossing . . . . . . . . . . .
2.3 Wicksells methode . . . . . . . . . . .
2.3.1 Voorbeeld . . . . . . . . . . . .
2.4 Moderne notatie en aanpak . . . . . .
2.5 Schatter toepassen . . . . . . . . . . .
2.5.1 Voorbeeld: F is exponentieel .
2.5.2 Voorbeeld: Andere verdelingen
2.6 Echte data . . . . . . . . . . . . . . . .
2.6.1 Wicksells schatting . . . . . . .
2.6.2 Nieuwe schatter . . . . . . . . .
2.6.3 Schatters vergelijken . . . . . .
.
.
.
.
.
.
.
.
.
.
.
.
4
4
4
6
8
8
11
12
13
14
15
16
19
3 Ellipso¨ıde lichamen
3.1 Ellipso¨ıden . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.2 Wicksells relatie . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.3 Stereologie . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
20
20
21
23
4 Discussie
25
A Berekening doorsnijding ellipso¨ıde met vlak
26
B R code
B.1 Verdelingsfunctie Wicksell (Figuur 2.7)
B.2 Schatten exponenti¨ele verdeling (Sectie
B.3 Rejection sampling (Sectie 2.5.2) . . .
B.4 Histogramschatter (Sectie 2.6.2) . . . .
28
28
28
29
30
Bibliografie
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
. . . .
2.5.1)
. . . .
. . . .
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
32
2
Hoofdstuk 1
Inleiding
In het begin van de jaren ’20 werd de Zweedse wiskundige Sven Dag Wicksell op de universiteit van Lund geconfronteerd met een interessant probleem. In verschillende organen in het
menselijk lichaam bevinden zich zogenaamde bolvormige lichaampjes, zoals de eilandjes van
Langerhans in de alvleesklier of de lichaampjes van Hassall in de thymus, maar Wicksell was
vooral ge¨ınteresseerd in de bolletjes in de milt. Dit omdat de anatoom T. Hellman experimenten
uitvoerde waarbij hij post mortem dwarsdoorsnedes van verschillende milten maakte. Hiermee
doorsneed hij ook de in de milt aanwezige bolletjes en op zijn doorsnedes waren dus cirkelvormige
profielen te zien. Hieruit concludeerde hij dat het aantal bolletjes per milt ontzettend verschilt
en dat ook de grootte van deze bolletjes sterk varieert. Omdat je milt echter niet doorzichtig
is, kon Hellman alleen de cirkelvormige profielen waarnemen en niet de bolletjes zelf. Hij wilde
echter wel graag weten hoe de bolletjes er uit zagen, dus hij noteerde de diameters van de door
hem waargenomen cirkels en schakelde de hulp van de wiskundige Wicksell in.
Het doel voor Wicksell was het vinden van een beschrijving van de bollen terwijl alleen de
diameters van de zichtbare cirkelprofielen bekend waren. Omdat dit praktische probleem voor
een wiskundige natuurlijk slecht gedefinieerd is, maakte Wicksell een model om de situatie te
beschrijven, waarover hij in 1925 zijn eerste artikel ”The corpuscle problem: A mathematical
study of a biometric problem” schreef. [1]
In het model gaat hij uit van R3 waarbij de punten volgens een homogeen Poisson proces
verdeeld zijn. Binnen deze omgeving bevinden zich bollen van verschillende groottes, waarbij
verondersteld wordt dat de bollen bij benadering uniform verdeeld zijn over het orgaan en dat
de diameters van deze bollen onafhankelijk van elkaar zijn. Verder nam Wicksell aan dat de
bollen ook netjes bolvormig zijn en ging hij er vanuit dat de diameters van de bollen verdeeld
zijn volgens een verdelingsfunctie F , dit is de functie die hij wilde schatten. Echter kunnen zoals
eerder genoemd deze diameters niet direct waargenomen worden, maar alleen de cirkelvormige
profielen in de dwarsdoorsnede, waarbij we een dwarsdoorsnede in dit model kunnen beschouwen
als het xy-vlak in R3 .
Het doel is dus het schatten van een verdelingsfunctie F voor de diameters van de bollen, aan
de hand van de diameters van de zichtbare cirkelprofielen.
3
Hoofdstuk 2
Bolvormige lichamen
2.1
Het probleem
Zoals in de inleiding beschreven hebben we in dit model te maken met twee soorten diameters,
de nog onbekende diameters van de bollen en de diameters van de waargenomen cirkelprofielen.
Het is duidelijk dat de diameters van deze cirkels over het algemeen niet gelijk zijn aan de
diameters van de bollen, dit is alleen het geval wanneer de doorsnede precies door het midden
van bol gaat. Omdat je de bollen echter voor de doorsnijding niet kunt waarnemen, is de kans
groot dat dit niet het geval is. Hoe verder het vlak van de doorsnede bij het midden van de bol
vandaan ligt, hoe kleiner de diameter van de cirkels zijn in vergelijking met de echte diameters
van de bollen. Omdat de afstand tussen het doorsnijdingsvlak en het midden van de bol over het
algemeen onbekend is, gaan we ervan uit dat de cirkeldiameters en de diameters van de bollen
verschillend zijn. De cirkeldiameters zullen we de zichtbare diameters noemen en diameters van
de bollen de echte diameters. Dus het probleem kan nu als volgt geformuleerd worden:
Het vinden van de verdeling van de echte diameters van de bollen uit de verdeling van de zichtbare diameters in de microtome doorsnijdingen.1
Dat deze twee verdelingen echt verschillend zijn volgt uit het hierboven beschreven feit dat de
zichtbare diameters over het algemeen kleiner zijn dan de echte diameters en dat een willekeurige doorsnede een relatief groter aantal grote bollen zal bevatten dan kleine bollen, omdat grote
bollen vaker in het bereik van een doorsnijding liggen.
Ten slotte moeten we nog opmerken dat als je series van opeenvolgende doorsnijdingen zou
nemen in plaats van enkele doorsnijdingen, dat je de echte diameters direct zou kunnen waarnemen. Echter is het maken van een doorsnijding over het algemeen erg omslachtig en omdat
er honderden milten onderzocht moeten worden is deze methode in de praktijk dus niet toepasbaar. Het wordt alleen gebruikt in speciale gevallen, bijvoorbeeld om de theoretische aannames
te controleren.
2.2
Wicksells oplossing
Het doel van Wicksell was dus het uitdrukken van de verdeling van de diameters van de bollen
in termen van de verdeling van de diameters van de cirkelprofielen. Hiertoe definieerde hij [1]
1
Een microtoom is een apparaat waarmee zeer nauwkeurig heel dunne plakjes kunnen worden gesneden voor
het maken van microscopische preparaten.
4
met F (r) de relatieve frequentiefunctie van de diameters van de bollen, waarbij r de grootte van
de diameters is. F (r) is hier dus de gezochte functie.
Als we dan aannemen dat R de grootste diameter is die is waargenomen, dan hebben we:
Z R
F (r)dr = 1
0
Vervolgens definieerde Wicksell met N het aantal bollen per volume-eenheid waaruit volgt dat
het aantal bollen per volume-eenheid met hun diameter in het interval r tot r + dr dat binnen
het bereik van een willekeurige doorsnijding valt, gelijk is aan
N rF (r)dr
Zij r0 nu de gemiddelde diameter van de bollen, dan zal N r0 het totaal aantal cirkelvormige
profielen in de doorsnijding zijn en daaruit volgt dan:
N r0 f (r)dr = N rF (r)dr,
waarbij f (r)dr de relatieve frequentie is van de bollen met diameter in het interval r tot r+dr die
in de doorsnede te zien zijn. Want het aantal bollen binnen het bereik van een doorsnijding moet
gelijk zijn aan het aantal cirkelvormige profielen dat in de doorsnijding wordt waargenomen.
Hieruit volgt nu dat
f (r) =
rF (r)
r0
(2.1)
De kans dat een bol in de doorsnijding zijn centrum op een afstand y tot y + dy van het
doorsnijdingsvlak heeft is gelijk aan
2dy
r
En dus is het relatieve aantal bollen met diameter in het interval r tot r + dr en hun centrum
op een afstand y tot y + dy van het doorsnijdingsvlak gelijk aan
2
f (r)dydr
r
Vervolgens noteerde Wicksell de zichtbare diameters van de waargenomen cirkels met x waaruit
hij concludeerde dat
x2 = r2 − 4y 2
Hieruit volgt direct dat
1p 2
r − x2
2
Nu geldt dat het relatieve aantal bollen met echte diameter in het interval r tot r + dr en de
zichtbare diameter in het interval x tot x + dx gelijk is aan:
y=
2
2
dy
f (r)drdy = f (r)dr dx
r
r
dx
2
x
= f (r)dr √
dx
2
r
2 r − x2
1
x
= f (r)dr √
dx
2
r
r − x2
1
x
= F (r)dr √
dx
2
r0
r − x2
5
Als we nu ten slotte φ(x) defini¨eren als de relatieve frequentie functie van de zichtbare diameters,
dan vinden we
Z
1
x R
F (r) √
φ(x) =
dr
(2.2)
r0 x
r 2 − x2
En omdat φ(x) de functie is die bekend is, hebben we nu een uitdrukking gevonden voor de
gezochte functie F (r).
2.3
Wicksells methode
Omdat Wicksell zijn artikel echt schreef voor een praktisch probleem, concludeerde hij dat
eerder beschreven afleidingen wel degelijk heel belangrijk waren, maar vanwege de benodigde
wiskundige vaardigheden in de praktijk waarschijnlijk niet toepasbaar zouden zijn. Daarom
ontwikkelde hij een kant en klare methode die voor de anatomen stapsgwijs te volgen zou zijn.
Hij ging ervan uit dat de cirkels die de anatomen waarnamen onderverdeeld werden in klassen,
afhankelijk van de grootte van hun diameters, de zogenoemde zichtbare diameters. Hiervoor
gebruikte hij de klassen 1.5 − 2; 2.5 − 3.5; 3.5 − 4.5 mm enz en de frequenties in deze klassen
noteerde hij als f2 , f3 , f4 , . . . , fR .
De diameters in de klassen 0 − 0.5 en 0.5 − 1.5 zijn over het algemeen te klein om waar te
nemen, waardoor de bijbehorende frequenties f 1 en f1 in de praktijk onbekend zijn. Voor de
4
R1
volledigheid nemen we ze toch mee en er geldt f 1 = kr0 02 φ(x)dx.
4
In het algemeen volgt uit vergelijking 2.2
i+ 21
Z
fi = kr0
Z
=k
φ(x)dx
i− 12
i+ 12
i− 12
R
Z
F (r) √
x
x
1
drdx
r 2 − x2
De constante k is hier een proportionele constante en zal in de afleidingen voor de relatieve
frequenties uiteindelijk verdwijnen.
En als we de integratievolgorde veranderen en integreren over de variabele x dan volgt:
Z
fi = k
i+ 21
Z
p
2
2
drF (r) r − (i − 0.5) + k
R
p
p
drF (r)( r2 − (i − 0.5)2 − r2 − (i + 0.5)2
i+ 12
i− 12
(2.3)
√
R1
RR
En voor f 1 geldt f 1 = k 02 drF (r)r + k 1 drF (r)(r − r2 − 0.25)
4
4
2
Omdat er maar sporadisch bollen met een diameter groter dan 15mm werden waargenomen,
stelde Wicksell de bovengrens op R = 15.5.
Vervolgens definieerde hij
Z
P1 = k
4
1
2
Z
drF (r); Pn =
0
6
n+ 12
n− 21
drF (r)
En om notatie van uitdrukking 2.3 wat compacter te maken definieerde Wicksell ten slotte:
f 1 = P 1 a00 +
4
4
15
X
Pj a0j ; fi =
j=1
15−i
X
Pi+j aij
(2.4)
j=0
waarbij de co¨effici¨enten aij dus gegeven worden door
p
p
a00 = 0.25; a0j = j − j 2 − 0.25; ai0 = i2 − (i − 0.5)2
p
p
aij = ( (i + j)2 − (i − 0.5)2 − ( (i + j)2 − (i + 0.5)2
Het systeem 2.4 is een systeem van geschatte vergelijkingen maar met een duidelijke geometrische betekenis. De hoeveelheden P 1 , P1 , P2 , . . . , P15 zijn proportioneel met het aantal bollen
4
waarvan de diameters exact gelijk zijn aan 41 , 1, 2, . . . , 15 mm respectievelijk. En als doorsneden
met een willekeurig vlak geeft dat de waargenomen verdeling f 1 , f1 , f2 , . . . , f15 van de zichtbare
4
cirkeldiameters.
Wanneer we de numerieke waarden berekenen, wordt het syteem 2.4 gegeven door:
Figuur 2.1: Numerieke waarden van verdeling van de zichtbare cirkeldiameters
Omdat niet de f 0 jes maar juist de hoeveelheden P onbekend zijn, moeten we bovenstaand systeem van lineaire vergelijkingen oplossen. Hierbij worden dus de P 0 s uitgedrukt in f 1 , f1 , f2 , . . . , f15
4
en deze oplossing wordt gegeven door:
Figuur 2.2: Hoeveelheden Pi uitgedrukt in de waargenomen fi
De toepassing van het systeem in Figuur 2.2 is als volgt: Door het invullen van de waargenomen
waarden van fi verkrijg je de bijbehorende waarden van Pi . Hoewel we in de praktijk de waarden
van f 1 en f1 vaak niet weten, zien we dat even goed de vergelijkingen van Figuur 2.2 kunnen
4
gebruiken, behalve de eerste twee. Wanneer je echter de waarden van P2 , P3 , . . . , P15 berekend
7
hebt, is het mogelijk de waarden van P 1 en P1 te vinden door grafische exrapolatie.
4
Definieer nu ten slotte
Fi =
Pi
,
S
waarbij S de som van de hoeveelheden P is.
Dan wordt met F 1 , F1 , F2 , . . . , F15 de relatieve frequenties van de echte boldiameters in de klas4
sen 0 − 0.5; 0.5 − 1.5; 1.5 − 2.5; . . . ; 14.5 − 15.5 genoteerd.
2.3.1
Voorbeeld
Om de toepassing van de vergelijkingen in Figuur 2.2 te verduidelijken, volgt nu een kort voorbeeld. Hierbij nemen we de volgende verdeling:
φ(x) =
x −x2
e 50
25
Daaruit volgt
Z
fi =
i+ 12
i− 21
2
2
−(i− 1
−(i+ 1
x −x2
2)
2)
e 50 dx = e 50 − e 50
25
(2.5)
En als we de numerieke waarden berekenen en vermenigvuldigen met 1000 vinden we:
i
fi
1/4
5
1
39
2
73
3
100
4
116
5
121
6
117
7
105
8
89
9
71
10
54
11
39
12
27
13
18
14
11
15
7
Tabel 2.1: De numerieke waarden in promille van de frequenties van de zichtbare diameters.
Wanneer we bovenstaande waarden invullen in het systeem in Figuur 2.2, dan vinden we:
i
Pi
1/4
0.8
1
6.9
2
11.9
3
16.3
4
19.0
5
19.7
6
19.0
7
16.9
8
14.3
9
11.3
10
8.6
11
6.2
12
4.2
13
2.9
14
1.7
Tabel 2.2: De met de vergelijkingen in Figuur 2.2 berekende hoeveelheden Pi
Wanneer we deze waarden optellen vinden we S = 161.5 en dus volgen ten slotte de volgende
aantallen voor de relatieve frequenties van echte boldiameters, weer uitgedrukt in promille.
i
Fi
1/4
5
1
43
2
74
3
101
4
118
5
122
6
118
7
105
8
89
9
70
10
53
11
38
12
26
13
18
14
11
15
11
Tabel 2.3: De gevonden relatieve frequenties van de echte boldiameters.
2.4
Moderne notatie en aanpak
In deze sectie zullen we Wicksells resultaten omschrijven in meer moderne notatie, waarbij we
de relaties niet langer zullen uitdrukken in termen van de boldiameters maar in termen van de
gekwadrateerde bolstralen. Hierbij nemen we aan dat de stralen onderling onafhankelijk zijn en
dat ze verdeeld zijn volgens een verdelingsfunctie F .
8
15
1.8
Verder zullen we nu ook een expliciete uitdrukking afleiden voor de gezochte functie F , in termen
van de bekende onderliggende verdeling van de cirkelstralen, in plaats van alleen de omgekeerde
relatie die Wicksell afleidde.
Voor de volledigheid volgt hieronder een overzicht van de notatie die Wicksell hanteerde.
N
x
r
r0
F (r)
f (r)
φ(x)
Aantal bollen per volume-eenheid
Zichtbare diameters van de cirkels
Echte boldiameters
Gemiddelde diameter van de bollen
Relatieve frequentie functie echte boldiameters
Relatieve frequentie functie boldiameters in doorsnijding
Relatieve frequentie functie van de zichtbare cirkeldiameters
Tabel 2.4: Notatie van Wicksell
Om in de nieuwe notatie verwarring met Wicksells notatie te voorkomen, zullen we aan alle
grootheden een andere variabele toewijzen. Zo staat B vanaf nu voor de echte boldiameters
en D voor de zichtbare cirkeldiameters. Met X noteren we de gekwadrateerde bolstralen dus
X = (B/2)2 en X heeft dichtheid fX en met Z noteren we de gekwadrateerde cirkelstralen dus
Z = (D/2)2 en Z heeft dichtheid gZ .
Hieronder volgt een volledig overzicht van de nieuwe notatie.
D
B
µB
fB
gD
X
Z
fX
gZ
Zichtbare diameters van de cirkels
Echte boldiameters
Gemiddelde diameter van de bollen
Dichtheid van de echte boldiameters
Dichtheid van de zichtbare cirkeldiameters
Gekwadrateerde bolstralen
Gekwadrateerde cirkelstralen
Dichtheid van de gekwadrateerde bolstralen
Dichtheid van de gekwadrateerde cirkelstralen
Tabel 2.5: Moderne notatie
Zoals we in sectie 2.2 hebben afgeleid, kwam Wicksell uit op de volgende uitdrukking voor de
verdeling van de diameters van de cirkels in termen van de boldiameters.
x
φ(x) =
r0
Z
R
F (r) √
x
1
dr
r 2 − x2
(2.6)
In de moderne notatie wordt de door Wicksell verkregen relatie 2.6 nu:
d
gD (d) =
µB
Z
∞
d
In de huidige notatie volgt hieruit:
9
f (b)
√ B
db
b2 − d2
(2.7)
d
P (Z ≤ z)
dz
d
D2
≤ z)
=
P(
dz
2
√
d
=
P (D ≤ 2 z)
dz
Z 2√z
d
gD (d)dd
=
dz 0
√
= gD (2 z)z −1/2
gZ (z) =
√ Z ∞
2 z
f (b)
√ B
=√
db
√
zµB 2 z b2 − 4z
(2.8)
Nu geldt voor fB (b) dat
fB (b) =
d
b2
d
B2
b2
b2
d
P (B ≤ b) = P (
≤ ) = P (X ≤ ) = fX ( )b
db
db
2
2
db
2
2
Als we dit nu invullen in de zojuist verkregen vergelijking 2.8, dan volgt
Z ∞
2
bfX (b/2)2
√
gZ (z) =
db
µB 2√z b2 − 4z
Z ∞ √
2 yfX (y) 1
2
√
=
√ dy
µB z
2 y−z y
Z ∞
2
fX (y)
√
=
dy
µB z
y−z
Hierbij hebben we de transformatie y = ( 2b )2 gebruikt. Wanneer we nu om de gewenste notatie
te verkrijgen de transformatie y = x toepassen, dan volgt:
1
gZ (z) =
2mF
waarbij mF =
R∞
0
Z
∞
z
f (x)
√X
dx,
x−z
(2.9)
√1 dF (y).
y
We hebben nu de dichtheid van de zichtbare cirkelstralen uitgedrukt in de dichtheid van de
bolstralen, maar omdat we meer ge¨ınteresseerd zijn in de bolstralen willen we liever de inverse
relatie hebben. Deze wordt gegeven door: [3]
R∞
(z − x)−1/2 g(z)dz
F (x) = 1 − xR ∞ −1/2
(2.10)
g(z)dz
0 (z)
We gaan nu laten zien dat deze relatie inderdaad overeenkomt met de gevonden uitdrukking
voor gZ (z) 2.9.
R∞
Definieer nu ψ(x) = x (z − x)−1/2 g(z)dz, dan de inverse relatie is uit te drukken als
F (x) = 1 −
10
ψ(x)
ψ(0)
(2.11)
Verder geldt
Z
ψ(x) =
x
∞
−1/2
(z − x)
Z ∞ Z
1
=
2mF
y=x
y
1
2mF
Z
∞
z
fX (y)
√
dydz
y−z
(z − x)−1/2 (y − z)−1/2 dzfX (y)dy
z=x
Voor de binnenste integraal geldt nu
Z 1
Z y
−1/2
−1/2
((y − x)u)−1/2 ((1 − u)(y − x))−1/2 (y − x)du
(z − x)
(y − z)
dz =
u=0
z=x
Z 1
u−1/2 (1 − u)−1/2 du
=
0
1 1
= B( , )
2 2
=π
Hierbij hebben we de transformatie z = x(1 − u) + yu gebruikt. Vullen we dit resultaat nu in in
bovenstaande vergelijking voor ψ(x), dan volgt:
Z ∞
1
ψ(x) =
πfX (y)dy
2mF y=x
π
(1 − F (x))
=
2mF
Als we dit resultaat nu invullen in vergelijking 2.11, dan volgt
1−
ψ(x)
=1−
ψ(0)
π
2mF
(1 − F (x))
π
2mF
= 1 − (1 − F (x))
= F (x)
Dus de inverse relatie in vergelijking 2.10 klopt.
2.5
Schatter toepassen
In sectie 2.4 hebben we een expliciete relatie afgeleid tussen de dichtheid van de zichtbare
cirkelstralen gz en de dichtheid van de bolstralen F . Hierbij werd de verdeling van de cirkelstralen
uitgedrukt in termen van onbekende verdeling van de bolstralen, dit is wat we het directe
probleem noemen. Vervolgens hebben we een inverse relatie gecontroleerd, waarbij de gezochte
verdeling F uitgedrukt wordt in termen van de bekende verdeling van de zichtbare cirkelstralen
gz . Dit noemen we het inverse probleem.
In deze sectie gaan we F schatten op basis van data gegenereerd door gz en dit noemen we het
statistische inverse probleem.
Het inverse probleem werd in vergelijking 2.10 gegeven door
R∞
(z − x)−1/2 g(z)dz
F (x) = 1 − xR ∞ −1/2
g(z)dz
0 (z)
11
Dit kunnen we schrijven als
T (x)
F (x) = 1 −
waarbij T (x) =
T (0)
Z
x
∞
g(z)
√
dz
z−x
De schatter voor T is nu gedefinieerd als
Z
Tn (x) =
x
∞
n
dGn (z)
1X
√
(Zi − x)−1/2 1(x,∞) (Zi )
=
n
z−x
i=1
(2.12)
Deze laatste uitdrukking kunnen we implementeren in R en hiermee kunnen we dus een schatting
maken van de verdeling van de bolstralen als de verdeling van de zichtbare cirkelstralen bekend
is.2
2.5.1
Voorbeeld: F is exponentieel
We zullen nu het gebruik van deze schatter illustreren en zijn kwaliteit beoordelen.
Hiertoe nemen we even aan dat F een exponenti¨ele verdeling is, dus stel
F (x) = 1 − e−λx
Dan weten we dat f (x) = λe−λx .
Hieruit volgt via vergelijking 2.9 dat voor de onderliggende verdeling van de cirkelstralen dan
geldt
Z ∞
1
f (x)
√X
g(z) =
dx
2mF z
x−z
r Z
2λ ∞
=
λ(x − z)−1/2 e−λx dx
π z
= λe−λz ,
waarbij we de transformatie x = z + y(1 − z) gebruikt hebben.
Dus als F (x) exponentieel verdeeld is, dan is g(z) dat ook. En dit is precies de reden waarom
wij dit als voorbeeld gebruiken.
Figuur 2.3 toont de schatting voor F en de echte verdelingsfunctie van F , waarbij we hebben
aangenomen dat λ = 2.
2
R is een softwarepakket en programmeertaal ontwikkeld voor statistiek en data-analysedoeleinden
12
Figuur 2.3: De grillige, zwarte, lijn is de gevonden schatter voor F. De rechte, rode, lijn is de
echte verdeling van F.
Zoals te zien in de figuur komt de schatter aardig in de buurt van de echte verdelingsfunctie,
maar is hij nog wel wat grillig en niet monotoon. Een monotone schatter kan worden verkregen
via ’isotonization’: [3]
Neem de rechter afgeleide van kleinste concave majorant van de functie Un gegeven door
Z x
Un (x) =
Tn (y)dy
0
De schatter voor F wordt nu gegeven door
U 0 (z)
F˜ = 1 − n0
Un (0)
Maar hier wordt in deze scriptie niet in verder detail op in gegaan.
2.5.2
Voorbeeld: Andere verdelingen
In het voorgaande voorbeeld hebben we aangenomen dat F exponentieel verdeeld was waardoor we de kwaliteit van onze schatter konden beoordelen. Wanneer F echter niet exponentieel
verdeeld is, kan niet op algebra¨ısche wijze een nette verdeling voor de bijbehorende g worden
gevonden. In het algemeen is het juist g die bekend is en dan weten we dus niet welke verdeling
voor F hierbij hoort. Maar natuurlijk kunnen we op iedere verdeling van g wel onze schatter
toepassen om een geschatte verdeling voor F te vinden, dit is precies waar een schatter voor is.
Omdat R echter niet elke dichtheid kent, kun je hier niet zo automatisch uit trekken als uit
de exponenti¨ele verdeling. Dit is op te lossen door rejection sampling te gebruiken, hiermee
genereer je observaties uit een zelf gekozen verdeling. Deze methode werkt voor alle verdelingen
in R3 met een dichtheid.
Ter illustratie zullen we nu met rejection sampling trekken uit de volgende verdeling voor g(z):
g(z) =
3√
1 − z1[0,1] (z)
2
13
In Figuur 2.4 is te zien hoe goed de dichtheid waar uit getrokken wordt met rejection sampling
de echte dichtheid benadert.
Figuur 2.4: Dichtheid van g(z)
Wanneer we onze schatter toepassen op g(z) door met rejection sampling uit deze dichtheid te
trekken, krijgen we onderstaande schatting voor de bijbehorende verdeling voor F .
Figuur 2.5: Geschatte verdelingsfunctie van F
2.6
Echte data
In deze sectie zullen we de data beschouwen die Wicksell ook tot zijn beschikking had, dit zijn
de resultaten van de experimenten die anatoom T. Hellman uitvoerde. Hij noteerde de grootte
van de cirkels die hij waarnam na post mortem dwarsdoorsnedes van verschillende milten en
deelde deze cirkeldiameters in in klassen.
14
We zullen eerst bekijken hoe Wicksell aan de hand van de beschikbare data, zie tabel 2.6, via
zijn kant en klare methode een schatting maakte voor de verdeling van de echte boldiameters.
Vervolgens zullen we met een nieuwe schatter, in termen van de gekwadrateerde bolstralen, een
schatting maken voor de verdelingsfunctie.
i
frequentie
2
52
3
146
4
197
5
210
6
184
7
143
8
95
9
57
10
31
11
15
12
7
13
4
14
2
15
1
Tabel 2.6: Frequenties waargenomen cirkeldiameters in mm
Zoals eerder genoemd werden er in de klassen 0 − 0.5 en 0.5 − 1.5 geen cirkelprofielen waargenomen, daarom komen deze niet in Tabel 2.6 voor.
2.6.1
Wicksells schatting
Voor de klasse 1.5 − 2.5 werd een frequentie van f2 = 52 gevonden, maar deze waarde levert
een negatieve waarde van P2 op, dit komt waarschijnlijk doordat de waarnemingen ook bij deze
diametergrootte niet compleet zijn.
Daarom gebruikte Wicksell alleen de waarden van f3 tot en met f1 5 en vulde deze in in het
systeem van Figuur 2.2. Hierbij vond hij de volgende waarden voor de Pi .
i
Pi
3
15.4
4
32.1
5
41.5
6
36.7
7
30.4
8
20.4
9
12.1
10
6.6
11
3.1
12
1.3
13
0.6
14
0.4
15
0.3
Tabel 2.7: De met de vergelijkingen in Figuur 2.2 berekende hoeveelheden Pi
Na het plotten van de gevonden waarden en extrapolatie van deze grafiek, zoals te zien in Figuur
2.6, vond hij de volgende waarden:
P 1 = 1; P1 = 3; P2 = 8
4
Figuur 2.6: Extrapolatie van de berekende waarden voor Pi
Nu alle waarden voor Pi gevonden zijn, is ook de waarde van de som bekend en deze is gelijk aan
S = 212.9. Hiermee vindt Wicksell ten slotte de volgende waarden voor de relatieve frequenties
15
van de echte boldiameters, uitgedrukt in promille.
i
Fi
1/4
5
1
14
2
38
3
72
4
151
5
194
6
172
7
143
8
96
9
57
10
31
11
15
12
6
13
3
14
1
15
1
Tabel 2.8: De gevonden relatieve frequenties van de echte boldiameters.
De gevonden frequenties leiden tot de volgende schatting voor de verdelingsfunctie van de echte
boldiameters.
Figuur 2.7: Verdelingsfunctie van de boldiameters.
2.6.2
Nieuwe schatter
Omdat van de zichtbare cirkeldiameters in dit geval geen verdeling bekend is, maar alleen in
klassen ingedeelde groottes, kunnen we op de echte data niet de schatter toepassen van sectie
2.5. Daarom zullen we nu een zogenoemde histogramschatter afleiden en toepassen.
Deze schatter willen we weer uitdrukken in termen van de gekwadrateerde cirkelstralen, maar
omdat de originele data in tabel 2.6 de frequenties van de cirkeldiameters weergeeft, zullen we
deze data eerst transformeren. Dit resulteert in de histogram in Figuur 2.8.
16
Figuur 2.8: Frequentiehistogram van de
gekwadrateerde bolstralen
Figuur 2.9: Kanshistogram van de gekwadrateerde bolstralen
Wanneer we nu de grenzen van het histogram noteren met ti (i = 1, 2, . . . , m) en aannemen dat
t0 = 0 en tm+1 = ∞, dan wordt de histogramdichtheid van Figuur 2.9 gerepresenteerd door [3]
g(x) =
m−1
X
gi 1[ti ,ti+1 ) (x), x ≥ 0
i=1
Waarbij geldt dat
gi =
1
n #{j
: Zj ∈ [ti , ti+1 )}
ti+1 − ti
Merk op dat g0 = gm = 0.
De histogramschatter voor F wordt nu gegeven door
Vn (z)
Fn (z) = 1 −
waarbij Vn (z) =
Vn (0)
Z
z
∞
g(y)
√
dy
y−z
Dus voor z ∈ [ti , ti+1 ) geldt nu
P
√
√
√
gi ti+1 − z + m−1
j=i+1 gj ( tj+1 − z − tj − z)
Fn =
Pm−1 √
√
i=1 gi ( ti+1 − ti
17
(2.13)
Wanneer we deze schatter nu toepassen op de histogramdichtheid van figuur 2.9, dan vinden we
de volgende schatting voor de verdelingsfunctie van de gekwadrateerde bolstralen.
Figuur 2.10: Geschatte verdelingsfunctie van de gekwadrateerde bolstralen
Merk op dat deze schatting geen nette verdelingsfunctie beschrijft, want hij neemt waarden
onder nul aan en begint dalend.
Evenals in sectie 2.5 kan een monotone schatter worden verkregen via ’isotonization’: [3]
Neem de rechter afgeleide van kleinste concave majorant van de functie Un gegeven door
Z z
Un (z) =
Vn (y)dy
0
De schatter voor F wordt nu gegeven door
U 0 (z)
F˜ = 1 − n0
Un (0)
Maar ook hier wordt niet in verder detail op dit proces in gegaan.
18
2.6.3
Schatters vergelijken
Om de zojuist gevonden geschatte verdelingsfunctie (Figuur 2.10) te kunnen vergelijken met de
verdelingsfunctie die Wicksell vond (Figuur 2.7), transformeren we Wicksells functie in diameters eerst naar een functie uitgedrukt in de gekwadrateerde bolstralen. Wanneer we vervolgens
beide verdelingsfunctie plotten, krijgen we als resultaat Figuur 2.11.
Figuur 2.11: De zwarte doorgetrokken lijn is de verdelingsfunctie van Wicksell, de rode lijn met
sterretjes is de geschatte verdelingsfunctie met de nieuwe schatter.
In de figuur is te zien dat onze schatter de verdelingsfunctie van Wicksell vrij aardig benadert.
Omdat onze schatter echter geen nette monotone verdelingfunctie oplevert, kunnen we hier geen
heel sterke conclusies aan verbinden. We kunnen ons echter niet aan de indruk onttrekken dat
wanneer onze schatter ook netjes in de oorsprong zou beginnen, de twee verdelingsfuncties wel
heel sterk op elkaar zouden lijken. We mogen dus wel concluderen dat Wicksell voor in zijn tijd
en met zijn kant en klare methode een heel aardige schatting voor de verdelingsfunctie van de
echte boldiameters heeft weten te produceren.
19
Hoofdstuk 3
Ellipso¨ıde lichamen
In het vorige hoofdstuk zijn we ervan uitgegaan dat de lichamen netjes bolvormig zijn. Dit
maakte het rekenwerk een stuk eenvoudiger, echter hoeft dit in de praktijk zeker niet zo te zijn.
Veel realistischer is het om te stellen dat de lichamen ellipso¨ıden zijn. Waar je voor de definitie
van een bol echter alleen een straal en een middelpunt nodig hebt, ligt dit bij een ellipso¨ıde een
stuk ingewikkelder.
3.1
Ellipso¨ıden
Een ellipso¨ıde wordt vastgelegd door zijn vorm, zijn ligging en zijn ori¨entatie, want anders dan
bij een bol heeft rotatie op een ellipso¨ıde wel invloed. De vorm wordt uitgedrukt door de lengtes van de drie hoofdassen, s1 , s2 en s3 . De ligging wordt bepaald door het centrum van de
ellipso¨ıde, (u1, u2, u3) en de ori¨entatie wil zeggen de stand van de hoofdassen ten opzichte van
de driedimensionale x, y en z as ( Zie Figuur 3.1)
Figuur 3.1: Ellipso¨ıde
20
In zijn meest algemene vorm is een ellipso¨ıde in wiskundige termen als volgt te beschrijven:




x1 − u1


E = (x1 , x2 , x3 ) ∈ R3 : x1 − u1 x2 − u2 x3 − u3 A x2 − u2  = c ,
(3.1)


x3 − u3
waarbij A een symmetrische en positief definiete 3x3 matrix is.
In deze notatie bepaalt de constante c samen met de diagonaalelementen van de matrix A de
diameters van de ellipso¨ıde. Hierbij geldt bijvoorbeeld dat de diameter in de richting van s1
gelijk is aan c/a11 .
De niet-diagonaaltermen van de matrix A bepalen de ori¨entatie van de ellipso¨ıde, dus wanneer
A een diagonaalmatrix is dan wijzen de hoofdassen van de ellipso¨ıde in dezelfde richtings als de
x, y en z as. (Zie Figuur 3.2). Merk op dat als de diagonaalelementen ook nog aan elkaar gelijk
zijn, dan beschrijft vergelijking 3.1 een bol.
Zoals eerder genoemd bepalen de constanten u1 , u2 en u3 het centrum van de ellipso¨ıde, dus
wanneer deze constanten gelijk zijn aan nul, is het centrum van de ellipso¨ıde gelijk aan de oorsprong. De ellipso¨ıde in Figuur 3.3 wordt dus beschreven met de vergelijking

x1 x2 x3
 
a11 0
0
x1
 0 a22 0  x2  = c
0
0 a33
x3
Figuur 3.2: A is een diagonaalmatrix
3.2
Figuur 3.3: (u1 , u2 , u3 ) = 0
Wicksells relatie
In zijn tweede artikel ”The corpuscle problem: Second memoir: Case of ellipsoidal corpuscles”[2]
bestudeert Wicksell de situatie waarin de lichamen niet bolvormig maar ellipso¨ıden zijn. Omdat dit echter vooral bestaat uit het beschouwen van verschillende gevallen - langwerpige en
afgevlakte ellipso¨ıden vragen bijvoorbeeld al om een verschillende aanpak - en niet leidt tot een
algemene uitdrukking voor de verdeling van de ellipso¨ıden, zal dit in deze scriptie niet beschreven worden.
21
Wat wel interessant is is dat Wicksell beweert dat de relatie in vergelijking 2.2 die hij voor de
bollen heeft afgeleid, ook voor ellipso¨ıde lichamen geldig is. Hierbij hoeven de ellipso¨ıden niet
van dezelfde vorm te zijn en zelfs hun ori¨entatie mag verschillen zijn.
Beschouw dus een willekeurige ellips waarvan twee doorsnedes gemaakt worden, ´e´en door het
centrum van de ellipso¨ıde en ´e´en op een afstand y van het centrum.
Definieer nu de zichtbare diameter x als het meetkundig gemiddelde van de grootste en de kleinste diameter van de ellips te zien in de willekeurige doorsnede en r als het meetkundig gemiddelde
van de ellips te zien in de doorsnede door het centrum van de ellipso¨ıde. De oppervlakte van
deze twee ellipsen noteren we met Ay en A0 respectievelijk en h is de afstand tussen het centrum
van de ellipso¨ıde en het horizontale raakvlak. Dan geldt er altijd:
Ay = A0 (1 −
y2
)
h2
Dit omdat voor ellipso¨ıden geldt dat de oppervlaktes van twee parallelle doorsnedes proportioneel
zijn met de oppervlaktes van de dezelfde doorsnedes in een bol.[4]
Nu geldt dat
πr2 = 4A0 en πx2 = 4Ay
En hieruit volgt
y2
y2
πx2 = 4A0 (1 − 2 ) = πr2 (1 − 2 )
h
h
r
2
y
⇔x=r 1− 2
h
x
y2
⇔ ( )2 = 1 − 2
r √
h
y
r 2 − x2
⇔ =
h
r
Hieruit volgt nu
f (r)dr|
dy
xdx
|dx = f (r)dr √
hdx
r r 2 − x2
En daaruit volgt nu ten slotte, net als bij de bollen,
Z
φ(x) = x
x
R
1
f (r) √
dr
r r2 − x2
Gebruik makend van vergelijking 2.1 is dit dus exact dezelfde uitdrukking voor F als degene die
Wicksell voor bolvormige lichamen vond (vergelijking 2.2).
Omdat in de afleiding van de schatters wel vereist is dat de lichamen bolvormig zijn, zullen de
in hoofdstuk 2 gevonden schatters voor ellipso¨ıde lichamen niet geldig zijn en kunnen wij deze
dus niet toepassen. Wel zullen we in de volgende sectie laten zien dat wanneer men ellipso¨ıde
lichamen doorsnijdt, er op de doorsnede ellipsvormige profielen te zien zijn.
22
3.3
Stereologie
Net als bij de bollen zullen we de ellipso¨ıde doorsnijden met een willekeurig vlak, zie Figuur 3.4.
Maar waar het bij de doorsnijding van een bol vrij duidelijk is dat er dan een cirkel ontstaat, is
dit bij een ellipso¨ıde nog niet zo makkelijk. Figuur 3.5 doet wel vermoeden dat er op het snijvlak
een ellips ontstaat, maar omdat het uit de figuur niet helemaal helder is dat dit inderdaad een
ellips is gaan we dit algebra¨ısch aantonen.
Figuur 3.4: Doorsneden ellipsoide
Figuur 3.5: Ellips
Hiertoe nemen we een algemene vorm van de ellipso¨ıde en doorsnijden deze met het vlak z = h,
dit omdat het resultaat dan netjes in het xy-vlak ligt en de interpretatie dan duidelijker is.
Natuurlijk heeft dit hetzelfde resultaat als doorsnijden met een vlak in een andere richting,
evenzo is doorsnijden op hoogte h hetzelfe als doorsnijden op een andere hoogte wanneer we de
ellipso¨ıde diezelfde afstand omhoog schuiven. Ten slotte hoeven we om dezelfde reden ook geen
translatie co¨effici¨enten in de x en y richting aan de ellipso¨ıde toe te kennen, omdat een opzij
geschoven ellipso¨ıde op exact dezelfde
 manier door het vlak z = h doorsneden

  wordt.
x1


We nemen dus de ellipso¨ıde E = (x1 , x2 , x3 ) ∈ R3 : x1 x2 x3 A x2  = c


x3
En doorsnijden deze met het vlak (x1 , x2 , x3 ) ∈ R3 : x3 = h .
Merk op dat het vlak de ellipso¨ıde alleen snijdt wanneer |h| < 2ac33 , anders loopt het vlak boven
of onder de ellipso¨ıde langs. Wanneer we ervan uitgaan dat dit het geval is, dan volgt nu
x1 x2

 
a
a
a
x1
11
12
13



h a12 a22 a23
x2  = c
a13 a23 a33
h
Als we dit uitwerken dan krijgen we
a11 x21 + a12 x1 x2 + a13 hx1 + a12 x1 x2 + a22 x22 + a23 hx2 + a13 hx1 + a23 hx2 + a33 h2 = c
⇔ a11 x21 + 2a12 x1 x2 + 2a13 hx1 + 2a23 hx2 + a22 x22 = c − a33 h2
0 x1 − u1
Dit is te schijven als x1 − u1 x2 − u2 A
= c0 ,
x2 − u2
23
En dit voldoet aan de algemene vorm van de vergelijking voor een ellips, waarbij
a12 a23 h − a13 a22 h
,
a11 a22 − a212
a12 a13 h − a11 a23 h
u2 =
,
a11 a22 − a212
a11 a223 h2 − 2a13 a12 a23 h2 + a213 a22 h2
,
c0 = c − a33 h2 +
a11 a22 − a212
a11 a12
0
A =
a12 a22
u1 =
Voor de gedetailleerde berekeningen, zie appendix A
24
Hoofdstuk 4
Discussie
We hebben kunnen concluderen dat Wicksell in zijn tijd een heel aardige schatting heeft gemaakt
van de verdeling van de diameters van de bollen. Omdat de monotone schatters in deze scriptie
niet aan bod komen, kunnen we echter geen heel sterke uitspraken doen over de kwaliteit van
de gevonden schatters. Een goed vervolg zou zijn om deze monotone schatters wel te berekenen
en toe te passen om zo betrouwbaardere resultaten te kunnen presenteren.
Verder hebben we gezien dat de doorsnijding van een ellipso¨ıde altijd een ellips oplevert, maar
hebben we geen methode kunnen ontwikkelen voor het vinden van een verdelingsfunctie van de
ellipso¨ıden wanneer alleen de verdeling van de ellipsvormige profielen bekend is.
Ik had mij meer willen verdiepen in de situatie waarin de lichaampjes ellipso¨ıden zijn, vooral
ook omdat bij TATA Steel, waarmee de statistiekgroep gezamelijke onderzoeksprojecten heeft,
de indruk bestaat dat de elliptische variant van het model van Wicksell toegepast zou kunnen
worden op structuren in staal. Ik had graag meer onderzocht over schatters voor de ellipso¨ıden
en of deze wellicht toepasbaar waren op data van TATA Steel.
Een andere interessante uitbreiding zou het zogenaamde ”Globular cluster problem”zijn. Wicksell schreef in zijn eerste artikel [1] al over de gelijkenis tussen de anatoom achter zijn microscoop en de sterrenkundige achter zijn telescoop. Een sterrenkundige ziet namelijk ook alleen
een projectie van de sterren. Wicksell beweerde zelfs dat de gevonden relatie tussen de zichtbare cirkeldiameters en de boldiameters hetzelfde is als de relatie tussen de echte dichtheid van
sterren en de waargenomen dichtheid. Het had leuk geweest om dit nog verder uit te zoeken.
Maar helaas, de tijd was eindig.
25
Bijlage A
Berekening doorsnijding ellipso¨ıde
met vlak


Doorsnijden van de ellipso¨ıde E = (x1 , x2 , x3 ) ∈ R3


a11
3

(x1 , x2 , x3 ) ∈ R : x3 = h levert x1 x2 h a12
a13
Dit uitwerken levert
: x1
a12
a22
a23

 
x1

x2 x3 A x2  = c met het vlak

x
  3
a13
x1


a23
x2  = c
a33
h
a11 x21 + a12 x1 x2 + a13 hx1 + a12 x1 x2 + a22 x22 + a23 hx2 + a13 hx1 + a23 hx2 + a33 h2 = c
(A.1)
⇔ a11 x21 + 2a12 x1 x2 + 2a13 hx1 + 2a23 hx2 + a22 x22 = c − a33 h2
(A.2)
x1 − u1
Als dit te schrijven is in de vorm x1 − u1 x2 − u2 B
= d, waarbij B symmetrisch
x2 − u2
en positief definiet is, dan is de figuur op het doorsnijdingsvlak een ellips.
Uitwerken van deze laatste uitdrukking levert:
b11 x21 − b11 u1 x1 + b12 x1 x2 − b12 u2 x1 − b11 u1 x1 + b11 u21 − b12 x2 u1 + b12 u1 u2 + b12 x1 x2
− b12 u1 x2 + b22 x22 − b22 u2 x2 − u2 b12 x1 + b12 u1 u2 − b22 x2 u2 + b22 u22 = d
⇔ b11 x21 − 2(b11 u1 + b12 u2 )x1 + 2b12 x1 x2 − 2(b12 u1 + b22 u2 )x2 + b22 x22 = d − (b11 u21 + 2b12 u1 u2 + b22 u22 )
(A.3)
Dus als onze doorsnijding in deze vorm te schrijven is, dan moet gelden
a11 = b11 , 2a12 = 2b12 en a22 = b22, dus
a11 a12
B=
a12 a22
Ook volgt er uit de vergelijkingen (A.2) en (A.3 dat
2a13 h = −2(b11 u1 + b12 u2 )
(A.4)
2a23 h = −2(b12 u1 + b22 u2 )
(A.5)
26
Uit (A.4) volgt u2 =
−b11 u1 −a13 h
b12
en na substitutie in (A.5) volgt
−b11 u1 − a13 h
= 2a23 h
b12
2b11 b22 − 2b212
2a23 hb12 − 2a13 hb22
u1 =
b12
b12
a23 hb12 − a13 hb22
u1 =
b11 b22 − b212
−2b12 u1 − 2b22
Substitueren we dit resultaat nu in de uitdrukking voor u2 dan volgt
−b11 a23 hb12 + b12 a13 hb22 a13 h
−
b12
b12 (b11 b22 − b212 )
a13 hb12 − a23 hb11
=
b11 b22 − b212
u2 =
Ten slotte moet er gelden dat c − a33 h2 = d − (b11 u21 + 2b12 u1 u2 + b22 u22 )
Na substitutie van de zojuist gevonden waarden voor u1 en u2 en uitwerking volgt:
b11 u21 + 2b12 u1 u2 + b22 u22 =
a223 h2 b11 − 2a13 a23 h2 b12 + a213 h2 b22
b11 b22 − b212
a2 h2 b −2a a h2 b +a2 h2 b
Dus volgt nu d = c − a33 h2 + 23 11 b 13b 23−b2 12 13 22
11 22
12
En hiermee zijn alle in paragraaf (3.1) genoemde resultaten aangetoond.
27
Bijlage B
R code
B.1
Verdelingsfunctie Wicksell (Figuur 2.7)
%construeren verdelingsfunctie van wicksell
wr<-c(5,14,38,72,151,194,172,143,96,57,31,15,6,3,2,1,0)%wicksells resultaten
grenzen<-c(0,0.5,1.5,2.5,3.5,4.5,5.5,6.5,7.5,8.5,9.5,10.5,11.5,12.5,13.5,14.5,15.5)
m=sum(wr)
prob<-numeric(17)
for (i in 1:16)
prob[i]<-(1/m*wr[i])/(grenzen[i+1]-grenzen[i])
som<-numeric(16)
for (i in 1:16){
for (j in 1:i){
som[i]<-som[i]+prob[j]}
}
%plotten verdelingsfunctie
plot(som,type=’l’,xlab="boldiameters in mm", ylab="")
B.2
Schatten exponenti¨
ele verdeling (Sectie 2.5.1)
%Rejection sampling
sample.x = runif(100000,0,1)
accept = c()
for(i in 1:length(sample.x)){
U = runif(1, 0, 5)
if(dunif(sample.x[i], 0, 1)*U <= 2*exp(-2*sample.x[i])) {
accept[i] = ’Yes’
}
else if(dunif(sample.x[i],0,1)*U > 2*exp(-2*sample.x[i])) {
accept[i] = ’No’
}
}
T = data.frame(sample.x, accept = factor(accept, levels= c(’Yes’,’No’)))
D = T[,1][T$accept==’Yes’]
28
%Schatter construeren
n=100
k=101
est<-numeric(k)
F<-numeric(k)
x<-seq(0,1,by=0.01)
z<-sample(D, size=100)
for (i in 1:k)
for (j in 1:n)
if (x[i]<z[j])
est[i]<-est[i]+(z[j]-x[i])^(-1/2)/n
for (l in 1:k)
F[l]<- 1-est[l]/est[1]
%Geschatte verdeling plotten
plot(x,F,type=’l’)
lines(x,1-exp(-2*x), col=2)
B.3
Rejection sampling (Sectie 2.5.2)
%Rejection sampling
sample.x = runif(100000,0,1)
accept = c()
for(i in 1:length(sample.x)){
U = runif(1, 0, 1)
if(dunif(sample.x[i], 0, 1)*3/2*U <= 3/2*sqrt(1-sample.x[i])) {
accept[i] = ’Yes’
}
else if(dunif(sample.x[i],0,1)*3/2*U > 3/2*sqrt(1-sample.x[i])) {
accept[i] = ’No’
}
}
%Histogram van de dichtheid
T = data.frame(sample.x, accept = factor(accept, levels= c(’Yes’,’No’)))
hist(T[,1][T$accept==’Yes’], breaks = seq(0,1,0.01), freq = FALSE, main = ’Histogram of X’
lines(x, 3/2*sqrt(1-x))
D = T[,1][T$accept==’Yes’]
sample(D, size=10, replace=TRUE)
%Schatter construeren
n=100
k=101
est<-numeric(k)
29
F<-numeric(k)
x<-seq(0,1,by=0.01)
z<-sample(D, size=100)
for (i in 1:k)
for (j in 1:n)
if (x[i]<z[j])
est[i]<-est[i]+(z[j]-x[i])^(-1/2)/n
for (l in 1:k)
F[l]<- 1-est[l]/est[1]
%Geschatte verdeling plotten
plot(x,F,type=’l’)
B.4
Histogramschatter (Sectie 2.6.2)
%construeren frequentiehistogram en histogramdichtheid
eg<- seq(0.5,15.5,by=1) %echte grenzen
ng<-(eg/2)^2 %nieuwe grenzen
data<-c(0,52,146,197,210,184,143,95,57,31,15,7,4,2,1,0)
n<-sum(data)
g<-numeric(16)
for (i in 1:15)
g[i]<-(1/n*data[i])/(ng[i+1]-ng[i])
%plotten frequentie histogram
plot(ng,data,type=’s’,xlab="gekwadrateerde cirkelstralen in mm", ylab="frequentie")
lines(ng,data,type=’h’)
%plotten histogramdichtheid
plot(ng,g,type=’s’,xlab="gekwadrateerde cirkelstralen in mm",ylab="")
lines(ng,g,type=’h’)
%Construeren histogramschatter
eg<- seq(0.5,15.5,by=1) % echte grenzen
ng<-(eg/2)^2 %nieuwe grenzen
n=1144
data<-c(0,52,146,197,210,184,143,95,57,31,15,7,4,2,1,0)
g<-numeric(16)
for (i in 1:15)
g[i]<-(1/n*data[i])/(ng[i+1]-ng[i])
%construeren van Vn
vnaive<-function(z){
%bepalen in welk interval z zit
bel<-numeric(15)
som=0
for (i in 1:15){
30
if (ng[i]<z){
bel[i]=1}}
for (i in 1:15){
if (bel[i]==1){
som=som+1}}
%bouwen van Vn
v<-g[som]*sqrt(ng[som+1]-z)
for (j in (som+1):15) {
v<-v+g[j]*(sqrt(ng[j+1]-z)-sqrt(ng[j]-z))
}
return(v)
}
%construeren van schatter voor F
F<-function(z){
noemer=0
for (i in 1:15){
noemer<-noemer+g[i]*(sqrt(ng[i+1])-sqrt(ng[i]))}
f<-1-vnaive(z)/noemer
return(f)
}
%plotten van verdelingsfunctie van F
grafiek<-numeric(52)
for (i in 1:52)
grafiek[i]<-F(i)
plot(grafiek, type=’l’,xlab="gekwadrateerde bolstralen in mm", ylab="")
31
Bibliografie
[1] Sven Dag Wicksell. The Corpuscle Problem: A Mathematical Study of a Biometric Problem,
Biometrika Trust, Vol. 17, No. 1/2, 1925.
[2] Sven Dag Wicksell. The Corpuscle Problem: Second Memoir: Case of Ellipsoidal Corpuscles,
Biometrika Trust, Vol. 18 No 1/2, 1926.
[3] P. Groeneboom en G. Jongbloed. Nonparametric Estimation under Shape Constraints: Estimators, Algorithms and Asymptotics., Cambridge Series in Statistical and Probabilistic Mathematics, eind 2014
[4] J.B. Fraleigh en R.A. Beauregard. Linear Algebra, University of Rhode Island, 3rd edition,
1995.
32

Download Report