VPFPM ZS2018/19 (Cvičení)

NMFM308 Výpočetní prostředky finanční a pojistné matematiky

NMFM308 Computational Tools for Finance and Insurance

$LicenseExpirationDate

Literatura

[1] Dupačová, J., Hurt, J., Štěpán, J.: Stochastic Modeling in Economics and Finance. Kluwer Academic Publishers. Dordrecht 2002.

[2] Cipra, T.: Finanční ekonometrie. Ekopress, Praha 2008.

[3] Cipra, T.: Matematika cenných papírů. Kamil Mařík Professional Publishing. Praha 2013.

[4] Hurt, J.: Yield curves with Mathematica 6.0. In: Wolfram Technology Conference 2007. http://library.wolfram.com/infocenter/Conferences/6956/ . Champaign (IL) 2007.

[5] Hurt, J.: Risk measures in finance. In: 2008 International Mathematica User Conference. http://library.wolfram.com/infocenter/Conferences/7230/.  Champaign (IL) 2008.

[6] Hurt, J.: Risk measures in finance revisited. In: Wolfram Technology Conference 2010. http://library.wolfram.com/infocenter/Conferences/7861/  . Champaign (IL) 2010.

[7] Sciencedirect:  http://www.sciencedirect.com/

[8] Mathematica Cookbook

[9] Luenberger, D.: Investment Science. Oxford University Press. New York 1998.

[10] Cipra, T.: Financial and Insurance Formulas. Springer-Verlag. Berlin Heidelberg 2010.

[11] Vencálek, Ondřej: Základy analýzy dat v softwaru Mathematica. Univerzita Palackého v Olomouci. Olomouc, 2015. Text je dostupný v knihovně Karlín a Malá Strana. Rovněž je přístupný elektronicky na adrese http://mathscience.upol.cz/materialy.html a to jak ve formátu pdf, tak i ve formě spustitelného “sešitu” Mathematica (.nb).

[12] An Elementary Introduction to the Wolfram Language by Stephen Wolfram: http://wolfr.am/eiwl

[13] Hands-on Start to Wolfram Mathematica and Programming with the Wolfram Language by Cliff Hastings, Kelvin Mischo, and Michael Morrison (Chapter 7 “Creating Interactive Models with a Single Command” available free online): http://www.handsonstart.com

[14] The Wolfram Language: Fast Introduction for Programmers
http://www.wolfram.com/language/fast-introduction-for-programmers/en
Zahrnuje též cvičení s odpovědmi.

[15] What is new 
http://www.wolfram.com/mathematica/quick-revision-history.html

[16] Mathematica & Wolfram Language
Fast Introduction for Math Students

http://www.wolfram.com/language/fast-introduction-for-math-students/en/

http://www.wolfram.com/language/fast-introduction-for-math-students/en/probability/

http://www.wolfram.com/language/fast-introduction-for-math-students/en/statistics/

http://www.wolfram.com/language/fast-introduction-for-math-students/en/data-plots-and-best-fit-curves/

[17] Podobné jako [12]:
http://www.wolfram.com/language/elementary-introduction/2nd-ed/?source=desktop

[18] NMFM301 Statistika pro finanční matematiku
http://www.karlin.mff.cuni.cz/~pesta/NMFM301/statistika_fm.pdf

[19] NMFM301 Statistika pro finanční matematiku Cvičení 
http://www.karlin.mff.cuni.cz/~maciak/nmfm301_1819.php

2018-10-08 Cvičení 2

FAQ

Diophantine equations and Fermats’s Last Theorem

Reference na příbuzné problémy:

http://mathworld.wolfram.com/FermatsLastTheorem.html

Tam je mj. velmi zajímavý poslední odstavec týkající se řešení pro n=12, pokud nepočítáme v přesnější aritmetice, tak ta rovnost platí. Vyskytuje se v seriálu Simpsons, vím, že bylo konsultováno s Wolfram Research.

Při pokusu o řešení odpovídající rovnice pomocí systému Mathematica je patrné, že v systému jsou zabudovány současné poznatky o konkrétních problémech, takže Mathematica to najde v nějaké databázi, neřešitelné problémy se pak ani nesnaží řešit. Viz dále.

Solve[(a^n)+(b^n)==(c^n),{a,b,c,n},Integers]

solveFermat=Reduce[(a^n)+(b^n)==(c^n),{a,b,c,n},Integers]

Snaha najít aspoň nějaké řešení:

findInstance=FindInstance[(a^n)+(b^n)==(c^n),{a,b,c,n},Integers]

Pro specifikované hodnoty n však nějaké výsledky dostaneme:

(* n=1, 10 řešení *)(findInstance1=FindInstance[((a^n)+(b^n)==(c^n))/.n->1,{a,b,c},Integers,10])//TableForm

(* n=2, 10 řešení *)(findInstance2=FindInstance[((a^n)+(b^n)==(c^n))/.n->2,{a,b,c},Integers,10])//TableForm

Pro n=3 připojíme podmínky a≠0,b≠0,c≠0, jinak FindInstance najde řešení triviální. V tomto případě je řešením prázdná množina { }.

FindInstance[((a^n)+(b^n)==(c^n))/.n->3,{a,b,c},Integers,2]

findInstance3=FindInstance[And[((a^n)+(b^n)==(c^n))/.n->3,a≠0,b≠0,c≠0],{a,b,c},Integers,2]

Cvičení: Jak napsat předchozí příkaz, máme-li podmínky v seznamu:
{((a^n)+(b^n)==(c^n))/.n->3,a≠0,b≠0,c≠0}

FindInstance[{((a^n)+(b^n)==(c^n))/.n->3,a≠0,b≠0,c≠0},{a,b,c},Integers,2]

Apply[And,{((a^n)+(b^n)==(c^n))/.n->3,a≠0,b≠0,c≠0}]

FindInstance[Apply[And,{((a^n)+(b^n)==(c^n))/.n->3,a≠0,b≠0,c≠0}],{a,b,c},Integers,2]

Podobně pro n=4:

findInstance3=FindInstance[And[((a^n)+(b^n)==(c^n))/.n->4,a≠0,b≠0,c≠0],{a,b,c},Integers,2]

list1={((a^n)+(b^n)==(c^n))/.n->2,a≠0,b≠0,c≠0}

findInstanceGubas=FindInstance[Apply[And,list1],{a,b,c},Integers,2]

Následující příkay je součin!!! Druhý řádek obsahuje mezeru interpretovanou jakou krát!! Za první řádek nutné vložit ENTER!

FindInstance[(((a^n)+(b^n)==(c^n))/.n->2)&&a≠0&&b≠0&&c≠0,{a,b,c},(Integers),2]FindInstance[(a^2)+(b^2)==(c^2)&&a≠0&&b≠0&&c≠0,{a,b,c},(Integers),2]
%//FullForm

FindInstance[(((a^n)+(b^n)==(c^n))/.n->2)&&a≠0&&b≠0&&c≠0,{a,b,c},(Integers),2]
FindInstance[(a^2)+(b^2)==(c^2)&&a≠0&&b≠0&&c≠0,{a,b,c},(Integers),2]

(((a^n)+(b^n)==(c^n))/.n->2)//FullForm
(a^2)+(b^2)==(c^2)//FullForm

FindInstance[(a^2)+(b^2)==(c^2)&&a≠0&&b≠0&&c≠0,{a,b,c},(Integers),2, RandomSeeding->1239]

2018-10-29 Cvičení 5

FAQ

Elementární výpočty

Jednoduché funkce

Compute 7+6+5 using the function Plus.

Compute 5×4×3×2 using Times.

Spočítejte maximum ze sta náhodných celých čísel mezi 1 a 100.

SeedRandom[29]
Max[RandomInteger[{1,100},100]]

Generování náhodných čísel (smíšená kongruenční metoda)

(Subscript[X, n+1])≡(a,(Subscript[X, n])+c)mod m

(Subscript[X, 0]),a,c,m vhodně zvolená přirozená čísla (c může být i nula, pak se jedná o multiplikativní kongruenční generátor).

a=134775813;m=(2^32);c=1;

xn1[xn_]:=Mod[a xn+c,m]

{xn1[0],xn1[xn1[0]],xn1[xn1[xn1[0]]]}

NestList[xn1,0,3]

Maximální možné generované číslo je m-1, tj. (2^32)-1:

(2^32)-1

nahodnacisla=NestList[xn1,0,1000000];

