Bonusopgaven Stochastische OR 3.2 2007/2008 1 December 2007 Cliff Voetelink (1554506) TR3
1.1: Het Metropolis‐Hastings Algoritme Probleembeschrijving Gegeven de kansdichtheid: f n , a ( x) = x
a 1 − n − 2x 2
e
x >0
(Eq. 1.1.1)
(Eq. 1.1.2)
1 2π . 8
(Eq. 1.1.3)
Ook kan er worden aangetoond dat E ( X ) = 4 en Var ( X ) = ∞ .
(Eq. 1.1.4)
Neem n = 5 en a = 4. −
2 5 − x 2
Dan geldt: π ( x) := f5,4 ( x) = x e
x > 0 ∞
∫
In dit geval kan er worden aangetoond dat: π ( x) dx = 0
De kansdichtheid is dus bekend op de constante 1/( 18 2π ) na, deze constante is echter niet nodig voor het algoritme, wat overigens een groot voordeel is als de constante niet analytisch te bepalen is. Situatie 1: We nemen q(t | s) ~ U (0,100) . Situatie 2: We nemen q (t | s ) ~ χ12 Situatie 3: We nemen q (t | s ) ~
p⎛b⎞ ⎜ ⎟ b⎝t⎠
p +1
t > b p = 1.5, b = 1
N.B. Situatie 3 heeft betrekking op een Paretoverdeling. N.B. Alle dichtheden zijn onafhankelijk van s . We gaan onderzoeken hoe goed de mixing van de Markovketen verloopt bij verschillende q(t | s) . De mixing van de Markovketen is van belang omdat de Markov‐keten niet te lang in de buurt van dezelfde toestand mag blijven, alle mogelijke toestanden moeten voldoende vaak bezocht worden. Het idee is dat we dus 500 toestanden genereren van de Markovketen { X n , n = 1,.., 500} . Daarna plotten we deze waarnemingen en bekijken of dit inderdaad het geval is. Onder een goede mixing verstaan we dat alle toestanden vaak genoeg bezocht worden en dat men niet te lang rond een bepaalde toestand blijft hangen. Indien de mixing goed verloopt dan hebben we een betrouwbare schatting voor de unieke evenwichtsverdeling (dichtheid) Eq. 1.1.2. Het voordeel hiervan is dat men op simpele wijze random trekkingen doen uit Eq.1.1.2 en daarom op simpele wijze bijvoorbeeld momenten kan schatten. Belangrijk is dat de kansdichtheid q( x) de staarten van de uit te simuleren kansdichtheid π ( x) domineert. We zullen de verwachting en varianties schatten met de volgende schatters:
1 m Eˆ ( X ) = ∑ X k m k =1
m ˆ ( X ) = ⎛ 1 ∑ X 2 ⎞ − Eˆ 2 ( X ) Var k ⎟ ⎜m ⎝ k =1 ⎠
(Eq. 1.1.5)
(Eq. 1.1.6)
Resultaten Fig. 1.1A: Een plot van de dichtheid π ( x) (Eq. 1.1.2) Theoretische Resultaten Situatie 1: q(t | s) ~ U (0,100)
q(t | s) heeft alleen waarden > 0 voor x ∈ (0,100] en zal daarom staarten van π ( x) nooit domineren.
Met andere woorden: de kansdichtheid is alleen groter dan nul op een bepaald interval. Dit leidt tot het feit dat je nooit een random trekking X zal doen die een waarde groter dan 100 heeft, dus dat je gewoon een deel van de kansdichtheid π ( x) verwaarloost (niet alle toestanden zullen bezocht worden) en dus nooit random trekking
X ∈ (100, ∞) zult krijgen terwijl π ( x) daar wel daadwerkelijk kansmassa heeft. Dit leidt weer tot onbetrouwbare schattingen voor de kansdichtheid π ( x) en tot onbetrouwbare schattingen voor de momenten van π ( x) . De mixing van de Markovketen { X n } verloopt dus beroerd! Theoretische Resultaten Situatie 2: q (t | s ) ~ χ12
2 In tegenstelling tot de homogene verdeling beslaat het domein van de χ1 kansdichtheid wel het gehele domein
van π ( x) . De bedoeling van de schrijver en deze opdracht was dan ook om aan te tonen dat we hier wel betrouwbare schattingen voor de verdeling π ( x) konden bepalen. Helaas blijkt dat niet zo te zijn. Het is namelijk zo dat de χ1 kansdichtheid veel dunnere staarten heeft dan de kansdichtheid π ( x) . Hierdoor wordt de 2
rechterstaart van π ( x) alsnog niet gedomineerd door de χ1 kansdichtheid. Waarom is dit een ramp? Vanaf 2
een bepaalde x geldt dat de χ1 kansdichtheid geëvalueerd in x kleiner is dan π ( x) . Dit zal betekenen dat sommige toestanden niet vaak genoeg bezocht worden, namelijk die toestanden waarvoor geldt dat 2
f χ 2 ( x ) < π ( x ) ! Met andere woorden: dit zal leiden tot onbetrouwbare schattingen voor π ( x) en de 1
momenten van π ( x) .
Theoretische Resultaten Situatie 3: q (t | s ) ~
3 ⎛1⎞ ⎜ ⎟ 2⎝t ⎠
5/ 2
t >1
We gebruiken nu een q (t | s ) die hetzelfde staartgedrag vertoont als π ( x) en deze ook domineert wanneer deze kansdichtheid π ( x) op de correcte manier verschoven is. Daarom zullen we alle toestanden wel in voldoende vaak bezoeken. Het staartgedrag is van de orde t −5 / 2 . Hoe een random trekking te doen uit een pareto verdeling? Dit kan via de inversiemethode. Los op p
⎛b⎞ F ( x) = 1 − ⎜ ⎟ ⎝x⎠
Los op: F ( x ) = u met u een random trekking uit (0,1)
We vinden op die manier dat: t =
b u
1/ p
dus t =
1 u2/3
(Eq. 1.1.7) (Eq. 1.1.8)
We vinden op deze manier random trekkingen (toestanden): Yt = X t + b . Met b = 1 uiteraard. Hier is dus overal één bij opgeteld! Omdat de waardeverzameling van de pareto(b,p) verdeling (b , ∞ ) is, en we trekkingen vinden voor Yt in plaats van X t , zullen we ook π ( x) in een ander punt moeten evalueren, nl: π ( yt − 1) = π ( xt ) . We vinden zo het kanshistogram voor de stochast Yt = X t + 1 . Dus we kunnen de verwachting en variantie van de stochast X vinden via:
E( X ) = E(Y ) −1 = E( X + 1) −1
2 2 2 2 Var ( X ) = E (Y ) − 2E (Y − 1) − 1 − E (Y ) = E (( X + 1) ) − 2E ( X ) − 1 − E ( X + 1)
(Eq. 1.1.9) (Eq. 1.1.10)
Het schatten van de momenten gaat weer met behulp van (Eq. 1.1.5) en (Eq. 1.1.6) toegepast op (Eq. 1.1.9) en (Eq. 1.1.10).
Praktische Resultaten Situatie 1: q(t | s) ~ U (0,100)
Procesverloop Xt, q(t|s)~Unif(0,100)
Geschaald Histogram Xt, q(t|s)~Unif(0,100)
100
0.7
90 0.6
80 0.5
Relatieve Frequentie)
60
50
40
30
0.4
0.3
0.2
20 0.1 10
0
0
100
200
300 Xt
400
500
0
600
0
20
40
60
80
100
x
Fig. 1.1B: Het procesverloop van Xt bij 500 trekkingen
Fig. 1.1C: Het histogram van Xt bij 500 trekkingen
We simuleren de verwachting en de variantie een aantal maal. We zullen zien dat de schattingen sterk variëren.
Eˆ ( X ) = 4.29
Eˆ ( X ) = 4.38
Eˆ ( X ) = 3.25
Eˆ ( X ) = 5.93
ˆ ( X ) = 77.40 Var
ˆ ( X ) = 44.5 Var
ˆ ( X ) = 62.59 Var
ˆ ( X ) = 110.15 Var
De histogrammen bij een steekproeven van 500 variëren ook behoorlijk veel. Al met al is te concluderen dat bij een steekproef van 500 random trekkingen de schattingen niet betrouwbaar zijn, dit geldt zowel voor de schattingen voor de kansverdeling als voor de momenten. Laten we nu eens 50.000 trekkingen doen. Geschaald Histogram Xt, q(t|s)~Unif(0,100) 0.8 0.7
ˆ ( X ) = 40.78 Var
Eˆ ( X ) = 3.53
0.6
Relatieve Frequentie)
Toestand // pi(Xt)
70
Eˆ ( X ) = 3.23
Eˆ ( X ) = 3.46
Eˆ ( X ) = 3.33
ˆ ( X ) = 34.53 Var
ˆ ( X ) = 35.78 Var
0.5
0.4
0.3
ˆ ( X ) = 37.61 Var
0.2
0.1
0
0
20
40
60
80
100
x
Fig. 1.1D: Het histogram van Xt bij 50.000 trekkingen
De histogrammen blijven nu vrijwel hetzelfde als we meerdere steekproeven simuleren. Ook de verwachtingen en de varianties blijken nu een stuk constanter. Het zijn alsnog echter zeer onbetrouwbare schattingen omdat men alleen toestanden uit een interval [0,100] toelaat. De mixing van de Markovketen verloopt dus slecht. Het geeft ook een vertekend beeld: we zien namelijk dat de geschatte verwachtingen een stuk lager liggen dan de theoretische verwachting (Zie Eq.1.1.4), idem voor de variantie, die nu indruk lijkt te wekken ongeveer constant te zijn. Praktische Resultaten Situatie 2: q (t | s ) ~ χ12
Procesverloop Xt, q(t|s)~ChiSq(1) 8
7
6
5 Toestand
4
3
2
1
0
Fig. 1.1E: De χ1 kansdichtheid 2
0
100
200
300 Xt
400
500
600
Fig. 1.1E: Het procesverloop van Xt bij 500 trekkingen
We merken op dat het bereik van de toestanden heel erg beperkt is: (0,8] lijkt het wel. Dit is ook niet zo gek omdat de χ1 kansdichtheid enorm veel kansmassa daar heeft. Bijna nooit wordt er een trekking voor X t met 2
een zeer grote waarde getrokken. Dit was bij de homogene verdeling wel anders. De χ1 kansdichtheid gaat dan 2
ook sneller naar nul dan π ( x) , zie ook Fig. 1.1A. Hier geldt dat alle toestanden wel kunnen worden bezocht maar dat niet alle toestanden voldoende vaak worden bezocht. De mixing van de Markovketen is dus nog steeds niet geweldig. In Fig. 1.1F zien het histogram van een steekproef van 500. Ook hier is het weer zo dat histogram significant wijzigt bij een andere steekproef omdat we maar 500 trekkingen doen. Ook schatten we de verwachting en de variantie een aantal keer.
Geschaald Histogram Xt, q(t|s)~ChiSq(1) 0.16
Eˆ ( X ) = 3.12
Relatieve Frequentie
0.14
Eˆ ( X ) = 2.82
0.12
Eˆ ( X ) = 1.97
0.1
Eˆ ( X ) = 2.20
ˆ ( X ) = 2.27 Var ˆ ( X ) = 4.23 Var
We merken op dat zowel de geschatte varianties als de geschatte verwachtingen veel variëren. De verwachting wordt ook hier weer te laag geschat en de geschatte varianties zijn laag omdat men nooit een trekking heeft met een hoge waarde. Dit geeft een enorm vertekend beeld!
0.08
0.06
0.04
0.02
0
ˆ ( X ) = 2.89 Var ˆ ( X ) = 8.01 Var
0
1
2
3
4 x
5
6
7
8
Fig. 1.1F: Het histogram van Xt bij 500 trekkingen
Al met al is te concluderen dat bij een steekproef van 500 random trekkingen de schattingen niet betrouwbaar zijn, dit geldt zowel voor de schattingen voor de kansverdeling als voor de momenten. Laten we nu eens 50.000 trekkingen doen.
Geschaald Histogram Xt, q(t|s)~ChiSq(1) 0.35
Eˆ ( X ) = 2.73 0.3
Eˆ ( X ) = 2.83
0.25
Eˆ ( X ) = 2.52
Eˆ ( X ) = 3.28
Relatieve Frequentie
0.2
ˆ ( X ) = 8.41 Var ˆ ( X ) = 10.48 Var ˆ ( X ) = 6.30 Var ˆ ( X ) = 18.12 Var
Merk op dat het bereik van toestanden wel iets groter is
0.15
geworden. Dat is ook niet zo vreemd aangezien we veel 0.1
meer trekkingen doen. De histogrammen blijven nu ongeveer hetzelfde. De geschatte varianties en
0.05
verwachtingen zijn nog steeds niet constant en zijn 0
0
5
10 x
15
20
Fig. 1.1G: Het histogram van Xt bij 50.000 trekkingen
.
bovendien ook onbetrouwbaar en komen niet overeen met de theoretische waarden (zie Eq. 1.1.4).
3 ⎛1⎞
5/ 2
Geschaald Histogram Xt+1, q(t|s)~Pareto(b,p) 0.35
0.3
Relatieve Frequentie)
0.25
0.2
0.15
0.1
0.05
0
Procesverloop Xt+1, q(t|s)~Pareto(b,p)
1000 900 800
700
600
500
400 300 200 100
0
0
2
4
6 Xt+1
8
10
0
10
20
30 Xt+1
40
50
Fig. 1.1H: Histogram / Geschatte kansdichtheid voor de stochast X+1 bij n= 100.000 simulaties
Fig. 1.1G: De paretodichtheid verschoven naar 1. Parameters b = 1, p =1.5.
Toestand
t >1
Praktische Resultaten Situatie 3: q (t | s ) ~ ⎜ ⎟ 2⎝t ⎠
12 4
x 10
Fig. 1.1I: Het procesverloop van Xt+1 bij 100.000 trekkingen
60
Het moge duidelijk zijn dat er enorm veel simulaties nodig zullen zijn om een betrouwbare schatting te kunnen doen voor de verwachting en de variantie (vanwege de dikke staart van π ( x ) ) en dus grote variatie in de toestanden zoals in Fig. 1.1H te zien is, in dit geval worden alle toestanden wel vaak genoeg bezocht en verloopt de mixing van de Markovketen goed. Om één miljoen keer te simuleren maken we twee hulpfuncties, eentje die m keer 50.000 waarnemingen genereert en vervolgens de geschatte totale verwachting en de totale geschatte variantie van de steekproef ter lengte 50.000*m retourneert en een functie die alleen de geschatte verwachting en geschatte variantie bepaalt van zo’n steekproef van 50.000 waarnemingen. We gebruiken m = 20.
Eˆ100.000 ( X ) = 3.903 Eˆ100.000 ( X ) = 3.801 Eˆ100.000 ( X ) = 4.043 Eˆ100.000 ( X ) = 3.981
ˆ Var 100.000 ( X ) = 1342 ˆ Var 100.000 ( X ) = 216 ˆ Var 100.000 ( X ) = 590
ˆ Var 100.000 ( X ) = 1119
Eˆ1.000.000 ( X ) = 3.91 Eˆ1.000.000 ( X ) = 3.92 Eˆ1.000.000 ( X ) = 3.93 Eˆ1.000.000 ( X ) = 3.92
ˆ Var 1.000.000 ( X ) = 539 ˆ Var 1.000.000 ( X ) = 495 ˆ Var 1.000.000 ( X ) = 814 ˆ Var 1.000.000 ( X ) = 451
Conclusies Zowel met de χ1 kansdichtheid als met de U (0,100) dichtheid kunnen er geen betrouwbare schattingen 2
worden gedaan wat betreft de kansverdeling en de momenten van π ( x) . Men zal op zoek moeten naar een dichtheid die de staarten van π ( x) wel volledig domineert anders zal het Metropolis‐Hastings algoritme niet goed werken en houden we slechts onbetrouwbare schattingen over. Deze kansdichtheid is nu gevonden en blijkt een Pareto( b = 1, p = 1.5 ) dichtheid te zijn. Bij 100.000 waarnemingen blijkt al een goede benadering van de verwachting plaats te vinden, we merken op dat de varianties hoog zijn (wat in verhouding staat tot de theorie). Als we echter een miljoen keer simuleren dan zien we dat de schattingen veel constanter zijn voor de verwachting. Het zijn goede benaderingen voor de theoretische verwachting: 4.
1.2 De Gibbs‐Sampler In deze paragraaf gaan we kijken naar een toepassing van de zogenaamde Gibbs‐sampler. De Gibbs sampler is eigenlijk niks anders dan een speciaal geval van het Metropolis‐Hastings‐algoritme, hier is de acceptatiekans om naar een bepaalde toestand te gaan namelijk gelijk aan 1. We kunnen het algoritme toepassen om trekkingen uit een multivariate kansdichtheid te doen waarvan de univariate conditionele dichtheden bekend zijn. In deze paragraaf passen we een variant op het algoritme toe, namelijk de variant waar iteratief elk van de componenten van de toestand wordt aangepast.
Probleembeschrijving We bekijken het volgende actuariele model:
X = # claims dat resulteert van de polissen in de portefeuille Y = kans dat een gegeven polis tot een claim leidt N = # polissen in een portefeuille
De stochastische vector ( X , Y , N ) heeft simultane kansdichtheid: n ⎛ n ⎞ x +α −1 n − x + β −1 − λ λ (1 − y) e π ( x, y, n) = ⎜ ⎟ y n! ⎝ x⎠ x ∈{0,1, 2,.., n} 0 ≤ y ≤ 1 n∈{0,1, 2,3,....}
( Eq. 1.2.1)
We gaan er vanuit dat elke polis dezelfde claimkans heeft en dat de polissen zich onafhankelijk van elkaar gedragen. π ( x, y, n) is een mengsel van discrete en continue stochastische variabelen. Er moet dus gelden dat: 1 ∞
n
∫ ∑∑ π ( x, y, n)dy = 1
( Eq. 1.2.2)
0 n =0 x =0
Idee: door middel van univariate conditionele kansdichtheden random trekkingen te doen genereert de Gibbs sampler een rij opeenvolgende toestanden x
(i )
= ( xi , yi , ni ) uit een Markovketen die π ( x, y, n) als
evenwichtsverdeling heeft. De conditionele kansdichtheden moeten echter expliciet kunnen worden bepaald, dat is hier het geval. Deze worden gegeven door:
π 1 ( x | y , n) ~ bin(n, y )
⎛n⎞ ⎝ x⎠
Æ π1 ( x | y, n) = ⎜ ⎟ y (1 − y) x
n− x
( Eq. 1.2.3)
1 x x +α (1 − x)n − x + β −1 ( Eq. 1.2.4a) Β( x + α , n − x + β ) met 0 < y < 1 en
Β( x + α , n − x + β ) = y x +α −1 (1 − y ) n − x + β −1 dy
π 2 ( y | x, n) ~ beta ( x + α , n − x + β ) Æ π 2 ( y | x, n) =
1
∫ 0
Dit is de normalisatieconstante.
( Eq. 1.2.4b)
π 3 ( n | x, y ) ~ Poisson (λ (1 − y )) verschoven naar n
Æ π 3 (n | x, y ) = e
− λ (1− y )
[λ (1 − y )]n − x n = x, x + 1,.... ( Eq. 1.2.5) (n − x)!
Men gaat random trekkingen uit deze conditionele verdelingen doen. Een random trekking voor N wordt gedaan door middel van een random trekking uit de Poisson( λ (1 − y) ) verdeling en vervolgens hier x bij op te tellen; x en y zijn gegeven.
We zijn uiteindelijk geïnteresseerd in de marginale dichtheid: π 1 ( x ) , die een goed beeld geeft over de kansverdeling van het uiteindelijke aantal claims in een willekeurige portefeuille. Deze is niet expliciet te bepalen, daarom gaan we aan de slag met de Gibbs‐sampler. Ook zijn we geinteresseerd in de verwachting en variantie van X .
Resultaten We nemen α = 2 , β = 8 en λ = 12 nu vast. Via de variant van het Gibbs algoritme kunnen we de verdeling van X op een naieve manier schatten: we plotten het histogram van de gevonden x1( k ) ' s . Echter passen we eerst inbrandingsperiode toe, dat wil (k )
zeggen: we laten de eerste 100 x ’s (toestanden) weg om de effecten van de begintoestand te verminderen. Vervolgens passen we nog een uitdunningsperiode toe door alleen elke 5de toestand op te slaan. We vinden op deze manier: Geschaald Histogram X Naief
Geschaald Histogram X Naief
0.25
0.45
0.4
0.2
0.35
0.3 Relatieve Frequentie)
Relatieve Frequentie)
0.15
0.1
0.25
0.2
0.15
0.1
0.05
0.05
0
0
2
4
6
8
10
12
14
x
Fig. 1.2.A: Het histogram van X via de naieve schatter n = 5000 simulaties Æ Uiteindelijk 980 waarnemingen gebruikt
0
0
5
10
15
x
Fig. 1.2.B: Het histogram van X via de naieve schatter n = 50000 simulaties Æ Uiteindelijk 9800 waarnemingen gebruikt
Geschaald Histogram X Naief
Geschaald Histogram X Naief
0.25
0.25
0.2
0.2
Relatieve Frequentie)
Relatieve Frequentie)
0.15
0.1
0.15
0.1
0.05
0.05
0 -5
0
5
10
15
0 -5
20
0
x
ˆ (X ) = ⎛ 1 Var ⎜
k =1
(k ) 2
15
20
Fig. 1.2.D: Het histogram van X via de naieve schatter n = 50000 simulaties Æ Uiteindelijk 9800 waarnemingen gebruikt Æ Staafbreedte ter lengte 1
( Eq. 1.2.6)
Æ Eˆ5000 ( X ) = 2.39
2 ⎞ ˆ ⎟ − E ( X ) ( Eq. 1.2.7) ⎠
ˆ ˆ Æ Var 5000 ( X ) = 4.47 Var50.000 ( X ) = 4.60
∑( x ) ⎝m m
10 x
Fig. 1.2.C: Het histogram van X via de naieve schatter n = 5000 simulaties Æ Uiteindelijk 980 waarnemingen gebruiktÆ Staafbreedte ter lengte 1
1 m Eˆ ( X ) = ∑ x ( k ) m k =1
5
(
)
Eˆ50.000 ( X ) = 2.44
We kunnen de dichtheid beter schatten door gebruik te maken van de univariate conditionele kansdichtheid van X .
πˆ1 ( x) =
1 m 1 m ⎛ n( k ) ⎞ ( k ) x (k ) (k ) ( k ) n( k ) − x = π ( | , ) x y n ⎜ ⎟ ( y ) (1 − y ) ∑ ∑ 1 m k =1 m k =1 ⎝ x ⎠
( Eq. 1.2.8)
We maken nu een vector x = (0,1, 2,...., 20) en we passen hierop de geschatte kansverdeling (die volgens Eq. 1.2.8 gevonden is) toe. We vullen één voor één de discrete waarden van x in Eq. 1.2.8 in. Zo verkrijgen we een benadering voor de kansverdeling van X.
Dit leidt tot de volgende grafieken:
Fig. 1.2.E: De geschatte continue kansdichtheid πˆ1 ( x)
Fig. 1.2.F: De geschatte kansmassafunctie πˆ1 ( x) van X
van X via conditioneren. n = 5000 simulaties Æ Uiteindelijk 980 waarnemingen gebruikt
via conditioneren. n = 5000 simulaties Æ Uiteindelijk 980 waarnemingen gebruikt
Als we n = 50.000 simulaties zouden nemen, dan blijkt dat de geschatte dichtheid er nog steeds hetzelfde uit te zien. Daarom is deze hier niet vermeld. We bepalen de verwachting en de variantie nu met de volgende formules: Eˆ ( X ) =
20
∑ kπˆ (k ) k =1
ˆ (X ) = Var
1
20
Æ Eˆ ( X ) = 2.34
Eˆ50.000 ( X ) = 2.41
( Eq. 1.2.10)
ˆ ( X ) = 4.23 Æ Var
ˆ Var 50.000 ( X ) = 4.55
20
∑ k 2πˆ1 (k ) −∑ kπˆ1 (k ) k =1
( Eq. 1.2.9)
k =1
Conclusies Interessant is het feit dat de geschatte kansverdelingen redelijk overeen komen. Hoewel het lijkt alsof Fig. 1.2B significant afwijkt van Fig. 1.2D, zijn het in feite dezelfde histogrammen maar dan anders geschaald op de x‐as. Men kan Fig. 1.2C en Fig 1.2D intepreteren als schattingen voor de kansverdeling van X op de naieve manier, met de hoogte van zo’n staaf gelijk aan een schatting voor P( X = k ) . We merken op dat de vorm van de histogrammen Fig. 1.2C en Fig 1.2D heel erg overeenkomt met de geschatte kansmassafunctie (Fig. 1.2F) die we via het sommeren van de conditionele kansdichtheden hebben gevonden. De vorm van de verdeling lijkt wat weg te hebben van een afgeknotte normale verdeling maar we moeten natuurlijk opmerken dat X eigenlijk een discrete verdeling heeft. Ook is de gamma verdeling een mogelijkheid.
Opvallend is dat de verwachtingen en varianties van X bij beide schattingsmethoden weinig verschillen maar dat de methode met betrekking tot conditioneren, d.w.z. via: πˆ1 ( x ) , toch altijd schatting oplevert met een lagere variantie voor X . We kunnen in het algemeen zeggen dat het aantal claims van polissen in het algemeen behoorlijk laag is, het gemiddelde bedraagt bij benadering 2.4 claims en de geschatte variantie rond de 4.5. De kans op hoogstens 10 claims (d.w.z.: P( X ≤ k ) ) is bij benadering zelfs groter dan 0.9.
Appendix 1.1A function p = probDensP(x) % Calculates de waarde van de density pi(x) op x % Op een multiplicatieve constante na % Die wordt immers toch weggedeeld later p = x^(-1/2*5)*exp(-1/2*4/x); end function mh = methahom(n, a, b, bins) % simuleert n maal Xt --> X0 niet meegeteld, dus eigenlijk n+1 keer % bins = # bins/staven bij histogram % q(t|s) ~ Uniform(a,b) aantalSimulaties = n; Xt = zeros(1,n); Xt(1) = a+(b-a)*rand; %X0 % Simuleren uit een Uniform(a,b) verdeling s(1) = Xt(1) %S0 n = 1; % n wordt nu een teller voor hoeveel simulaties we zitten for k = 2:(aantalSimulaties+1) %Stap 1 t(k) = a+(b-a)*rand; u = rand; %Stap 2 alfa = (probDensP(t(k))*1/(b-a))/(probDensP(s(k-1))*1/(b-a)); if( u < alfa) s(k) = t(k); else s(k) = s(k-1); end %Stap 3 Xt(k) = s(k); n = n+1; % Herhaal Vanaf Stap 1 end s; Xt; % Histogram plotten [y,b] = hist(Xt,bins); Frequentie = y/aantalSimulaties; figure('Position', [200 200 500 500 ]); bar(b,Frequentie,1,'k') title('Geschaald Histogram Xt, q(t|s)~Unif(0,100)') xlabel('x') ylabel('Relatieve Frequentie)')
%X1, ... , Xn plotten om het mixing proces van de MarkovKeten te zien figure('Position', [200 200 500 500 ]); plot(Xt); title('Procesverloop Xt, q(t|s)~Unif(0,100)') xlabel('Xt') ylabel('Toestand // pi(Xt)') % daarnaast theoretische dichtheid x = linspace(0.0001,30,10000); f = x.^(-1/2.*5).*exp(-1/2.*4./x); figure('Position', [200 200 500 500 ]); subplot(1,2,2); plot(x,f,'k'); title('Ware Dichtheid'); xlabel('x'); ylabel('Dichtheid pi(x)'); VerwachtingGeschatNaief = 1/length(Xt)*sum(Xt) VariantieGeschatNaief = (1/length(Xt)*sum(Xt.^2))-VerwachtingGeschatNaief^2 end
Appendix 1.1B function p = probDensP(x) % Calculates de waarde van de density pi(x) op x % Op een multiplicatieve constante na % Die wordt immers toch weggedeeld later p = x^(-1/2*5)*exp(-1/2*4/x); end
function mh = methachi(n, df, bins) % simuleert n maal Xt --> X0 niet meegeteld, dus eigenlijk n+1 keer % bins = # bins/staven bij histogram % q(s,t) ~ Chisquared(df) aantalSimulaties = n; Xt = zeros(1,n); Xt(1) = chi2rnd(df); %X0, %Random trekking uit ChiKwadr(df) verdeling s(1) = Xt(1); %S0 n = 1; % n wordt nu een teller voor hoeveel simulaties we zitten for k = 2:(aantalSimulaties+1) %Stap 1 t(k) = chi2rnd(df); u = rand; %Stap 2 alfa = (probDensP(t(k))*chi2pdf(s(k-1),df))/(probDensP(s(k1))*chi2pdf(t(k),df)); %Let op de argumenten! Chi2pdf(s(k-1)) is niks anders dan q(tk,sk-1) en %Chi2pdf(t(k)) is q(sk-1,t(k)) want q(s,t) is onafhankelijk van s ! if( u < alfa) s(k) = t(k); else s(k) = s(k-1); end %Stap 3 Xt(k) = s(k); n = n+1; % Herhaal Vanaf Stap 1 end s; max(Xt) % Histogram plotten [y,b] = hist(Xt,bins); Frequentie = y/aantalSimulaties; figure('Position', [200 200 500 500 ]);
bar(b,Frequentie,1,'k') title('Geschaald Histogram Xt, q(t|s)~ChiSq(1)') xlabel('x') ylabel('Relatieve Frequentie') %X1, ... , Xn plotten om het mixing proces van de MarkovKeten te zien figure('Position', [300 300 500 500]); plot(Xt); title('Procesverloop Xt, q(t|s)~ChiSq(1)') xlabel('Xt') ylabel('Toestand') %ChiSq(df) verdeling plotten x = linspace(0,30,100); f = chi2pdf(x,df); figure('Position', [200 200 500 500 ]); subplot(1,2,2); plot(x,f,'k'); title('Ware Dichtheid Chisq(df)'); xlabel('x'); ylabel('Dichtheid Chisq(x)'); VerwachtingGeschatNaief = 1/length(Xt)*sum(Xt) VariantieGeschatNaief = (1/length(Xt)*sum(Xt.^2))-VerwachtingGeschatNaief^2 end
Appendix 1.1C function mh = methapareto(n, b, p, bins) % % % %
simuleert n maal Xt --> X0 niet meegeteld, dus eigenlijk n+1 keer bins = # bins/staven bij histogram q(t|s) ~ Pareto (b, p) 1 verschoven naar rechts We gaan simuleren voor b = 1, p = 1.5
aantalSimulaties = n; Yt = zeros(1,n); Yt(1) = b/((rand)^(1/p)); %X0, %Random trekking uit Pareto(b,p) verdeling. N.B. Dit is een random trekking voor Yt = Xt+1 ! s(1) = Yt(1); %S0 n = 1; % n wordt nu een teller voor hoeveel simulaties we zitten for k = 2:(aantalSimulaties+1) %Stap 1 t(k) = b/((rand)^(1/p)); %N.B. Dit is een random trekking voor X+1 ! u = rand; %Stap 2 alfa = (probDensP(t(k)-1)*((p/b)*(b/(s(k-1)))^(p+1)))/(probDensP(s(k-1)1)*((p/b)*(b/(t(k)))^(p+1))); %N.B. We moeten X invullen, niet X+1 daarom bij (Yt-1) invullen in %pi(x) if( u < alfa) s(k) = t(k); else s(k) = s(k-1); end %Stap 3 Yt(k) = s(k); n = n+1; % Herhaal Vanaf Stap 1 end Yt; max(Yt); % Histogram plotten % We krijgen nu een histogram van de naar 1 verschoven dichtheid pi [y,g] = hist(Yt,bins); Frequentie = y/aantalSimulaties; figure('Position', [200 200 500 500 ]); bar(g,Frequentie,1,'k') title('Geschaald Histogram Xt+1, q(t|s)~Pareto(b,p)') xlabel('x') ylabel('Relatieve Frequentie') %X1, ... , Xn plotten om het mixing proces van de MarkovKeten te zien figure('Position', [300 300 500 500]); plot(Yt);
title('Procesverloop Xt+1, q(t|s)~Pareto(b,p)') xlabel('Xt+1') ylabel('Toestand') %%%%%%% Histogram maken bij benadering edges = [1:50]; n = histc(Yt,edges); Frequentie = n/sum(n); figure('Position', [200 200 500 500 ]); bar(edges,Frequentie,1,'k') title('Geschaald Histogram Xt+1, q(t|s)~Pareto(b,p)') xlabel('Xt+1') ylabel('Relatieve Frequentie)')
%Pareto(b=1 ,p = 1.5) verdeling plotten x = linspace(b,30,100); f = ((p./b).*(b./x).^(p+1)) figure('Position', [200 200 500 500 ]); subplot(1,2,2); plot(x,f,'k'); title('Ware Dichtheid Pareto 1 Rechts Verschoven(b=1,p=1.5)'); xlabel('x'); ylabel('f(x)'); %N.B. Dit zou de verwachting van X+1 zijn dus we moeten er nog 1 vanaf %trekken, VerwachtingGeschat = 1/length(Yt)*sum(Yt)-1 VariantieGeschat = (1/length(Yt)*sum(Yt.^2))-2*VerwachtingGeschat-1VerwachtingGeschat^2 end
Appendix 1.1D function [VerwachtingGeschat,VariantieGeschat] = methapareto2(b, p) % % % %
simuleert n maal Xt --> X0 niet meegeteld, dus eigenlijk n+1 keer bins = # bins/staven bij histogram q(t|s) ~ Pareto (b, p) 1 verschoven naar rechts We gaan simuleren voor b = 1, p = 1.5
n = 50000; aantalSimulaties = n; Yt = zeros(1,n); Yt(1) = b/((rand)^(1/p)); %X0, %Random trekking uit Pareto(b,p) verdeling, N.B. Dit is een random trekking voor X+1 ! s(1) = Yt(1); %S0 n = 1; % n wordt nu een teller voor hoeveel simulaties we zitten for k = 2:(aantalSimulaties+1) %Stap 1 t(k) = b/((rand)^(1/p)); u = rand; %Stap 2 alfa = (probDensP(t(k)-1)*((p/b)*(b/(s(k-1)))^(p+1)))/(probDensP(s(k-1)1)*((p/b)*(b/(t(k)))^(p+1))); %N.B. We moeten X invullen, niet X+1 daarom Yt-1 invullen in %pi(x) , N.B. Voor de pareto vullen we gewoon Yt in want deze heeft wel %gewoon domein (1, inf) ! Deze is verschoven! if( u < alfa) s(k) = t(k); else s(k) = s(k-1); end %Stap 3 Yt(k) = s(k); n = n+1; % Herhaal Vanaf Stap 1 end
%N.B. Dit zou de verwachting van X+1 zijn dus we moeten er nog 1 vanaf %trekken VerwachtingGeschat = 1/length(Yt)*sum(Yt)-1; VariantieGeschat =(1/length(Yt)*sum(Yt.^2))-2*VerwachtingGeschat-1VerwachtingGeschat^2;
end
Appendix 1.1E function vector = hulpfunctpareto(m,b,p) %Simuleert m keer 50.000 simulaties %m = # keer dat je 50.000 simulaties gaat uitvoeren %Print uiteindelijk de geschatte verwachting en variantie van al deze %steekproeven samen for k = 1:m [verw(k),var(k)] = methapareto2(b,p); end verw; var; verwachting = sum(verw)/m; variantie = sum(var)/m; vector(1) = verwachting; vector(2) = variantie; vector; end
Appendix 1.2A function gibbs = gibbsVariant(alfa, beta, lambda, n, bins) % % % %
EIS: n > 100 Genereert een rij van opeenvolgende toestanden uit een Markovketen die als evenwichtsverdeling pi(x,y,z) heeft. bins = # bins/staven bij histogram
aantalSimulaties = n; aantalStochasten = 3; A = zeros(n,aantalStochasten); %Matrix A: elke rij stelt een toestand voor, elke kolom is de waarde vd stoch variabele % de eerste kolom: X, tweede kolom Y, derde kolom N % Stap A(1,1) A(1,2) A(1,3)
0: Initialisatie = 3; = 0.4; = 7;
for i = 1:n %Random trekking uit univariate conditionele verdeling van x: pi1 v(1) = binornd(A(i,3),A(i,2)); v(2) = betarnd(v(1)+alfa,A(i,3)-v(1)+beta); v(3) = poissrnd(lambda*(1-v(2)))+v(1); A(i+1,1) = v(1); A(i+1,2) = v(2); A(i+1,3) = v(3); end A; %Inbrandingsperiode --> Eerste zoveel waarnemingen weggooien --> zeg 100 m = n; A = A(100:m,1:3); m = m-100; %Uitdunningsprocedure %We gebruiken alleen iedere 5de waarneming m = floor(m/5); for i = 1:m A(i,1) = A(5*i,1); A(i,2) = A(5*i,2); A(i,3) = A(5*i,3); end A = A(1:m,1:3); A m
% Nu de naieve schatting voor pi1(x) bepalen voor de stochast X1 v = A(:,1)
%Histogram plotten van v [y,b] = hist(v,bins); Frequentie = y/m; figure('Position', [200 200 500 500 ]); bar(b,Frequentie,1,'k') title('Geschaald Histogram X Naief') xlabel('x') ylabel('Relatieve Frequentie)') % Naieve verwachting & Variantie VerwachtingGeschatNaief = 1/m*sum(v) VariantieGeschatNaief = 1/m*sum(v.^2)-VerwachtingGeschatNaief^2 % % % % %
Nu met de schatter de expliciete formule pi1(x|y,n) Eerst alle pi1(x|y,N) termen bepalen We maken een vector x met discrete waarden 0.... 20 (in dit geval, we zullen zien dat de kansen als x > 20 vrijwel nul zijn. x is discreet omdat pi1(x|y,n) een discrete verdeling heeft: binomiaal
x = [0:20]; x; % Nu sum pi(x|y (k ), z (k) ) bepalen for j = 1:length(x) f(j) = 0; for k = 1:m f(j) = f(j) + binopdf(x(j),A(k,3),A(k,2)); end end % nu pi1-hoed(x) bepalen f = 1/m.*f; f % We merken op dat de kansen vrijwel nul zijn vanaf x > 20 % Geschatte Kansdichtheid tekenen figure('Position', [0 200 300 300]); subplot(1,2,2); plot(x,f,'k'); title('Schatting Dichtheid pi1hoed(x)'); xlabel('x'); ylabel('pi1hoed(x)'); %f(j) = P(X=j), j = 0, ... , 20 VerwachtingSchatting = x*f' VariantieSchatting = (x.^2)*f'-VerwachtingSchatting^2 end
Appendix 1.2B: Geschaalde Histogrammen function gibbs = gibbsVariant(alfa, beta, lambda, n) % EIS: n > 100 % Genereert een rij van opeenvolgende toestanden uit een Markovketen die % als evenwichtsverdeling pi(x,y,z) heeft. aantalSimulaties = n; aantalStochasten = 3; A = zeros(n,aantalStochasten); %Matrix A: elke rij stelt een toestand voor, elke kolom is de waarde vd stoch variabele % de eerste kolom: X, tweede kolom Y, derde kolom N % Stap A(1,1) A(1,2) A(1,3)
0: Initialisatie = 3; = 0.4; = 7;
for i = 1:n %Random trekking uit univariate conditionele verdeling van x: pi1 v(1) = binornd(A(i,3),A(i,2)); v(2) = betarnd(v(1)+alfa,A(i,3)-v(1)+beta); v(3) = poissrnd(lambda*(1-v(2)))+v(1); A(i+1,1) = v(1); A(i+1,2) = v(2); A(i+1,3) = v(3); end A; %Inbrandingsperiode --> Eerste zoveel waarnemingen weggooien --> zeg 100 m = n; A = A(100:m,1:3); m = m-100; %Uitdunningsprocedure %We gebruiken alleen iedere 5de waarneming m = floor(m/5); for i = 1:m A(i,1) = A(5*i,1); A(i,2) = A(5*i,2); A(i,3) = A(5*i,3); end A = A(1:m,1:3); A m
% Nu de naieve schatting voor pi1(x) bepalen voor de stochast X1 v = A(:,1) %Histogram plotten van v
edges = [0:18]; n = histc(v,edges); Frequentie = n/sum(n); figure('Position', [200 200 500 500 ]); bar(edges,Frequentie,1,'k') title('Geschaald Histogram X Naief') xlabel('x') ylabel('Relatieve Frequentie)') % Naieve verwachting & Variantie VerwachtingGeschatNaief = 1/m*sum(v) VariantieGeschatNaief = 1/m*sum(v.^2)-VerwachtingGeschatNaief^2 % % % % %
Nu met de schatter de expliciete formule pi1(x|y,n) Eerst alle pi1(x|y,N) termen bepalen We maken een vector x met discrete waarden 0.... 20 (in dit geval, we zullen zien dat de kansen als x > 20 vrijwel nul zijn. x is discreet omdat pi1(x|y,n) een discrete verdeling heeft: binomiaal
x = [0:20]; x; % Nu sum pi(x|y (k ), z (k) ) bepalen for j = 1:length(x) f(j) = 0; for k = 1:m f(j) = f(j) + binopdf(x(j),A(k,3),A(k,2)); end end % nu pi1-hoed(x) bepalen f = 1/m.*f; f % We merken op dat de kansen vrijwel nul zijn vanaf x > 20 % Geschatte Kansmassafunctie tekenen figure('Position', [0 200 300 300]); subplot(1,2,2); bar(x,f,1,'k') title('Schatting Kansmassafunctie pi1hoed(x)'); xlabel('x'); ylabel('pi1hoed(x)'); %f(j) = P(X=j), j = 0, ... , 20 VerwachtingSchatting = x*f' VariantieSchatting = (x.^2)*f'-VerwachtingSchatting^2 end
Appendix 1.3A function r=binornd(n,p,mm,nn) % BINORND Random matrices from a binomial distribution. % R = BINORND(N,P,MM,NN) is an MM-by-NN matrix of random % numbers chosen from a binomial distribution with parameters N and P. % % The size of R is the common size of N and P if both are matrices. % If either parameter is a scalar, the size of R is the size of the other % parameter. % Alternatively, R = BINORND(N,P,MM,NN) returns an MM by NN matrix. % % The method is direct using the definition of the binomial % distribution as a sum of Bernoulli random variables. % % % %
References: [1] L. Devroye, "Non-Uniform Random Variate Generation", Springer-Verlag, 1986 See Lemma 4.1 on page 428, method on page 524.
% % %
B.A. Jones 1-12-93 Copyright (c) 1993-98 by The MathWorks, Inc. $Revision: 2.6 $ $Date: 1997/11/29 01:44:57 $
if nargin < 2, error('Requires at least two input arguments.'); end if nargin == 2 [errorcode rows columns] = rndcheck(2,2,n,p); if max(size(p)) == 1 & max(size(n)) > 1 p = p(ones(size(n))); end if max(size(n)) == 1 & max(size(p)) > 1 n = n(ones(size(p))); end end if nargin == 3 [errorcode rows columns] = rndcheck(3,2,n,p,mm); n = n(ones(mm(1),1),ones(mm(2),1)); p = p(ones(mm(1),1),ones(mm(2),1)); end if nargin == 4 [errorcode rows columns] = rndcheck(4,2,n,p,mm,nn); n = n(ones(mm,1),ones(nn,1)); p = p(ones(mm,1),ones(nn,1)); end if errorcode > 0 error('Size information is inconsistent.'); end
r = zeros(rows,columns); for i = 1:max(max(n)) u = rand(rows,columns); k1 = find(n >= i); k2 = find(u(k1) < p(k1)); r(k1(k2)) = r(k1(k2)) + 1; end k= find(p < 0 | p > 1 | n < 0); if any(k) tmp = NaN; r(k) = tmp(ones(size(k))); end
Appendix 1.3B function r = betarnd(a,b,m,n); %BETARND Random matrices from beta distribution. % R = BETARND(A,B) returns a matrix of random numbers chosen % from the beta distribution with parameters A and B. % The size of R is the common size of A and B if both are matrices. % If either parameter is a scalar, the size of R is the size of the other % parameter. Alternatively, R = BETARND(A,B,M,N) returns an M by N matrix. % % %
Reference: [1] L. Devroye, "Non-Uniform Random Variate Generation", Springer-Verlag, 1986
% %
Copyright (c) 1993-98 by The MathWorks, Inc. $Revision: 2.5 $ $Date: 1997/11/29 01:44:53 $
if nargin < 2, error('Requires at least two input arguments'); end if nargin == 2 [errorcode rows columns] = rndcheck(2,2,a,b); end if nargin == 3 [errorcode rows columns] = rndcheck(3,2,a,b,m); end if nargin == 4 [errorcode rows columns] = rndcheck(4,2,a,b,m,n); end if errorcode > 0 error('Size information is inconsistent.'); end r = zeros(rows,columns); % Use Theorem 4.1, case A (Devroye, page 430) to derive beta % random numbers as a ratio of gamma random numbers. if prod(size(a)) == 1 a1 = a(ones(rows,1),ones(columns,1)); g1 = gamrnd(a1,1); else g1 = gamrnd(a,1); end if prod(size(b)) == 1 b1 = b(ones(rows,1),ones(columns,1)); g2 = gamrnd(b1,1); else g2 = gamrnd(b,1); end r = g1 ./ (g1 + g2);
% Return NaN if b is not positive. if any(any(b <= 0)); if prod(size(b) == 1) tmp = NaN; r = tmp(ones(rows,columns)); else k = find(b <= 0); tmp = NaN; r(k) = tmp(ones(size(k))); end end % Return NaN if a is not positive. if any(any(a <= 0)); if prod(size(a) == 1) tmp = NaN; r = tmp(ones(rows,columns)); else k = find(a <= 0); tmp = NaN; r(k) = tmp(ones(size(k))); end end
Appendix 1.3C function r = poissrnd(lambda,m,n) %POISSRND Random matrices from Poisson distribution. % R = POISSRND(LAMBDA) returns a matrix of random numbers chosen % from the Poisson distribution with parameter LAMBDA. % % The size of R is the size of LAMBDA. Alternatively, % R = POISSRND(LAMBDA,M,N) returns an M by N matrix. % % POISSRND uses a waiting time method for small values of LAMBDA, % and the method of Ahrens and Dieter for larger values of LAMBDA. % % %
References: [1] L. Devroye, "Non-Uniform Random Variate Generation", Springer-Verlag, 1986 page 504.
% %
Copyright 1993-2000 The MathWorks, Inc. $Revision: 2.9 $ $Date: 2000/07/28 19:33:02 $
if nargin < 1, error('Requires at least one input argument.'); end if nargin == 1 [errorcode rows columns] = rndcheck(1,1,lambda); end if nargin == 2 [errorcode rows columns] = rndcheck(2,1,lambda,m); end if nargin == 3 [errorcode rows columns] = rndcheck(3,1,lambda,m,n); end if errorcode > 0 error('Size information is inconsistent.'); end if (prod(size(lambda)) == 1) lambda = lambda(ones(rows*columns,1)); else lambda = lambda(:); end %Initialize r to zero. r = zeros(rows, columns); j = (1:(rows*columns))';
% indices remaining to generate
% For large lambda, use the method of Ahrens and Dieter as % described in Knuth, Volume 2, 1998 edition. k = find(lambda >= 15);
if ~isempty(k) alpha = 7/8; lk = lambda(k); m = floor(alpha * lk); % Generate m waiting times, all at once x = gamrnd(m,1); t = x <= lk; % If we did not overshoot, then the number of additional times % has a Poisson distribution with a smaller mean. r(k(t)) = m(t) + poissrnd(lk(t)-x(t)); % If we did overshoot, then the times up to m-1 are uniformly % distributed on the interval to x, so the count of times less % than lambda has a binomial distribution. r(k(~t)) = binornd(m(~t)-1, lk(~t)./x(~t)); j(k) = []; end % For small lambda, generate and count waiting times. p = zeros(length(j),1); while ~isempty(j) p = p - log(rand(length(j),1)); kc = [1:length(k)]'; t = (p < lambda(j)); j = j(t); p = p(t); r(j) = r(j) + 1; end % Return NaN if LAMBDA is negative. r(lambda < 0) = NaN;
Appendix 1.3D function [errorcode, rows, columns] = rndcheck(nargs,nparms,arg1,arg2,arg3,arg4,arg5) %RNDCHECK error checks the argument list for the random number generators. % % %
B.A. Jones 1-22-93 Copyright (c) 1993-98 by The MathWorks, Inc. $Revision: 2.5 $ $Date: 1997/11/29 01:46:40 $
sizeinfo = nargs - nparms; errorcode = 0; if nparms == 3 [r1 c1] = size(arg1); [r2 c2] = size(arg2); [r3 c3] = size(arg3); end if nparms == 2 [r1 c1] = size(arg1); [r2 c2] = size(arg2); end if sizeinfo == 0 if nparms == 1 [rows columns] = size(arg1); end if nparms == 2 scalararg1 = (prod(size(arg1)) == 1); scalararg2 = (prod(size(arg2)) == 1); if ~scalararg1 & ~scalararg2 if r1 ~= r2 | c1 ~= c2 errorcode = 1; return; end end if ~scalararg1 [rows columns] = size(arg1); elseif ~scalararg2 [rows columns] = size(arg2); else [rows columns] = size(arg1); end end if nparms == 3 scalararg1 = (prod(size(arg1)) == 1); scalararg2 = (prod(size(arg2)) == 1); scalararg3 = (prod(size(arg3)) == 1); if ~scalararg1 & ~scalararg2 if r1 ~= r2 | c1 ~= c2 errorcode = 1;
return; end end if ~scalararg1 & ~scalararg3 if r1 ~= r3 | c1 ~= c3 errorcode = 1; return; end end if ~scalararg3 & ~scalararg2 if r3 ~= r2 | c3 ~= c2 errorcode = 1; return; end end if ~scalararg1 [rows columns] = size(arg1); elseif ~scalararg2 [rows columns] = size(arg2); else [rows columns] = size(arg3); end end end if sizeinfo == 1 scalararg1 = (prod(size(arg1)) == 1); if nparms == 1 if prod(size(arg2)) ~= 2 errorcode = 2; return; end if ~scalararg1 & arg2 ~= size(arg1) errorcode = 3; return; end if (arg2(1) < 0 | arg2(2) < 0 | arg2(1) ~= round(arg2(1)) | arg2(2) ~= round(arg2(2))), errorcode = 4; return; end rows = arg2(1); columns = arg2(2); end if nparms == 2 if prod(size(arg3)) ~= 2 errorcode = 2; return; end scalararg2 = (prod(size(arg2)) == 1); if ~scalararg1 & ~scalararg2 if r1 ~= r2 | c1 ~= c2 errorcode = 1; return;
end end if (arg3(1) < 0 | arg3(2) < 0 | arg3(1) ~= round(arg3(1)) | arg3(2) ~= round(arg3(2))), errorcode = 4; return; end if ~scalararg1 if any(arg3 ~= size(arg1)) errorcode = 3; return; end [rows columns] = size(arg1); elseif ~scalararg2 if any(arg3 ~= size(arg2)) errorcode = 3; return; end [rows columns] = size(arg2); else rows = arg3(1); columns = arg3(2); end end if nparms == 3 if prod(size(arg4)) ~= 2 errorcode = 2; return; end scalararg1 = (prod(size(arg1)) == 1); scalararg2 = (prod(size(arg2)) == 1); scalararg3 = (prod(size(arg3)) == 1); if (arg4(1) < 0 | arg4(2) < 0 | arg4(1) ~= round(arg4(1)) | arg4(2) ~= round(arg4(2))), errorcode = 4; return; end if ~scalararg1 & ~scalararg2 if r1 ~= r2 | c1 ~= c2 errorcode = 1; return; end end if ~scalararg1 & ~scalararg3 if r1 ~= r3 | c1 ~= c3 errorcode = 1; return; end end if ~scalararg3 & ~scalararg2 if r3 ~= r2 | c3 ~= c2 errorcode = 1;
return; end end if ~scalararg1 if any(arg4 ~= size(arg1)) errorcode = 3; return; end [rows columns] = size(arg1); elseif ~scalararg2 if any(arg4 ~= size(arg2)) errorcode = 3; return; end [rows columns] = size(arg2); elseif ~scalararg3 if any(arg4 ~= size(arg3)) errorcode = 3; return; end [rows columns] = size(arg3); else rows = arg4(1); columns = arg4(2); end end end if sizeinfo == 2 if nparms == 1 scalararg1 = (prod(size(arg1)) == 1); if ~scalararg1 [rows columns] = size(arg1); if rows ~= arg2 | columns ~= arg3 errorcode = 3; return; end end if (arg2 < 0 | arg3 < 0 | arg2 ~= round(arg2) | arg3 ~= round(arg3)), errorcode = 4; return; end rows = arg2; columns = arg3; end if nparms == 2 scalararg1 = (prod(size(arg1)) == 1); scalararg2 = (prod(size(arg2)) == 1); if ~scalararg1 & ~scalararg2 if r1 ~= r2 | c1 ~= c2 errorcode = 1; return; end end if ~scalararg1 [rows columns] = size(arg1); if rows ~= arg3 | columns ~= arg4
errorcode = 3; return; end elseif ~scalararg2 [rows columns] = size(arg2); if rows ~= arg3 | columns ~= arg4 errorcode = 3; return; end else if (arg3 < 0 | arg4 < 0 | arg3 ~= round(arg3) | arg4 ~= round(arg4)), errorcode = 4; return; end rows = arg3; columns = arg4; end end if nparms == 3 scalararg1 = (prod(size(arg1)) == 1); scalararg2 = (prod(size(arg2)) == 1); scalararg3 = (prod(size(arg3)) == 1); if ~scalararg1 & ~scalararg2 if r1 ~= r2 | c1 ~= c2 errorcode = 1; return; end end if ~scalararg1 & ~scalararg3 if r1 ~= r3 | c1 ~= c3 errorcode = 1; return; end end if ~scalararg3 & ~scalararg2 if r3 ~= r2 | c3 ~= c2 errorcode = 1; return; end end if ~scalararg1 [rows columns] = size(arg1); if rows ~= arg4 | columns ~= arg5 errorcode = 3; return; end elseif ~scalararg2 [rows columns] = size(arg2); if rows ~= arg4 | columns ~= arg5 errorcode = 3; return; end
elseif ~scalararg3 [rows columns] = size(arg3); if rows ~= arg4 | columns ~= arg5 errorcode = 3; return; end else if (arg4 < 0 | arg5 < 0 | arg4 ~= round(arg4) | arg5 ~= round(arg5)), errorcode = 4; return; end rows = arg4; columns = arg5; end end end