Predavanje 19

FlipView[{Text[Style[" Branko Urošević, imqfbrankourosevic@gmail.com ",20,"Text",Italic]],Text[Style[" Faculty of Economics, University of Belgrade, NBS and CESIfo Munich",20,"Text", Italic]]}]

Teorija verovatnoće u WL

Kako modelirati šuteve na koš jednog košarkaša?

U podaci19.xlsx, u drugom Sheetu imamo podatke o 100 šuteva na koš. Naime, imamo informacije o rednom broju šuta kao i da li je igrač pogodio ili nije. Unesimo ih na standardni način

sutevi=Import["C:\\Users\\Prof Branko Urosevic\\Dropbox\\Branko Urosevic\\operaciona2\\downloads\\podaci19.xlsx"][[2]]

Radi dalje analize ovih podataka, transformišimo ih tako da “DA” -> 1, “NE”->0, i izvadimo samo drugu kolonu (bez naslova). U tom slučaju

sutevi1=Rest[sutevi/.{"DA"->1,"NE"->0}][[All,2]]

Koliko je dobar bio ovaj košarkaš? Najjednostavniji način da se o tome nešto zaključi jeste da se nađe srednja vrednost druge kolone, odnosno ovih jedinica i nula. Taj broj daje informaciju o tome koliko je puta košarkaš pogodio, u odnosu na ukupan broj puta.

Total[sutevi1]/Length[sutevi1]

Na osnovu ovu uzorka vidimo da je, u proseku, košarkaš pogađao 50 posto šuteva. To smo jednostavnije mogli da izračunamo na osnovu komande Mean

Mean[sutevi1]

Naravno, u pitanju je samo jedan relativno mali uzorak. Ono što bi bilo mnogo bolje jeste ukoliko bismo znali raspodelu verovatnoće koja najbolje odgovara ovim podacima. Ispostavljla se da je to Bernoulli raspodela.

Šta je Bernoulli raspodela

Bernoulli raspodela, kao i naš primer, ima moguće ishode 1 i 0, sa nekom konstantnom verovatnoćom p.  To je ujedno i jedini parametar koji definiše ovu raspodelu. To je najjednostavniji primer diskretne raspodele verovatnoće. Pretpostavimo da je verovatnoća da ćemo pogoditi u koš p=1/2. U tom slučaju, raspodela verovatnoće je:

𝒟=BernoulliDistribution[1/2];

Definisanje raspodele verovatnoće omogućava da koristimo ogroman aparat koji u WL postoji za rad sa njima. Ukoliko želimo da izračunamo verovatnoću realizacije x=0 ili x=1, to jest da dobijemo matematički izraz za raspodelu verovatnoće koristimo PDF (dolazi od reči Probability Density Function). Sintaksa je PDF[raspodela, x]:

PDF[𝒟,x]

Ako nas zanima kumulativna raspodela, to jest verovatnoća da je realizacija slučajne promenljive manje ili jednaka x, dobijamo:

CDF[𝒟,x]

Primetimo da je CDF[𝒟,x] suma PDF[𝒟,y] po y, gde y ide od početne realizacije do realizacije x. To važi za diskretne raspodele. U slučaju neprekidnih raspodela, ta suma postaje integral.

Očekivana vrednost Bernoulli slučajne promenljive

Pošto smo definisali funkciju raspodele nije problem da to umesto nas uradi WL. Ima više načina za to. Jedan način je da koristimo funkciju Mean. Ona izračunava ne samo uzoračku sredinu nego i teorijsku očekivanu vrednost za datu raspodele verovatnoće. Sintaksa je u tom slučaju Mean[ raspodela ]:

Mean[𝒟]

Srednja vrednost realizacije slučaje promenljives, stoga, kod Bernoulli raspodele je jednak verovatnoći da će pojedinačni šut ući u koš. Drugim rečima, iz procenjene srednje vrednosti dobijamo ocenu za jedini parametar koji treba da znamo da bismo definisali ovu raspodelu u potpunosti.

Drugi način da izračunamo očekivanu vrednost neke slučajne promenljive je da iskoristimo funkciju Expectation. Znak  \[Distributed]  (raspodeljeno) kuca se kao Esc dist Esc. Sintaksa je

Expectation[ x, x \[Distributed] raspodela] ili, u opštem slučaju, Expectation[ f[x], x \[Distributed] raspodela]

ocekivanaVrednost=Expectation[x,x\[Distributed]𝒟]

Varijansa i standardna devijacija

Varijansa ili na engleskom Variance je mera kvadratnog odstupanja od srednje vrednosti. Za zadatu raspodelu 𝒟 ona se najjednostavnije izračunava ovako:

Variance[𝒟]

Evo i drugog načina za izračunavanja varijanse:

varijansa=Expectation[(((x-ocekivanaVrednost))^2),x\[Distributed]𝒟]