nahodnacisla

Max[nahodnacisla]

SeedRandom[29]
Max[RandomInteger[{0,(2^32)-1},1000000]]

2018-11-05 Cvičení 6

FAQ

Náhodný výběr z množiny prvků

seznam1=Partition[Range[20],2]

S navracením (s opakováním)

SeedRandom[13]
RandomChoice[seznam1,5]

Bez navracení (bez opakování)

SeedRandom[13]
RandomSample[seznam1,5]
(* V případě RandomSample[množina,Length[množina]] dostaneme náhodnou permutaci *)

Grafické znázorňování

Grafické  znázornění nejprve pro Bi(n,p)

Clear[rozdeleniBi]
rozdeleniBi[n_,p_]:=Table[{k,PDF[BinomialDistribution[n,p],k]},{k,0,n}]
ListPlot[rozdeleniBi[5,1/10],PlotStyle->{Red,Thick}]

Pospojování bodů

ListLinePlot[rozdeleniBi[5,1/10],PlotStyle->{Red,Thick}]

ListPlot[rozdeleniBi[5,1/10],PlotStyle->{Red,Thick},Joined->True]

ListPlot[rozdeleniBi[40,1/10],PlotStyle->{Red,Thick}]

ListLinePlot[rozdeleniBi[40,1/10],PlotStyle->{Red,Thick}]

ListLinePlot[rozdeleniBi[40,1/10],PlotStyle->{Red,Thick},PlotRange->All]

ListLinePlot[rozdeleniBi[40,1/10],PlotStyle->{Red,Thick},PlotRange->{0,0.000001}]

ListLogPlot[rozdeleniBi[40,1/10],PlotStyle->{Red,Thick},Joined->True]

Histogram

SeedRandom[13]
nahodnyvyberBi40jednadesetina10000=RandomVariate[BinomialDistribution[40,1/10],10000];
Max[%]
Text@Row@{"Průměr, rozptyl, směr. odchylka:   ",Through[{Mean,Variance,StandardDeviation}[%%]]//N}
Histogram[%%%]

Through[{Mean,Variance,StandardDeviation}[nahodnyvyberBi40jednadesetina10000]]//N

Empirická hustota a porovnání s normálním rozdělením
Histogram ⟶ F1 ⟶ Ctrl+F⟶pdf ⟶ Find Next ⟶ F3 ⟶ F3 ...

p1=Histogram[nahodnyvyberBi40jednadesetina10000,Automatic,"PDF"]

{Mean[BinomialDistribution[n,p]],StandardDeviation[BinomialDistribution[n,p]]}
%/.{n->40,p->(1/10)}

p2=Plot[PDF[NormalDistribution[4,3 (Sqrt[(2/5)])],x],{x,0,11},PlotStyle->{Red,Thick}]

Show[p1,p2]

***Různé

?*Genera*

WordList

2018-11-12 Cvičení 7

FAQ

Metody Monte Carlo (Simulační metody)

Princip

MC metody: odvětví experimentální matematiky, které se zabývá experimenty s náhodnými čísly. Numerická technika provádění hromadných experimentů s modely pomocí počítače umožňující generovat velké množství náhodných výběrů (umělých náhodných výběrů). Většinou jsou efektivnější v situacích, kdy sama podstata problému je náhodná. Při řešení deterministických problémů tomu jsou většinou úspěšnější klasické metody numerické matematiky.

Generování náhodných čísel

*Náhodná čísla z rovnoměrného rozdělení

Náhodná čísla by měly mít vlastnosti náhodného výběru, nezávislé stejně rozdělené náhodné veličiny (iid, independent identically distributed r. v.))

Pro úspěšné použití metod MC je nutné mít kvalitní generátory náhodných čísel.

V současnosti nejúspěšnější jsou generátory posloupností náhodných čísel založeny na rafinovaných deterministických algoritmech, které produkují tzv. pseudonáhodná čísla. Při znalosti vstupních parametrů je možné příslušnou posloupnost přesně zrekonstruovat.  Ukazuje se, že generátory založené na deterministických algoritmech jsou podstatně lepší než ty, které jsou založené na “opravdové náhodě” (tahání očíslovaných lístků z urny, výběr očíslovaných pingpongových míčků, ruleta, házení kostkou, šum elektronek, házení mincí aj.)

Na generátorech náhodných resp. pseudonáhodných čísel se v zásadě požaduje: dlouhá perioda (nesmí se začít opakovat příliš brzy, vyžaduje se řádově (10^8) členů posloupnosti bez opakování), nezávislost nebo nekorelovanost, konečně to, aby se řídila požadovaným rozdělením, sdruženým rozdělením dvojic, trojic, atd. Teorie čísel (v poslední době teorie celulárních automatů aj.) poskytuje fundamentální podklad pro konstrukci generátorů, nicméně je stejně dobrá praxe podrobit nový generátor baterii statistických testů. Takové testy by měly být specielně zaměřeny na to, k jakému účelu mají být náhodná čísla použita.

⚠ V drtivé většině případů se primárně generují náhodná čísla z rovnoměrného rozdělení a ta se vhodnými postupy transformují na náhodná čísla z požadovaného rozdělení. Na počítači je přirozeně možné generovat náhodná čísla z diskrétního rovnoměrného rozdělení, řekněme R(0,...,M-1), kde M je “velké číslo”, např. (2^31). Po normování je možné takto generovaná náhodná čísla považovat za náhodná čísla ze spojitého rovnoměrného rozdělení R(0,1).

Podrobnosti mj. Hurt, J. (1982): Simulační metody. Skripta. SPN, Praha.

Hyperlink[Button["Random Number Generation"],"paclet:tutorial/RandomNumberGeneration"]

**Kongruenční generátory

Kongruence zde i dále znamena zbytek po deleni. Generuji se cisla 0,1,...,M-1. Generator typu (Subscript[X, n+1])≡a (Subscript[X, n]) mod M se nazyva multiplikativni, 
(Subscript[X, n+1])≡a (Subscript[X, n])+b  mod M smisena kongruencni metoda. Při vhodně zvolených konstantách a,b,M mají dlouhou periodu (začnou se opakovat až za “dlouho”), někdy dokonce maximální možnou, tj. M.

Pro ilustraci uzijeme zabudovanou funkci Mod:

(illustrationModulo=Prepend[Flatten[Table[{x,y,Mod[x,y]},{x,2,5},{y,2,x}],1],{" x "," M ","Mod[x,M]"}])//Grid

Následující výstup je mírně modifikovaný výstup s využitím příkladu z Help k funkci Grid:

Text@Grid[illustrationModulo,Background->{None,{Lighter[Yellow,.9],{White,Lighter[Blend[{Blue,Green}],.8]}}},Dividers->{{Darker[Gray,.6],{Lighter[Gray,.5]},Darker[Gray,.6]},{Darker[Gray,.6],Darker[Gray,.6],{False},Darker[Gray,.6]}},Alignment->{{Center,Center,Center}},ItemSize->{{2,2,5}},Frame->Darker[Gray,.6],ItemStyle->14,Spacings->{Automatic,.8}]

**Příklad 1 (Perioda)

Generovani posloupnosti 17 nahodnych cisel     (Subscript[X, n+1])≡5(Subscript[X, n])+1 mod 16. Perioda je maximalni rovna 16.

gen16=NestList[Mod[5 # +1,16]&,0,16]

**Příklad 2 (Perioda) 
(Vetsi problem, rozdeleni dvojic. Stejny generator (Subscript[X, n+1])=5(Subscript[X, n])+1 mod 16.)

Hypoteticky má každá dvojice pst 1/(16^2)=1/256. Vzhledem k rekurenci však není možné generovat více než 16-1 dvojic. Generovaná čísla leží na dvou systémech přímek.

Zkoumáme rozdělení dvojic ((Subscript[X, n]),(Subscript[X, n+1])).

gen16shift=RotateLeft[gen16]

ListPlot[Transpose[{gen16,gen16shift}],PlotStyle-> PointSize[0.02]]

**Příklad 3 (Generator kapesniho kalkulatoru TI-59)

Nasledujici funkce generuje nahodna cisla smisenou kongruencni metodou delky n+1. Perioda by mela byt 199017.

Clear[gen59];
(* Délka generované posloupnosti *)
gen59[n_]:=NestList[Mod[24298 # +9991,199017]&,0,n]

gen59[10]

Počet různých generovaných náhodných čísel v posloupnosti délky 200000:

{Length[gen59[200000]],Length[Union[gen59[200000]]]}

Po sobe jdouci dvojice

Clear[gen592];
gen592[n_]:=Module[{shift,g},
shift=RotateLeft[g=gen59[n]];
Transpose[{g,shift}]
]

gen592[50000]//ListPlot

V detailu:

ListPlot[gen59compl=Select[gen592[199017],(#[[1]]<1000)&]]

ListPlot[gen59compl=Select[gen592[199017],(#[[1]]<2000∧#[[2]]<2000)&]]

**Příklad 4 (Velmi špatný kongruenční generátor IBM ze šedesátých let RANDU)

http://demonstrations.wolfram.com/PoorStatisticalQualitiesForTheRANDURandomNumberGenerator/

*Náhodná čísla z jiných rozdělení

Universální metoda

Věta (Metoda inverzni transformace)

Nechť náhodná veličina X má distribuční funkci F a kvantilovou funkci (F^(-1))(y)=inf{x|F(x)≥y}, 0<y<1. Nechť U∼R(0,1). Potom (F^(-1))(U) má rozdělení s ditribuční funkcí F, tj., má rozdělení jako X. Platí pro libovolné rozdělení, spojité, diskrétní a libovolné jiné!!!

Priklad (exponencialni rozdeleni se stredni hodnotou 1)

dFexp[x_]:=1-(E^(-x)) (* Pozor, to není obyčejné e, ale E! *)
(* Může být též psáno jako E^(-x). *)
Plot[dFexp[x],{x,0,10}]

sol1=Solve[dFexp[x]==y,x,Reals]//Simplify[#,Assumptions->0<y<1]&
invF[y_]=First[x/.sol1]
invF[z]

alternativne

InverseCDF[ExponentialDistribution[1],y]
Quantile[ExponentialDistribution[1],y]

Attributes[Log]

nahexp=invF[RandomVariate[UniformDistribution[{0,1}],{100000}]];
{Mean[nahexp],Variance[nahexp],Max[nahexp]}
Show[{Histogram[nahexp,100,"PDF"],Plot[PDF[ExponentialDistribution[1],x],{x,0,5},PlotStyle->{Red,Thick}]}]
(* *** Histogram[nahexp,100,"PDF",Epilog->First@Plot[PDF[ExponentialDistribution[1],x],{x,0,5},PlotStyle->{Red,Thick}]]
*)

(* ** *)DistributionFitTest[nahexp,ExponentialDistribution[1],"HypothesisTestData"]["TestDataTable",All]

RandomVariate[ExponentialDistribution[λ]]

Box-Mullerova transformace

http://en.wikipedia.org/wiki/Box% E2 %80 %93 Muller_transform

***Příklad 5 Generovani nahodnych cisel z velmi zvlastnich rozdeleni (prevzato z MathGroup Archive: December 1998 [00243], Paul Abbott)

Pripominka: Metoda inversni transformace je zalozena na znalosti kvantilove funkce.

Clear[pdf];
pdf[A_,s_][r_]:=(r/s^2*BesselI[0,(r*A)/s^2])/E^((A^2+r^2)/(2*s^2));

pdf[1,1][x]

After some 93 (version 7) 25 (version 8, 9) seconds we get that pdf is really a PDF.

Integrate[pdf[A,s][r],{r,0,∞},Assumptions->Re[(s^2)]>0&&(Re[(A/(s^2))]>0||(Re[(A/(s^2))]==0&&Im[(A/(s^2))]>0))]//Simplify

Print["Mean =  ",e11=Integrate[x pdf[1,1][x],{x,0,∞}]//N,"   Variance =   ",Integrate[(x^2) pdf[1,1][x],{x,0,∞}]-(e11^2)//N]

Plot[pdf[1,1][x],{x,0,5}]

Remark: Use NDSolve (and dynamic programming) instead of NIntegrate to determine the CDF:

Clear[cdf,r];
cdf[A_,s_]:=cdf[A,s]=Module[{y},
First[y/.NDSolve[{y'[r]==pdf[A,s][r],y[0]==0},y,{r,0,10}]]]

Porovnejte s vypoctem CDF pomoci integrace.

From the plot,

cdfplot=Plot[cdf[1,1][x],{x,0,10},PlotRange->All]

Short[InputForm[cdfplot],10]

extract the plotpoints (which takes advantage of the adaptive code built into Plot):

Clear[x];
pts=Cases[cdfplot,Line[{x__}]->x,Infinity];

ListPlot[pts]

and then numerically compute the inverse function using interpolation:

Map[Reverse,pts]//ListPlot

inverseCDF[1,1]=Interpolation[Map[Reverse,pts]];

Plot[inverseCDF[1,1][x],{x,0.0000001,0.9999999999}]

Table[inverseCDF[1,1][Random[]],{10}]

Uziti vymozenosti verze 8, 9, 10

𝒟velmizvlastni=ProbabilityDistribution[(E^((1/2) (-1-(x^2)))) x BesselI[0,x],{x,0,∞}]

{Mean[𝒟velmizvlastni],Variance[𝒟velmizvlastni]}//N

nahodnacislazvelmizvlastnihorozdeleni=RandomVariate[𝒟velmizvlastni,{1000}];
Histogram[%]

Porovnání binomického a Poissonova rozdělení

Motivace (Aproximace binomického rozdělení rozdělením Poissonovým):

Nechť n,(Subscript[p, n]) ⟶ λ>0 pro n ⟶ ∞. Potom Bi(n,(Subscript[p, n])) ⟶ Po(λ) v distribuci. Jinými slovy, Pro “velká” n a “malá” p je Bi(n,p) přibližně Poissonovo s parametrem λ.

https://www.youtube.com/watch?v=ceOwlHnVCqo

n=100;p=0.01;np=n p;
Show[{ListPlot[Table[{k,PDF[BinomialDistribution[n,p],k]},{k,0,n/10}],PlotStyle->{Red,PointSize[0.02]},PlotRange->All,PlotLabel->"Probabilities"],ListLinePlot[Table[{k,PDF[PoissonDistribution[1],k]},{k,0,n/10}],PlotStyle->{Thick},PlotRange->All,PlotLabel->"Probabilities"]}]

nebo

Clear[comparisonBinomialPoisson]
comparisonBinomialPoisson[n_,p_]:=
Module[{max,np,pD,bD},
np=n p;
pD=PoissonDistribution[np];
bD=BinomialDistribution[n,p];
(* Upper bound for the plots: *)
max=Ceiling[Max[Quantile[pD,.999],Quantile[bD,.999]]];
GraphicsColumn@{Show[ListPlot[Table[{k,PDF[pD,k]},{k,0,max}],PlotStyle->{Red,PointSize[0.02]}],ListLinePlot[Table[{k,PDF[bD,k]},{k,0,max}],PlotStyle->Thick],PlotLabel->"Probabilities of Poisson (red points) and binomial (blue line) distribution",PlotRange->All],ListLinePlot[Table[{k,PDF[bD,k]-PDF[pD,k]},{k,0,max}],PlotStyle->Thick, PlotLabel->"Differencies between binomial and Poisson probabilities",PlotRange->All,PlotTheme->"Scientific"]}
]

comparisonBinomialPoisson[100,.01]

comparisonBinomialPoisson[100,.5]

2018-11-19 Cvičení 8

FAQ

Poznámky k 2018-11-12 Cvičení 7

Porovnání binomického a Poissonova rozdělení

Zde jsou jen některá mnou objevená nedopatření. U zkoušky nemusíte používat diakritiku (zvláště, jestli Vás to zdržuje jako mne - specielně přepínají-li se klávesnice QWERTY a QWERTZ).

Zadání připouští celou řadu řešení. Nejjednodušší je grafické porovnání rozdělení pstí obou rozdělení a graf rozdílu odpovídajících pstí, viz zde.

1

Cviceni 5

Tato buňka má Style Input. Očividně se má jednat o odstavci, sekci, kapitole ... . Toho se dosáhne: Menu ⟶ Format ⟶  Subsubsection nebo přímo Alt+6:

Cviceni 5

2

Table[PDF[BinomialDistribution[50,p],k],{p,{0.1}}]

DiscretePlot[Evaluate[%],{k,1,50}]

Tabulka není tabulka. Navíc pro DiscretePlot potřebujeme jenom hodnoty, obejdeme se bez Table. Odvolání na předchozí výstup je v jiné buňce - nebezpečná konstrukce. Možná konstrukce:

n=50;p=0.1;
DiscretePlot[PDF[BinomialDistribution[n,p],k],{k,0,n}]

Takový graf ale toho mnoho nevypovídá, vyšší hodnoty jsou velmi málo pravděpodobné. Možná vylepšení:

n=50;p=0.1;
GraphicsRow@{DiscretePlot[PDF[BinomialDistribution[n,p],k],{k,0,n/5},PlotLabel->"Bi(n,p)",PlotStyle->Thick,AxesLabel->{k,Probability}],ListLogPlot[Table[{k,PDF[BinomialDistribution[n,p],k]},{k,0,n}],PlotLabel->"Log(Bi(n,p))",PlotStyle->{Thick,PointSize[0.015]},AxesLabel->{k,Probability}]}

3

Clear[n,p]
rozdil[n_,p_]=Table[PDF[PoissonDistribution[n p],k]-PDF[BinomialDistribution[n,p],k],{k,0,n}]

Definice funkce rozdil musí být s odloženým přiřazením (:=):

Clear[rozdilDelayed]
rozdilDelayed[n_,p_]:=Table[PDF[PoissonDistribution[n p],k]-PDF[BinomialDistribution[n,p],k],{k,0,n}]
rozdilDelayed[10,0.5]
%//ListPlot

Ještě lepší je vytvářet tabulku i s hodnotami:

Clear[rozdilDelayedTabulka]
rozdilDelayedTabulka[n_,p_]:=Table[{k,PDF[PoissonDistribution[n p],k]-PDF[BinomialDistribution[n,p],k]},{k,0,n}]
TableForm[rozdilDelayedTabulka[10,0.5]]
%//ListPlot
(* ale *)
ListPlot[TableForm[rozdilDelayedTabulka[10,0.5]]]

4

Jeden student si vzpomněl na to, že pro velká n a malá p je binomické rozdělení přibližně Poissonovo s parametrem n,p. V technickém provedení se opět vyskytuje nebezpečná konstrukce, odkazy na předchozí výstup jsou v různých buňkách.

n,(Subscript[p, n]) ⟶ λ>0 ⟹ Bi(n,(Subscript[p, n])) ⟶ Po(λ) pro n ⟶ ∞.

V technickém provedení se opět vyskytuje nebezpečná konstrukce, odkazy na předchozí výstup jsou v různých buňkách.

Limit[PDF[BinomialDistribution[n,μ/n],k],n->∞]

PDF[PoissonDistribution[μ],k]

FullSimplify[%-%%,k∈Integers]

Správně má být:

Limit[PDF[BinomialDistribution[n,μ/n],k],n->∞]
PDF[PoissonDistribution[μ],k]
FullSimplify[%-%%,k∈Integers]

5

n1=100;p1=0.01;Bi[n_,p_]:=Table[{k,PDF[BinomialDistribution[n,p],k]},{k,0,n}]
p1=Plot[PDF[Bi[n1,p1],x],{x,0,11}]

Chybná konstrukce. Plot[f[x],{x,(Subscript[x, min]),(Subscript[x, max])}] zobrazuje funkci na odpovídajícím intervalu. Grafické znázornění tabulky pomocí funkce ListPlot a jejími variacemi. Navíc se použilo označení grafu p1, což dříve označovalo pst.

Bi[n1,p1]

Správný začátek může být např.

ListPlot[Bi[100,0.01]]

a postupným vylepšováním

ListPlot[Take[Bi[100,0.01],10]]

6

Všechny zabudované funkce začínají velkým písmenem. Zde má být Show.

show[p1,p2]

7

dFexp[x_] := 1-(e^(-x))
Plot[dFexp[x], {x,0,10}]

“e” není známý symbol. Má to být základ přirozených logaritmů. Správně e (z palety nebo Escee Esc) nebo Exp[-x] nebo E.

dFexp[x]/.x->1
?e

Clear[dFexp]
dFexp[x_] := 1-(E^(-x))
GraphicsRow@{Plot[dFexp[x], {x,0,10}],Plot[dFexp[x], {x,0,10},PlotRange->{0,1}]}

sol1 = Solve[dFexp[x] =y,x, Reals] // Simplify[#,Assumptions -> 0<y<1]&
inv[y_]=First[x/.sol1]
invF[z]

Argumentem Solve není rovnice, ale přiřazovací příkaz. Má tam být  dFexp[x] ==y. Navíc se tím poničí definice funkce dFexp:

?dFexp

Plot[dFexp[x],{x,0,10}]

Další chyby:

IniformDistribution[{0,1}] (* Uniform ... *)

plot[ExponentialDistribution[x],{x,0,1}] 
(* plot s malým p, ExponentialDistribution[x] je exp. rozdělení s parametrem x, a ne žádná jeho hustota, distr. funkce atd., což by byly přípustné argumenty funkce Plot *)

8

binomicke[k_]:=Binomial[100,k]*0.01^k*(1-0.01)^(100-k)

poissonovo[k_]:=(1/Factorial[k])*Exp[-1]

TableForm[Table[{u,binomicke[u],N[poissonovo[u]]},{u,0,6}]]

Je pohodlnější sdružit 3 příkazy do jedné buňky a místo vlastních definicí použít zabudované funkce:

(tabBinPoi=Table[{k,PDF[BinomialDistribution[100,0.01],k],PDF[PoissonDistribution[1.],k]},{k,0,6}])//TableForm

9

Mean[BinomialDistribution[100,0.01]]

Mean[PoissonDistribution[1]]

Variance[BinomialDistribution[100,0.01]]

Variance[BinomialDistribution[100,0.01]]

Sdružit 3 příkazy do jedné buňky nebo do seznamu:

Mean[BinomialDistribution[100,0.01]]
Mean[PoissonDistribution[1]]
Variance[BinomialDistribution[100,0.01]]
Variance[PoissonDistribution[1]](* Variance ? *)

{Mean[BinomialDistribution[100,0.01]],
Mean[PoissonDistribution[1]],
Variance[BinomialDistribution[100,0.01]],
Mean[PoissonDistribution[1]],
Variance[PoissonDistribution[1]]}

SmoothHistogram

SmoothHistogram (obecně *Histogram*) fungují pouze pro náhodné výběry. Navíc je má smysl použít pouze pro náhodné výběry ze spojitých rozdělení.

Birozdeleni[n_,p_]:=Table[{k, PDF[BinomialDistribution[n,p],k]},{k,0,n}]
Porozdeleni[x_]:=Table[{k, PDF[PoissonDistribution[x],k]},{k,0,100}]
SmoothHistogram[{Birozdelenie[100,0.01],Porozdelenie[1]}]

SmoothHistogram pro náhodný výběr z diskrétního rozdělení nedává smysl:

SmoothHistogram[{SeedRandom[13];RandomVariate[BinomialDistribution[100,.01],10000],SeedRandom[13];RandomVariate[PoissonDistribution[1],10000]}]

SmoothHistogram[SeedRandom[13];RandomVariate[BinomialDistribution[100,.01],10000],PlotRange->All]

Histogram a SmoothHistogram pro náhodný výběr ze spojitého rozdělení ve tvaru hustoty

GraphicsRow@{Histogram[SeedRandom[13];RandomVariate[ChiSquareDistribution[2.5],10000],Automatic,"PDF"],SmoothHistogram[SeedRandom[13];RandomVariate[ChiSquareDistribution[2.5],10000],Automatic,"PDF"],Plot[PDF[ChiSquareDistribution[2.5],x],{x,0,20},PlotStyle->Red]}

Statistika 1

Princip testování statistických hypotéz

tabulka, chyba 1. druhu, chyba 2. druhu, síla testu

SmoothHistogram

SmoothHistogram,SmoothKernelDistribution, SmoothHistogram3D, SmoothDensityHistogram

CLV, pravidlo 2σ, pravidlo 3σ

Standardní normální roydělení N(0,1)

distrNormal=NormalDistribution[] 
{Probability[Abs[x]>2.,x\[Distributed]distrNormal],Probability[Abs[x]>3.,x\[Distributed]distrNormal]}

Clear[μ,σ]
{Probability[Abs[x-μ]>2.σ,x\[Distributed]NormalDistribution[μ,σ]],Probability[Abs[x-μ]>3.σ,x\[Distributed]NormalDistribution[μ,σ]]}

Uniform, Exponential, (χ^2), Student

Rozdělení výběrového průměru pro různá rozdělení

Generovat N n

2018-11-26 Cvičení 9

FAQ

NMaximize

Statistika 2

CLV, pravidlo 2σ, pravidlo 3σ

Počítání pravděpodobnosti, že se náhodná veličina liší od své střední hodnoty v absolutní hodnotě o více než 2σ resp. 3σ.

Standardní normální rozdělení N(0,1)

distrNormal=NormalDistribution[] 
{Probability[Abs[x]>2.0,x\[Distributed]distrNormal],Probability[Abs[x]>3.0,x\[Distributed]distrNormal]}

Clear[μ,σ]
{Probability[Abs[x-μ]>2.0×σ,x\[Distributed]NormalDistribution[μ,σ]],Probability[Abs[x-μ]>3.0×σ,x\[Distributed]NormalDistribution[μ,σ]]}

Uniform distribution

Clear[a,b]
distrUniform=UniformDistribution[{a,b}];
{Mean[distrUniform],StandardDeviation[distrUniform]}

{((a+b)/2),((-a+b)/(2 (Sqrt[3])))}/.{a->-1,b->1.}

(χ^2)- distribution

Student t- distribution

Clear[distrStudent]
distrStudent[n_]:=StudentTDistribution[μ,σ,n (* = degrees of freedom *)]

{Mean[distrStudent[n]],StandardDeviation[distrStudent[n]]}

Volte n=1, 5, 10, 50. Výsledek ve formě tabulky

Cauchy distribution (≡ t(1))

PDF[CauchyDistribution[],x]==PDF[StudentTDistribution[1],x]
{PDF[CauchyDistribution[],x],PDF[StudentTDistribution[1],x]}

Exponential distribution

Rozdělení výběrového průměru pro různá rozdělení

Je-li (OverBar[X]) výběrový průměr náhodného výběru z libovolného rozdělení s konečným rozptylem (σ^2) a střední hodnotou μ, pak E,(OverBar[X])=μ a var (OverBar[X])=(σ^2)/n.

Úloha: Generovat N náhodných výběrů o rozsahu n, z každého výběru spočítat výběrový průměr. Dostaneme tedy N hodnot, vlastně náhodný výběr výběrových průměrů. Graficky znázornit empirickou hustotu (histogram), “vyhlazený” histogram (SmoothHistogram) a do stejného grafu odpovídající hustotu normálního rozdělení, pozor na stejné parametry! Použijte rozdělení z předchozího odstavce. Efektivní postup spočívá ve vytvoření obecné procedury (Module), kde vstup je N, n, rozdělení  a volitelné semínko pro generování náhodných čísel (Optional, ), celkem 4 vstupní parametry.

Normal and Other Distributions

pocetVyberu=1000;n=1000;
SeedRandom[13]
tab10000=Table[xx=RandomVariate[NormalDistribution[],{n}];Mean[xx],{pocetVyberu}];
%//Dimensions
Histogram[%%,Automatic,"PDF"]

pocetVyberu=10000;n=10;
SeedRandom[13]
tab10000=Table[xx=RandomVariate[NormalDistribution[],{n}];Mean[xx],pocetVyberu];
Histogram[%]

Obecná procedura

Clear[rozdeleniPrumeru]
rozdeleniPrumeru[pocetVyberuN_,rozsahVyberun_,rozdeleni_,seminko_:13]:=
Module[{n=rozsahVyberun,nahodneVybery,xx,histogram,smoothHistogram,normalHustota,mean,sd},
SeedRandom[seminko];
mean=Mean[rozdeleni];sd=StandardDeviation[rozdeleni];
nahodneVybery=Table[xx=RandomVariate[rozdeleni,{n}];Mean[xx],{pocetVyberuN}];
histogram=Histogram[nahodneVybery,Automatic,"PDF",PlotLabel->rozdeleni];
smoothHistogram=SmoothHistogram[nahodneVybery,PlotLabel->rozdeleni,PlotStyle->Thick];
normalHustota=Plot[PDF[NormalDistribution[mean,sd/(Sqrt[n])],x],{x,mean-4sd/(Sqrt[n]),mean+4sd/(Sqrt[n])},PlotStyle->{Thick,Red}];
Show[{histogram,smoothHistogram,normalHustota}]
]

rozdeleniPrumeru[10000,3,NormalDistribution[]]

GraphicsRow@{rozdeleniPrumeru[10000,10,StudentTDistribution[3]],rozdeleniPrumeru[10000,10,StudentTDistribution[30]]}

GraphicsRow@{rozdeleniPrumeru[2000,3,ExponentialDistribution[1]],rozdeleniPrumeru[2000,30,ExponentialDistribution[1]]}

2018-12-03 Cvičení 10

FAQ

Poznámky k 2018-11-26 Cvičení 9

Rozdělení výběrového průměru pro různá rozdělení

Je-li (OverBar[X]) výběrový průměr náhodného výběru z libovolného rozdělení s konečným rozptylem (σ^2) a střední hodnotou μ, pak E,(OverBar[X])=μ a var (OverBar[X])=(σ^2)/n.

Úloha: Generovat N náhodných výběrů o rozsahu n, z každého výběru spočítat výběrový průměr. Dostaneme tedy N hodnot, vlastně náhodný výběr výběrových průměrů. Graficky znázornit empirickou hustotu (histogram), “vyhlazený” histogram (SmoothHistogram) a do stejného grafu odpovídající hustotu normálního rozdělení, pozor na stejné parametry! Použijte rozdělení z předchozího odstavce. Efektivní postup spočívá ve vytvoření obecné procedury (Module), kde vstup je N, n, rozdělení  a volitelné semínko pro generování náhodných čísel (Optional, ), celkem 4 vstupní parametry.

Obecná procedura

Clear[rozdeleniPrumeru]
rozdeleniPrumeru[pocetVyberuN_,rozsahVyberun_,rozdeleni_,seminko_:13]:=
Module[{n=rozsahVyberun,nahodneVybery,xx,histogram,smoothHistogram,normalHustota,mean,sd},
SeedRandom[seminko];
mean=Mean[rozdeleni];sd=StandardDeviation[rozdeleni];
nahodneVybery=Table[xx=RandomVariate[rozdeleni,{n}];Mean[xx],{pocetVyberuN}];
histogram=Histogram[nahodneVybery,Automatic,"PDF",PlotLabel->rozdeleni];
smoothHistogram=SmoothHistogram[nahodneVybery,PlotLabel->rozdeleni,PlotStyle->Thick];
normalHustota=Plot[PDF[NormalDistribution[mean,sd/(Sqrt[n])],x],{x,mean-4sd/(Sqrt[n]),mean+4sd/(Sqrt[n])},PlotStyle->{Thick,Red}];
Show[histogram,smoothHistogram,normalHustota]
]

rozdeleniPrumeru[10000,3,NormalDistribution[]]

GraphicsRow@{rozdeleniPrumeru[10000,10,StudentTDistribution[3]],rozdeleniPrumeru[10000,10,StudentTDistribution[30]]}

GraphicsRow@{rozdeleniPrumeru[2000,3,ExponentialDistribution[1]],rozdeleniPrumeru[2000,30,ExponentialDistribution[1]]}

rozdeleniPrumeru[2000,80,ExponentialDistribution[1],333]

1

Ve slozenych prikazech se nesmi zapominat na strednik

Clear[nahvyber2]
nahvyber2[NNN_,nnn_]:=Module[{xxx,tab10000},
tab10000=Table[xxx=RandomVariate[NormalDistribution[],{nnn}]Mean[xxx],{NNN}];
{Histogram[tab10000]}]
nahvyber2[1000,10]

Clear[nahvyber2]
nahvyber2[NNN_,nnn_]:=Module[{xxx,tab10000},
tab10000=Table[xxx=RandomVariate[NormalDistribution[],{nnn}] ;Mean[xxx],{NNN}];
{Histogram[tab10000]}]
nahvyber2[1000,10]

2

Tady se nasobi, vyjde nesmysl; prikazy musi vlevo licovat:

SeedRandom[seed]
data=Table[RandomVariate[dist,rozsah],n];

StandardDeviation jej jiz smerodatna odchylka. Sqrt je navic:

Plot[PDF[NormalDistribution[Mean[dist],Sqrt[StandardDeviation[dist]/n]],x]

3

Iterator nedava smysl:

Module[{n,N},{n,10,100,10}; N=100000/n ]

4

Tady zustavaji predchozi hodnoty μ, σ:

distrStudent[n_]:=StudentTDistribution[μ,σ,n(*degrees of freedom*)]
distrStudent[1]

Kdyz uz s parametry, tak korektne:

distrStudent[n_,μ_,σ_]:=...

5

Promenna distrUniform byla vyse definovana bez argumentu (distrUniform = UniformDistribution[]). Vyvolani s argumentem vede k chybe:

pocetVyberu =10000;
n=10;
SeedRandom[13]
tab10000=Table[xx=RandomVariate[distrUniform[{a,b}],{n}];
Mean[xx], {pocetVyberu}];
Histogram[%]

Misto

Table[xx=RandomVariate[distrUniform[{a,b}]],{n}];

stacilo pouzit a parametry a,b drive specifikovat:

Table[xx=RandomVariate[UniformDistribution[{a,b}]],{n}];

6

vybery - chybi strednik za ... {pocet}]

Histogram[vybery,”Probability”] ... ”Probability” muze byt az treti argument, radeji “PDF”. Zde a ve vsech dalsich funkcich se nespecifikovane argumenty musi zadat jako Automatic. U voleb (Options) je to jedno.  Zde melo byt Histogram[vybery,Automatic,”PDF”]

Misto Variance ma byt StandardDeviation (viz definice NormalDistribution) :
NormalDistribution[Mean[vybery],Variance[vybery]]

pocet=10000;n=10;SeedRandom[5];rozdeleni=NormalDistribution[];
vybery=Table [x=RandomVariate[rozdeleni,{n}];Mean[x],{pocet}]
Show[Histogram[vybery,"Probability"],SmoothHistogram[vybery],Plot[PDF[NormalDistribution[Mean[vybery],Variance[vybery]]][x],{x,-3,3}]]

7

Pokus1[N0_,n0_,dist0_,rand0_]:= Module[{N=N0,n=n0,dist=dist0,rand=rand0},
											SeedRandom[rand];
											tab =Table[xlist=RandomVariate[dist,n];Mean[xlist],N];
								Show[SmoothHistogram[tab],					Histogram[tab]]]

Pokus1[10000,10,NormalDistribution[],13]

Nazvy uzivatelskych  promennych a procedur by mely zacinat malym pismenem, nazev procedury pro jistotu predem vycisten:
Clear[pokus1]

SeedRandom se ve slozenem prikaze modulu muze primo pouzit se skutecnou hodnotou parametru: SeedRandom[rand0].

Podivne formatovani.

Zbytecne prirazovani hodnot lokalnich promennych argumentu procedury. Zabudovana funkce N se neda v modulu pouzit. Chybou ale je, ze jako lokalni promenne nejsou deklarovany tab, xlist. Ma tedy byt

Module[{tab,xlist,...}, telo procedury]

8

Clear[PocetVyberu];
Clear[n];
Clear[rozd];
Clear[sem];
generator[PocetVyberu_,n_,rozd_,sem_]:=Module[{PocetVyberu,n,rozd,sem},SeedRandom[sem];Table[Mean[RandomVariate[rozd,{n}]],{PocetVyberu}]]

tabnew=generator[10000,10,NormalDistribution[],13]

Argumenty modulu neni mozne deklarovat jako lokalni promenne modulu. Nemusi se vynulovat vzory, ktere figuruji v definici funkce generator, ale je bezpecne vynulovat tu definovanou funkci, nejlepe ClearAll:

ClearAll[generator];
generator[pocetVyberu_,n_,rozd_,sem_]:=Module[{pouze promenne, ktere se skutecne jako lokalni promenne v nasledujicim slozenem prikaze vyskytuji},
SeedRandom[sem];Table[Mean[RandomVariate[rozd,{n}]],{pocetVyberu}]]

9

Table[PDF[StudentTDistribution[n (*= stupne volnosti*)]], {n,{1,5,10,50}}]

Lepe pouzit i argument hustoty, teba x:

Clear[x]
Table[PDF[StudentTDistribution[n (*= stupne volnosti*)],x], {n,{1,5,10,50}}]

Jak to vypada pro jina rozdeleni nez normalni. N=? n=? z CLV, když existuje druhý moment, tak ... ma normalni rozdeleni 0,1

Toto je ve formatu Input, ma to byt ve formatu Text (z Menu nebo ALT+7):

Jak to vypada pro jina rozdeleni nez normalni. N = ? n=? z CLV, když existuje druhý moment, tak ... ma normalni rozdeleni 0,1

clear[n];clear[N];
SeedRandom[13]
tab[n_]=Table[xx= RandomVariate[NormalDistribution[],{n}];
Mean[xx], {N}];

Histogram[%]

clear[N], svit modre symbolu clear neni prirazena zadna hodnota. Vsechny zabudovane funkce zacinaji velkym pismenem!!!
N je zabudovana funkce, takze

Clear[N]

Dale misto okamziteho prirazeni ma byt uzito odlozeneho prirazeni a funkci definovat i pro ostatni argumenty (nejenom n).

10

NN= 2000; n=50;

tabl2000=Table[xx=RandomVariate[NormalDistribution[],{n}];
Mean[xx],{NN}];

sh=SmoothHistogram[%]

Vsechno by melo byt ve stejne vstupni bunce, uz s ohledem na referenci na posledni vystup pomocoi %:

NN= 2000; n=50;
tabl2000=Table[xx=RandomVariate[NormalDistribution[],{n}];
Mean[xx],{NN}];
sh=SmoothHistogram[%]

Dalsi je velmi casta chyba, ze se zapomene na argument u funkci typu hustota, distribucni funkce, ... . Dalsi obrazek tudiz nezobrazi hustotu N(0,1):

Show[sh,Plot[PDF[NormalDistribution[]],{x,-1,1}]]

Lepe s nekterymi vylepsenimi:

Show[sh,Plot[PDF[NormalDistribution[],x],{x,-1,1},PlotStyle->{Thick,Red}]]

Proc to neni to, co cekame? Protoze vyberovy prumer ma rozptyl (σ^2)/n a ne (σ^2)! Musime tedy do grafu zahrnout srovnavatelnou hustotu N(0,1/(Sqrt[n])), v nasem pripade je n=50. Dohromady:

Show[sh,Plot[PDF[NormalDistribution[0,1/(Sqrt[n])],x],{x,-1,1},PlotStyle->{Thick,Red}],PlotRange->All]

11

Pekne graficke upravy, ale chybi to podstatne, uvedomit si, ze smerodatna odchylka vyberoveho prumeru je σ/(Sqrt[n]), a ze histogram musi mit podobu hustoty (treti argument funkce Histogram je “PDF”.
Inicializace SeedRandom musi byt pred kazdym volanim RandomVariate!

Clear[nd,Nd]

Nd=1000; nd=100;
SeedRandom[5]
table1=Table[x=RandomVariate[NormalDistribution[],{nd}];Mean[x],{Nd}];

hist=Histogram[Table[x=RandomVariate[NormalDistribution[],{nd}];Mean[x],{Nd}],ChartElementFunction->"FadingRectangle",ChartStyle->{Blue}]

Vysledek je pak (graf funkce normalniho rozdeleni je prakticky jenom cervena cara temer totozna s vodorovnou osou:

Show[hist,Plot[PDF[NormalDistribution[0,1],x],{x,-4,4},PlotStyle->Red]]

Spravne reseni:

histHustota=Histogram[Table[x=RandomVariate[NormalDistribution[],{nd}];Mean[x],{Nd}],Automatic,"PDF",ChartElementFunction->"FadingRectangle",ChartStyle->{Blue}]

nd

Plot[PDF[NormalDistribution[0,1/(Sqrt[nd])],x],{x,-4/(Sqrt[nd]),4/(Sqrt[nd])},PlotStyle->Red]

Show[histHustota,Plot[PDF[NormalDistribution[0,1/(Sqrt[nd])],x],{x,-4/(Sqrt[nd]),4/(Sqrt[nd])},PlotStyle->Red]]

12

Module[{N=1000,n=100,rozdel=StudentTDistribution[5],sem=SeedRandom[13],xx,vyber},sem;vyber=Table[xx=RandomVariate[rozdel,{n}];Mean[xx],{N}];SmoothHistogram[vyber]]

Takovou lokalizaci je mozne provest pomoci prikazu With, zde je lepsi primo definovat funkci uzivajici Module :

Clear[function12]
function12[iN_,in_,iRozdeleni_,iSeminko_:13]:=
	Module[{xx,vyber},
	SeedRandom[iSeminko];
vyber=Table[xx=RandomVariate[iRozdeleni,{in}];Mean[xx],{iN}];
SmoothHistogram[vyber]]

function12[1000,100,StudentTDistribution[5]]

13

nahodnyVyber := Module[{N,n,rozdelenie, seminko},]

Module musi mit dva argumenty: seznam lokalnich promennych a slozeny prikaz.

Řešení

Clear[rozdeleniPrumeru]
rozdeleniPrumeru[pocetVyberuN_,rozsahVyberun_,rozdeleni_,seminko_:13]:=
Module[{n=rozsahVyberun,nahodneVybery,xx,histogram,smoothHistogram,normalHustota,mean,sd},
SeedRandom[seminko];
mean=Mean[rozdeleni];sd=StandardDeviation[rozdeleni];
nahodneVybery=Table[xx=RandomVariate[rozdeleni,{n}];Mean[xx],{pocetVyberuN}];
histogram=Histogram[nahodneVybery,Automatic,"PDF",PlotLabel->rozdeleni];
smoothHistogram=SmoothHistogram[nahodneVybery,PlotLabel->rozdeleni,PlotStyle->Thick];
normalHustota=Plot[PDF[NormalDistribution[mean,sd/(Sqrt[n])],x],{x,mean-4sd/(Sqrt[n]),mean+4sd/(Sqrt[n])},PlotStyle->{Thick,Red}];
Show[{histogram,smoothHistogram,normalHustota}]
]

Úloha

Najděte empirické rozdělení výnosů společnosti Apple (Ticker AAPL) v období 1. 1. 2015 až 20. 11. 2018.

2018-12-10 Cvičení 11

FAQ

Úloha

Najděte empirické rozdělení výnosů společnosti Apple (Ticker AAPL) v období 1. 1. 2015 až 20. 11. 2018.

dataAAPL=FinancialData["AAPL","Return",{{2015,1,1},{2018,11,20}}];
%//DateListPlot[#,PlotRange->All]&

hodnotyAAPL=Map[Last,dataAAPL];

ClearAll[grafyJH]
grafyJH[data_?VectorQ]:=Module[{n=Length[data],histogram,smoothHistogram,normalHustota,hustotaFitDistribution,fitDistribution,mean,sd},
mean=Mean[data];sd=StandardDeviation[data];
histogram=Histogram[data,Automatic,"PDF"];
smoothHistogram=SmoothHistogram[data,PlotStyle->Thick];
normalHustota=Plot[PDF[NormalDistribution[mean,sd],x],{x,mean-5sd,mean+5sd},PlotStyle->{Thick,Red}];
fitDistribution=FindDistribution[data];
hustotaFitDistribution=Plot[PDF[fitDistribution,x],{x,mean-5sd,mean+5sd},PlotStyle->{Thick,Green},PlotRange->All];
Print["Nalezené rozdělení:   ",fitDistribution];
Print["p-hodnota testu dobré shody s normálním rozdělením = ",DistributionFitTest[hodnotyAAPL]];
Print["Použitý test:   ",DistributionFitTest[hodnotyAAPL,Automatic,"AutomaticTest"]];
Print[{Mean,StandardDeviation,Skewness,Kurtosis},":"];
Print[Through[{Mean,StandardDeviation,Skewness,Kurtosis}[hodnotyAAPL]]];
Show[{normalHustota,histogram,smoothHistogram,hustotaFitDistribution},PlotRange->All]]

Apple

grafyJH[hodnotyAAPL]
With[{m=3},Print["První ",ToString[m]," nalezená rozdělení:\n" , FindDistribution[hodnotyAAPL,3]]]

EUR/CZK

cenyEURCZK=FinancialData["EUR/CZK",{{2016,1,1},{2018,11,20}}];
cenyEURCZKWithoutDates=Map[Last,cenyEURCZK];
vynosyEURCZK=((Rest[cenyEURCZKWithoutDates]-Most[cenyEURCZKWithoutDates])/(Most[cenyEURCZKWithoutDates]));
grafyJH[vynosyEURCZK]

FindDistribution[vynosyEURCZK,3]

Směs rozdělení

Směs tří normálních rozdělení

Clear[mixture]
mixture[μ1_,μ2_,μ3_,σ1_,σ2_,σ3_,ϵ_]:=
ProbabilityDistribution[ϵ PDF[NormalDistribution[μ1,σ1],x]+(1-2ϵ)PDF[NormalDistribution[μ2,σ2],x]+ϵ PDF[NormalDistribution[μ3,σ3],x],{x,-∞,∞}]

mixture[-2,0,1,1,1,1,0.1]

Plot[PDF[mixture[-5,0,4,1,1,1,0.05],x],{x,-7,7},PlotStyle->Thick]

Statistika 2 - Jádrové odhady

Vyhlazené jádrové rozdělení (SmoothKernelDistribution)

(Subscript[X, 1]),...,(Subscript[X, n ]) náhodný výběr, k symetrická hustota zvaná jádro (kernel). Úkolem je najít neparametrický odhad hustoty f (neznáme její analytický tvar).

Jádrový odhad je  (OverHat[f])(x)=(1/(n h))(∑_(i=1)^n)k(((x-(Subscript[X, i]))/h)) pro vyhlazovací jádro k(x) a šířku okénka (bandwidth) h.

The probability density function for (RefLink[SmoothKernelDistribution, paclet:ref/SmoothKernelDistribution]) for a value x is given by a linearly interpolated version of (1/(n h))(∑_(i=1)^n)k(((x-(Subscript[x, i]))/h)) for a smoothing kernel k(x) and bandwidth parameter h.

Odpovídající zabudovaná funkce je

SmoothKernelDistribution

(*Jádrový odhad rozdělení*)
Clear[data];
SeedRandom[13];
data=RandomReal[NormalDistribution[10,3 ],1000];
distribuce=SmoothKernelDistribution[data]
GraphicsRow[{Plot[PDF[distribuce,x],{x,-10,20}],
Plot[CDF[distribuce,x],{x,-10,20}]}]

(* Ceny Volkswagen *)cenyVW=FinancialData["F:VOW",{{2014,9,22},{2018,(* 8 *)11,30}}];
%//DateListPlot

(* Výnosy Volkswagen *)returnsVW=FinancialData["F:VOW","Return",{{2014,9,22},{2015,(* 8 *)11,30}},"Value"];
grafyJH[returnsVW]

smooothHistogramVW=SmoothHistogram[returnsVW]
Text@Row@{"Šířka okénka:   ",SmoothKernelDistribution[returnsVW]["Bandwidth"]}

Manipulate[SmoothHistogram[returnsVW,bw,"PDF", Filling->Axis,PlotLabel->"Odhadnutá hustota výnosů VW",PlotRange->All],{{bw,.005},0.0000001,0.1,Appearance->"Labeled"},SaveDefinitions->True]

Apple

SmoothHistogram[hodnotyAAPL]

Text@Row@{"Šířka okénka:   ",SmoothKernelDistribution[hodnotyAAPL]["Bandwidth"]}

Manipulate[SmoothHistogram[hodnotyAAPL,bw,"PDF", Filling->Axis,PlotLabel->"Odhadnutá hustota výnosů Apple",PlotRange->All],{{bw,.005},0.0000001,0.1,Appearance->"Labeled"},SaveDefinitions->True]

Konec 2018-1210

2018-12-17 Cvičení 12

FAQ

20171211_Nonparametric_Regression, Kernel_Regression

https://en.wikipedia.org/wiki/Kernel_regression

https://stats.stackexchange.com/questions/6670/

Cambell, J. Y. et al.: Econometrics of Financial Markets

We wish to estimate the relation between two variables X and Y based on the observations ((Subscript[X, 1]),(Subscript[Y, 1])),...,((Subscript[X, N]),(Subscript[Y, N])) satisfying the relationship

(Subscript[Y, n])=m((Subscript[X, n]))+(Subscript[ε, n]),   n=1,...,N.

where m is an unknown function, (Subscript[ε, n])’s are iid with zero mean. A crucial role in kernel regression plays a kernel, here denoted K which is simply a symmetric probability density function. In practise we use a rescaling kernel

(Subscript[K, h])(x)=(1/h)K((x/h))

for some h>0 called bandwidth or window width. Note that (Subscript[K, h]) is also a symmetric density. The Nadaraya-Watson kernel estimator of m is

(OverHat[m])(x)=(((∑_(n=1)^N)(Subscript[K, h])(x-(Subscript[X, n]))(Subscript[Y, n]))/((∑_(n=1)^N)(Subscript[K, h])(x-(Subscript[X, n])))).

There are plenty of kernels and a vast number of recommendations in various situations. We will restrict our illustrations to two frequently used kernels: Gauss and Epanechnikov. These are also often used in kernel density estimation.

Epanechnikov kernel (Subscript[K, h])(u)=(1/h)(3/4)(SuperPlus[(1-(u^2))])

kernelGauss[x_]=PDF[NormalDistribution[],x]

Plot[{kernelEpanechnikov[x],kernelGauss[x]},{x,-3,3}]

plotEpaGau=Plot[{kernelEstimator[ir3,400,kernelEpanechnikov,x],kernelEstimator[ir3,400,kernelGauss,x]},{x,7,1000},PlotStyle->Thick]

Show[plotEpaGau,listPlotir3,PlotStyle->Thick]

ir3={{7,4.74`},{14,4.72833333333333`},{30,4.705`},{60,4.595`},{90,4.53`},{180,4.42`},{270,4.415`},{360,4.41`},{730,4.51`},{1095,4.68`},{1460,4.85`},{1825,4.99`},{2190,5.09`},{2555,5.19`},{2920,5.27`},{3285,5.33`},{3650,5.38`},{4015,5.42`},{4380,5.46`},{4745,5.48666666666667`},{5110,5.51333333333333`},{5475,5.54`}};pts=ir3;
listPlotir3=ListPlot[pts]

2019-01-07 Cvičení 13

FAQ

Obligace

FinancialBond[{"FaceValue"->1000,"Coupon"->.035,"Maturity"->10,"CouponInterval"->1/2},{"InterestRate"->{1->.02,2->.03,3->.04,5->.045,10->.06},"Settlement"->9.5,"DayCountBasis"->"ActualISDA"}]

Information[FinancialBond]

FinancialBond[

Download dat z Wolfram|Alpha

yields={WolframAlpha["yield curve",{{"TreasuryYieldCurve:EconomicData",1},"ComputableData"}],WolframAlpha["yield curve",{{"TreasuryYieldCurve:EconomicData",2},"ComputableData"}],WolframAlpha["yield curve",{{"TreasuryYieldCurve:EconomicData",3},"ComputableData"}]}

yields1=DeleteCases[yields,Missing["NotAvailable"]]//First

(QuantityPostfix[2.37`, "%", percent, "Percent"])//FullForm
(QuantityPostfix[2.37`, "%", percent, "Percent"])//Normal

yields2=yields1[[All,2]]
yields3=yields1[[All,2,1]]

Lépe:

yields3better=yields2//Normal

100 yields2/100

Quantity[2.37`,"Percent"]/100

Interpolace

yields1=DeleteCases[yields,Missing["NotAvailable"]]//First

yieldsyields1//Normal

timesAndYields={{1/4,0.024100000000000003`},{1,0.025`},{2,0.0239`},{5,0.023700000000000002`},{10,0.0256`},{30,0.0292`}};
plot1=ListPlot[%];
plot2=Plot[Interpolation[timesAndYields][t],{t,0,30},PlotStyle->{Thick,Red}];
plot3=Plot[Interpolation[timesAndYields,InterpolationOrder->1][t],{t,0,30},PlotStyle->{Thick,Green}];
plot4=Plot[Interpolation[timesAndYields,InterpolationOrder->4][t],{t,0,30},PlotStyle->{Thick,Black}];

Show[plot1,plot2,plot3,plot4,PlotRange->All]

Novoroční úloha

Nedostatky:

1. Formátování

Špatně je měnit fonty vstupních buněk:

Vstupní data

Kurzy

4+5

Citace zadání:

“Vypracujte protokol (slohové cvičení) podle následující osnovy (užívejte odstavce, tj. Section, Subsection apod.):”

Správně (umožňuje uzavírat a otvírat odstavce apod., viz Menu Format):

Vstupní data (Alt+5)

Kurzy (Alt+6)

4+5

2. Kontrola definicí

Před definitivním odevzdáním práce notes uložte, Mathematicu uzavřete a znovu otevřete a výpočty zopakujte. Nemůže se Vám pak stát toto:

sazbyA={1,2,3,4}
Total[sazbyB]

3. Maturity

Parametr Maturity je v našem příkladě vždy T=30. Mění se hodnotící okamžik t (Settlement). Doba do splatnosti je tedy T-t. Graficky možno znázornit buď v proměnné t nebo v inverzním čase n=T-t. Viz

[1] Dupačová, J., Hurt, J., Štěpán, J.: Stochastic Modeling in Economics and Finance. Kluwer Academic Publishers. Dordrecht 2002. str. 48 a dále.

V každém případě má funkce znázorňující spravedlivou hodnoty pilovitý tvar.

4. Kalendářová konvence

Při použití abstraktních časů (ne datumů) je irelevantní. Nevede však na chybovou zprávu. V zadání jistá past.

5. Znovu: Úspornost vstupu !!!

Místo

T=30

t=4

použít

T=30;
t=4;

6. Využívat zabudovaných funkcí, v tomto případě Thread

yields={{"1\[Hyphen]month treasury bill",(QuantityPostfix[2.42`, "%", percent, "Percent"])},{"3\[Hyphen]month treasury bill",(QuantityPostfix[2.41`, "%", percent, "Percent"])},{"6\[Hyphen]month treasury bill",(QuantityPostfix[2.47`, "%", percent, "Percent"])},{"1\[Hyphen]year treasury bill",(QuantityPostfix[2.5`, "%", percent, "Percent"])}};

Místo

intersetRates1={1/12->yields[[1,2]],1/4->yields[[2,2]],1/2->yields[[3,2]],1->yields[[4,2]]}

použít

intersetRates2=Thread[{1/12,1/4,1/2,1}->yields[[All,2]]]

7. Formální úprava

Místo

Row@{"F=",nominalValue=1000}
Row@{"c=",couponRate=0.027}

použít např.

Grid[{{"F",nominalValue=1000},{"c",couponRate=0.027}},Frame->All]

8. Nesmyslná konstrukce výnosové křivky

V konečném čase následující (byť fakticky správná) konstrukce nemůže vést k cíli pro dejme tomu 300 pozorování. Existuje např. zabudovaná funkce Interpolation. Ostatně i v jiných systémech existují podobné skripty.

interest[x_]=If[x<1/12,
0.0243,
If[x<1/4, 
0.0243+(0.0241-0.0243)/(1/6)*(x-1/12),
If[x<1/2,
 0.0241+(0.0249-0.0241)/(1/4)*(x-1/4),...

9. Interpolace

Velmi nebezpečná věc pokud se s ní nezachází s největší opatrností. Pokud se použije pouze funkce Interpolation bez zadání InterpolationOrder, dostaneme někdy těžko použitelný výsledek.

Úlohy

Úloha 1

Graficky znázorněte profil současné hodnoty a zjistěte vnitřní míru výnosnosti (IRR) peněžního toku

{{{2010,1,1},-600},{{2011,1,1},300},{{2012,1,1},400}}

(* Vstupní forma *){{{2010,1,1},-600},{{2011,1,1},300},{{2012,1,1},400}}

Nápověda: Cashflow, Timevalue

Úloha 2

Do investičního fondu dne 31. 1. 2019 vložíte PV=CZK 100000. Po dobu 36 měsíců Vám budou k aktuální částce připisovány výnosy ke konci měsíce v náhodné výši s rozdělením N((1.5/12)[%],rozptyl=((((5/12)))^2)[((%)^2)]).  Použijte SeedRandom[12345].

1. Jaký je stav fondu na konci uvažovaného období?
2. Graficky znázorněte průběh současné hodnoty odpovídajícího peněžního toku v závisloti na hodnotící úrokové míře.

Návod:
1. Pro generování datumů použijte funkci DatePlus s volbou step = ”Month” namapovanou na Range[37]-1.
2. Pro generování výše částek na konci měsíce použijte funkci NestList s počáteční hodnotou 100000.
3. Použijte funkci Cashflow.
4. Pro grafické znázornění použijte funkci DateListPlot po vhodném sloučení datumů a částek.