Varijansa meri stepen varijabilnosti reaiizacije slučajne promenljive. Međutim, pošto sadrži kvadrat promenljive, dimenziono se ne može uporediti sa očekivanom vrednošću. U tu svrhu najčešće se koristi standardna devijacija. To je kvadratni koren od varijanse. Ona se može izračunati, opet, na više načina:

StandardDeviation[𝒟]

ili preko varijanse:

standardnaDevijacija=(Sqrt[varijansa])

Ocena očekivane vrednosti, varijanse i standardne devijacije iz uzorka

Da li na osnovu uzorka od 100 šuteva možemo biti sigurni da je verovatnoća da će igrač pogoditi u jednom šutu biti 1/2? Razume se da ne možemo u to biti sigurni. Pretpostavimo, naime, da je zaista verovatnoća pogotka 0.5 i napravimo novi slučajni uzorak dužine 100.

uzorak=RandomVariate[𝒟,100]

RandomVariate[raspodela, veličina uzorka] generiše za zadatu raspodelu (prvi argument) jedan broj slučajnih realizacija te slučajne promenljive (dimenzije uzorka su date drugim argumentom). Evo kako možemo da bacimo fer novčić 100 puta koristeći RandomVariate:

Uzoračka sredina i standardna devijacija izračunavaju se pomoću istih funkcija kao i teorijske vrednosti:

Mean[uzorak]//N

StandardDeviation[uzorak]//N

Dakle, kada napravimo novi uzorak od veličine 100, nećemo obavezno dobiti istu ocenu za srednju vrednost (odnosno za verovatnoću). To možemo znatno popraviti ukoliko povećamo broj šuteva koje “snimamo” u uzorku.

Model prinosa na akciju koja ima dve podjedako verovatne realizacije: 1/2 ili -1/3

Pretpostavimo da ako uložimo novac u neku akciju možemo da kroz godinu dana ili zaradimo 50 posto, ili izgubimo jednu trećinu svog novca, i to sa podjednakom verovatnoćom:

TreePlot[{{0->(1/2),"(1/2)"},{0->-(1/3),"(1/2)"}},Left,DirectedEdges->True,VertexLabeling->True,VertexCoordinateRules->{{0,0},{1,0.5},{1,-0.5}},ImageSize->250,PlotLabel-> Style["Godišnji prinos na akciju 
sa dva jednako verovatna stanja",14,Bold,Blue]]

Kako modelirati ovu slučajnu promenljvu? Za to nam je potrebna takozvana transformisana Bernoulli promenljiva. To se ostvaruje preko TransformedDistribution funkcije. Ovo je veoma korisna i veoma moćna funkcija za rad sa raspodelama verovatnoća. Njena sintaksa je:

NovaRaspodela = TransformedDistribution[ f [x], x \[Distributed] Stara Raspodela]

Naći transformaciju Bernoulli promenljive

Originalna Bernoulli promenljiva x ima realizacije 1 i 0. Mi želimo da je transformišemo tako da nova promenljiva ima realizaciju x1=1/2 ukoliko je x=1, kao i vrednost x1= -1/3 kada je x=0. Pretpostavimo da je veza između stare i nove slučajne promenljive linearna, to jest da, za neke nepoznate a i b, važi da je:

 x1=a x + b 
 
Tražimo takve vrednosti za a i b da se zadovolje prethodni uslovi:

Solve[{1*a  +b==1/2,0*a  +b==-1/3},{a,b}]

Novu raspodelu dobijamo na sledeći način:

𝒟1=Module[{𝒟=BernoulliDistribution[1/2]},TransformedDistribution[(5/6) x-(1/3),x\[Distributed]𝒟]];

Da bismo proverili da zaista generišemo odgovarajuće vrednosti, napravimo uzorak nove promenljive

RandomVariate[𝒟1,100]

Sa novom funkcijom raspodele možemo da radimo sve što smo radili i sa starom

Primera radi, očekivani prinos na ovu akciju jednak:

1/2 *(1/2) + 1/2*(-1/3) =1/12

Razume se, on se može izračunati i ovako

Mean[𝒟1]

Expectation[r,r\[Distributed]𝒟1]

Varijansa i standardna devijacija su:

{Variance[𝒟1],StandardDeviation[𝒟1]}

Numeričke vrednosti ovih veličina su:

{Mean[𝒟1],Variance[𝒟1],StandardDeviation[𝒟1]}//N

Možemo da generišemo slučajni uzorak i izračunamo uzoračku sredinu i standardnu devijaciju:

Module[{sample1=RandomVariate[𝒟1,100000]},{Mean[sample1],StandardDeviation[sample1]}//N]

Najzad, evo kako možemo da izračunamo očekivanu vrednost, recimo, kuba prinosa E[ (r^3)]:

Expectation[(r^3),r\[Distributed]𝒟1]

Opšti slučaj Bernoulli raspodele

Sve što smo radili za konkretnu Bernoulli raspodelu može da se radi i u opštem slučaju. Ukoliko je verovatnoća da ćemo pogoditi koš  u jednom bacanju jednaka p a da ćemo promašiti jednaka 1-p, Bernoulli raspodela se definiše na sledeći način:

𝒟=BernoulliDistribution[p];

Ovako dobijamo očekivanu vrednost i varijansu

{Mean[𝒟],Variance[𝒟]}

Evo kako linearno transformišemo Bernoulli promenljivu u opštem slučaju:

𝒟1=TransformedDistribution[a x +b,x\[Distributed]𝒟];

Očekivana vrednost i standardna devijacija transformisane promenljive su:

{Mean[𝒟1],Variance[𝒟1]}//Simplify

2. Proizvod Bernoulli raspodela: vežbanje 3 slobodna bacanja

Pretpostavimo da neki drugi košarkaš vežba slobodna bacanja tako što baca puno puta (100) seriju od po 3 bacanja. Rezultati se mogu naći u podaci19.xlsx fajlu, Sheet 3.

podacitriBacanja=Import["C:\\Users\\Branko Urosevic\\Dropbox\\Branko Urosevic\\operaciona2\\downloads\\podaci19.xlsx"][[3]]

Vidimo da za svaku seriju od po 3 šuta imamo podatak o tome da li je koš pogođen ili ne. Možemo da transformišemo ove podatke tako da umesto DA i NE imamo 1 i 0, respektivno. Ako, pri tome izdvojimo samo podatke, vidimo oni izgledaju ovako:

podaci3=Rest[podacitriBacanja/.{"DA"->1,"NE"->0}][[All,2;;4]]

Dakle, u prvoj seriji, igrač je promašio prvo a pogodio druga dva slobodna bacanja. U drugoj seriji, promašio je prva dva a pogodio treće, itd. Pošto je svako bacanje nezavisno od drugih, slučajna promenljiva “ishod tri slobodna bacanja” može se predstaviti kao proizvod 3 nezavisne identične Bernoulli raspodele. Da bismo procenili verovatnoću jednog pogodka, koristimo uzoračku sredinu ukupnog broja pogodata u odnosu na celokupan broj šuteva.

Mean[Flatten[podaci3]]//N

Definišimo funkciju raspodele koristeći ProductDistribution funkciju za 3 identične Bernoulli asphodel sa parametrom p

𝒟p=Module[{𝒟=BernoulliDistribution[0.39]},ProductDistribution[{𝒟,3}]];

uzorak=RandomVariate[𝒟p,100]

Ideja binomne raspodele

Vidimo da proizvod Bernoulli raspodela daje detaljnu informaciju o tome koji šut neke  serije jeste a koji nije završio u košu. Često to je previše detaljna informacija za naš ukus. Naime, želimo da znamo samo koliko je ukupno koševa postignuto u svakoj od serija. Očigledno taj broj može biti 0,1,2, ili 3, pri čemu se neki od tih ishoda verovatniji nego drugi. U slučaju našeg “snimka” treninga imali bismo sledeći rezultat:

Total/@podaci3

Histogram[%]

Ukoliko je procena p=0.39 uzeta kao tačna, teorijska raspodela verovatnoće koja odgorava ovakvoj realizaciji bi bila

𝒟ukup=TransformedDistribution[x1+x2+x3,{x1,x2,x3}\[Distributed]𝒟p]

U pitanju je, dakle, binomna raspodela. To je raspodela sume od 3 (ili u opštem slučaju n) istovetno raspodeljenih Bernoulli slučajnih promenljivih

Binomna raspodela: opšti slučaj

Vidimo da je rezultujuća raspodela binomna raspodela.  U opštem slučaju, binomna raspodela se definiše preko dva parameta, ukupnog broja bacanja (n) i verovatnoće da će u jednom bacanju izaći pismo (p). Sintaksa je BinomialDistribution[n,p]. PDF binomne raspodele je:

PDF[BinomialDistribution[n,p],k]//Refine[#,k≥ 0&&k≤ n]&

Funkcija Refine služi da uprostimo izraze prilikom nekih izračunavanja. Binomial je broj kombinacija n nad k. To je broj načina na koji se može dogoditi da iz n bacanja imamo tačno k pogodaka (pisama).

Binomial[n,k]//TraditionalForm

Na primer:

Binomial[7,3]

Podsetimo se da se ovo izračunava ovako:

(((7!))/(3! (7-3)!))

Očekivana vrednost i varijansa binomne raspodele su u opštem slučaju date poznatim izrazima:

Module[{𝒟=BinomialDistribution[n,p]},{Mean[𝒟],Variance[𝒟]}]

Tri slobodna bacanja

Isti igrač koji je trenirao prethodno bacanje tri slobodna bacanja sada učestvuje na utakmici. Ima tri slobodna bacanja u odlučujućem momentu. Verovatnoća da pogodi jedno bacanje je 0.39. Nađite verovatnoću da će pogoditi sva 3. Verovatnoću nekog događaja, u slučaju diskretne promenljive, možete da računate preko funkcije Probability. Sintaksa je Probability[ f[x] ==nesto, x \[Distributed]raspodela ]. Odgovajuća raspodela je:

𝒟3=BinomialDistribution[3,0.39];

Probability[pogotka==3,pogotka\[Distributed]𝒟3]

To možemo da proverimo i ovako. Prvo, da bismo imali 3 od 3, moramo da pogodimo sva tri. To će se desiti sa verovatnoćom

(0.39^3)

To se, takođe, može izračunati i ovako:

PDF[𝒟3,3]

Koliko je šansa da se spetljate i ne pogodite nijedno?

Probability[pogotka==0,pogotka\[Distributed]𝒟3]

(((1-0.39))^3)

Nije teško naći i verovatnoće da ćemo pogoditi jedno odnosno dva od tri slobodna bacanja

{PDF[𝒟3,1],PDF[𝒟3,2]}

{0.39 (((1-0.39))^2)Binomial[3,1], (0.39^2) (1-0.39)Binomial[3,2]}

Nastavak prethodnog primera

Pretpostavimo da je dodatno vežbao i 200 puta bacio 3 slobodna bacanja. Evo moguće realizacije (opet koristimo RandomVariate)

uzorak=RandomVariate[𝒟3,200]

Histogram[uzorak]

Vidimo da raspodela počinje da liči pomalo na normalnu raspodelu. To nije uopšte slučajno. Takođe vidimo da je 2 najverovatniji ishod. Nađimo očekivanu vrednost i standardnu devijaciju za uzorak i teorijski:

{Mean[uzorak],StandardDeviation[uzorak]}//N

Sa druge strane, teorijske vrednosti su

{Mean[BinomialDistribution[3,0.39]],StandardDeviation[BinomialDistribution[3,0.39]]}

3. Uniformna raspodela

Šta je uniformna raspodela

Razmotrimo sada slučaj jedne od najjednostavnijih neprekidnih raspodela. Pretpostavimo da promenljiva x može imati bilo kakvu vrednost od 0 do 1, i da su sve vrednosti podjednako verovatne. Mi bismo i bez znanja funkcije raspodele mogli da vršimo uzorkovanje tako što bismo iskoristili RandomReal:

RandomReal[1,{50}]

Histogram[RandomReal[1,{100000}]]

Međutim, zgodno je definisati direktno slučajnu promeljivu koja ima takvu raspodelu. Raspodela je UniformDistribution[{a,b}]. U slučaju kada je donja granica a=0, a gornja b=1, imamo

𝒟=UniformDistribution[{0,1}];

Kako izgledaju gustina verovatnoće PDF i kumulativna funkcija raspodele za standardnu uniformnu raspodelu? One su date sledećim izrazima:

PDF[𝒟,x]

CDF[𝒟,x]

Izbor tačke na lenjiru dužine 1 metar

Pretpostavimo da imamo lenjir veličine 1 metar. Zatvorili su nam oči i rekli da nasumično odredimo neku tačku na tom metru. Naći verovatnoću da tačka neće biti veća od 60 cm ili 0.6 metara. Ova verovatnoća je:

CDF[𝒟,0.6]

Pošto je uniformna raspodela neprekidna, umesto sume verovatnoća sada imamo integral od gustine verovatnoće PDF. Stoga je CDF[𝒟,0.6] isto što i:

Integrate[PDF[𝒟,x],{x,0,0.6}]

Verovatnoća je površina ispod prave linije (osenčena žutim). Pošto je u pitanju pravougaonik, površina je 0.6.

Plot[PDF[𝒟,x],{x,0,1},Mesh->{{0}},MeshShading->{Directive@{Thick,Blue}},Filling->Axis,ColorFunction->(If[0.≤ #<0.6,Yellow,White]&),ColorFunctionScaling->False]

Izbor tačke na delu metra

Razmotrimo sada situaciju ako nas zanimaju samo tačke koje nisu manje od 0.2 i ne veće od 0.8 (na primer, zalepilli smo selotejp na tački 0.2 i 0.8, pa kad naslepo napipavamo znamo da ne treba da idemo levo od 0.2 ili desno od 0.8). U tom slučaju, raspodela verovatnoće je data sa:

𝒟=UniformDistribution[{0.2,0.8}];

Sada imamo:

PDF[𝒟,x]

Plot[PDF[𝒟,x],{x,0,2},Filling->Axis]

Nađimo sada verovatnoću da je izabrana tačka manja ili jednaka 0.6

CDF[𝒟,0.6]

Zašto je ova verovatnoća veća nego kada su moguće vrednosti od 0 do 1? Pa zato što znamo da smo u užem intervalu koji sadrži 0.6. Ova verovatnoća može da se izračuna i ovako:

Integrate[PDF[𝒟,x],{x,0.2,0.6}]

Plot[PDF[𝒟,x],{x,0,1},Mesh->{{0}},MeshShading->{Directive@{Thick,Blue}},Filling->Axis,ColorFunction->(If[0.≤ #<0.6,Yellow,White]&),ColorFunctionScaling->False]

Opšti slučaj uniformne raspodele, očekivana vrednost i standardna devijacija

U opštem slučaju, uniformna raspodela u jednoj dimenziji je:

𝒟=UniformDistribution[{a,b}];

Gustina verovatnoće i cumulativna funkcija raspodele su, u tom slučaju, date izrazima:

PDF[𝒟,x]

CDF[𝒟,x]

Očekivana vrednost ove promenljive je:

ocekivanaVrednost=Mean[𝒟]

Integrate[x PDF[𝒟,x],{x,a,b}]//Refine[#,a<b]&//Quiet

a varijansa

Variance[𝒟]

Integrate[(((x-ocekivanaVrednost))^2) PDF[𝒟,x],{x,a,b}]//Refine[#, a<b]&//Quiet

4. Normalna raspodela (NormalDistribution)

Uvod u normalnu promenljivu u jednoj dimenziji

Normalna raspodela je najčešće korišćena funkcija raspodele. U finansijama, na primer, često se pretpostavlja da su prinosi na akcije normalno raspodeljeni.

𝒟=NormalDistribution[μ,σ];

PDF[𝒟,x]

{Mean[𝒟],StandardDeviation[𝒟]}

Vidimo da je μ očekivana vrednost a σ je standardna devijacija normalne slučajne promenljive. U  slučaju prinosa, μ bi bio očekivani prinos a σ bi bila mera rizika koja odgovara raspodeli prinosa.

Manipulate[Show[{Plot[{((Exp[-((((r-μ))^2)/(2 (σ^2)))])/(Sqrt[2 π (σ^2)]))},{r,-10,10},PlotRange-> {0,0.8},ImageSize->400,Epilog-> Inset[Grid[{{"μ-2 σ",μ-2 σ},{"μ+2 σ",μ+2 σ},{"Shaded area:",CDF[NormalDistribution[μ,σ],μ+2 σ]-CDF[NormalDistribution[μ,σ],μ-2 σ]}},Frame-> True, Alignment-> Center, Spacings-> {2,1},ItemStyle->Directive[Bold,15]],{5,0.6}]],Plot[((Exp[-((((r-μ))^2)/(2 (σ^2)))])/(Sqrt[2 π (σ^2)])),{r,μ-2 σ,μ+2 σ},PlotRange-> {0,0.8},ImageSize->600,Filling->Bottom],Graphics[{PointSize[0.02],Green,Point[{μ-2 σ,0}]}],Graphics[{PointSize[0.02],Green,Point[{μ+2 σ,0}]}]}],{μ,-2,2,Appearance-> "Labeled"},{σ,0.5,2,Appearance-> "Labeled"},SaveDefinitions->True]

Linearna transformacija normalne slučajne promenljive je normalna

Važno svojstvo normalnih promenljivih jeste da je inearna transformacija normalne promenljive opet normalna promenljiva.

Module[{𝒟=NormalDistribution[μ,σ]},TransformedDistribution[a x+b,x\[Distributed]𝒟]]

sa očekivanom vrednošću a μ +b i varijansom jednakom  (σ^2) (a^2).

Ukoliko je x normalna promenljiva sa očekivanom vrednošću μ i standardnom devijacijom jednakom σ, onda je z =( x -μ)/σ  normalna promenljiva sa očekivanom vrednošću 0 i standardnom devijacijom 1. To je takozvana standardna normalna promenljiva.

𝒟s=TransformedDistribution[(x-μ)/σ,x\[Distributed]NormalDistribution[μ,σ]]

U WL, standardna normalna raspodela se označava prosto kao NormalDistribution[]. Gustina verovatnoće i kumulativna funkcija raspodele za standardnu normalnu promenljivu su date izrazima.

{PDF[𝒟s,x],CDF[𝒟s,x]}

Linearna kombinacija nezavisnih normalnih promenljivih je normalna

𝒟1=NormalDistribution[(Subscript[μ, 1]),(Subscript[σ, 1])];𝒟2=NormalDistribution[(Subscript[μ, 2]),(Subscript[σ, 2])];TransformedDistribution[a1 x1+a2 x2,{x1\[Distributed]𝒟1,x2\[Distributed]𝒟2}]

Viši momenti normalne raspodele

Mada su оčekivana vrednost i varijansa najčešće korišćeni momenti normalne raspodele, korisno je znati i kako se izračunavaju i momenti višeg reda, bilo da su centrirani (oko očekivane vrednosti) ili necentrirani (oko 0). Ovo pokazuje da je normalna raspodela simetrična i ima debele repove.

{Skewness[NormalDistribution[μ,σ]],Kurtosis[NormalDistribution[μ,σ]]}

Tabela necentriranih momenata normalne raspodele NormalDistribution[μ,σ]:

TableForm[Table[{k,Expand@Moment[NormalDistribution[μ,σ],k]},{k,8}],TableHeadings->{None,{"k","E[(x^k)]"}}]

A evo i centriranih momenata te iste raspodele

TableForm[Table[{k,CentralMoment[NormalDistribution[μ,σ],k]},{k,8}],TableHeadings->{None,{"k","E[(x-E[x])^k]"}}]

Vidimo da je u WL spljoštenost (kurtosis) definisan kao četvrti centralni moment podeljen sa kvadratom varijanse, odnosno četvrtim stepenom standardne devijacije (σ^4). U specijalnom slučaju standardne normalne raspodele, prvih 8 necentriranih momenata su:

Moment[NormalDistribution[0,1],#]&/@Range[8]

Kako bismo to izračunali iz prvih principa:

PDF[NormalDistribution[],x]

Integrate[(x^(#)) PDF[NormalDistribution[],x],{x,-∞,∞}]&/@Range[8]

Intervali poverenja

Naći verovatnoće da vrednost normalne promenljive date raspodelom NormalDistribution[μ,σ] leži u invervalu od 2 standardne devijacije u odnosu na srednju vrednost:

Module[{𝒟=NormalDistribution[μ,σ]},Probability[μ-2σ<x<μ+2σ,x\[Distributed] 𝒟]]

N[%]

To se može predstaviti i ovako:

CDF[𝒟,μ+2σ]-CDF[𝒟,μ-2σ]//N

Ponovimo to isto ali za interval od 3 standardne devijacije od srednje vrednosti:

Module[{𝒟=NormalDistribution[μ,σ]},Probability[μ-3 σ<x<μ+3σ,x\[Distributed] NormalDistribution[μ,σ]]]

N[%]

Pretpostavimo da je životni vek sijalice normalno raspodeljen tako da je očekivani životni vek 12 a standardna devijacija 2 meseca. Naći verovatnoću da će sijalica raditi najmanje 18 meseci pre nego što pregori.

Module[{𝒟=NormalDistribution[12,2]},Probability[x≥ 18,x\[Distributed]𝒟]]

N[%]

Vidimo da je ta verovatnoća oko 0.1 posto, to jest veoma mala.

Verovatnoća gubitka ili dobitka novca prilikom investiranja u akciju

Pretpostavimo da naša investicija ima očekivani prinos od 15 posto godišnje kao i rizik, meren preko standardne devijacije, od 25 posto godišnje. Prinos na investiciju je normalna promenljiva. Nađi verovatnoću da ćemo kroz godinu dana biti u gubitku kao i verovatnoću da ćemo ostvariti dobitak od najmanje 10 posto na kraju godine.

Probability[r<0,r\[Distributed]NormalDistribution[0.15,0.25]]

Probability[r>0.1,r\[Distributed]NormalDistribution[0.15,0.25]]

Kako možemo da ove verovatnoće da izračunamo iz prvih principa? Ukoliko uvedemo oznaku:

𝒟=NormalDistribution[0.15,0.25];

Ove verovatnoće su, onda, integrali

{NIntegrate[PDF[𝒟,x],{x,-∞,0}],NIntegrate[PDF[𝒟,x],{x,0.1,∞}]}

Grafički to su površine ispod krive

Plot[PDF[𝒟,x],{x,-1.5,1.5},Mesh->{{0}},MeshShading->{Directive@{Thick,Blue}},Filling->Axis,ColorFunction->(If[-1.5≤ #<0.,Yellow,White]&),ColorFunctionScaling->False,PlotRange-> All]

Plot[PDF[𝒟,x],{x,-1.5,1.5},Mesh->{{0}},MeshShading->{Directive@{Thick,Blue}},Filling->Axis,ColorFunction->(If[0.1≤ #<1.5,Yellow,White]&),ColorFunctionScaling->False,PlotRange-> All]

Kako na osnovu podataka možemo da ocenimo o kojoj raspodeli se radi?

Razmorimo uzorak

uzorak=RandomVariate[NormalDistribution[0.7,0.3],10^3];

Ukoliko pretpostavimo da je u pitanju uzorak iz normalne raspodele, parametre možemo da probamo da ocenimo preko funkcije EstimatedDistribution (rada na osvovu metoda maksimalne verodostojnosti).

edist=EstimatedDistribution[uzorak, NormalDistribution[m,sigma]]

Proverimo sada koliko dobro fituje ova raspodela. U principu, ne možemo da je odbacimo ukoliko je P vrednost dovoljno velika (bar 0.1). Zaista:

DistributionFitTest[uzorak,edist,{"TestDataTable",All}]

Dakle, to radi odlično. E sada, pretpostavimo da uopšte ne znamo o kojoj se raspodeli radi (ovo za sada radi samo za jednodimenzione raspodele):

FindDistribution[uzorak]

Ovo daje najbolju distribuciju po defoltu. Međutim, možemo da uzmemo i prve 3 po BIC i AIC kriterijumima (to su takozvani informacioni kriterijumi)

FindDistribution[uzorak,3,{"BIC", "AIC"}]

Ovako možemo da nađemo detaljniji izveštaj

report=FindDistribution[uzorak,3,All]

5. Multinormalna raspodela (MultinormalDistribution)

Multinormalna raspodela je raspodela vektorske promenljive zadate vektorom očekivanih vrednosti μ i matricom kovariansi Σ. Sintaksa je MultinormalDistribution[ μ, Σ].

Ako znamo raspodelu svih slučajnih promenljivih u skupu, iz nje možemo naći raspodele pojedinačnih promenljivih pomoću MarginalDistribution funkcije. Ukoliko imamo funkciju raspodele od dve (ili više) promenljivih, marginalnu raspodelu prve promenljive dobijamo integraljenjem  po svim ostalim promenljivim

g1[x]==∫f[x,y]dy

Marginalna raspodela za drugu drugu promenljivu:

g2[y]==∫f[x,y]dx

Ukoliko je raspodela promenljivih multinormalna sa vektorom očekivanih prinosa μ i matricom kovarijansi Σ

μ={μ1,μ2};Σ={{((σ1)^2),σ1 σ2 ρ},{σ1 σ2 ρ,((σ2)^2)}};

ovo su marginalne raspodelaeza prvu i drugu promenljivu:

Refine[MarginalDistribution[MultinormalDistribution[μ,Σ],1],σ1>0]

Refine[MarginalDistribution[MultinormalDistribution[μ,Σ],2],σ2>0]

Ako je raspodela multinormalna, sve marginalne raspodele su normalne.

Očekivane vrednost i varijansa MultinormalDistribution raspodele

Pošto je vektorska promenljiva neprekidno raspodeljena, očekivana vrednost vektorske promenljive je vektor očekivanih vrednosti čije su komponente date sa (u drugoj jednačini koristimo marginalne raspodele):

𝔼[x]==∫x f[x,y]dxdy ==∫x g1[x]dx

𝔼[y]==∫y f[x,y]dxdy==∫y g2[y]dy

Mean[MultinormalDistribution[μ,Σ]]

Prilikom izračunavanja očekivanih vrednosti mi treba da znamo samo odgovarajuće marginalne raspodele. Slična je situacija sa varijansom. Po definiciji, varijansa za prvu promenljivu je:

𝕍[x]==∫(((x-𝔼[x]))^2)f[x,y]dxdy ==∫(((x-𝔼[x]))^2)g1[x]dx

Na sličan način, za varijansu druge promenljive dobijamo:

𝕍[y]==∫(((y-𝔼[y]))^2)f[x,y]dxdy ==∫(((y-𝔼[y]))^2)g2[y]dy

U slučaju multinormalne promenljive imamo da je:

Variance[MultinormalDistribution[μ,Σ]]

Kako očekivane vrednosti tako i varijansa (a samim tim i standardne devijacije kao kvadratni koren iz varijanse), su veličine koje su karakteristike pojedinačnih komponenti vektorske slučajne promenljive. To nije slučaj sa kovarijansom.

Matrica kovarijanse u slučaju MultiNormal raspodele

Kovarijansa meri u kojoj meri su vrednosti prve promenljive u vezi sa vrednostima druge promenljive. Definiše se kao očekivana vrednost otklona od srednjih vrednosti te dve promeljive. Ukoliko x, y imaju zajedničku raspodelu, kovarijansa je:

cov[x,y]==𝔼[(x-𝔼[x])(y-𝔼[y])] ==

∫(x-𝔼[x])(y-𝔼[y])f[x,y]dxdy

U naše slučaju, matrica kovarianse je:

Covariance[MultinormalDistribution[μ,Σ]]

Na dijagonali su varijanse. Kovarijansa je σ1 σ2 ρ. Ovde je ρ koeficijent korelacije, a σ1 i σ2 su standardne devijacije. Koeficijent korelacije je kovarijansa normirana da bude između -1 i 1. Ako su u pitanju prinosi na akcije, korelacija jednaka 0 znači da ne postoji linearna statistička veza između prinosa na dve akcije. Naime, to šta se dešama sa prvom akcijom nema veza sa tim šta se dešava sa drugom. Ako je koeficijent korelacije jednak 1, dve akcije su potpuno korelalisane, to jest ako prva akcija raste, raste i druga i obrnuto. Sa druge strane, ako je koeficijent korelacije -1, ukoliko prva akcija raste druga opada i obrnuto.

Refine[Correlation[MultinormalDistribution[μ,Σ]],σ1>0&&σ2>0]//MatrixForm

Evo i centalnih momenata 3. i 4. reda:

Skewness[MultinormalDistribution[μ,Σ]]

Kurtosis[MultinormalDistribution[μ,Σ]]

Proceniti vrednosti parametara multinormalne raspodele

Generišemo neke podatke. Prvo to uradimo na malom uzorku. Dobijamo uzorke prinosa na prvu i drugu akcija sa zadatim parametrima.

RandomVariate[MultinormalDistribution[{1,4},{{2,-1/4},{-1/4,2/3}}],15]

Vidimo da je svaka simulacija par brojeva (tačke) pošto je u pitanju multi normalna promenljiva u 2 dimenzije. Sada generišemo uzorak od 1000 tačaka:

sample=RandomVariate[MultinormalDistribution[{1,4},{{2,-1/4},{-1/4,2/3}}],10^3];

Procenimo vrednosti parametara multinormalne raspodele:

Clear[μ,σ]

edist=EstimatedDistribution[sample, MultinormalDistribution[Array[μ,2],Array[σ,{2,2}]]]

Proverimo kako to sada radi na podacima (to jest, da li su P vrednosti dovoljno velike da možemo da odbacimo nultu hipotezu):

DistributionFitTest[sample,edist,{"TestDataTable",All}]

Prinos na portfolio koji se sastoji od dve akcije

Prepostavimo da je vektor očekivanih prinosa na dve akcije dat vektorom:

μp={0.1,0.15};

Definišimo matricu korelacije i vektor standardnih devijacija ovako (0.4 je koeficijent korelacije između prinosa na dve akcije):

corr={{1,0.4},{0.4,1}};sigma={0.25,0.35};

Prvo izračunajmo matricu kovarijanse (ovakva dekompozicija se u praksi često koristi):

Σp=DiagonalMatrix[sigma].corr.DiagonalMatrix[sigma]

Pretpostavimo da su prinosi multinormalno raspodeljeni. Onda bi funkcija raspodele bila:

𝒟=MultinormalDistribution[μp,Σp];

Nađimo prvo očekivane prinose i standardne devijacije za ove dve akcije. Dobijamo da su jednaki:

{Mean[𝒟],StandardDeviation[𝒟]}

A evo i matrice kovarijansi i matrice korelacije:

Covariance[𝒟]

Correlation[𝒟]

Simulacija prinosa na akcije za različite nivoe korelacije

Da bismo bolje sagledali vezu između korelacije i toga kako se ponašaju prinosi na akcije, razmotrimo uzorak od 100 prinosa uzorkovanih iz raspodele 𝒟. Prvo to uradimo u slučaju kadaje koeficijent korelacije 0.4

uzorak=RandomVariate[𝒟,100];

Graphics[{Red,PointSize[0.01],Point[uzorak]},Axes->True,AxesLabel->{"Prinos 1","Prinos 2"}]

Razmotrimo sada situaciju kada je sve isto kao i pre samo što je koeficijent korelacije 0.9. U tom slučaju

Σp1={{(0.25^2),0.25 0.35 0.9},{0.25 0.35 0.9,(0.35^2)}};𝒟1=MultinormalDistribution[μp,Σp1];

uzorak1=RandomVariate[𝒟1,100];

Graphics[{Red,PointSize[0.01],Point[uzorak1]},Axes->True,AxesLabel->{"Prinos 1","Prinos 2"}]

Vidimo da je sada trend mnogo izraženiji!

Naći verovatnoću gubitka od 15 posto na godišnjem nivou za pojedinačne akcije i za portfolio akcija

Vratimo se na situaciju sa korelacijom 0.4. Ukoliko bismo sav novac uložili u prvu akciju, verovatnoća gubitka od 15 posto bila bi jednaka:

Probability[r1<-0.15,{r1,r2}\[Distributed]𝒟]

a ukoliko bismo sve uložili u drugu, verovatnoće gubitka od 15 posto bila bi jednaka:

Probability[r2<-0.15,{r1,r2}\[Distributed]𝒟]

Kada napravimo portfolio od dve hartije tako da po 50 posto para uložimo u obe akcije, prinos portfolia bi bio: 0.5 r1+0.5 r2. U tom slučaju, verovatnoća gubitka od 15 posto bi bila:

Module[{𝒟port=TransformedDistribution[0.5 r1+0.5 r2,{r1,r2}\[Distributed]𝒟]},Probability[rp<-0.15,rp\[Distributed]𝒟port]]

Vidimo da je verovatnoća takvog gubitka manja kada investiramo u portfolio nego u pojedinačne akcije. Možemo da proverimo kakva će biti verovatnoća gubitka od 15 posto za različite kompozicije portfolio, odnosno za različite pondere. Ukoliko je w udeo novca koji smo uložili u prvu akciju, a 1-w udeo uložen u drugu akciju:

Table[{w,Module[{𝒟port=TransformedDistribution[w r1+(1-w) r2,{r1,r2}\[Distributed]𝒟]},Probability[rp<-0.15,rp\[Distributed]𝒟port]]},{w,0.1,0.9,0.05}]

ili, grafički

ListLinePlot[%,AxesLabel->{"w","Verovatnoća gubitka>15 posto"}]

Muzika iz filma “Dečko koji obećava (Beograd, 1981)”

https://www.youtube.com/watch?v=-UZ_2m6LAcI