VPFPM 2009

Doporučený metodický postup pro studium systému Mathematica  (součást přednášky NFAP008)

Úplný začátek

Podívejte se na Help ⟶ Welcome Screen ...  ⟶ GETTING STARTED. Z přednášky si čiňte poznámky. Záznamy z přednášek na mé domovské stránce. Při zvětšujícím se rozsahu dostupná bez výstupů.

Pokročilejší přístup

Help ⟶ VIRTUAL BOOK  ⟶ Introduction ⟶ Getting Started
Užitečné a z matematického hlediska zajímavé informace jsou  v záložce "Demonstrations".

Důkladný průvodce

Help ⟶  Documentation Center. Projít letmo všechny kapitoly. Všimnout si též záložky How To Topics vpravo dole.

Přednášející radí následující postup pro porozumění systému:

(Nečiní si nicméně nárok na to, že jeho doporučení jsou suboptimální, natož optimální. Veškeré konstruktivní připomínky, podněty a náměty k přednášce VPFPM sdělené přednášejícímu jsou vítány. Je možné je podávat ústně, písemně i anonymně a e-mailem. Na připomínky podané e-mailem bude přednášejícím reagováno.)

(i)

Projít celou The Mathematica Book (záložka Help, Virtual Book). Tím se získá přehled, co všechno je možné pomocí systému realizovat. V této fázi není rozumné se věnovat detailům, navíc jest v prvním čtení nemožné porozumět všemu. Cílem je učinit si mlhavou představu o více než 2000 zabudovaných (build-in) funkcí, případně funkcí z Add-On Packages. Pak pomocí "Help" je v případě potřeby syntakticky správně zabudovat do uživatelského projektu.

(ii)

Jako jeden z efektivních postupů se jeví tento: Help ⟶ přibližný název funkce ⟶ přesný název funkce ⟶ (vybrat příbuznou funkci) ⟶ Further Examples ⟶ příslušnou část označit a  Copy+Paste do vlastního notebooku. U některých funkcí je k dispozici "Advanced Documentation". Doporučuji její přečtení. Čtenář se zde srozumitelným způsobem seznámi se současnými metodami, které obvykle nejsou v syllabech přednášek konaných na MFF.

(iii)

Není možné se naučit, natož pochopit všechno ihned. K nejmocnějším funkcím systému Mathematica (Apply, Map, Thread, případně Fold, Nest a funkcím podobným) je nezbytné se neustále vracet a zkoušet je v konkrétních situacích.

(iv)

Jádro systému Mathematica bylo naprogramováno v jazyku C. Půjdeme-li zpátky, pak C byl naprogramován v Assembleru, a ten přímo ve strojovém kódu. Pro efektivní užití systému Mathematica však není vhodná zakořeněná znalost ideologie procedurálních jazyků jako jsou Pascal, FORTRAN, C++ a další. Znalosti objektově orientovaného programování jsou naopak vítány. V systému Mathematica je sice možné prakticky doslova přepsat program napsaný v některém z těchto procedurálních jazyků, ale tím by byla zcela negována hlavní zbraň systému Mathematica, tj. funkcionalita, přehlednost, rychlost, kontrolovatelná přesnost a verifikovatelnost při modelování. Zejména z hlediska přehlednosti konstrukce modelů (a ve vyšší fázi použití nadstavby WorkBench pro vytváření projektů) jsou funkce typu "cyklus" naprosto nevhodné. Výhodou pro porozumnění systému Mathematica  je znalost logiky jazyků Prolog, LISP (byť velmi starých, které v době svého vzniku neměly podporu v potřebné počítačové technice) .

(v)

Pro úspěšné zvládnutí přednášené látky stačí zhruba 100 až 200 hodin samostudia. Aktivní zvládnutí celého systému je prakticky nemožné.

Lekce 1 (26. 2. 2009)

Motivace, odkazy

Zakladni stranky:

http://www.wolfram.com

http://www.wolfram.cz

Vyhradni distributor v CR a SR: http://www.elkan.cz/

Ukazky:

http://demonstrations.wolfram.com/

Podivejte se na Latest, pak Search: Finance, Insurance, Economics, Probability, Statistics. Demonstrace je mozne stahnout a studovat. Zacnete studovat ty jednodussi. Cim hezci grafika, tim komplikovanejsi kod a pro zacatecnika hure srozumitelny.

Klavesove zkratky

http://reference.wolfram.com/mathematica/guide/SelectingAndTypingInNotebooks.html

Velmi uzitecne: F1, Ctrl+K, Ctrl+Shift+K, Ctrl+Shift+B

Klavesove zkratky (vyhodna prevence neustaleho prehmatavani mezi mysi a klavesnici, da se vytisknout na dvou strankach podelne): http://media.wolfram.com/brochures/mathematicashortcuts.pdf

Mathematica book v papirove podobe v knihovne (strasi verze)

Novinky verze 7

http://www.wolfram.com/products/mathematica/newin7/

Mathematica pro zacatecniky (Antonin Slavik)

http://www.karlin.mff.cuni.cz/~slavik/info.html

The Mathematica Journal (TMJ) Pro pristup pouzijte licencni cislo L******** (viz Help -> About Mathematica):

The Mathematica Journal
Vyvolejte 
http://www.wolfram.com/services/premiersupport/

a jako licencni cislo zadejte L4697 - 4335 resp. L4667 - 9660.  Tam najdete vyznaceny The Mathematica Journal. Je pravdepodobne, ze bude fungovat pouze ze site na fakulte.

Specializovana cisla, statistika viz  Vol. 9(4)

4 pilire

Numericke vypocty

Vypocet pomoci Shift + Enter nebo Enter na numericke klavesnici

1+1

(2^1000)

N[(2^100000)]

π

npi=N[π,20]

N[π,{15,10}]

Sdruzovani vice vypoctu do jedne bunky

1+1
2+2
2+
3
2
+3

{π,N[π],"log 2 numericky = ",Log[2.]," log 2 presne = ",Log[2],N[Log[2]]}

Symbolicke vypocty

e1=Expand[(((x+y))^7)]

e2=Factor[e1]

TeXForm[e2/e1]

TraditionalForm@Expand[(((x-1))^5)]
TraditionalForm[Expand[(((x-1))^5)]]
Expand[(((x-1))^5)]//TraditionalForm

Clear[a,b,c,x]
reseni=Solve[a (x^2)+b x +c==0,x]

?x

?dbknsdfashgjfsfge

?reseni

koreny=x/.reseni

Integrate[Sin[(x^2)],x]

i1=(∫_0^∞)(E^(a (x^2)/2))\[DifferentialD]x

Assuming[a>0,(∫_0^∞)(E^(-a (x^2)/2))\[DifferentialD]x]

Integrate[(E^(-((a (x^2))/2))),{x,0,∞},Assumptions->{a>0}]

i1/.If[a_,b_,c___]->b

Grafika

{-Range[5],0,Range[5]}
{-Range[5],0,Range[5]}//Flatten
Plot[Sin[x+x^2],{x,0,3π},PlotRange->{-5,5},Ticks->{{0,Pi,2Pi,3Pi},{-5,0,1,4,5}}]
Plot[Sin[x+x^2],{x,0,π},PlotRange->{-5,5},Ticks->{{0,(Pi/2),Pi},{-Range[5],0,Range[5]}//Flatten}]
Plot3D[Sin[x+y^2],{x,0,π},{y,0,π}]

Options[Plot3D]

FinancialData["EUR/CZK"]

(* Dynamic[Refresh[FinancialData["USD/CZK","LatestTrade"],UpdateInterval->1]]
Dynamic[Refresh[FinancialData["EUR/CZK","LatestTrade"],UpdateInterval->1]] *)

Programovani

e1=Expand[(((x+y))^7)]

hu[x_]:=x^2

"="je okamzite prirazeni, ":=" je odlozene prirazeni

{hu[2],hu[Log[2]],hu[Log[2.]],hu[e1]}

Expand[hu[e1]]

{(Sin^2)[1.],((Sin[1.])^2)}

Poznamka k sazbe: (x^2) se vysazi jako x Ctrl ^ 2, (Subscript[x, 3]) jako xCtrl _  (Sqrt[x]) Ctrl2 x nebo z palety (Menu File ⟶ Palettes nebo pomoci Esc sekvenci

Π ⟶⟶⟵↔↘↗⥮\[ControlKey]\[EscapeKey]⌘

hu[x_,y_]:=x^y

{hu[2.],hu[2,3],hu[2,3,4]}

?hu

hu[2]=0
{hu[3],hu[2],hu[2.]}

?hu

{Unset[hu[2]],hu[2]}

?hu

hu[x_]:=(x^2)

hu[2]

hu1[x_]:=(x/2)

hu1[y]

?hu1

Clear[hu1] (* Totez je hu1=. *)

?hu1

Remove[hu1]
?hu1

?Global`*

Odstraneni casti definice

hu[2]=0;
hu[x_]:=(x^2)
hu[x_,y_]:=(x^y)

{hu[2],hu[3]}
hu[2]=.
{hu[2],hu[3]}

Clear[hil1];
hil1[n_]:=Table[(1/(i+j-1)),{i,n},{j,n}]

hil1[3]

h3m=hil1[3]//MatrixForm

Poznamka: f[x] je totez, co x//f resp. f@x

Poznamka: h3m neni totez co hil1[3]!!!

FullForm[h3m]

FullForm[hil1[3]]

Jak odstranit MatrixForm?
1. Nedefinovat symbol jako MatrixForm, ale MatrixForm pouzit pouze pro znazorneni:

(h4m=hil1[3])//MatrixForm

2. Aplikovat funkci Identity. Ta vymeni hlavicku vyrazu, tj. MatrixForm, za Identity.

Apply[Identity,h3m]//FullForm

Inverse[hil1[3]]//MatrixForm
Inverse[h3m]

h20=hil1[20];
First[h20inv=Inverse[h20]]
dh20=Det[h20]
N[dh20]
Det[h20inv]
Det[h20]*Det[h20inv]
First[h20.h20inv]

?Global`*

Funkce hj (libovolny pocet argumentu)

Clear[hj];
hj[x___]:=
Which[Length[{x}]==0,Print[x,"     No arguments"],
	  Length[{x}]==1,Print[x,"    Just one argument"],
          Length[{x}]>1,Print[{x},"   More than one argument so that ","\nThe sum of the elements = ",Total[{x}]]		]
{hj[],hj[3],hj[{1,2,3}],hj[1,2,34,5,6]}

hj[{1,2}⋂{3,4}]

Lekce 2

FAQ

Obycejny soucin * nebo mezera. Maticovy soucin . nebo funkce Dot

?a

aa=Table[a[i,j],{i,2},{j,3}]
(bb=Table[b[i,j],{i,3},{j,2}])//MatrixForm
cc={c[1],c[2],c[3]};
dd={d[1],d[2]};

Maticove nasobeni

aa.bb

bb.aa

* u matic znamena nasobeni po slozkach (Hadamarduv soucin):

aa*bb

aa*Transpose[bb]

cc cc

Definice funkce s jednim a vice argumenty

Pri vyvolani procedury se nejdrive vyhodnoti argument. Nevyhovuje-li predepsanemu vzoru (v nasem pripade dvojite __ znamena jeden a vice aergumentu), procedura se vubec nespusti. Viz priklad nize hj2[].

Clear[hj2];
hj2[x__]:=
Which[Length[{x}]==0,Print[x,"     No arguments"],
	  Length[{x}]==1,Print[x,"    Just one argument"],
          Length[{x}]>1,Print[{x},"   More than one argument so that ","\nThe sum of the elements = ",Total[{x}]]		]
{hj2[],hj2[3],hj2[{1,2,3}],hj2[1,2,34,5,6]}

Vnitrni struktura systemu Mathematica a jeho jazyk  I (Alt+5)

Vyrazy (Expressions)

Hlavicky

{Head[7],Head[(1/7)],Head[7.0],Head[2+I]}

Atomy = Symbol ∨ Cislo ∨ Retezec (String)

h=Head;
{h[f],h[Integrate],h[π],h[EulerGamma],h[hj2]}

{Head[7],Head[(1/7)],Head[7.0],Head[2+I]} (* cisla *)

h["Toto je retezec."]

Struktura vyrazu h[(Subscript[e, 1]),...,(Subscript[e, n])], kde h je hlavicka a (Subscript[e, 1]),...,(Subscript[e, n]) jsou vyrazy pripadne atomy.

Head[a+b+c]

FullForm[a+b+c]

{Head[{a,b,c}],FullForm[{a,b,c}]}

Length[{a,b,c}]
Length[{{a,b,c}}]
Length[a+b+c]

expr=Sin[x](a (x^2)+b x +c)

Head[expr]
FullForm[expr]

FullForm[(5/7)]

FullForm[HoldForm[(5./7)]]

Casti vyrazu( [[ se dosahne escape sekvenci Esc[[Esc)

soucet=a+b;
FullForm[soucet]
{Part[soucet,1],Part[soucet,2],Part[soucet,3],Part[soucet,0]}
{soucet[[1]],soucet[[2]],soucet[[3]],soucet[[0]]}

(a+b)[[1]]

{1,4,7}[[2]]

{"Co dostanu ted?"[[1]],"Co dostanu ted?"[[0]]}

{π[[1]],π[[0]]}

expr[[2,1]]

Table[expr[[j]],{j,0,3}]

expr[[1,3,2,2]]

(c+b x+a (x^2)) Sin[x]//TreeForm

aa[[1,2]]
aa[[1]]
aa[[All,2]]

Definice

Okamzite a odlozene prirazeni (globalni pravidlo, Immediate and delayed assignment): "=" a ":="

r1=Random[]
SeedRandom[13]
r2:=Random[]

{r1,r1,r1}
{r2,r2,r2}
SeedRandom[13]
{Random[],Random[],Random[]}

?r1

?r2

Funkce definovane po castech

Clear[g];
g[x_]:=(x^3)/;(*za podminky *)x≤0
g[x_]:=x/;0<x≤1
g[x_]:=Sin[x]/;x>1
Plot[g[x],{x,-1,2}]

D[g[x],x]

Odstraneni svisle carky

Plot[g[x],{x,-1,1,2}]

Plot[g[x],{x,-1,2},Exclusions->{x==1}]

Plot[g[x],{x,-1,2},Exclusions->True]

Plot[Tan[x],{x,-π,π},Exclusions->{x==π/2}]

Konec 5. 3. 2009

Zahada

Plot[Tan[x],{x,-π,π},Exclusions->{Tan[x]==π/2}]

{FindRoot[Tan[x]==x,{x,1}],N[(π/2)]}

Clear[g1];
g1[x_]:=Piecewise[{
{(x^3),x≤0},
{x,0<x≤1},
{Sin[x],x>1}
				}]
Plot[g1[x],{x,-1,2}]

Clear[g2];
g2[x_]:=Which[x<0,(x^3),
			0≤x<1,x,
			x≥1,Sin[x]
			]
Plot[g2[x],{x,-1,2}]

D[g2[x],x]

Plot[D[g2[x],x],{x,-1,2}]

Zde je videt dalsi zaludnost systemu: bez aplikace funkce Evaluate se funkce nezobrazi a objevi se chybova zprava. Viz Help k funkci Evaluate. Souvisi to s atributy funkce Plot. Prilis to nemusite studovat, spokojte se s tim, ze kdyz nektera funkce se nechce zobrazit, aplikujte na ni Evaluate.

Attributes[Plot]

Nyni jiz funguje :

Plot[D[g2[x],x]//Evaluate,{x,-1,2}]

Pouziti Evaluate

der2[x_]:=D[g2[x],x]
Plot[der2[x],{x,-1,2}] (* s odlozenym prirazenim nefunguje. *)

der2i[x_]=D[g2[x],x]
Plot[der2i[x],{x,-1,2}] (* s okamzitym prirazenim funguje. *)

Dalsi mozny, ale velmi spatny zpusob je pomoci funkce If.

g3[x_]:=If[x<0,(x^3),If[0≤x<1,x,If[x≥1,Sin[x]]]]
Plot[g2[x],{x,-1,2}]

g4[x_]:=If[x<0,(x^3),If[0≤x<1,x,If[True,Sin[x]]]]

Lekce 3

FAQ

$ImportFormats

Vnitrni struktura systemu Mathematica a jeho jazyk  II (Alt+5)

Funkce s vice definicemi

1+2=4

Clear[f,x,y,z];
f[x_]:=1+x+(x^2)
f[x_,y_]:=x+y
f[x_,y_,z_]:=(1/(x+y-z
))
f[x_,y_,z_]:=If[x==y==z,1,(1/(x+y-z
))]
f[1]=2;

?f

fz[x_]:=x
fz[x_]:=3x
?fz

fzz[x_]:=3x
fzz[x_]:=x
fzz[x_]:=x+2
?fzz

Cviceni 1

SeedRandom[13]
r[1][n_]:=Table[Random[],{n}]
r[2][n_]:=(x:=Random[];Table[x,{n}])
r[3][n_]:=(x=Random[];Table[x,{n}])
r[4][n_]=Table[Random[],{n}]

Table[r[i][3],{i,4}]//TableForm
r[4][3]
r[4][3]

Cviceni 2

fsum[x_]=Sum[(1/(((1+x))^j)),{j,n}]
gsum[x_]:=Sum[(1/(((1+x))^j)),{j,n}]

fsum[i]
gsum[i]

Clear[y];
f1[y_]=Sum[(1/(((1+y))^j)),{j,n}]
f1[i]

gsum je pochopitelne casove narocnejsi nez f1:

Timing[Table[f1[i],{i,10000}];]
Timing[Table[gsum[i],{i,10000}];]

Predikaty a booleovske operace

Predikaty

?*Q

PrimeQ[144]

{OddQ[21],EvenQ[21]}

IntegerQ[5/9]

{NumericQ[π],NumericQ[∞]}

{NumberQ[5/9],NumberQ[π],NumberQ[1+4.5I]}

AtomQ["asbhjo[dfuhdb\. n"]

{ListQ[{a,{b,c}}],VectorQ[{a,{b,c}}],MatrixQ[{{a,b},{b,c}}]}

IntervalMemberQ[Interval[{3,4}],π]

Booleovske funkce a logicke operace

7<5

Equal[3,7-4,6/2]

3==7-4==6/2
HoldForm[3==7-4==6/2]//FullForm

expr=((x^2)-1==(((x^4)-1)/((x^2)+1)))

Tan[z]==((Sin[z])/(Cos[z]))

Nekdy je nutne pouzit funkci Simplify.

{6≠12/6,Unequal[6,12/2]}

Truth Table

Example of Distribute and a Table

TruthTable[vars_List,expr_]:=Distribute[Table[{True,False},{Length[vars]}],List,List,List,{##,Function[vars,expr][##]}&]
DisplayTruthTable[vars_List,expr_]:=With[{t=ReplaceAll[TruthTable[vars,expr],{True->"T",False->"F"}]},Grid[Prepend[t,Append[vars,expr]],GridFrame->True,RowLines->Prepend[Table[0,{Length[t]-1}],2],ColumnLines->Append[Table[0,{Length[vars]-1}],2]]]

DisplayTruthTable[{A,B,C},Implies[Or[A,B],C]]

DisplayTruthTable[{A,B,C,D},Implies[Or[A,B],Or[C,D]]]

DisplayTruthTable[{A,B,C,D},Implies[(A∨B),(C∨D)]]

DisplayTruthTable[{A,B},A∧B]

DisplayTruthTable[{A,B},Xor[A,B]]

BaseForm[14,2]

2^^1110

BaseForm[BitOr[2^^10101,2^^10011],2]

Atributy

Attributes[Plus]

Plus[Plus[a,b]]//FullForm

Attributes[Sin]

Range[10]

Sin[{1,2,3}]
Sin[(π/4)(Range[10]-1)]

Clear[druhemocniny];
druhemocniny[x_]:=(x^2)+Length[x]
druhemocniny[Range[3]]
{Length[3],druhemocniny[3]}
Attributes[druhemocniny]
SetAttributes[druhemocniny,Listable]
Attributes[druhemocniny]

Remove[plus];
plus[a,plus[b,c]]
SetAttributes[plus,Flat]
plus[a,plus[c,b]]
SetAttributes[plus,Orderless]
plus[a,plus[c,b]]

Ryzi, anonymni funkce (Pure functions)

Ryzi funkce ma formu Function[promenna, predpis]

Function[x,(x^2)]

f=Function[x,(x^2)]
{f[x],f[3],f[I]}

Function[x,(x^2)][y]

v uspornem zapise  # znaci promennou ryzi funkce a & znaci ryzi funkci. Ekvivalentni zapis druhe mocniny

((#)^2)&

{((#)^2)&[x],((#)^2)&[3],((#)^2)&[I]}

Function[{x,y},x+y][cc,bbbb]

Function[x,x.x][{aa,bb,cc}]

Function[x,x.x][qq]
Function[x,x.x][{qq}]
Function[x,x.x][(({{a, b, c}, {d, e, f}, {g, h, i}}))]
Function[x,x.x][(({{a, b, c}, {d, e, f[g]}, {g, h, i}}))]

Definice ryzich funkci vice promennych

Function[{x,y},x+(y^2)][a,b]
(#1+((#2)^2))&[a,b]

Co se stane, kdyz nasobime matice, ktere nemaji stejne stycne rozmery?

(({{a, b, c}, {a1, a2, a3}, {d, e, ff}, {g, h, i}})).(({{a, b, c}, {a1, a2, a3}, {d, e, ff}, {g, h, i}}))
(({{a, b, c}, {a1, a2, a3}, {d, e, ff}, {g, h, i}}))(({{a, b, c}, {a1, a2, a3}, {d, e, ff}, {g, h, i}}))

Off[Dot::"dotsh"]
(({{a, b, c}, {a1, a2, a3}, {d, e, ff}, {g, h, i}})).(({{a, b, c}, {a1, a2, a3}, {d, e, ff}, {g, h, i}}))
On[Dot::"dotsh"]

Lekce 4 (19. 3. 2009)

Viz novy odkaz na klavesove zkratky
Klavesove zkratky (vyhodna prevence neustaleho prehmatavani mezi mysi a klavesnici): http://media.wolfram.com/brochures/mathematicashortcuts.pdf

FAQ

KJ questions

Q1 Jsou v systemu Mathematica zabudovany funkce IRR, "vnitřní míra výnosnosti" (nebo i např. jiné finanční funkce jako je durace apod.), neboť se mi nic takového nepodařilo najít, ale mám pocit, že to existuje. Našla jsem pár demonstrací a tam byla durace i IRR počítána "klasicky", dle známé rovnice..... Samozřejmě, že se to dá obejít, pokud si vyvtoříme funkci NPV(i) a pak hledáme řešení rovné 0, ale umí to Mathematica nějak sama pouhým zadáním jednotlivých cash flow (a popřípadě časovým určením)?

A1 Nejsou , ale

 (i) velmi jednoduche zpusoby programovani charakteristik peneznich toku budou uvedeny na prednasce, 

(ii) financni funkce jsou v komercnich balicich (ty jsou sice drahe, ale porad podstatne levnejsi a vykonejsi nez treba v systemu MATLAB) napr.

Unrisk
http://www.unrisk.com/
dokumentace viz
http://www.unrisk.com/download/pe/UnRiskPE30Docu.pdf

Financial Essentials (v pristi verzi bude soucasti systemu)
http://www.wolfram.com/products/applications/finance/

Derivatives Expert IV
http://www.wolfram.com/products/applications/derivatives/
dokumentace viz
http://www.wolfram.com/products/applications/derivatives/manual.pdf

V teto souvislosti jest vhodne tez upozornit na balik 

Time Series
http://www.wolfram.com/products/applications/timeseries/
inspirujici dokumentace viz
http://media.wolfram.com/documents/TimeSeriesDocumentation.pdf

Q2 Existuje nějaká klávesová zkratka na přepínání mezi jednotlivými notebooky?  
A2  Mozna, ze ano, ale nejjednodussi je: Menu ⟶ Window ⟶ spodni okno vybrat zadouci notebook. Ja sam pouziji organizaci na desktopu  s vhodnymi prekryvy, musim take pocitat s velmi castym volanim Help resp. F1.

Graficke znazorneni stromu

TreePlot

Ryzi, anonymni funkce (Pure functions)

Prvni ukazka aplikace ryzi funkce

data=Range[20];
(* Nasledujicimi dvema prikazy vytvorime velmi odlehla pozorovani (outliers) *)
data[[5]]=1000;
data[[15]]=-1000;
ListPlot[data,PlotStyle->PointSize[0.02]]

Vnitrni algoritmus pro grafiku extremy nezobrazi. Muzeme vsak pouzit nasledujici Option, pak se nam vsak s nejvetsi pravdepodobnosti ztrati to, co jsme chteli videt.

ListPlot[data,PlotStyle->PointSize[0.02],PlotRange->All]

Pouzitim ryzi fce v zabudovane funkci Select vyradime odlehla pozorovani. Za odlehla pozorovani budeme uvazovat ta, ktera jsou v avbsolutni hodnote ≥ 30. Podivejte se do Help na ListPlot, pro zduvodneni nasledujiciho vystupu.

withoutoutliers=Select[data,Abs[#]<30&]
ListPlot[withoutoutliers,PlotStyle->PointSize[0.02]]

Pridame x-souradnice

data2={Range[20],((Range[20])^2)}//Transpose
ListPlot[data2,PlotStyle->PointSize[0.01]]
data2[[5,2]]=-10000
ListPlot[data2]

withoutoutliers2=Select[data2,Abs[#[[2]]]<1000&]
ListPlot[withoutoutliers2]

Map, Apply, Thread (provleceni?, vlakno?) I (Obecny princip)

Map preskoci hlavicku druheho argumentu a apluje prvni argument na to, co je za hlavickou.
Apply vymeni hlavicku druhueho vyrazu za prvni.

Clear[f,g,a,b,c];
Map[f,g[a,b,c]]
Apply[f,g[a,b,c]]

Thread[f[{a,b,c}]]

Thread[f[{a,b,c},x]]

Thread[f[{a,b,c},{x,y,z}]]

Convert equations for lists to lists of equations:

Thread[{a,b,c}=={x,y,z}]

Apply a function to both sides of an equation:

Thread[Log[x==y],Equal]

Map a Apply

Clear[f,g,a,b,c];
li={a,b,c};
soucin=a b c;
{Head[li],Head[soucin]}
{FullForm[li],FullForm[soucin]}

Map preskoci hlavicku druheho argumentu a applikuje prvni argument na to, co je za hlavickou druheho argumentu

Map[Plus,soucin]

Apply[Plus,soucin]

Apply[Plus,{a,b,c}]

Plus[{a,b,c}]

Total[{a,b,c}]

Apply je funkce, ktera vymeni hlavicku druheho argumentu za prvni argument

Apply[Times,a+b+c]

Apply[Times,li]
Apply[List,soucin]
Apply[Plus,soucin]

Priklad vypoctu aritmetickeho prumeru s pojistkou na to, aby argument funkce byl vektor delky alespon jedna.

prumer[x_]:=Apply[Plus,x]/Length[x]/;VectorQ[x]&&Length[x]>0
{prumer[{a,b,c}],prumer[{{a},b,c}],prumer[a,b,c],prumer[2]}

{{{a}},b,c}//Total

{{{a}},b,c}//Total//Flatten

{{{a}},b,c}//Total//Flatten//First

Pouziti zabudovane funkce Total (pozor na formu vstupu!!!

Clear[pr];
pr[data_]:=((Total[dt=data])/(Length[dt]))
pr[{{{a,d}},b,c}]

Clear[pr];
pr[data_]:=((Total[dt=data//Flatten])/(Length[dt]))
pr[{{{a,d}},b,c}]

Priklad (vypocet rozptylu)

rozptyl[x_?VectorQ]:=prumer[(((x-prumer[x]))^2)]
{rozptyl[Range[5]],rozptyl[{Range[5]}]}

Jsou vsechny prvky seznamu suda cisla?

{EvenQ[3],EvenQ[6]}
Map[EvenQ,{1,2,3,4,5}]
And[{False,True,False,True,False}]
And[False,True,True]

jsouvsechnyprvkyseznamusude[list_]:=Apply[And,Map[EvenQ,list]]
jsouvsechnyprvkyseznamusude[Range[5]]
jsouvsechnyprvkyseznamusude[2Range[5]]

Priklad na postupna maxima (rekordy)

{First[r5=Range[5]],Last[r5],Most[r5],Rest[r5]}

Max[Range[5]]
FoldList[f,x,{a,b,c,d}]
FoldList[Plus,0,{a,b,c,d}]
r6=FoldList[Max,0,{5,3,6,3,2,77,4,888}]
r7=Rest[r6]
Union[r7]

Clear[rekordy];
rekordy[x_]:=Union[Rest[FoldList[Max,-1+Min[x],x]]]
rekordy[{5,3,6,3,2,77,4,888}]
rekordy[-{5,3,6,3,2,77,4,888}]

Modelovani ceny akcie

Clear[cena];
cena[s0_?NumericQ,delka_?IntegerQ]:=FoldList[Plus,s0,RandomReal[{-s0/20,s0/20},delka]]

cena[100,10]

ListPlot[Table[cena[100,1000],{5}],Joined->True]

Seznamy (Lists)

{b,c,a,aa}⋃{d,e,f,e,a}
{b,c,a,aa}⋂{d,e,f,e,a}
Complement[{a,b,c,d,e},{a,c},{d}]

r8=RandomInteger[{-4,4},15]
Sort[r8]
Select[r8,Positive]

Table[(i^2),{i,1,5}] (* spatne *)
Table[(i^2),{i,5}] (* spatne *)
((Range[5])^2)(* spravne *)

Range[4,8]
CharacterRange["a","z"]
CharacterRange["a","z"]//FullForm
{a,b,c}
{a,b,c}//FullForm

t1=Table[3k,{k,1,10,2}]
Length[t1]
Range[3,27,6]

Table[{i,Expand[(((1+α))^i)]},{i,3}]//TableForm

{Range[3],Expand[(((1+α))^(Range[3]))]}//Transpose//TableForm

Table[i+j,{i,5},{j,i}]//TableForm

(∑_(i=1)^n)(∑_(j=1)^i)1

Sum[1,{i,n},{j,i}]

Integrate[x y,{x,0,1},{y,0,x}]

Length[{Range[4],((Range[4])^2)}]
Length[{Range[4],((Range[4])^2)}//Transpose]

Dimensions[r10={Range[4],((Range[4])^2)}]
Dimensions[{Range[3],((Range[4])^2)}]
Dimensions[{Range[3],((Range[4])^2)},Infinity]

Position[{1+x^2,1+x,5,x^4,a+(1+x^2)^2},x^_] (* x sammotne neni vwe forme x^ !!!*)

{r10[[0]],r10[[2]],r10[[2,1]],r10[[All,2]],r10[[2,All]],r10[[6]]}
{Sin[x][[0]],Sin[x][[1]]}

mat=Table[(Subscript[a, i,j]),{i,3},{j,4}]

Table[(Subscript[a, ij]),{i,3},{j,4}]

Table[(Subscript[a, i j]),{i,3},{j,4}]

Map[FullForm,{(Subscript[a, 1,2]),a12,a[1,2]}]

{(Subscript[a, 1,2]),a12,a[1,2]}={3,3,3}

{(Subscript[a, 1,2]),a12,a[1,2]}

Clear[(Subscript[a, 1,2]),a12,a[1,2]]

{(Subscript[a, 1,2]),a12,a[1,2]}

{(Subscript[a, 1,2]),a12,a[1,2]}={3,3,3}
Remove[(Subscript[a, 1,2]),a12,a[1,2]]
{(Subscript[a, 1,2]),a12,a[1,2]}

?a[1,2]

?(Subscript[a, 1,2])

(Subscript[a, 1,2])=.

Casti seznamu

r11=Range[7];
Take[r11,2]
Take[r11,-2]

Drop[r11,2]
Drop[r11,-2]

CharacterRange["1","4"]
CharacterRange["1","14"]

indices=Outer[StringJoin,CharacterRange["1","4"],CharacterRange["1","8"]]//FullForm

Lekce 5 (26. 3. 2009)

FAQ

Dotazy tykajici se /. a /;

/.

{a,b,c}/.x_->{x,d}

Map[f,g[a]]
Map[f,g[a,b,c]]
Map[((#)^2)&,g[a,b,c]]

Map[ReplaceAll[#,x_->{x,d}]&,{a,b,c}]

Map[#/.x_->{x,d}&,{a,b,c}]

{a,b,c}/.Rule[x_,Hold[Thread[{x,y}]]]//ReleaseHold

/;

fac[n_/;n>0]:=n!
fc[n_]:=n!/;n>0
fcc[n_?Positive]:=n!
{fac[1.5]+fac[-2],fc[1.5]+fc[-2],fcc[1.5]+fcc[-2]}

(1/2)!

Indexovani matic bez carek mezi indexy

indices=Outer[StringJoin,CharacterRange["1","4"],CharacterRange["1","8"]]//Flatten

Clear[aa];
(aa=Partition[Table[(Subscript[a, ij]),{ij,indices}],8])//MatrixForm

Spravny odkaz (reference)

aa[[3,4]]

Spatny odkaz

(Subscript[a, 34])

Porovnejte testem a podivejte se na uplnou formu

(Subscript[a, 34])===aa[[3,4]]

{(Subscript[a, 34]),aa[[3,4]]}//FullForm

x==y
x===y

Pozdeji bude zminka o asociovanych polich (velmi mohutny nastroj indexovani a odkazovani)

Seznamy (Lists) pokracovani I

Timing[Table[Sin[π k/4],{k,0,50000}];]
Timing[Sin[(Range[50001]-1)π/4];]

Table[Sin[π k/4],{k,0,5,1/2}]

Table[Sin[π k/4],{k,{0,5,,3}}]

Vektory a matice

vektor={a,b};
xx={(Subscript[x, 1]),(Subscript[x, 2])};
yy={(Subscript[y, 1]),(Subscript[y, 2])};
matice={{x,y},{z,v}};
vektor.matice
matice.vektor
{xx.yy,xx*yy,xx yy}

Transpose[{xx}]
Grid[Transpose[{xx}].{yy},Frame->All]
Transpose[{xx}].{yy}//Grid[#,Frame->All]&

Transpose[{xx}].yy
yy.Transpose[{xx}]

Grid[Outer[Times,xx,yy],Frame->All]

Cross[{x1,x2,x3},{y1,y2,y3}]//TableForm

(aa=Array[a,{2,3}])//MatrixForm

aa.Transpose[aa]//MatrixForm

IdentityMatrix[3]//MatrixForm

(dm=DiagonalMatrix[{x1,x2,x3}])//MatrixForm

dm+1

Diagonal[dm+1]

Priklad

Vytvoreni positivne definitni matice. Vytvorime nahodnou ctvercovou matici C. Ta je s. j. regularni Potom C.(C^T) je s. j. positvne definitni..

c[n_]:=RandomReal[{0,1},{n,n}]
c[2]
posdef[n_]:=Module[{pom},pom=c[n];pom.Transpose[pom]]
Det[posdef[50]]
Det[(p=posdef[50]).Inverse[p]]

?pom

Eigenvalues[p]

cisla=HilbertMatrix[2]//Eigenvalues
vektory=HilbertMatrix[2]//Eigenvectors
eig1=HilbertMatrix[5]//Eigenvalues
N[eig1]

Map[Norm,vektory]//N
First[vektory].Last[vektory]

cisla50=Eigenvalues[p];
vektory50=Eigenvectors[p];

Map[Norm,vektory50]

Outer[Dot,vektory50,vektory50]

D{{-0.5212491759206568`,-0.5185970200880448`,-0.4156332345853547`,-0.2965023701634475`,-0.4457485683133595`},{-0.5212491759206568`,-0.5185970200880448`,-0.4156332345853547`,-0.2965023701634475`,-0.4457485683133595`}}

p5=posdef[5];
cisla5=Eigenvalues[p5];
vektory5=Eigenvectors[p5];
tup5=Tuples[vektory5,2];
tup5//First
Map[Dot[#[[1]],#[[2]]]&,tup5]//Chop

Transpose[vektory5].DiagonalMatrix[cisla5].vektory5-p5//Chop

Odmocninova matice

Transpose[vektory5].DiagonalMatrix[(Sqrt[cisla5])].vektory5-MatrixPower[p5,(1/2)]//Chop

Timing[HilbertMatrix[30]//Eigenvalues//N;]
Timing[HilbertMatrix[60]//N//Eigenvalues;]

Ukazka programu viz napr. C:\Program Files\Wolfram Research\Mathematica\7.0\AddOns\Packages\Calendar

Testovani a hledani v listech

{Position[li={a,b,c,d,e},c],Position[{a,b,c,d,e},cc]}

Count[li2={a,b,a,c,a,d,c,b,b,b,a,e,e},a]

Union[li2]

{Union[li2],Map[Count[li2,#]&,Union[li2]]}//Transpose//TableForm

{MemberQ[li2,a],MemberQ[li2,f]}

Lekce 6 (2. 4. 2009)

FAQ

Dodatek k /. (K listu {a,b,c} pridat d ke kazdeme elementu.)

Map[ReplaceAll[#,x_->{x,d}]&,{a,b,c}]

Map[#/.x_->{x,d}&,{a,b,c}]

{a,b,c}/.Rule[x_,Hold[Thread[{x,y}]]]//ReleaseHold

Uplne nejjednodussi je ale :

Map[{#,d}&,{a,b,c}]

c[n_]:=RandomReal[{0,1},{n,n}]
posdef[n_]:=Module[{pom},pom=c[n];pom.Transpose[pom]]

Test, zda matice je ortogonalni matice (Klara Jelenova)

p5=posdef[5];
cisla5=Eigenvalues[p5];
vektory5=Eigenvectors[p5];
vektory5.Transpose[vektory5]//Chop

Aplikace funkci na radky nebo sloupce matice

mat={{a,b,c},{d,e,f}};
{Total[mat],Total[Transpose[mat]],Map[Total,mat]}

Analogicky Mean:

{Mean[mat],Mean[Transpose[mat]],Map[Mean,mat]}

Generovani LDAP hesla

LDAP je heslo pro autentizacni system UK, 16 znaku, 1. znak velke pismeno, dalsich 14 znaku mala pismena, 16. znak cislice (asi s vyloucenim nuly), z pismen jsou vylouceny "y" a "z".

capitalswithoutYandZ=Complement[CharacterRange["A","Z"],{"Y","Z"}]
lowercasewithoutyandz=Complement[CharacterRange["a","z"],{"y","z"}]
digitsascharacters=Map[ToString,Range[9]]
StringJoin[RandomChoice[capitalswithoutYandZ],
RandomChoice[lowercasewithoutyandz,14],
RandomChoice[digitsascharacters]]

Seznamy (Lists) pokracovani II a dokonceni

m=IdentityMatrix[3];
{FreeQ[mat,0],FreeQ[m,0]}
Position[m,0]

Pridavani, odstranovani a modifikovani prvku seznamu

Prepend[m={a,b,c},aa]
Append[m,aa]
AppendTo[m,aa]
m
AppendTo[m,aa]

Riffle[m={a,b,c},x]

🐺

ReplacePart[m,2->x]
mat🐺=ReplacePart[mat={{a,b,c},{d,e,f}},{1,2}->🐺]

{mat🐺[[2]]=🐺🐺🐺,mat🐺}

mat[[All,1]]={x,y};
mat

Kombinace listu

Join[{a,b,c},{x,y},{t,u},{b,c},a]
Join[{a,b,c},{x,y},{t,u},{b,c},{a}]

Union[{a,b,c},{x,y},{t,u},{b,c},{a}]

Riffle[{a,b,c},{x,y,z}]

Konstrukce listu

Range, Table, Array

NestList[f,x,5]

NestList[2#&,1,5]
NestList[Function[x,2x],1,5]

Sparse arrays (Ridke matice)

sa1=SparseArray[{3->x,4->y},15]

sa1.sa1
sa1*sa1

Normal[sa1]
Normal[sa1*sa1]

SparseArray[{i_->(i^2)},5]//Normal

Select[{1,3,7,4,10,2},#<5&]

TakeWhile[{1,3,7,4,10,2},#<5&]

{a,b,c,d}[[{2,1,4}]]

Pick[{a,b,c,d,e},{1,1,0,1,0},1]
Pick[{a,b,c,d,e},{1,1,0,1,0},0]

Jak zasejes, tak sklidis. (Sow, Reap)

NestList[((#)^2)&,2,6]
Nest[((#)^2)&,2,6]

Reap[Nest[(If[#>1000,Sow[#]];((#)^2))&,2,6]]

Poznamky k indexovani

mat={{x,b,c},{y,e,f},{h,i,j}};
mat[[1,2]]
mat[[1]][[2]]
mat[[{1,2}]]
mat[[{1,2},1]]

mat[[2;;3]]

Take[mat,2]
Take[mat,-2]
Drop[mat,1]
Drop[mat,-1]

Extract[mat,{{1,3},{1,2}}]
mat[[{{1,3},{1,2}}]]

Velmi dulezity priklad

f[x_,y_]:=x (y^2)
{f[2,1],f[{2,1}],Apply[f,{2,1}]}

Viz take Sequence!

ff[{x_,y_}]:=x (y^2)
{ff[2,1],ff[{2,1}]}

fff[x_,y_]:=x (y^2)
fff[{x_,y_}]:=x (y^2)
?fff

{fff[2,1],fff[{2,1}]}

{Sin[2.,1],Sin[{2.,1}]}

Clear[sin];
sin[y__]:=Sin[{y}]//Flatten
{sin[1.],sin[1.,2],sin[1.,2,3],{sin[1.,2],sin[1.,2,3]}}

Dalsi konstrukce

Outer[f,mat]

Deleni a vyplnovani seznamu (partitioning and padding)

Partition[{a,b,c,d,e,f,g},3]

Deleni s prekrytim (offset)

Partition[{a,b,c,d,e,f,g},3,1]
Partition[{a,b,c,d,e,f,g},3,3,1]
Partition[{a,b,c,d,e,f,g},3,3,1,x]
Partition[{a,b,c,d,e,f,g},3,3,1,{x,y,z,v,w}]
Partition[{a,b,c,d,e,f,g},3,1,{1,1},x]

Sparse arrays

SparseArray[{2->a,5->b}]//Normal

SparseArray[{{1,2}->a,{3,2}->b,{3,3}->c}]//Normal

(tab1=Table[{i,j},{i,4},{j,i}])//Length
Flatten[tab1,1]
(SparseArray[(indexy1=Flatten[tab1,1])->Range[Length[indexy1]]]//Normal)//MatrixForm

SparseArray[{Band[{1,1}]->x,Band[{2,1}]->y},{5,5}]//MatrixForm

SparseArray[{i_,i_}->1,{4,4}]//MatrixForm

SparseArray[{{i_,i_}->1,{_,3}->a,{2,_}->b},{4,4}]//MatrixForm

SparseArray[{i_,i_}->RandomReal[],{4,4}]//MatrixForm
SparseArray[{i_,i_}:> RandomReal[],{4,4}]//MatrixForm

SparseArray[i_->(i^2),{10}]//Normal

(* Vektor, ktery bude mit na sudych pozicich p a na ostatnich druhou mocninu *)SparseArray[{_?EvenQ->p,i_-> (i^2)},{15}]//Normal

SparseArray[{1|3,2|4}->a,{4,4}]//MatrixForm

SparseArray[i_/;3<i<7->p,{10}]//Normal

SparseArray[{i_,j_}/;Abs[i-j]<2->i+j,{5,5}]//MatrixForm

SparseArray[{Band[{1,1}]->2Range[5],Band[{1,2}]->2Range[4]+1,Band[{2,1}]->2Range[4]+1}]//MatrixForm

SparseArray[{Band[{1,1}]->2Range[5],Band[{1,2}]->2Range[4]+1,Band[{2,1}]->2Range[10]+1}]//MatrixForm

Nespecifikovana pole jsou misto nul nahrazeny karetnimi symboly

SparseArray[{i_,j_}/;Abs[i-j]<2->i+j,{5,5},♡♢♣♠]//MatrixForm

Konec 2. 4. 2009

Lekce 7 (9. 4. 2009)

FAQ

Principy programovani

Definice (funkce logaritmus)

Definice je prepisovaci pravidlo, ktere je automaticky pouzito na libovolny vyraz. Priklad na definici logaritmu:

Clear[log];
log[x_*y_]:=log[x]+log[y]
log[2x (y^2)]

Prestoze z matematickeho hlediska je (y^2) totez, co y*y, potrebujeme dalsi prepisovaci pravidlo, protoze porovnani vzoru (pattern matching) je v zasade syntakticke, a ne semanticke!

log[((a_)^(n_))]:=n*log[a]
log[2x (y^2)]

Clear[klaralog];
klaralog[((a_)^(n_))*b_]:=n klaralog[a]+klaralog[b]
klaralog[b_*((a_)^(n_))]:=n klaralog[a]+klaralog[b]
klaralog[x*2*y^3*z]

Definice (linearni funkce, M220, p. 63)

NumericQ[1,π,Sin[(Sqrt[E])],a,Exp[x]]

NumericQ[{1,π,Sin[(Sqrt[E])],a,Exp[x]}]

Map[NumericQ,{1,π,Sin[(Sqrt[E])],a,Exp[x]}]
Map[NumericQ[#]&,{1,π,Sin[(Sqrt[E])],a,Exp[x]}]
Map[Function[x,NumericQ[x]],{1,π,Sin[(Sqrt[E])],a,Exp[x]}]

Pridani listovatelnosti k NumericQ

Clear[numericQ];
numericQ[x_]:=NumericQ[x];
SetAttributes[numericQ,Listable]
numericQ[{1,π,Sin[(Sqrt[E])],a,Exp[x]}]

Clear[f,x,y,a];
f[c_?NumericQ*x_]:=c*f[x]
f[x_+y_]:=f[x]+f[y]
f[π+(Sqrt[2])x+3 a (x^2)]

Ale zcela v poradku to neni:

{f[2π],f[π 2]}

Pridame

f[c_?NumericQ]:=c f[1]/;c≠1
{f[2π],f[π 2]}

Ruzne programovaci styly

Nasledujici funkce zjistuji cetnosti.

Proceduralni styl

SeedRandom[130];
(* data=Table[Random[Integer,{1,10}],{100}]; *)
data=RandomInteger[{1,10},{100}];

proc[x_]:=Module[{bins=Table[0,{10}]},
Do[bins[[x[[i]]]]++,{i,1,Length[x]}];bins]
pr=proc[data]
Total[pr]

Funkcionalni styl

Clear[fun];
fun[x_?VectorQ]:=fun[Map[{#,1}&,x]]
fun[{1,2,3,2}]

Poznamka: _ prave jeden argument, __ jeden a vice argumentu, ___ nula a vice argumentu, x___

fun[{a___,{b_,n_},c___,{b_,m_},d___}]:=fun[{a,{b,n+m},c,d}]

fun[{1,2,3,2}]

fun[a_List]:=Map[Last,a];
data
fun[data]
fun[data//Sort]

fun[a_List]:=a//Sort;
cetnosti=fun[data]

Pochazeji tato cisla skutecne z diskretniho rovnomerneho rozdeleni?

(χ^2)-test dobre shody pro tento konkretni pripad

cetnosti2=Map[Last,cetnosti]
pocetruznychcisel=Length[cetnosti2]
pocetpozorovani=Total[cetnosti2]
((((cetnosti2-pocetpozorovani/pocetruznychcisel))^2)/(pocetpozorovani/pocetruznychcisel))
chi2=Total[((((cetnosti2-pocetpozorovani/pocetruznychcisel))^2)/(pocetpozorovani/pocetruznychcisel))]//N

Print["Hodnota (χ^2) na ",pocetruznychcisel-1," stupnich volnosti = ",chi2,"\nOdpovidajici p - hodnota = ",1-CDF[ChiSquareDistribution[pocetruznychcisel-1],chi2]]

Complement[data,{1}]//Length

DeleteCases[data,1]//Length
fun[DeleteCases[data,1]]

Mathematica styl

Table[Count[data,i],{i,10}]

Map[Count[data,#]&,Union[data]]

Dalsi pouziti funkce Tally

Tally[{a,a,b,a,c,b,a}]

Use test argument to count expressions with the same Head:

Tally[{{a,b},{w,x,y,z},E,a,vvv,11111,{w,x,y,z},2,3,4,5,E},Head[#1]===Head[#2]&]

Tally[data]//Sort

ExampleData["Text"]

text=StringTake[ExampleData[{"Text","Hamlet"}],10000];

thamlet=Tally[Characters[text]]
Sort[Map[Reverse,thamlet]]//Reverse

Substitucni pravidla (replacement rules)

{1,2,3,2,3,1,2,4}/. 3->"a"/.x_Integer->0

v1={1,2,3,2,3,1,2,4}/. 3->"a"/.x_Integer->0/."a"->1

Total[v1]

pravidla[data_List,n_]:=
Total[data/.n->"a"/.x_Integer->0/."a"->1]
pravidla[{1,2,3,2,3,1,2,4},2]

Map[pravidla[data,#]&,Range[10]]
Map[pravidla[data,#]&,Union[data]]

char2=Characters["Map[pravidla[data,#]&,Range[10]] Map[pravidla[data,#]&,Union[data]]'"]

pravidla2[data_List,n_]:=
Count[data,n]
pravidla2[{"a","p","[","p","r","a","v"},"p"]

{Union[char2],Map[pravidla2[char2,#]&,Union[char2]]}//Transpose//TableForm

Map[pravidla[C,#]&,Union[data]]

Zakladni financni funkce s ohledem na zmeny zahrnuty do prednasky 16. 4. 2009

Lekce 8 (16. 4. 2009)

FAQ

Graficke zobrazeni minima nekolika funkci

Clear[f1,f2,fmin];
f1[x_]:=(x^2)
f2[x_]:=Sin[2x]
fmin[x_]:=Min[f1[x],f2[x]]
Plot[fmin[x],{x,-5,5}]

Nalezeni bodu zlomu

FullSimplify[Min[f1[x],f2[x]],x∈Reals]

FullSimplify[f1[x]-fmin[x]]

Reduce[{f1[x]==fmin[x],f2[x]==fmin[x]},x]

fminexpand[x_]=PiecewiseExpand[fmin[x]]

Plot[fminexpand[x],{x,-5,5}]

i1=Reduce[{f1[x]≤fminexpand[x],-5≤x<5},x]
i2=Reduce[{f2[x]≤fminexpand[x],-5≤x<5},x]

forder[x_]:=Ordering[{f1[x],f2[x]}]

Position[forder[2],2]

Reseni pro linearni funkce viz Lekce 9 -> FAQ -> Nalezeni bodu zlomu pro linearni funkce

Nalezeni bodu zlomu pro linearni funkce

Clear[g1,g2,g3];
g1[x_]:=2x+3;g2[x_]:=3x+4;g3[x_]:=-x+5;
Plot[{g1[x],g2[x],g3[x]},{x,-3,6}]
gmin[x_]:=Min[g1[x],g2[x],g3[x]]
gminexpand[x_]=PiecewiseExpand[gmin[x]]
Plot[gmin[x],{x,-3,6}]

(regions=Map[Reduce[(gminexpand[x]==#),x]&,{g1[x],g2[x],g3[x]}])//FullForm
(regions2=regions//.{Inequality->List,LessEqual->List,GreaterEqual->List}//Flatten)//FullForm
crossings=Select[regions2,NumericQ]
Union[crossings]

Problem z bakalarske prace (Celkove doby setrvani ve stavech markovskeho retezce)

Data ve tvaru {{(Subscript[dobasetrvani, 1]),(Subscript[stav, 1])},...,{(Subscript[dobastervani, n]),(Subscript[stav, n])}}

SeedRandom[13];
markcasy=RandomReal[ExponentialDistribution[1],1000];
markstavy=RandomInteger[{1,10},{1000}];
moznestavy=Union[markstavy]
markdata={markcasy,markstavy}//Transpose;

Vybereme vsechna data, kdy je retezec ve stavu 1

(s1=Select[markdata,#[[2]]==1&])//Length

Z toho secteme prvni slozky (nejdrive vybereme prvni slozky)

Map[First,s1] (* ekvivalentne *)
Transpose[s1][[1]]

Total[Map[First,s1]]

Short[(ss=Map[Function[x,Select[markdata,#[[2]]==x&]],moznestavy]),10]

Total[Map[First,ss[[1]]]]

Map[{#,Total[Map[First,ss[[#]]]]}&,moznestavy]//TableForm

Reseni uloh celocoselneho linearniho programovani

FindMinimum[{x+y,x+2y≥3&&x≥0&&y≥0&&y∈Integers},{x,y}]

Prohlednete si funkci FindInstance.

Viz tez tutorial/ConstrainedOptimizationExact #683534434.

Zakladni financni funkce I

Korekce k prednasce 9. 4. 2009 (vyznaceny zelene)

PV: ekvidistatni casove intervaly mezi platbami

Uvodni poznamky

x=.;
NestList[Function[x,(x/(1+r))],c,4]
FoldList[Function[{x,r},(x/(1+r))],c,{(Subscript[r, 1]),(Subscript[r, 2]),(Subscript[r, 3])}]

Predpoklada se, ze cf je penezni tok, cf={(Subscript[cf, 0]),(Subscript[cf, 1]),...,(Subscript[cf, T])}, i hodnotici urokova mira, ii={(Subscript[i, 1]),...,(Subscript[i, T])} promenne urokove miry pro obdobi 1,...,T, rr={(Subscript[R, 1]),...,(Subscript[R, T])} urokove miry pro investice se splatnostmi 1,...,T, (vynosova krivka).

Soucasna hodnota I

Clear[pv];
Options[pv]={interestrates->"constant"};
(* dalsi mozne volby "variable" a "yieldcurve" *)
pv[cf_?VectorQ,i_,opts___?OptionQ]:=Module[{n=Length[cf],interests},
If[n==0,Return[Print["Penezni tok ma nulovou delku."]]];
interests=interestrates/.Flatten[{opts,Options[pv]}];
Which[interests=="constant",cf.(1/(((1+i))^(Range[n]-1))),
	interests== "variable",cf.(1/(FoldList[#1 * (1+#2)&,1,i])),
	interests=="yieldcurve",First[cf]+Rest[cf].(1/(((1+i))^(Range[n-1])))]
	 ]

Options[pv]

cf=.;ii=.;rr=.;i1=.;i2=.;i3=.;
cf={cf0,cf1,cf2,cf3};
ii={i1,i2,i3};
rr={r1,r2,r3};
{pv[1,ii],pv[{},ii],pv[,ii]}

pv[cf,i]//TraditionalForm
pv[cf,ii,interestrates->"variable"]
pv[cf,rr,interestrates->"yieldcurve"]
pv[cf,i,in] (* in nemuze byt uvazovano jako option, vrati tedy pv[... *)

Nejsou sladeny argumenty a volby:

pv[cf,rr,interestrate->"yieldcurve"]

{a,b,c}/.{a->s,a->d}
{a,b,c}//.{a->s,a->d}

Soucasna hodnota II (alternativni postup pri zpracovani voleb)

Clear[pv2];
Options[pv2]={interestrates->"constant"};
(* dalsi mozne volby "variable" a "yieldcurve" *)
pv2[cf_?VectorQ,i_,opts___?OptionQ]:=Module[{n=Length[cf],interests},
If[n==0,Return[Print["Penezni tok ma nulovou delku."]]];
interests=interestrates/.{opts}/.Options[pv2];
Which[interests=="constant",cf.(1/(((1+i))^(Range[n]-1))),
	interests== "variable",cf.(1/(FoldList[#1 * (1+#2)&,1,i])),
	interests=="yieldcurve",First[cf]+Rest[cf].(1/(((1+i))^(Range[n-1])))]
	 ]

pv2[cf,i]//TraditionalForm
pv2[cf,ii,interestrates->"variable"]
pv2[cf,rr,interestrates->"yieldcurve"]
pv2[cf,i,in] (* in nemuze byt uvazovano jako option, vrati tedy pv[... *)
pv2[cf,ii,inter->"variable"]

IRR: ekvidistatni casove intervaly mezi platbami

Nejjednodussi pristup (vyhledani vsech korenu) spociva v pouziti funkce NSolve. Odpovidajici algebraicka rovnice je v sice v promenne v, ale Mathematica rozpozna, ze se jedna o algebraickou rovnici i kdyz je promennou urokova mira i.

Nejprve ukazka projektu, ktery asi zadny potencialni investor nebude povazovat za realny:

i=.
cfnerealny={-1000,3600,-4310,1716};
Plot[pv[cfnerealny,i],{i,0,0.4}]

Vnitrni mira vynosnosti se ziska takto:

i=.
irrsolv=NSolve[pv[cfnerealny,i]==0,i]

i/.irrsolv

i=.
irrsolv=NSolve[pv[cfnerealny,i]==0`40,i]

0.3000000000000002`

Zde je podstatne, ze se pouzije numericka varianta Solve, tj. NSolve. Vsimneme si, co udela obecna procedura pro vyhledani korenu FindRoot (zde je nutne zadat pocatecni reseni x0)

i0=.;i=.;
Table[{i0,FindRoot[pv[cfnerealny,i],{i,i0}]},{i0,0,0.5,0.05}]//TableForm

Nyni priklad nahodneho penezniho toku. Uvadime tez prislusne statisticke charakeristiky

i=.;
SeedRandom[313131];
{Mean[cf13=RandomReal[{-25,33},{30}]],Variance[cf13],"   Smerodatna odchylka = ",StandardDeviation[cf13],Quantile[cf13,0.05],Quantile[cf13,0.95]}
{Min[cf13],Max[cf13]}
ListPlot[cf13,Joined->True]
Plot[pv[cf13,i],{i,-0.18,.5},PlotRange->All]
NSolve[pv[cf13,i]==0,i]

Samozrejme, ze nas nezajimaji imaginarni, a z hlediska ekonomickeho ani zaporna vnitrni vynosova procenta. Separace v nasledujici obecne procedure:

Clear[irr];
(* Predpoklada vyse definovanou funkci pv *)
irr[cf_]:=Module[{solutions,i},
solutions=i/.NSolve[pv[cf,i]==0,i];
Select[solutions,(Im[#]==0&&Re[#]>0)&]
				]

irr[cf13]

Alternativne

Cases[i/.NSolve[pv[cf13,i]==0,i],i_Real]
Cases[i/.NSolve[pv[cf13,i]==0,i],i_Real/;i≥0]

Cviceni

Vyzkousejte dalsi zabudovane funkce k ziskani stejneho vysledku.

Konec 16. 4. 2009

Lekce 9 (23. 4. 2009)

Uzitecny odkaz

http://demonstrations.wolfram.com/BinomialTree/

Obecne poznamky k presnosti vypoctu a FAQ

Obecna poznamka k presnosti I - ilustrace na vypoctu IRR

Clear[irrsimple,i];
(* Predpoklada vyse definovanou funkci pv *)
irrsimple[cf_]:=Module[{solutions,i},
solutions=i/.NSolve[cf.(1/(((1+i))^(Range[Length[cf]]-1)))==0,i];
Select[solutions,(Im[#]==0&&Re[#]>0)&]
				]
irrsimple[cfnerealny={-1000,3600,-4310,1716}]

Zkopirovana hodnota 0.3:

0.3000000000000002`

irrsimple[cfnerealny2={-1000.05,3600,-4310,1716}]

Zkopirovana prvni hodnota

0.294064372011615`

Avsak presnost vypoctu na 30 mist dosahneme

irrsimple[cfnerealny3={-1000.05`30,3600,-4310,1716}]

Zkopirovana prvni hodnota

0.294064372011614836801956765663300974059618719095`29.69897000433602

Obecna poznamka k presnosti II - vypocet funkce sin komplexniho argumentu

Sin[4.5`100+I]

Obecna poznamka k presnosti III - hledani pevneho bodu zobrazeni (porovnani s programem Excel)

3 Try in MS Excel: 0.2->C2. Then copy formula 11*C2-2 15 times. This results (dependent on the used computer) in:

0.20.20.20.20.20.20.20.20.20.20.20.2000050.2000510.2005570.2061320.2674570.9420288.36231289.98543987.839710864.24

4 The same in Mathematica is solved as follows:

fjo[x_]:=11*x-2
NestList[fjo,0.2,30]

Obviously, in rationals it works:

NestList[fjo,(1/5),30]

Look at

BaseForm[1/5,2]
BaseForm[0.2,2]
BaseForm[0.2`30,2]

But (e.g.):

fjo30=NestList[fjo,0.2`30,30]

Nalezeni bodu zlomu pro linearni funkce (z minule lekce)

Priklad 1

Clear[g1,g2,g3];
g1[x_]:=2x+3;g2[x_]:=3x+4;g3[x_]:=-x+5;
Plot[{g1[x],g2[x],g3[x]},{x,-3,6}]
gmin[x_]:=Min[g1[x],g2[x],g3[x]]
gminexpand[x_]=PiecewiseExpand[gmin[x]]
Plot[gmin[x],{x,-3,6}]

(regions=Map[Reduce[(gminexpand[x]==#),x]&,{g1[x],g2[x],g3[x]}])//FullForm
regions
(regions2=regions//.{Inequality->List,LessEqual->List,GreaterEqual->List}//Flatten)//FullForm
crossings=Select[regions2,NumericQ]
Union[crossings]

Priklad 2 (nahodne linearni funkce)

Clear[garray, x]; (* generuje pole n linearnich funkci *)
garray[n_]:=RandomInteger[{-150,150},n] x+RandomInteger[{-150,50},n]
garray[3]

SeedRandom[1313131]
g23=garray[23]
Plot[g23,{x,-30,30}]
Plot[g23,{x,-3,3}]

garrayminimum[x_]=Min[g23]
garrayminimumpiecewise[x_]=PiecewiseExpand[garrayminimum[x]]
Plot[garrayminimumpiecewise[x],{x,-10,10},PlotPoints->500]

Graphics[Line[{{-19/5,0},{-19/5,(-81+115(-19/5))}}]]

Show[Plot[garrayminimumpiecewise[x],{x,-10,10},PlotPoints->500],Graphics[Line[{{-19/5,0},{-19/5,(-81+115x)/.x->(-19/5)}}]]]

(regionsarray23=Map[Reduce[(garrayminimumpiecewise[x]==#),x]&,g23])//FullForm

(regions2array23=regionsarray23//.{Inequality->List,LessEqual->List,GreaterEqual->List}//Flatten)//FullForm

crossingsarray23=Select[regions2array23,NumericQ]
Union[crossingsarray23]

Soucty pri podminkach kladenych na indexy

rr=.;k=.;aa=.;
rr=RandomChoice[Range[10],10]
k=4;
aa=Array[a,{10}]
Sum[a[i]*Boole[rr[[i]]≥k],{i,Length[aa]}]

Boole[rr≥k]
Thread[rr≥k]
Boole[Thread[rr≥k]]
Total[aa*Boole[Thread[rr≥k]]]

aa[[2]]

Nahodny vyber s navracenim

SeedRandom[13];
rwith=RandomChoice[Range[12],10]
Length[rwith]==Length[Union[rwith]]

Nahodny vyber bez navraceni

SeedRandom[13];
rwithout=RandomSample[Range[12],10]
Length[rwithout]==Length[Union[rwithout]]

Zjistit, zda dana mnozina je permutace

Zabudovana funkce PermutationQ

Needs["Combinatorica`"]
{PermutationQ[{3,1,2}],PermutationQ[{3,1,1}]}

perQ[x_?VectorQ]:=Sort[x]==Range[Length[x]]
{perQ[{2,1,3}],perQ[{2,2,1}],perQ[{{1},2,3}]}

per2Q[x_,y_]:=Sort[x]===Union[y]
{per2Q[{a,c,b},{c,b,a}],per2Q[{a,c,b},{c,c,b,a}],per2Q[{a,a,b},{a,b,c}]}

Attributes[Plus]

Attributes[Thread]

Opakovane pouziti Thread

Thread[aa==(aa==3)]
Thread[Thread[aa==3]]

Thread[aa==(aa==3)]
Thread[Thread[aa==(aa==3)]]

Nalezeni zakonitosti v posloupnostech

FindSequenceFunction[{1,1,2,3},n]
FindSequenceFunction[RandomInteger[{0,1},5],n]

Table[(1/2) (1+(((-1))^n)),{n,5}]

FindSequenceFunction[{1,1,2,4,3,9,4,16},n]

FindSequenceFunction[{1,3,4,7,11,18},n]

Zvyrazneni bodu a spojovani v DateListPlot

Show[DateListPlot[{1,1,2,3,5,8,11},{2000,8},PlotStyle->PointSize[0.04]],DateListPlot[{1,1,2,3,5,8,11},{2000,8},Joined->True]]

DateListPlot[{1,1,2,3,5,8,11},{2000,8},PlotStyle->PointSize[0.04],Joined->True,Mesh->All]

Zakladni financni funkce II

Durace (ekvidistantni casove intervaly, bez overovani vstupu)

Clear[pv];
Options[pv]={interestrates->"constant"};
(* dalsi mozne volby "variable" a "yieldcurve" *)
pv[cf_?VectorQ,i_,opts___?OptionQ]:=Module[{n=Length[cf],interests},
If[n==0,Return[Print["Penezni tok ma nulovou delku."]]];
interests=interestrates/.Flatten[{opts,Options[pv]}];
Which[interests=="constant",cf.(1/(((1+i))^(Range[n]-1))),
	interests== "variable",cf.(1/(FoldList[#1 * (1+#2)&,1,i])),
	interests=="yieldcurve",First[cf]+Rest[cf].(1/(((1+i))^(Range[n-1])))]
	 ]

Clear[dur];
Options[dur]={interestrates->"constant"};
(* dalsi mozne volby "variable" a "yieldcurve" *)
dur[cf_?VectorQ,i_,opts___?OptionQ]:=Module[{n=Length[cf],interests},
If[n==0,Return[Print["Penezni tok ma nulovou delku."]]];
interests=interestrates/.Flatten[{opts,Options[pv]}];
Which[interests=="constant",((cf*(Range[n]-1)).(1/(((1+i))^(Range[n]-1))))/pv[cf,i,opts],
	interests== "variable",((cf*(Range[n]-1)).(1/(FoldList[#1 * (1+#2)&,1,i])))/pv[cf,i,opts],
	interests=="yieldcurve",((Rest[cf]*Range[n-1]).(1/(((1+i))^(Range[n-1]))))/pv[cf,i,opts]]
	 ]

cf=.;ii=.;rr=.;i1=.;i2=.;i3=.;
cf={cf0,cf1,cf2,cf3};
ii={i1,i2,i3};
rr={r1,r2,r3};
{dur[1,ii],dur[{},ii],dur[,ii]}

dur[cf,i]//TraditionalForm
dur[cf,ii,interestrates->"variable"]
dur[cf,rr,interestrates->"yieldcurve"]
dur[cf,i,in] (* in nemuze byt uvazovano jako option, vrati tedy dur[... *)

Zakladni financni funkce III: Fair values of bonds and durations of bonds

Fair = dirty = hruba = spravedliva,
Pure = cista 
c...kuponova sazba!!!
n...pocet obdobi do splatnosti

Present value of a coupon bond

(∑_(j=1)^n)(v^j)//TraditionalForm

Clear[dirtyvalue,purevalue,v];
dirtyvalue[c_,f_,n_,i_]:=c f (v^(n-Floor[n])) (∑_(j=0)^(Floor[n]))(v^j)+f (v^n)/.v->(1/(1+i))
purevalue[c_,f_,n_,i_]:=f ((v^n)+((c v ((v^n)-1))/(v-1)))/.v->(1/(1+i))

Kuponova obligace s kuponovou sazbou c=0.1, splatnost T=10, nominalni hodnota F=1, hodnotici urokova sazba i=0.15, tj. diskontni obligace.

p12=Plot[{dirtyvalue[0.1,1,10-n,0.15],purevalue[0.1`,1,10-n,0.15]},{n,0,10},PlotPoints->500,Exclusions->None,AxesLabel->{t,PV}]

Kuponova obligace s kuponovou sazbou c=0.1, splatnost T=10, nominalni hodnota F=1, hodnotici urokova sazba i=0.1.

p3=Plot[dirtyvalue[0.1`,1,10-n,0.1],{n,0,10},Exclusions->None,PlotStyle->{{Dashing[{0.01`,0.01`}],Black}},AxesLabel->{t,PV}]

Kuponova obligace s kuponovou sazbou c=0.1, splatnost T=10, nominalni hodnota F=1, hodnotici urokova sazba i=0.05, tj. premiova obligace.

p45=Plot[{dirtyvalue[0.1,1,10-n,0.05],purevalue[0.1,1,10-n,0.05]},{n,0,10},PlotPoints->500,Exclusions->None,AxesLabel->{t,PV}]

Show[p12,p3,p45,PlotRange->All]

Manipulate[Plot[{dirtyvalue[c,f,tT-t,i],purevalue[c,f,tT-t,i],purevalue[c,f,tT-t,c]},{t,0,tT},PlotPoints->500,Exclusions->None,AxesLabel->{t,PV},PlotRange->{{0,tT},{0,2f}}],
{{tT,10,"T"},5,30,1},
{{c,0.1,"c"},0.01,0.20,0.01},
{{f,1,"F"},1,1000,50},
{{i,0.15,"i"},0.01,0.20,0.01}
]

Duration of a coupon bond

((c f(∑_(j=1)^n)j (v^j)+f n (v^n))/(c f(∑_(j=1)^n) (v^j)+f  (v^n)))/.v->(1/(1+i))//PowerExpand//FullSimplify

Clear[pureduration,dirtyduration,v];
pureduration[c_,f_,n_,i_]:=1+(1/i)-((1+i+c n-i n)/(i+c (-1+(((1+i))^n))))
dirtyduration[c_,f_,n_,i_]:=(c f (∑_(j=0)^(Floor[n])) (n-Floor[n]+j) (v^(n-Floor[n]+j))+f n (v^n))/(c f (v^(n-Floor[n])) (∑_(j=0)^(Floor[n]))(v^j)+f (v^n))/.v->(1/(1+i))

Plot[{dirtyduration[0.1,1,10-n,0.15],pureduration[0.1,1,10-n,0.15]},{n,0,10},PlotPoints->500,Exclusions->None]

Plot[{dirtyduration[0.1,1,10-n,0.05],pureduration[0.1,1,10-n,0.05]},{n,0,10},PlotPoints->500,Exclusions->None]

PV, FV, D, C, IRR, nestejne intervaly mezi platbami

Measures of cash flows

Calendar

Calendar conventions play a crucial role in time value of money. There is a variety of calendar conventions for  financial instruments. Within Mathematica framework, they are included in Applications Packages like Financial Essentials or Derivatives Expert. They must be bought separately, however. In Mathematica Add-On Packages the only available calendar is Actual/Actual, i.e., actual number of days in a month and actual number of days in a year. For the illustration of the financial calculations we will make use of this one.

The actual date in the form {Year, Month, Day, Hour, Minute, Second}:

Date[]

We may extract just the date:

Date[]
Date[][[Range[3]]]

Needs["Calendar`"]

DayOfWeek[today={2002,6,28}]

{DaysBetween[today,{2010,3,28}],DaysBetween[{2010,3,28},today]}

DaysPlus[today,50]

Present Value, calendar Actual/360

For a cash flow CF = ((Subscript[CF, 0]),...,(Subscript[CF, T]))and a  valuation interest rate r the formula for the present value reads:

PV(CF,r)=(∑_(t=0)^T)((Subscript[CF, t])/(((1+r))^t))

In fact, this is the scalar product of the cash flow and the vector of discount factors. Implementation heavily uses listability:

cf={cf0,cf1,cf2,cf3};
len=Length[cf]
disc=1/(1+r)^(Range[len]-1)
cf.disc

Clear[pv1];
pv1[cf_,r_]:=cf.(1/(1+r)^(Range[Length[cf]]-1))
pv1[cf,r]

Modification for the payments going on given different dates and we evaluate the present value with respect to the startday:

dif=Map[DaysBetween[today,#]&,{{2004,3,28},{2003,1,13},{2003,1,1}}]
discounts=1/((1+r)^(dif/360))

Clear[pv2];
pv2[cf_,dates_,r_,startday_]:=Module[{diskonty,dif},
dif=Map[DaysBetween[startday,#]&,dates];
diskonty=1/((1+r)^(dif/360));
cf.diskonty]

dates={{2004,3,28},{2003,1,13},{2003,1,1}};
cf={1000,-900,-50};
pv2[cf,dates,0.1,today]

Plot[pv2[cf,dates,r,today],{r,0,0.1}]

IRR (Internal Rate of Return)

cf2={-1000,Table[100,{10}],1100}//Flatten;

NSolve[pv1[cf2,r]==0,r]

Cases[r/.NSolve[pv1[cf2,r]==0,r],z_Real/;z≥0][[1]]

dates2=Map[DaysPlus[today,#]&,Table[365*(j-1),{j,Length[cf2]}]]

Plot[pv2[cf2,dates2,r,today],{r,0.09,0.11}]

FindRoot[pv2[cf2,dates2,r,today]==0,{r,0.01}]

Konec 23. 4. 2009

Lekce 10 (30. 4.2009)

Manualy k versi 7.0

http://url.wolfram.com/6eVW8E

FAQ

Mnohorozmerne normalni rozdeleni

Needs["MultivariateStatistics`"];
pdf[xx_,μμ_,Σ_]:=PDF[MultinormalDistribution[μμ,Σ],xx]
pdf[Array[x,5],Array[μ,5],Array[σ,{5,5}]+Transpose[Array[σ,{5,5}]]]

Nacitani dat z Excelu

$ImportFormats

SetDirectory[NotebookDirectory[]]

Import["pokusnadataExcel.xls"]//TableForm

({{({{"a"}, {"b"}, {"a"}, {"b"}, {"a"}, {"b"}}), ({{"b"}, {1.`}, {2.`}, {3.`}, {4.`}, {5.`}}), ({{"a"}, {2.`}, {3.`}, {4.`}, {5.`}, {6.`}}), ({{"b"}, {3.`}, {4.`}, {5.`}, {6.`}, {7.`}}), ({{"a"}, {4.`}, {5.`}, {6.`}, {7.`}, {8.`}}), ({{"b"}, {5.`}, {6.`}, {7.`}, {8.`}, {9.`}}), ({{"a"}, {6.`}, {7.`}, {8.`}, {9.`}, {10.`}}), ({{"b"}, {7.`}, {8.`}, {9.`}, {10.`}, {11.`}}), ({{"a"}, {8.`}, {9.`}, {10.`}, {11.`}, {12.`}}), ({{"b"}, {9.`}, {10.`}, {11.`}, {12.`}, {13.`}}), ({{"a"}, {10.`}, {11.`}, {12.`}, {13.`}, {14.`}}), ({{"b"}, {11.`}, {12.`}, {13.`}, {14.`}, {15.`}})}})

Poznamky

Poznamka k definicim (globalni pravidla)

a=Sin
{a[π],a[3(π/2)],a[x],3a}
?a

a[x]+Sin[x]

Poznamky ke grafice

optionsPlotbutnotListPlot=Complement[Options[Plot],Options[ListPlot]]
optionsListPlotbutnotPlot=Complement[Options[ListPlot],Options[Plot]]

a=.;x=.;
FullForm[x->a]
x->a//Last

Map[First,optionsPlotbutnotListPlot]

Poznamky k analytickemu zjednoduseni

i1=Integrate[v1=(Sqrt[x]) ArcTan[x],x]//FullSimplify
d1=D[i1,x]
d1==v1
d1===v1
d1//Simplify
(d1//Simplify)==v1

Clear[k,n];
soucet=(∑_(k=1)^n)(1/(k^6))
Limit[soucet,n->∞]
(∑_(k=1)^∞)(1/(k^6))

soucin=(∏_(n=1)^5)Gamma[((2n)/5)]
soucin//Simplify
soucin//FullSimplify

Series[n!,{n,∞,3}]//PowerExpand//FullSimplify

Series[n!,{n,∞,3}]//PowerExpand//FullSimplify//Normal

((Cos[x])^4)//FullSimplify
TrigReduce[((Cos[x])^4)]

Diferencialni rovnice symbolicky

eqdif=y'[x]+x y[x]==x Exp[-(x^2)]
DSolve[eqdif,y[x],x]
DSolve[{eqdif,y[0]==0},y,x]

Zabudovani pocatecni podminky pridanim dalsi rovnice

dsol=DSolve[{eqdif,y[0]==0},y[x],x]

yx=y[x]/.dsol[[1]]
(y[x]/.dsol)//First

?yx

y1[x_]:=yx
y1[z]

y2[x_]=yx
y2[z]

yx/.x->z

y3[x_]=y[x]/.dsol[[1]]
y3[z]

Diferencialni rovnice numericky

ndsol=NDSolve[{x''[t]+((x[t])^3)==Sin[t],x[0]==x'[0]==0},x,{t,0,50}]
FullForm[ndsol]
ParametricPlot[{x[t],x'[t]}/.ndsol,{t,0,50}]

Kvantilova funkce normalniho rozdeleni

Quantile[NormalDistribution[0,1],y]//InputForm

Generovani nahodnych cisel z velmi zvlastnich rozdeleni (prevzato z MathGroup Archive: December 1998 [00243], Paul Abbott)

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

After some 93 seconds we get that pdf is real 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

{"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}]

Use NDSolve (and dynamic programming) instead of NIntegrate tp 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}]]]

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 adaptivecode built into Plot):

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

ListPlot[pts]

and then numerically compute the inverse function using interpolation:

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

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

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

Zvyrazneni bodu a spojovani v DateListPlot II (uzitecna myslenka i v jinych pripadech)

data1={{{2006,10,1,0,0,0},10},{{2006,10,15,0,0,0},12},{{2006,10,30,0,0,0},15}, {{2006,11,20,0,0,0},20}};
data2={{{2006,10,5,0,0,0},15},{{2006,10,20,0,0,0},8},{{2006,11,10,0,0,0},5}, {{2006,11,15,0,0,0},1}};
DateListPlot[{data1,data2},Joined->{True,False},PlotStyle->PointSize[0.04]]
(* Pouziti na jednu mnozinu dat *)
DateListPlot[{data1,data1},Joined->{True,False},PlotStyle->PointSize[0.04]]

Complement[Options[ListPlot],Options[DateListPlot]]

<<PlotLegends`
DateListPlot[{data1,data1},Joined->{True,False},PlotStyle->PointSize[0.04],Frame->False,AxesLabel->{"Sine","Cosine"}]

Options[DateListPlot]

Petrikuv problem

tabDP=Table[NMaxValue[{Exp[-n(n-1)/(2m)]-(∏_(k=1)^n)(1-(k-1)/m),m≥n},m],{n,2,100}]

ListPlot[tabDP]

Ocividne spatne!!! Chyba ve funkci NMaxvalue.

(∏_(k=1)^n)(1-(k-1)/m)

Clear[m,n,approxDP];approxDP[m_,n_]:=Exp[-n(n-1.`30)/(2m)]-(((-(1/m)))^n) Pochhammer[-m,n]

t40=Table[{m,approxDP[m,40]},{m,40,1000}];
maxt40=Max[t40[[All,2]]]
Position[t40[[All,2]],maxt40]
ListPlot[t40]

Popisy grafu

Needs["PlotLegends`"]
Plot[{Sin[x],Cos[x]},{x,0,2 π},BaseStyle->{FontFamily->"Tahoma",FontSize->14},PlotLabel->"Graf funkce Cos(x) a Sin(x)",PlotStyle->{Hue[0.3`],Hue[0.6`]},PlotLegend->{"Sin(x)","Cos(x)"},LegendPosition->{1.2`,-0.5`},LegendShadow->{0.03`,-0.03`},LegendLabel->"Legenda grafu",LegendOrientation->Horizontal,Epilog->{Arrow[{{3.5,0.8},{2.5,Sin[2.5]}}],Text["Sin(x)",{3.5`,0.8`},{0,-1}],Arrow[{{1,-0.75`},{(π/2),Cos[(π/2)]}}],Text["Cos(x)",{1,-0.9`},{0,-1}]}]

Balicek (Package): zabudovane balicky (Add-On Packages, StandardExtraPackages)

Cim novejsi verze, tim vice funkci je primo v jadru (kernel). Nektere jsou vsak mimo jadro a musi se vyvolat pomoci funkci Get (Unix verse <<), resp. Needs. Velmi citlive na syntaxi!!! V nejnovejsi versi viz (za cviceni zkuste vytvorit odkaz na Documentation Center s timto heslem)

guide/StandardExtraPackages

{$UserDocumentsDirectory,$Path} (* Poznamenejme, ze $ je povazovano za pismeno, v Mathemaice oznacuje systemove promenne a vnitrni promenne modulu. *)

Directory[]
SetDirectory["C:\\Program Files\\Wolfram Research\\Mathematica\\6.0\\AddOns\\Packages"];
Directory[]

FileNames[]

SetDirectory[$UserDocumentsDirectory]

SetDirectory[NotebookDirectory[]]

SetDirectory["C:\Program Files\Wolfram Research\Mathematica\7.0\AddOns\ExtraPackages"]
FileNames[]

SetDirectory["C:\Program Files\Wolfram Research\Mathematica\7.0\AddOns\Packages"]
FileNames[]

Help: UnitsPackage

Needs["Units`"]
units=Names["Units`*"]

Copied from Help (Units of volume): results in Input form. UnitsJH added.

unitsJH=({{Bag, Barrel}, {BoardFoot, Bucket}, {Bushel, Butt}, {Cord, Cup}, {Drop, Ephah}, {Fifth, Firkin}, {FluidDram, FluidOunce}, {Gallon, Gill}, {Hogshead, Jeroboam}, {Jigger, Last}, {Liter, Magnum}, {Minim, Noggin}, {Omer, Pint}, {Pony, Puncheon}, {Quart, RegisterTon}, {Seam, Shot}, {Stere, Tablespoon}, {Teaspoon, Tun}, {UKGallon, UKPint}, {WineBottle, \[Null]}})

There is the comma after "WineBottle". This command will eliminate it and, morevoer flatten the buttons in an one-dimensional list:

unitsJH=unitsJH//Flatten//Most

Methodology

unitsJH[[1]]//Head
unitsJH[[1]]//FullForm
f1=Apply[List,unitsJH[[1]]]
f1[[1]]//FullForm

unitsJH2=Map[First,Map[Apply[List,#]&,unitsJH]]

Conversion

Convert[Bag,Liter]
Convert[Liter,Pony]
Convert[Pony,Liter]

Transpose[{unitsJH2,Map[Convert[#,Liter]&,unitsJH2]}]//TableForm

Vetsina nasledujiciho vystupu je nesmyslna. Proto provedem selekci nasledujicim prikazem.

(konverzevsechjednoteknalitr={units//ToExpression,Map[Convert[#,Liter]&,units//ToExpression]}//Transpose)//Short

Select[konverzevsechjednoteknalitr,#[[1]]=!=#[[2]]&]//TableForm

Needs["ANOVA`"]
Names["ANOVA`*"]

?StudentNewmanKeuls

Needs["RegressionCommon`"]
Names["RegressionCommon`*"]

?AdjustedRSquared

Balicek (Package): struktura a tvorba

Pripona souboru Package je *.m. Muze byt matouci, MATLAB ma to same. Struktura viz priklad

Na siti vyhledejte Skeleton.m. Odkazy viz Maeder: Programming in Mathematica

Kde muze Mathematica najit balicek?

Tri moznosti: (i) tam, kam vede cesta, (ii) pridat cestu (iii) tam, odkud je vyvolana seance (spusten notebook). (iii) je nejjednodussi.

$Path (* (i) sem vedou cesty *)

AppendTo[$Path,"C:\\Documents and Settings\\hurt"] (* (ii) Pridatcestu *)

?*Directory*

(* iii *)
Directory[]
SetDirectory[NotebookDirectory[]]
Directory[]

Priklad IPM (An Introduction to Programming with Mathematica)

Autori:
Paul R. Wellin
Organization: Wolfram Research, Inc.

Richard J. Gaylord 
Organization: University of Illinois at Urbana-Champaign
Department: Department of Material Science and Engineering

Samuel N. Kamin

http://library.wolfram.com/infocenter/Books/5169/#downloads

SetDirectory[NotebookDirectory[]]
<<RandomWalks`

?ShowWalk

walk1D[n_] := NestList[#1 + (-1)^Random[Integer] & , 0, n]

walk1D[10]

RandomWalk[10,Dimension->2]

Package RandomWalks

BeginPackage["IPM3`RandomWalks`"]
Needs["Graphics`ArgColors`"]
Needs["Graphics`Polyhedra`"]
RandomWalk::"usage"="RandomWalk[n] generates an n-step walk in two dimensions. The default behavior gives a lattice walk with steps in one of the four compass directions. Setting the option LatticeWalk->False will generate an off-lattice walk. Setting the option Dimensions->3 will generate a 3-dimensional walk; Dimensions->1 gives a 1-dimensional walk.";

LatticeWalk::"usage"="LatticeWalk->val is an option to RandomWalk that determines whether the random walk will be a lattice walk or an off-lattice walk. Possible values are True and False.";

ShowWalk::"usage"="ShowWalk[walk] displays a two-dimensional random walk connecting each site with a line. Any Graphics options can be passed to ShowWalk. E.g., ShowWalk[walk,  Background->GrayLevel[0]] to produce a black background.";

AnimateWalk::"usage"="AnimateWalk[walk, opts] creates an animation of a two-dimensional random walk. A red ball will be seen to move to the current position in the walk to aid in visualizing the animation.";
RandomWalk::rwn="Argument `1` is not a positive integer.";
Options[RandomWalk]={LatticeWalk->True,Dimension->2}
Begin["`Private`"]
walk1D[n_] := NestList[#1 + (-1)^Random[Integer] & , 0, n]

walk1DOffLattice[n_]:=
  FoldList[Plus,0,Table[Random[Real,{-1,1}],{n}]]

walk2D[n_]:=
  Module[{NSEW={{0,1},{1,0},{0,-1},{-1,0}}},
    FoldList[Plus,{0,0},
      NSEW[[
        Table[Random[Integer,{1,4}],{n}]]]]]

walk2DOffLattice[n_]:=
  FoldList[Plus,{0,0},
    Map[{Cos[#],Sin[#]}&,Table[Random[Real,{0,2 π}],{n}]]]

walk3D[n_] := Module[{NSEW3 = Sqrt[2]*Vertices[Cube]}, 
   FoldList[Plus, {0, 0, 0}, 
    NSEW3[[Table[Random[Integer, {1, 8}], {n}]]]]]

walk3DOffLattice[n_] := Module[{x, y, z}, 
   x[ϕ_, θ_] := Sqrt[1 - Cos[ϕ]^2]*Cos[θ]; 
    y[ϕ_, θ_] := Sqrt[1 - Cos[ϕ]^2]*Sin[θ]; 
    z[ϕ_] := Cos[ϕ]; FoldList[Plus, {0, 0, 0}, 
     Apply[{x[#1, #2], y[#1, #2], z[#1]} & , 
      Table[Random[Real, {0, 2*Pi}], {n}, {2}], {1}]]]

RandomWalk[n_,opts___?OptionQ]:=Module[{dim,latticeQ},
    If[Not[IntegerQ[n]&&n>0],
      Message[RandomWalk::rwn,
        n],{latticeQ,
          dim}={LatticeWalk,Dimension}/.Flatten[{opts,Options[RandomWalk]}];
      Which[
        dim==1,If[latticeQ,walk1D[n],walk1DOffLattice[n]],
        dim==2,If[latticeQ,walk2D[n],walk2DOffLattice[n]],
        dim==3,If[latticeQ,walk3D[n],walk3DOffLattice[n]]
        ]]]

ShowWalk[coords_,opts___]:=Which[
    Length[Dimensions[coords]]==1,
    ListPlot[coords,opts,PlotJoined->True],
    Dimensions[coords][[2]]==2,
    Show[Graphics[Line[coords],opts,AspectRatio->Automatic]],
    Dimensions[coords][[2]]==3,
    Show[Graphics3D[Line[coords],opts,AspectRatio->Automatic]]]
AnimateWalk[coords_,opts___]:=
  Scan[Show[
        Graphics[{{RGBColor[1,0,0],PointSize[0.025],
              Point[coords[[#1]]]},
            Line[Take[coords,#1]]}],opts,AspectRatio->Automatic,
        PlotRange->({Min[#1]-0.2,Max[#1]+0.2}&)/@Transpose[coords]]&,
    Range[2,Length[coords]]]
End[]
EndPackage[]

Vicerozmerna nahodna prochazka (m-rozmerna)

Table[(((-1))^(RandomInteger[{0,1},4])),{5}]//TableForm

Clear[walkJH];
walkJH[m_(* dimenze *),n_(* delka *)]:=NestList[#1+(((-1))^(RandomInteger[{0,1},m]))&,Table[0,{m}],n]

walkJH[1,10]

ListPlot[walkJH[1,10000]//Flatten]

ListPlot[walkJH[2,10000],Joined->True]

ListPlot3D[walkJH[3,10000]]

Konec 30. 4. 2009

Lekce 11 (7. 5. 2009)

V pristi prednasce se bude hovorit o opcich, binomickem modelu hodnoceni opci a teorii optimalniho portfolia. Pripomente si prislusne partie z prednasek Financni management.

Pomocne odkazy

Obecne

http://reference.wolfram.com/mathematica/tutorial/StatisticalModelAnalysis.html

http://reference.wolfram.com/mathematica/ref/LinearModelFit.html

http://demonstrations.wolfram.com/ExpectedReturnsOfTheDowIndustrialsBetaModel/

http://demonstrations.wolfram.com/ExpectedReturnsOfTheDowIndustrialsFamaFrenchModel/

http://mba.tuck.dartmouth.edu/pages/faculty/ken.french/data_library.html

tutorial/UnconstrainedOptimizationOverview

Jan Hurt

2005 Discriminant Analysis and Multidimensional Scaling: http://library.wolfram.com/infocenter/Conferences/5777/
2006 The Kermack-McKendrick Model of Epidemics and Its Extensionshttp://library.wolfram.com/infocenter/Conferences/6462/
2007 Yield Curves: http://library.wolfram.com/infocenter/Conferences/6956/
2008 Risk Measures in Finance: http://library.wolfram.com/infocenter/Conferences/7230/

Poznamky

Poznamka 1. V lekci 10 opraven odstavec Popisy grafu

Poznamka 2. K vlastnim definicim

a=Sin;
{a==Sin,a===Sin}

Poznamka 3. Jak zaridit, abychom stale nemuseli psat napr. Joined -> True?

Options[ListPlot]

j=Joined->True;
ListPlot[Sin[Range[100]],Joined->False]

ListPlot[Sin[Range[100]],j]

ListPlot[{1,2,3}]

listPlotJH[x_]:=ListPlot[x,Joined->True]
listPlotJH[{1,2,3}]
ListPlot[{1,2,3}]

jj=Sequence[Joined->True,AxesLabel->{x,y},AxesStyle->Directive[Orange,12]]
pp=ListPlot[Sin[Range[100]],jj]

Poznamka 4. (Options v grafu)

U vetsiny, zvlaste grafickych funkci, je siroka moznost voleb. To, co je mozne volit se ziska uzitecnym prikazem (uvadime na konkretnim priklade):

Toto ocividne nedava kruznici:

Plot[{-(Sqrt[1-(x^2)]),(Sqrt[1-(x^2)])},{x,-1,1}]

Options[Plot]

Musime zmenit AspectRatio. Takova ponekud nasilna zmena vsak neni obvykla.

Plot[{-(Sqrt[1-(x^2)]),(Sqrt[1-(x^2)])},{x,-1,1},AspectRatio->Automatic]

FAQ

Misto Erf vyskytujici se ve vzorcich s normalnim rozdelenim pozadujeme Φ.

Φ[x_]=CDF[NormalDistribution[0,1],x]
Φ[x]
Remove[Φ]
{(1/2) (1+Erf[(x/(Sqrt[2]))]),(2^(Erf[y]))+Erf[π]}/.Erf[t_]->2Φ[t(Sqrt[2])]-1

Pokud vsak pozadujeme konkretni hodnoty, musime Φ aktivovat, a ta to zase prevede na Erf.

Vypocet funkcni hodnoty ve zretenych promennych

Clear[petra];
petra[cf_,i_,t_]:=(cf^i)+t
Map[Function[x,petra[cf,x,t]],{i1,i2,i3}]

Map[Function[y,Map[Function[x,petra[cf,x,y]],{i1,i2,i3}]],{t1,t2,t3,t4}]

Vystup vzorce v TeXu

PDF[NormalDistribution[μ,σ],x]//TeXForm

Jeste ke grafum

Show se ridi direktivami u prvniho grafu

{Show[pp,p1,p2],Show[p1,pp,p2]}

Implicitni hodnoty volitelnych promennych

Clear[f,a,b];
f[x_,y_:3]:={x,y}
{f[a,b],f[a]}

Testy v definici funkce

Clear[f];
f[x_?Positive]:=Log[x]
{f[-1],f[1]}

Clear[f]; (* f vraci x pokud je x delitelne tremi *)
f[x_?(Mod[#,3]==0&)]:=x
Map[f[#]&,Range[10]]
Map[Mod[#,3]==0&,Range[7]]

Statisticke funkce

TUTORIALS
Basic Statistics
Continuous Distributions
Discrete Distributions
Descriptive Statistics
Convolutions and Correlations
Curve Fitting
Statistical Model Analysis
Importing and Exporting Data
Pseudorandom Numbers
Random Number Generation

Statistika: ukazky

Dokumentace viz heslo "Statistics" (Help -> Documentation Center). Vetsina funkci, ktere byly puvodne pouze v Standard Add Ons Packages jsou nyni funkce zabudovane. Ilustrace na Poissonove rozdeleni. Pro dalsi kliknuti na modre polozky je nutne Menu: "Evaluation   ->   Dynamic Evaluation Enabled"!

?Poi*

PDF (Power Density Function) = hustota 
CDF (Cumulative Distribution Function) = distribucni funkce

λ=2;k=7;
λ=.;k=.;
PDF[distr=PoissonDistribution[λ],k]
distrGauss=NormalDistribution[μ,σ]

{CDF[distr,k],CDF[distrGauss,x],Quantile[distrGauss,q]}

Cviceni: Jaky je vztah mezi Erf a Φ?

Momenty

{Mean[distr],Variance[distr],Skewness[distr],Kurtosis[distr],Kurtosis[distrGauss]}

ExpectedValue[(x^4),distr,x]
ExpectedValue[((#)^4)&,distr]
ExpectedValue[Function[c,(c^4)],distr]
ExpectedValue[((#)^(4/3))&,distr]

Table[{k,ExpectedValue[(((x-λ))^k),distr,x]},{k,6}]//TableForm

Obecne podrobnosti: Help -> DiscreteDistributions a ContinuousDistributions

Nahodna cisla z Poissonova rozdeleni s parametrem λ

RandomInteger[distr,{4}]

Neni mozne generovat nahodna cisla s obecnym parametrem. Obecny parametr nutne speciikovat:

RandomInteger[distr/.λ->3,{4}]

To sice funguje, ale lepsi je prehlednejsi

tabPois1=RandomInteger[PoissonDistribution[3],{4}]

Jak videt, tato 4 nahodna cisla jsou jina nez v predchozim pripade. Po kazdem novem vyvolani dostaneme jedinecnou posloupnost. Pomoci nasledujici konstrukce je mozne zajistit reprodukovatelnost nahodnych posloupnosti. Vrele se doporucuje!!! Ten, kdo vystupy sleduje si muze nahodne experimenty overit.

SeedRandom[13];
tabPois2=RandomInteger[PoissonDistribution[3],{4}]
SeedRandom[13];
tabPois3=RandomInteger[PoissonDistribution[3],{5}]
tabPois4=Table[Random[PoissonDistribution[3]],{5}]

Matice nahodnych cisel

SeedRandom[1113];
(tabPois3=RandomInteger[PoissonDistribution[3],{3,5}])//MatrixForm

SeedRandom[1113];
(tabPois4=RandomInteger[PoissonDistribution[3],{5,3}])//MatrixForm

MaxMemoryUsed[]

SeedRandom[1113];tabPois5=RandomInteger[PoissonDistribution[3],{100000}];
MaxMemoryUsed[]

Table[ByteCount[(4^j)],{j,20}]

Compatibility/tutorial/Statistics/DescriptiveStatistics

Clear[LocationReport];LocationReport[data_]:={Mean->Mean[data],HarmonicMean->HarmonicMean[data],Median->Median[data]}

LocationReport[tabPois5]
%//N

DispersionReport[tabPois5]//N

{Min[tabPois5],m5=Max[tabPois5]}
Print["P{Po[3]>15} = ",1-CDF[PoissonDistribution[3.],15]]

fr5=Tally[tabPois5] (* zuctovani, tabulka *)

fr5//Sort

Sort[fr5,#1[[2]]<#2[[2]]&]

Map[Reverse,Sort[Map[Reverse,fr5]]]

Needs["BarCharts`"]
BarChart[fr5]

Histogram[tabPois5,PlotRange->All]

Priklad (vlastni distribucni funkce a generovani nahodnych cisel)

dFuniform[x_]:=Which[
x≤0,0,
0<x<1,x,
True,1]
Plot[dFuniform[x],{x,-1,2},AxesOrigin->{0,-1}]

Solve[dFuniform[x]==y,x]

Jiny priklad (exponencialni rozdeleni se stredni hodnotou 1)

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

sol1=Solve[dFexp[x]==y,x]
invF[y_]=First[x/.sol1]
invF[z]

nahexp=invF[RandomReal[UniformDistribution[{0,1}],{100000}]];
Mean[nahexp]

Reduce[dFexp[x]==y,x]

Vyrovnani, regrese

x^(Range[6]-1)

fr5

fit5=Fit[fr5,x^(Range[6]-1),x]

p1=ListPlot[fr5,PlotStyle->PointSize[0.02]];
p2=Plot[fit5,{x,0,25},PlotRange->All];
Show[p1,p2]

FindMinimum[fit5,{x,9}]

FindMinimum[fit5,{x,9.52618271607493`}]

9.526182675682644`

regr5=Regress[fr5,x^(Range[6]-1),x]

RSquared/.regr5

Mean[FRatioDistribution[n,m]]

RSquared is the ratio of the model sum of squares to the total sum of squares.

AdjustedRSquared is equal to 1-(n-1)/(n-p)(1-RSquared)
where n=Length[data]andp is the number of parameters in the model.

Regress[fr5,x^(Range[5]-1),x]

Vybrane statisticke metody Mathematica v. 7

Vsechny nove rysy

(Tento odkaz se spatne hleda.Ostatni odkazy nejsou uplne.)

http://reference.wolfram.com/mathematica/guide/SummaryOfNewFeaturesIn70.html

Statisticke funkce

Toto je vypis z vyse uvedene stranky:

Statistics & Data Analysis
FindCurvePath for reconstruction of sampled curves
Median, Quartiles, QuartileDeviation, ... added support for built - in distributions
LinearModelFit, NonlinearModelFit to construct regression models from data
LogitModelFit, ProbitModelFit, GeneralizedLinearModelFit added
Interpolation now allows immediate calculation of interpolated values
InverseGammaDistribution, InverseChiSquareDistribution added
LevyDistribution added
StudentTDistribution extended to have multiple parameters
ExpectedValue extended to support sampled data
Standardize to adjust constant shift and variance of the data
FittedModel for symbolic representation of fitted models
DesignMatrix to construct the design matrix for a linear model
Entropy function to measure of randomness in data added
ContraharmonicMean added

FindCurvePath

list=Table[{Cos[t],Sin[t]},{t,RandomReal[{0,2π},(* 1500 *)15]}];
ListLinePlot[list,AspectRatio->Automatic]
curve=FindCurvePath[list]
ListLinePlot[list[[curve[[1]]]],AspectRatio->Automatic]
ListCurvePathPlot[list]

Viz tez

FindShortestTour[{{4,3},{1,1},{2,3},{3,-5},{-1,2},{3,4}}]

Median, Quartiles, QuartileDeviation, ...

Nyni funguje i pro zabudovana rozdeleni

Clear[λ];
{Median[ExponentialDistribution[λ]],Quartiles[ExponentialDistribution[λ]]}

Quartiles[distrGauss]

NonlinearModelFit

data={{0,1},{1,0},{3,2},{5,4},{6,4},{7,5}};
nlm=NonlinearModelFit[data,Log[a+b x^2],{a,b},x]
Normal[nlm]

Show[ListPlot[data],Plot[nlm[x],{x,0,7}],Frame->True]

nlm["BestFitParameters"]

ProbitModelFit

V probitove analyze se predpoklada, ze pravdepodobnost P jefunkci nezavisle promennych (Subscript[x, 1]),...(Subscript[x, k]):

P=Φ((Subscript[β, 0])+(Subscript[β, 1])(Subscript[x, 1])+...+(Subscript[β, k])(Subscript[x, k])),

kde (Subscript[β, 0]),...,(Subscript[β, k]) jsou nezname parametry a Φ je distribucni funkce N(0,1). Pozorovany jsou hodnoty Y nabyvajici hodnot 0,1 a jim odpovidajici hodnoty x. V nasledujicim priklade je k=1 a data jsou ve standardnim tvaru {x,y}.

data={{1,0},{1,1},{1,1},{2,0},{2,0},{2,1},{2,1},{2,1},{2,1},{3,1},{3,1},{3,1},{3,1},{3,1}};
probit=ProbitModelFit[data,x,x]
Normal[probit]
Show[ListPlot[data],Plot[probit[x],{x,1,3}],PlotRange->All,AxesOrigin->{.9,-.1}]
pr1=probit["PredictedResponse"]

Podobna je logitova analyza. Misto Φ se vsak pouziva distribucni funkce logistickeho rozdeleni.

LogitModelFit

{PDF[LogisticDistribution[μ,β],x],CDF[LogisticDistribution[μ,β],x]}

logit=LogitModelFit[data,x,x]
Normal[logit]
Show[ListPlot[data],Plot[logit[x],{x,1,3}],PlotRange->All,AxesOrigin->{.9,-.1}]

lo1=logit["PredictedResponse"]

Porovnani predikced pomoci probitu a logitu:

Max[Abs[pr1-lo1]]

LevyDistribution

Rozdeleni s tezkym pravym chvostem (heavy right tail). Pozor, nosic rozdeleni je (μ,∞)!!!

PDF[LevyDistribution[μ,σ],x]
Mean[LevyDistribution[μ,σ]]
Plot[PDF[LevyDistribution[0,1],x],{x,0,5}]
Median[LevyDistribution[μ,σ]]

Integrate[Exp[-a (x^2)],{x,-∞,∞}]
Integrate[Exp[-a (x^2)],{x,-∞,∞},Assumptions->a>0]
Assuming[a>0,Integrate[Exp[-a (x^2)],{x,-∞,∞}]]

Assuming[σ>0,(∫_μ^∞)PDF[LevyDistribution[μ,σ],x]\[DifferentialD]x]

SeedRandom[13];rannorm=RandomReal[NormalDistribution[0,1],{10000}];
SeedRandom[13];ranlev=RandomReal[LevyDistribution[0,1],{10000}];
SeedRandom[13];rancau=RandomReal[CauchyDistribution[0,1],{10000}];

Map[Max,{rannorm,ranlev,rancau}]

Map[Min,{rannorm,ranlev,rancau}]

Map[Mean,{rannorm,ranlev,rancau}]
Map[Median,{rannorm,ranlev,rancau}]

Select[rancau,Abs[#]>5000&]

Mean[CauchyDistribution[0,1]]

StudentTDistribution

Nyni take s parametry polohy (location) μ a meritka (scale) σ.

PDF[StudentTDistribution[μ,σ,ν],x]
{Mean[StudentTDistribution[μ,σ,ν]],Variance[StudentTDistribution[μ,σ,ν]]}

Standardize

Prekvapujici moznosti

Bezne

strs1=Standardize[rs1=RandomReal[{0,10},20]];
{Mean[strs1],Variance[strs1]}

Pouze nulovy prumer:

strs2=Standardize[rs1,Mean,1&];
{Mean[strs2],Variance[strs2]}

strs3=Standardize[rs1,0&,StandardDeviation];
{Mean[strs3],Variance[strs3]}

Standardize[rs1,Median,InterquartileRange]

Standardize[list,(Subscript[f, 1]),(Subscript[f, 2])] is effectively (list-(Subscript[f, 1])[list])/(Subscript[f, 2])[list].

DesignMatrix

Matice planu pro kvadratickou regresi s absolutnim clenem a bez nej.

data=Table[{i,i^(3/2)+RandomReal[]},{i,5}]
DesignMatrix[data,{x,x^2},x]//MatrixForm
DesignMatrix[data,{x,x^2},x,IncludeConstantBasis->False]//MatrixForm

ContraharmonicMean

ContraharmonicMean[{a,b,c,d}]

(* the order p Lehmer contraharmonic mean *)ContraharmonicMean[{a,b,c,d},p]

Statistical Packages

Directory[]
SetDirectory["C:\\Program Files\\Wolfram Research\\Mathematica\\7.0\\AddOns\\Packages"];
Directory[]
FileNames[]
SetDirectory[NotebookDirectory[]]

Needs["MultivariateStatistics`"]
Names["MultivariateStatistics`*"]

MaxMemoryUsed[]
Needs["NumericalDifferentialEquationAnalysis`"]
Needs["PrimalityProving`"]
MaxMemoryUsed[]

ExampleData

ExampleData[]

ExampleData["Text"]

StringTake[ExampleData[{"Text","HomerOdysseyGreek"}],1000]

ExampleData[{"Text","TheRaven"}]

ExampleData[{"Geometry3D","Beethoven"}]

ExampleData[{"Geometry3D","SpaceShuttle"}]

ExampleData[{"Text","UNHumanRightsChinese"}]

Ruzne

Pocasi

WeatherData["Prague","Temperature"]

Text[Grid[{CityData[#,"Name"],Row[{Round@WeatherData[#,"Temperature"],"°C"}]}&/@Take[CityData[],6],Frame->All]]

WeatherData["Prague", #]& /@ {"Temperature", "Pressure", "WindSpeed"}

DateListPlot[WeatherData["Prague","MeanTemperature",{{2008,1,1},{2009,5,2},"Day"}],Joined->True]

DateListPlot[WeatherData["Prague","Temperature",{{2008,1,1},{2009,4,26},"Day"}],Joined->True]

stations=WeatherData[{"Prague",5},"NearestStation"]

WeatherData[{"Prague",5},"StationDistance"]

WeatherData[{"Prague",5},"StationDistance","Units"]

Slider (=posuvnik) bez poklepnuti mysi

Slider[0.5,AutoAction->True]

Speak

Speak["zmrzlina, zmrzlina, zmrzlina, zmrzlina. Once upon a midnight dreary, while I pondered, weak and weary, Over many a quaint and curious volume of forgotten lore-- While I nodded, nearly napping, suddenly there came a tapping, As of some one gently rapping--rapping at my chamber door."]

Konvoluce

Convolve[Exp[-x] UnitStep[x],Sin[x],x,y]

h=Exp[-x] UnitStep[x];
Convolve[h,UnitStep[x],x,y]

pdf=PDF[NormalDistribution[μ,σ],x];
Convolve[pdf,pdf,x,y]//Simplify[#,Assumptions->{σ>0}]&

Poloha pocitace nebo IP

FindGeoLocation[]

FinancialData a CountryData

FinancialData

FinancialData["GE"]

DateListPlot[fd=FinancialData["GE","Jan. 1, 2000"],Joined->True]

fd//First
{Min[fd[[All,2]]],Max[fd[[All,2]]]}

f1=FinancialData["EUR/USD",{2009,2,3}]

Country Data

CountryData["SA","Flag"]

CountryData["CZ","Population"]

ListPlot[Map[Last,CountryData["CZ",{{"GDP"},{1995,2008}}]]]

DateListPlot[CountryData["CZ",{{"GDP"},{1995,2008}}]]

Poznamka k praci s jednotkami, import atypickych formatu a jejich zpracovani

Import kursovniho listku CNB a jeho zpracovani

SetDirectory[NotebookDirectory[]]

Aktualni kursy

http ://www.cnb.cz/cs/financni_trhy/devizovy_trh/kurzy_devizoveho _trhu/denni_kurz.jsp

-> Stazeni v textovem formatu (nazev souboru denni_kurz.txt)

Rozsireni $Path:

(* AppendTo[$Path,"E:\Dokumenty\VPFPM\VPFPM2007 2008"] *)
SetDirectory[NotebookDirectory[]]

Cteni souboru

kursytxt=Import["denni_kurz20090429.txt"];
kursytxt//FullForm

Nejprve zamenit desetinnou carku (vymozenost z roku 1941 za Protektoratu Boehmen und Maehren, kdy CSN prevzala normu DIN, kterou se dosud v ceske matematicke literature psala desetinna tecka)

kursydesetinnatecka=StringReplace[kursytxt,","->"."]

Dale oddelovac "|" nahradit carkami:

(kursyoddelovace=StringReplace[kursydesetinnatecka,"|"->","])//FullForm

Jeednotlive zaznamy (radky) jsou zde oddeleny (jako v UNIX) "\n". Vyse uvedeny jeden kompaktni retezec "kursyoddelovace" tedy rozdelime na retezce podle radku:

(kursyporadcich=StringSplit[kursyoddelovace,"\n"])//FullForm

Prvni dva radky jsou pro dalsi praci nezadouci: odstranime je:

(kursyporadcichbezhlavicky=Drop[kursyporadcich,2])//FullForm

Vytvorime matici pocet men × 5:

(kursymatice5=Map[StringSplit[#,","]&,kursyporadcichbezhlavicky])//MatrixForm

Kazdy zaznam jej jeden zaznam tvoreny 5 polozkami. Vytvorime tedy novy soubor ve tvaru matice pocet men × 2:

Potrebujeme jen posledni dva sloupce: (Pozor, i cislo ma tvar String!!!)

(kursymatice2string=kursymatice5[[All,{4,5}]])//FullForm

(kursymatice2=kursymatice2string//ToExpression)//FullForm

kursymatice2

((kursymatice2[[1]])/(kursymatice2[[2]]))

(kursyAUD=Map[((#)/(kursymatice2[[1]]))&,kursymatice2])//TableForm

HoldForm[(AUD/AUD)]

Select[kursyAUD,#[[2]]>.15&]

Select[kursyAUD,0.63092≤ #[[2]]≤ 0.63093&]

Select[kursyAUD,#[[1]]==RUB/AUD&]

CZK=1;
EUR=25.30 CZK
15 EUR

Nahrazeni desetinne carky desetinnou teckou a ponechani carky v gramaticke vete.

Vychazi z pravidel ceskeho pravopisu: v cisle se za desetinnou carkou nepise mezera, zatimco v textu ano. Provizorne v textu nahradime retezec "carka, mezera" retezcem "XXX", v cislech nahradime carku teckou a pote v textech opet nahradime "XXX" carkou a mezerou.

retez={"245,345","Nevim, jestli to dokazu.", "5,34"}
r1=StringReplace[retez,", "->"XXX"]
(r2=StringReplace[r1,","->"."])
(r3=StringReplace[r2,"XXX"->", "])//FullForm

Ted si jeste postavime ulohu: Text ponechat textem a z cisel ve forme retezec vytvorime skutecne cisla.

ToExpression[r3]//FullForm

Takto jednodusse to vsak nejde. Ty elementy listu, ktere vrati $Failed, ponechame v nezmenene podobe.

Map[If[ToExpression[#]===$Failed,#,ToExpression[#]//Evaluate]&,r3]

Konvoluce Cauchyovych rozdeleni

Clear[μ,σ,x];
PDF[CauchyDistribution[μ,σ],x]

Clear[x,y];
PDF[CauchyDistribution[0,1],x]

Clear[pdf

Konvoluce dvou a ctyr

Clear[pdf2];
pdf2[x_]=Integrate[PDF[CauchyDistribution[0,1],y] PDF[CauchyDistribution[0,1],x-y],{y,-∞,∞},Assumptions->x∈Reals]
pdf4=Integrate[pdf2[y] pdf2[x-y],{y,-∞,∞},Assumptions->x∈Reals]

Rozdeleni aritmetickeho prumeru nahodneho vyberu z Cauchyova rozdeleni

Rozdeleni souctu n Caychy iid ma hustotu

Clear[n,pdfn];
pdfn[x_,n_]=PDF[CauchyDistribution[0,n],x]//Simplify
pdfn[x,4]

Prumer ma tedy hustotu:

Clear[pdfnprumer];
pdfnprumer[x_,n_]=n pdfn[n x,n]
Clear[pdfnprumer];
pdfnprumer[x_,n_]=n pdfn[n x,n]//Simplify

(∫_(-∞)^∞)pdfnprumer[x,n]\[DifferentialD]x

Plot[Table[pdfnprumer[z,n],{n,1,10}]//Evaluate,{z,-5,5},PlotRange->All,AxesOrigin->{0,0}]

Konec 7. 5. 2009

Lekce 12 posledni (14. 5. 2009)

Poznamky k posledni lekci

Testy v definici funkce

Clear[f];
f[x_?Positive]:=Log[x]
{f[-1],f[1]}

Clear[f]; (* f vraci x pokud je x delitelne tremi *)
f[x_?(Mod[#,3]==0&)]:=x
Map[f[#]&,Range[10]]
Map[Mod[#,3]==0&,Range[7]]

Nahrazeni desetinne carky desetinnou teckou a ponechani carky v gramaticke vete.

Vychazi z pravidel ceskeho pravopisu: v cisle se za desetinnou carkou nepise mezera, zatimco v textu ano. Provizorne v textu nahradime retezec "carka, mezera" retezcem "XXX", v cislech nahradime carku teckou a pote v textech opet nahradime "XXX" carkou a mezerou.

retez={"245,345","Nevim, jestli to dokazu.", "5,34"}
r1=StringReplace[retez,", "->"XXX"]
(r2=StringReplace[r1,","->"."])
(r3=StringReplace[r2,"XXX"->", "])//FullForm

Ted si jeste postavime ulohu: Text ponechat textem a z cisel ve forme retezec vytvorime skutecne cisla.

ToExpression[r3]//FullForm

Takto jednodusse to vsak nejde. Ty elementy listu, ktere vrati $Failed, ponechame v nezmenene podobe.

Map[If[ToExpression[#]===$Failed,#,ToExpression[#]//Evaluate]&,r3]

Konvoluce Cauchyovych rozdeleni

Clear[μ,σ,x];
PDF[CauchyDistribution[μ,σ],x]

Clear[x,y];
PDF[CauchyDistribution[0,1],x]

Clear[pdf

Konvoluce dvou a ctyr

Clear[pdf2];
pdf2[x_]=Integrate[PDF[CauchyDistribution[0,1],y] PDF[CauchyDistribution[0,1],x-y],{y,-∞,∞},Assumptions->x∈Reals]
pdf4=Integrate[pdf2[y] pdf2[x-y],{y,-∞,∞},Assumptions->x∈Reals]

Rozdeleni aritmetickeho prumeru nahodneho vyberu z Cauchyova rozdeleni

Rozdeleni souctu n Caychy iid ma hustotu

Clear[n,pdfn];
pdfn[x_,n_]=PDF[CauchyDistribution[0,n],x]//Simplify
pdfn[x,4]

Prumer ma tedy hustotu stejnou jako je zakladni rozdeleni:

Clear[pdfnprumer];
pdfnprumer[x_,n_]=n pdfn[n x,n]
Clear[pdfnprumer];
pdfnprumer[x_,n_]=n pdfn[n x,n]//Simplify

(∫_(-∞)^∞)pdfnprumer[x,n]\[DifferentialD]x

Plot[Table[pdfnprumer[z,n],{n,1,10}]//Evaluate,{z,-5,5},PlotRange->All,AxesOrigin->{0,0}]

Odhady parametru polohy v Cauchyove rozdeleni

Clear[distrCauchy];
distrCauchy[μ_]:=CauchyDistribution[μ,1]

Zde je simulace 1000 nahodnych vyberu z vyse uvedeneho rozdeleni rozsahu 100:

Needs["HypothesisTesting`"];
SeedRandom[131];
prumery=Mean[randomCauchy1=RandomReal[distrCauchy[1],{100,1000}]];
mediany=Median[randomCauchy1];
Print["Prumer ze simulovanych prumeru = ",Mean[prumery],"\nStandardni odchylka = ",  StandardDeviation[prumery],"\n95%-ni interval spolehlivosti",MeanCI[prumery],
"\n{Min, Max} =  ",{Min[prumery],Max[prumery]}]
Histogram[prumery]
Histogram[Select[prumery,Quantile[prumery,0.25]<#<Quantile[prumery,0.75]&]]
Histogram[Select[prumery,Quantile[prumery,0.05]<#<Quantile[prumery,0.95]&]]
Print["\n\nPrumer ze simulovanych medianu = ",Mean[mediany],"\nStandardni odchylka = ",  StandardDeviation[mediany],"\n95%-ni interval spolehlivosti",MeanCI[mediany],
"\n{Min, Max} =  ",{Min[mediany],Max[mediany]}]
Histogram[mediany]

Poznamka k asymptotice vyberovych kvantilu a ukazka v pripade Cauchyova rozdeleni

Teorie pravi, ze vyberovy p kvantil ma asymptoticky normalni rozdeleni. Necht (Subscript[ξ, p]) je teoreticky a     (OverHat[(Subscript[ξ, p])]) je vyberovy p kvantil, 0<p<1. Potom

(Sqrt[n])((OverHat[(Subscript[ξ, p])])-(Subscript[ξ, p]))  ⟶ N(0,((p(1-p))/(([f((Subscript[ξ, p]))])^2)))

v distribuci, kdyz n ⟶ ∞, kde f je hustota prislusneho rozdeleni, ktera musi byt spojita a nenulova v (Subscript[ξ, p]). Pozor, tato limitni veta se nevztahuje na extremni hodnoty, tj. minimum, maximum!!!

Teoreticky rozptyl a smerodatna odchylka v nasem pripade jsou

Clear[p,ξ];
{Quantile[CauchyDistribution[μ,1],p],
rozptylmedianuCauchyovarozdeleni=(1/n)((p(1-p))/((PDF[CauchyDistribution[μ,1],Quantile[CauchyDistribution[μ,1],p]])^2))/.{n->100,μ->1,p->(1/2)},
Print["Asymptoticka smerodatna odchylka = ",(Sqrt[rozptylmedianuCauchyovarozdeleni])//N]}

MeanTest[prumery,1,FullReport->True]
VarianceTest[prumery,1,TwoSided->True]

MeanTest[mediany,1,FullReport->True,TwoSided->True]
VarianceTest[mediany,((π^2)/400),TwoSided->True,FullReport->True]

FAQ

Pocitani exponenciely

Clear[t];
Log[Exp[t]]//FullSimplify
{Assuming[t∈Reals,Log[Exp[t]]//Simplify],
Simplify[Log[Exp[t]],Assumptions->t∈Reals]}

{(E^I),(E^I)}//TrigReduce//FullSimplify

{Log[(E^I)],Log[(((E+I))^I)]}//FullSimplify

Struktura notebooku

Rozmerny vystup

t=RandomReal[{0,1},1000000];
t//Mean
Mean@t

Specialni znaky

Viz Basic typesetting nebo klavesove zkratky Esc mezera esc Esc dava Esc

Cyklus

i=0;
Do[i++;Print[i],{i,10}]

n=1;While[n<4,Print[n];n++]

While

guide/CLanguageInterface

Baliky

?Agglomerate

Needs["Histograms`"]

Shlukova analyza (Cluster analysis)

Zakladni v jadru zabudovana funkce je FindClusters. To je velmi ucinna funkce pro nehierarchicke shlukovani. Ma zabudovane volby, ktere mimo jine umoznuji pracovat s analyzovanymi pozorovanimi databazovymi prostredky. Shlukovane objekty mohou byt libovolneho typu. Miry nepodobnosti (dissimilarities, nekdy tez ne zcela spravne nazyvane vzdalenosti) je mozne volit z pomerne velke mnoziny. Zatim jsem se nesetkal s ucebnici o shlukove analyze, ktera by pojednavala o vetsine mer systemu Mathematica. Miry nepodobnosti v systemu pouzivane jsou pro kvantitativni a kvalitativni promenne a pro retezce. Z novych mer ve verzi 7 upozorneme na

(* Similarity and Distance Measures *)
DamerauLevenshteinDistance 
SmithWatermanSimilarity
NeedlemanWunschSimilarity 
ChessboardDistance

Umistenim kursoru na prislusnou miru a stisknutim klavesy F1 ziskate podrobnou informaci. Doporucuji tez procist si pekne vypracovanou napovedu s dalsimi inspirativnimi priklady (graficke znazorneni shluku v rovine), ktere je mozne s minimalnimi upravami prevzit pro Vase konkretni problemy. Help -> Documentation Centre -> tutorial/PartitioningDataIntoClusters nebo jednoduse FindClusters.

Z nepochopitelnych duvodu u FindClusters chybi jakykoli odkaz na hierarchicke procedury. Ty jsou ve "Standard AddOn Packages". Nize jsou uvedeny i funkce teto procedury. I zde je poucne precist si inspirativni napovedu.

Needs["HierarchicalClustering`"]
Names["HierarchicalClustering`*"]

Balicek (Package) financnich funkci vytvorte za cviceni!

Vzory, predikaty, porovnani vzoru

Vzory (patterns) (viz tez Slavik)

Vzor popisuje tridu vyrazu, o kterych se rika, ze vyhovuji vzoru (they match the pattern). Vzor je vyraz, ktery muze obsahovat vzorove objekty.

Patterns

One of the unique strengths of Mathematica's core language is its powerful and succinct—yet highly readable—symbolic pattern language. Convenient both for immediate use in individual functions, and for systematic large-scale programming, Mathematica's pattern language generalizes concepts like regular expressions to describe general patterns for arbitrary symbolic structures.

Basic Pattern Objects

_  (Blank)— any expression (a "blank" to be filled in)

x_ — any expression, to be referred to as x

__ (BlankSequence)— any sequence of one or more expressions

___ (BlankNullSequence)— any sequence of zero or more expressions

Composite Patterns

p|p|p (Alternatives) — any of several alternatives

p.. (Repeated), p... (RepeatedNull) — a pattern to be repeated

x:p (Pattern) — an arbitrary pattern, to be referred to as x

Except — anything except a specified pattern

paclet:ref/Longest, Shortest — longest, shortest possible matches

Restrictions on Patterns

_h — pattern with a specified head h (e.g. _Integer)

Condition (/;) — condition on a pattern (e.g. x_/;x>7)

PatternTest (?) — pattern with a function test (e.g. _?NumberQ)

Pattern Defaults

_:e (Optional) — pattern that defaults to e if omitted

_. (Optional) — pattern with predefined default

Default — predefined default arguments for a function

Attributes Affecting Patterns

String Patterns »

StringExpression — symbolic string pattern

Pattern Matching Functions »

Predikaty (criterion, crit) (viz tez Slavik)

Predikat je funkce, ktera vraci True nebo False. Jestlize podmineny vzor ma tvar x_/;pred[x], muze byt zapsan v jednodussim tvaru _?pred.

Testing Expressions

Mathematica symbolic expressions can represent an immense range of types of objects. Mathematica provides a rich collection of functions to test expressions. Functions that "ask a question" have names that end in Q. They return True for an explicit true answer, and False otherwise.

Equality & Inequality

Equal (==), Unequal (!=) — numerical equality

SameQ (===), UnsameQ (=!=) — literal symbolic identity

Logical Combinations »

Numerical Properties »

NumberQ, MachineNumberQ — literal number, literal machine-precision number

NumericQ — numeric quantity (consisting of numbers, numeric constants and functions)

List-Oriented Tests

VectorQ, MatrixQ — vector, matrix, perhaps with elements of specified types

ArrayQ — full array of any depth

MemberQ — test whether a pattern occurs in a list

Expression Tests

MatchQ — test whether an expression matches a pattern

String Tests »

Porovnavani vzoru (Pattern matching): funkce

Pattern Matching Functions

Pattern matching makes possible some of the most succinct and elegant programs in the Mathematica language—immediately compressing large numbers of conditional cases into simple, readable and efficient pattern specifications.

Cases — a list of elements matching a pattern

Position — positions of elements matching a pattern

Count — a count of occurrences of a pattern

ReplaceAll (/.) — replace all occurrences of a pattern

DeleteCases — delete all occurrences of a pattern

Set (=), SetDelayed (:=) — define transformations for patterns

MemberQ — test if any member of a list matches a pattern

FreeQ — test if an expression is free of a pattern

MatchQ — test whether an expression matches a pattern

Patterns for Part Specifications

Porovnavani vzoru (Pattern matching): priklady

Vzory, pravidla, definice

(Patterns, rules, definitions)

_  nebo ekvivalentne Blank[] ...vzor, kteremu vyhovuje vse (any expression matches it)

e3={1,x,f[],f[a],f[a,b],g[e4,e5]};

MatchQ[f[]g[a,b]h[1,3,4]2[3,4],_]

{MatchQ[h[a,b],_h],MatchQ[h[a,b],_a]}

Vyhledavani vyhovujicich si vyrazu

Finding Matching Expressions

?Cases

Cases[e3,_] (* vsechny vyrazy *)

Cases[e3,_f] (* vyrazy s hlavickou f *)

Cases[e3,_[_,_]] (* libovolne funkce s prave dvema libovolnymi argumenty *)

Cases[e3,_[__]]  (* libovolne funkce s jednim a vice argumenty *)

Cases[e3,f[___]] (* funkce s hlavickou f s libovolnym, tedy i nulovym poctem argumentu *)

Cases[e3,f[_,_]] (* funkce s hlavickou f s prave dvema argumenty *)

?*Q

Podminene vzory (conditional patterns)

Tvar: predloha/; podminka. Nejdrive se vyhledaji vyrazy vyhovujici predloze a pote ty, ktere splnuji podminku.

Cases[e6={f[1],f[0],f[-1],f[a]},f[x_]/;x≥0]

Jiny mozny zapis:

Cases[e6,f[x_/;x≥0]]

e7=Table[fff[j],{j,9}];

(2≤x≤5||x>6)==Or[2≤x≤5,x>6]

Cases[e7,fff[x_]/;2≤x≤5||x>6]

Predikaty

Predikat je funkce, ktera vraci True nebo False. Nasledujici zapisy jsou ekvivalentni:

{Cases[e6,f[x_/;NonNegative[x]]],Cases[e6,f[x_?NonNegative]],
Cases[e6,f[_?NonNegative]]}

Vlastni map a nest

Zde ukazeme, jak snadne je naprogramovat zabudovane funkce Map a Nest

map (zde pouze pro seznamy, za cviceni sestrojte uplny ekvivalent zabudovane funkce Map)

Casta chyba pri definici vlastnich funkci je, ze se zapomina na pocatecni hodnoty, prazdne mnoziny apod.

Clear[map,f];
map[f_,{}]:={}
map[f_,{x_,y___}]:=Join[{f[x]},map[f,{y}]]
map[f,g[a]]
map[f,{1,2,3}]

nest

Clear[nest];
nest[f_,x_,0]:=x
nest[f_,x_,n_]:=f[nest[f,x,n-1]]
nest[Sin,€,4]

apply I

Clear[apply];
{list={a,b,c},list[[0]]=Plus,list,list[[0]]=Times,list}
sum1=a+b+c;sum1[[0]]=Times;sum1

Clear[apply];
apply[f_,g_]:=g[[0]]=f
apply[f[x],g[x]]
apply[Plus,{a,b,c}]

apply

f[x][[0]]

ap[li_]:=li[[0]]
ap[{a,b,c}//Plus]

Clear[apply];
apply[f_,g_]:=(g[[0]]=f;g)
apply[f[x],g[x]]
apply[Plus,{a,b,c}]

Nelinearni regrese a vynosove krivky

There is a number of parametric models used for the either fitting or interpolation of empirical data coming from markets. We will deal with just two of them which, under the author's opinion are most vital.

Realna data, dvojice ve tvaru {splatnost ve dnech, urokova mira}:

plot3=ListPlot[(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}}),PlotStyle->{Red,PointSize[0.01]}]

Nelson-Siegel

Clear[nelsonsiegelspot]; (* Spot curve or zero-coupon *)
nelsonsiegelspot[b0_,b1_,b2_,c_,t_]:=b0+(((b1+b2) (1-(E^(-(t/c)))))/(t/c))-b2 (E^(-(t/c)))

The parameters have the following meaning: the function starts at level b0+b1 at time 0 and has asymptote at level b0; a local extreme is determined by parameter c.

Clear[t];
Limit[nelsonsiegelspot[b0,b1,b2,c,t],t->0]
Limit[nelsonsiegelspot[b0,b1,b2,c,t],t->∞,Assumptions->{b0>0∧b1>0∧b2>0∧c>0}]

Manipulating spot:

Manipulate[Plot[nelsonsiegelspot[b0,b1,b2,c,t],{t,0,10000},PlotRange->{{0,10000},{3,6}},AxesOrigin->{0,3}],
{{b0,5.8,"b0"},5,6,0.1},
{{b1,-1,"b1"},2,-0.5,0.1},
{{b2,-3,"b2"},-5,-1,0.2},
{{c,400,"c"},100,1000,10}
]

(* Forward curve *)
Clear[nelsonsiegelforward];
nelsonsiegelforward[b0_,b1_,b2_,c_,t_]:=b0+b1 (E^(-t/c))+b2 (t/c) (E^(-t/c))

Manipulating forward:

Manipulate[Plot[nelsonsiegelforward[b0,b1,b2,c,t],{t,0,10000},PlotRange->{{0,10000},{3,6}},AxesOrigin->{0,3},PlotStyle->Orange],
{{b0,5.8,"b0"},5,6,0.1},
{{b1,-1,"b1"},2,-0.5,0.1},
{{b2,-3,"b2"},-5,-1,0.2},
{{c,400,"c"},100,1000,10}
]

Fitting Nelson-Siegel: on the graph, estimated spot is blue and estimated forward is orange.

fitnelsonsiegel3= FindFit[ir3,nelsonsiegelspot[b0,b1,b2,c,t],{b0,b1,b2,c},t]
curvenel3=nelsonsiegelspot[b0,b1,b2,c,t]/.fitnelsonsiegel3;
curvenelfc3[t_]:=nelsonsiegelspot[b0,b1,b2,c,t]/.fitnelsonsiegel3;
curvenel3forward=nelsonsiegelforward[b0,b1,b2,c,t]/.fitnelsonsiegel3;
Show[plot3,Plot[curvenel3forward,{t,0,10000},PlotRange->All,PlotStyle->Orange],Plot[curvenel3,{t,0,10000},PlotRange->All],PlotRange->All]

Swensson

By adding the 5th parameter in Nelson-Siegel model we get Swensson model. Note that possibly two extremes are determined by parameters c1 and c2.

Clear[swenssonforward];
swenssonforward[b0_,b1_,b2_,b3_,c1_,c2_,t_]:=b0+b1 (E^(-t/c1))+b2 (t/c1) (E^(-t/c1))+b3 (t/c2)(E^(-t/c2))
Clear[swenssonspot];
swenssonspot[b0_,b1_,b2_,b3_,c1_,c2_,t_]=(1/t)(∫_0^t)swenssonforward[b0,b1,b2,b3,c1,c2,τ]\[DifferentialD]τ//Simplify

Limit[swenssonspot[b0,b1,b2,b3,c1,c2,t],t->0,Assumptions->{b0>0&&b1>0&&b2>0&&b3>0&&c1>0&&c2>0}]
Limit[swenssonspot[b0,b1,b2,b3,c1,c2,t],t->∞,Assumptions->{b0>0&&b1>0&&b2>0&&b3>0&&c1>0&&c2>0}]

Manipulate[Plot[swenssonspot[b0,b1,b2,b3,c1,c2,t],{t,0,5000},PlotRange->{{0,5000},{3,6}},AxesOrigin->{0,3}],
{{b0,5.8,"b0"},5,6,0.1},
{{b1,-1,"b1"},2,-0.5,0.1},
{{b2,-3,"b2"},-5,-1,0.2},
{{b3,-3,"b3"},-5,-1,0.2},
{{c1,200,"c1"},100,1000,10},
{{c2,5000,"c2"},100,10000,10}
]

Manipulate[Plot[swenssonforward[b0,b1,b2,b3,c1,c2,t],{t,0,5000},PlotRange->{{0,5000},{3,6}},AxesOrigin->{0,3},PlotStyle->Orange],
{{b0,5.8,"b0"},5,6,0.1},
{{b1,-1,"b1"},2,-0.5,0.1},
{{b2,-3,"b2"},-5,-1,0.2},
{{b3,-3,"b3"},-5,-1,0.2},
{{c1,200,"c1"},100,1000,10},
{{c2,5000,"c2"},100,10000,10}
]

fitswensson3= FindFit[ir3,swenssonspot[b0,b1,b2,b3,c1,c2,t],{b0,b1,b2,b3,c1,c2},t]
(curveswe3=swenssonspot[b0,b1,b2,b3,c1,c2,t]/.fitswensson3)//Simplify;
curveswefc3[t_]:=swenssonspot[b0,b1,b2,b3,c1,c2,t]/.fitswensson3
(curveswe3forward=swenssonforward[b0,b1,b2,b3,c1,c2,t]/.fitswensson3)//Simplify;
Show[plot3,Plot[curveswe3,{t,0,10000},PlotRange->All],Plot[curveswe3forward,{t,0,10000},PlotRange->All,PlotStyle->Orange],PlotRange->All]

Interpolation by polynomials

Polynomials of  the degree n-1

Stupid but sometimes observable method: for data size n, use polynomial of the degree n-1. Since the maximum allowed order for interpolation is 21 we use it for the demonstration of the catastrophic behaviour of the interpolated function

intpol3=Interpolation[ir3,InterpolationOrder->Min[Length[ir3],21]];
plotpol3=Plot[intpol3[t],{t,Min[First[ir3][[1]]],Max[Last[ir3][[1]]]},PlotRange->All];
Show[plot3,plotpol3,PlotRange->All]

Manipulate[intpol3=Interpolation[ir3,InterpolationOrder->order];
Show[plotpol3=Plot[intpol3[t],{t,Min[First[ir3][[1]]],Max[Last[ir3][[1]]]},PlotRange->All],plot3,PlotRange->All],
{{order,4,"order"},1,21,1}]

Dalsi priklad na vynosove krivky (PRIBOR)

PRIBOR (Prague Interbank Offered Rate) je prumerna urokova sazba, za kterou si banky navzajem poskytuji uvery na mezibankovnim trhu. PRIBID (Prague Interbank Bid Rate) je prumerna urokova sazba, za kterou si ceske banky mohou ukladat penize u jinych ceskych bank.

http://www.cnb.cz/cs/financni_trhy/penezni_trh/pribor/denni.jsp

http://www.cnb.cz/cs/financni_trhy/penezni_trh/pribor/denni.txt?date=05.05.2009

Soubor ulozit jako Pribor20090505.txt

Zpracovani souboru pomoci Import

ˇ

SetDirectory[NotebookDirectory[]];(pribortxt1=Import["Pribor20090505.txt"])//FullForm
(* Import["Pribor20090505.txt","Text"] vraci to same *)

(* Desetinne carky nahradit teckami, oddelovace | carkami *)pribortxt2=StringReplace[pribortxt1,","->"."];
pribortxt3=StringReplace[pribortxt2,"|"->","];
(* Rozdelit podle \n *)
(pribortxt4=StringSplit[pribortxt3,"\n"])//FullForm
(* Prvni dva radky k nicemu nejsou *)
(pribortxt5=Drop[pribortxt4,2])//FullForm

Map[StringReplace[#," "->""]&,pribortxt5]//InputForm

"d"~~__->"d"

(* Rozdelime radky na retezce *)(pribortxt6=Map[StringSplit[#,","]&,pribortxt5])

Tato konstrukce je obtiznejsi. Misto jiraskovskeho popisu dodaneho poskytovatelem dat potrebuje strucne a jednoznacne vyjadreni. Povsimnete si, ze mezera pred "d", "m", "r" je pohlcena, a vse, co nasleduje za temito znaky je vynechano.

(pribortxt7=Map[StringReplace[#,{" d"~~__->"d"," m"~~__->"m"," r"~~__->"r"}]&,pribortxt6])//FullForm

Nyni vsechno vyjadrime ve dnech.

pribortxt7[[All,1]]//FullForm

(rules7=Thread[pribortxt7[[All,1]]-> Map[ToString,{1,7,14,30,60,90,180,270,360}]])

Zatim ale vsechno je ve forme retezcu.

(pribortxt8=Map[StringReplace[#,rules7]&,pribortxt7])//FullForm

Pripomenme: 1. sloupec termin(splatnost), 2. sloupec PRIBID, #. sloupec PRIBOR. Konecna podoba dat, se kterymi jiz muzeme snadno pracovat.

(pribortxt9=pribortxt8//ToExpression)//TableForm

Zajimat nas bude jenom PRIBOR.

pribortxt10=Map[Drop[#,{2}]&,pribortxt9]

plot10=ListPlot[pribortxt10,PlotStyle->{Red,PointSize[0.01]}]

Zpracovani souboru pomoci ReadList

SetDirectory[NotebookDirectory[]];
stream=OpenRead["Pribor20090505.txt"]
Read[stream,String]
Read[stream,String]
ReadList[stream,{Word,Word,Word},WordSeparators->"|",RecordSeparators->"\n"]
Close[stream]

Dale se postupuje stejne jako v pripade Import.

Odhad parametru vynosove krivky

Clear[nelsonsiegelspot]; (* Spot curve or zero-coupon *)
nelsonsiegelspot[b0_,b1_,b2_,c_,t_]:=b0+(((b1+b2) (1-(E^(-(t/c)))))/(t/c))-b2 (E^(-(t/c)))

fitnelsonsiegel10= FindFit[pribortxt10,nelsonsiegelspot[b0,b1,b2,c,t],{b0,b1,b2,c},t]
curvenel10=nelsonsiegelspot[b0,b1,b2,c,t]/.fitnelsonsiegel10;
curvenelfc10[t_]:=nelsonsiegelspot[b0,b1,b2,c,t]/.fitnelsonsiegel10;
Show[plot10,Plot[curvenel10,{t,0,360},PlotRange->All],PlotRange->All]

Cviceni: Opticky Nelson-Siegelova krivka prilis nesedi. Vyzkousejte Swenssonovu.

Automaticke stazeni aktualni tabulky PRIBOR

"http://www.cnb.cz/cs/financni_trhy/penezni_trh/pribor/denni.txt?date=11.05.2009"

Needs["Calendar`"]
Names["Calendar`*"]

dennipribor1="http://www.cnb.cz/cs/financni_trhy/penezni_trh/pribor/denni.txt?date=11.05.2009"
dennipribor2=StringDrop[dennipribor1,-10]
(actualdate1=DateString[Date[],{"Day",".","Month",".","Year"}])//FullForm
cestaaktualnipriborCNB=StringJoin[dennipribor2,actualdate1]

pribotaktualni=Import[cestaaktualnipriborCNB,"Text"]

Aktualni data jsou vsak k dispozici az po 13. hodine. Proto pro jistotu aktualni data budeme stahovat az po 14. hodine, jinak stahneme data z predchoziho dne. Pozor, nejsou osetreny neobchodni dny!!! Musi byt aktivovan kalendar.

DaysPlus[Date[],-1]

dennipribor1="http://www.cnb.cz/cs/financni_trhy/penezni_trh/pribor/denni.txt?date=11.05.2009"
dennipribor2=StringDrop[dennipribor1,-10]
actualdate2=DateString[If[Date[][[4]]<14,DaysPlus[Date[],-1],Date[]],{"Day",".","Month",".","Year"}]
cestaaktualnipriborCNB2=StringJoin[dennipribor2,actualdate2]

pribotaktualni2=Import[cestaaktualnipriborCNB2,"Text"]

Import[cestaaktualnipriborCNB2,"Table"]

$ImportFormats

Prace s databazi

Approach in elder versions

Simulation of a simple insurance database with fields {names, sex (Male, Female), age, insured amount}

c1=CharacterRange["A","Z"]
ToCharacterCode[c1]//Flatten

ToCharacterCode[{"M","F"}]//Flatten

genernames:=FromCharacterCode[RandomInteger[{65,90},5]]

sex2=Table[FromCharacterCode[RandomInteger[] (77-70)+70],{30}];
names2=Table[genernames,{30}];
ages2=RandomInteger[{20,85},30];
insured2=Table[RandomInteger[{20,100}] 1000,{30}];
database=Transpose[{names2,sex2,ages2,insured2}];

database

Find all men elder than 60, insured over 50000. Both of the following constructions give the same result.

Cases[database,_?(And[#[[2]]=="M",#[[3]]>60,#[[4]]>50000]&)]

sele1=Select[database,(And[#[[2]]=="M",#[[3]]>60,#[[4]]>50000]&)]

Sort these selecting records according to the insured amount:

Sort[sele1,#1[[4]]>#2[[4]]&]

Export the database in the Excel format:

Export["C:\Dokumenty\VPFPM\VPFPM2007\database.xls",database,"XLS"]

Modern approach

SetDirectory[NotebookDirectory[]]

Simulation of a simple insurance database with fields {names, sex (Male, Female), age, insured amount}

n=30; (* The size of the database *)

genername:=StringJoin[RandomChoice[CharacterRange["A","Z"],RandomInteger[{3,10}]]]
names=Table[genername,{n}];

sex=RandomChoice[{"M","F"},n]

ages=RandomInteger[{20,85},n]

insuredamount=Table[RandomInteger[{20,100}] 1000,{n}]

head={"Name","Sex ","Age ","Insured amount"}

database=Transpose[{names,sex,ages,insuredamount}]
(database=PrependTo[database,head])//TableForm

Find all men elder than 60, insured over 50000. Both of the following constructions give the same result.

Join[{head},Cases[database,_?(And[#[[2]]=="M",#[[3]]>60,#[[4]]>50000]&)]]//TableForm

The same

sele=Select[database,(And[#[[2]]=="M",#[[3]]>60,#[[4]]>50000]&)]

Sort these selecting records according to the insured amount:

Sort[sele1,#1[[4]]>#2[[4]]&]//TableForm

Export the database in the Excel format:

Directory[]

Export["database.xls",database]

Databaze: Asociovana pole (Associated Arrays)

mat=(({{"names", "BC", "EFG", "HIJ", "KLM"}, {"a", 1, 2, 3, 4}, {"b", 5, 6, 7, 8}, {"c", 9, 10, 11, 12}, {"d", 13, 14, 15, 16}}));

Clear[data,firstrow];
firstrow=mat[[1]]
fields=Map[data[#]&,firstrow]
matdrop=Drop[mat,1]//Transpose
MapThread[Set,{fields,matdrop}]

data["KLM"]

Select[mat,#[[2]]>7&]
Select[fields,#[[4]]==16&]

??fields

fields//Dimensions

fields

{fields[[1,All]],fields[[All,1]]}

firstrow

FullForm[firstrow]
ToExpression[firstrow]

Remove[symbols];
(symbols=ToExpression[firstrow])//FullForm
nnn=Length[symbols]
MapThread[Set,{symbols,Range[nnn]}]
?symbols
symbols

BC

Diferencni rovnice

Pro reseni diferencnich rovnic (nebo jejich soustav) je zabudovana mohutna funkce RSolve. Syntaxe zadani je prakticky stejna jako u obycejnych rovnic, diferencialnich rovnic, &c. Jak jest obecne z kalkulu znamo, spojite pripady jsou vetsinou daleko snadneji zvladnutelne nez pripady diskretni. Mathematica poskytuje neuveritelne rozsahly aparat pro symbolicke i numericke reseni diferencnich rovnic. Ilustraci vsak zacneme na banalnim prikladu vypoctu Fibonacciovych cisel. Pomineme na chvili skutecnost, ze funkce pro vypocet Fibonacciovych cisel jest funkce zabudovana, a z metodickych duvodu uvedeme nekolik algoritmu jejich vypoctu. Take se pritom poucime, jak efektivne vyuzit prikazy pro rekursivni procedury.

Fibonacci A

Spatna definice

Remove[fibS];
fibS[1]=1;
fibS[2]=1;
fibS[n_]:=fibS[n-1]+fibS[n-2]

Table[fibS[j],{j,6}]
Map[fibS,Range[6]]

Toto nefunguje

fibS[3.5]

fibS[270]

$RecursionLimit=5000;

fibS[270]

fibS[-3]

Trivialni definice:

Remove[fibA];
fibA[1]=1;
fibA[2]=1;
fibA[n_?IntegerQ]:=fibA[n-1]+fibA[n-2]/;Positive[n]

Tady uzivame test, zda n je cele cislo a navic, jestli je kladne.

{fibA[3],fibA[4],fibA[-6],fibA[3.4],fibA[10]}

Casova narocnost je velka (za cviceni zjistete řád a proč)

fibA[29]//Timing

(fibA[24]//Timing)[[1]]

ListPlot[Table[(fibA[j]//Timing)[[1]],{j,24}],Joined->True,PlotRange->All]

Fibonacci B (tzv. Dynamicke programovani)

Vyhodou je, kdyz si jiz vypocitana cisla muzeme zapamatovat. To Mathematica umoznuje nasledujici konstrukci:

$RecursionLimit=256;

Clear[fibB];
fibB[1]=fibB[2]=1;
fibB[n_?IntegerQ]:=fibB[n]=fibB[n-1]+fibB[n-2]

?fibB

Vsimneme si, jak se po vyvolani definice funkci A a B lisi:

{fibA[6],fibB[6]}

?fibA

?fibB

Vsimneme si tez casu potrebneho pro vypocet:

$RecursionLimit=500;

fibB[270]//Timing

?fibB

Ale i tato konstrukce selze pri vyvolani

fibB[300]

To lze odstranit zmenou hodnoty systemove promenne, uvidime, ze to bude nezbytne pri hodnoceni americkych opci:

$RecursionLimit

$RecursionLimit=5000;
fibB[500]//Timing

Na tomto miste je take mozne demonstrovat "listovatelnost" funkci, ktera jiz byla demonstrovana pri vyvolavani BuiltIn funkci. Pro uzivatelem definovane funkce pouzijeme nasledujici postup:

Nejprve zkusme aplikovat funkci fibA na list:

fibA[Range[6]]
Attributes[fibA]

Map[fibA,Range[6]]

Po zmene vlastnosti vsak jiz dosahneme kyzeneho vysledku:

SetAttributes[fibA,Listable];
fibA[Range[6]]
Attributes[fibA]

Fibonacci M

fibM[n_]:=First[MatrixPower[(({{1, 1}, {1, 0}})),n-2].{1,1}]/;n>2

Timing[fibM[1600];]

Fibonacci W

fibW[n_]:=Module[{f1=1,f2=1},Do[{f1,f2}={f2,f1+f2},{n-2}];f2]

Map[fibW,Range[5]]

Fibonacci C (Reseni diferencnich rovnic)

Vcetne pocatecnich podminek se jedna o soustavu 3 rovnic.

Clear[n];
reseniC=RSolve[{fibC[n]==fibC[n-1]+fibC[n-2],fibC[1]==fibC[2]==1},fibC[n],n]

To je obecne reseni. Neni vsak patrne, jak vypada v konkretnim pripade. Uzijeme

reseniC2=reseniC//FunctionExpand

Uvedomme si, ze cos(n π)=(((-1))^n) pro n cele. Kyzeny tvar dosahneme pouzitim predpokladu

Assuming[n∈Integers,Simplify[reseniC2]]

reseniC//FullSimplify

Konkretni hodnoty dostaneme napr. takto:

reseniC/.n->200

Fibonacci D (zabudovana funkce)

Fibonacci[200]//Timing

Fibonacci R (zalozene na lokalnich pravidlech)

Opakovana substituce (ReplaceRepeated)

Clear[fac];
fac[4]/.{fac[0]->1,fac[n_]:>n*fac[n-1]}

Nahrazovani se deje pouze pouze jednim pruchodem listu pravidel. Nasledujici "//." (ReplaceRepeated) aplikuje pravidla tak dlouho, az se vyraz nemeni.

fac[4]//.{fac[0]->1,fac[n_]:>n*fac[n-1]}

⚠ Varovny priklad

Clear[fact];
TimeConstraint[fact[-3]//.{fact[1]->1,fact[n_]:>n*fact[n-1]},3]

Duvod spociva v tom, ze prvni pravidlo se neaplikuje nikdy, a tudiz druhe se aplikuje do nekonecna. Odpomoc spociva v pojistce, dotazu, na typ argumentu. V pripade, ze skutecny parametr nevyhovuje, vypocet se neprovede.

fact[-3]//.{fact[1]->1,fact[n_Integer?Positive]:>n*fact[n-1]}

fact[5]//.{fact[1]->1,fact[n_Integer?Positive]:>n*fact[n-1]}

Fibonacci F (zabudovana funkce, IPM, p. 218)

Clear[fibF];
fibF[n_?(IntegerQ[#]&&Positive[#]&)]:=2n
{fibF[2.5],fibF[3],fibF[-1]}

Diference vyssich radu

Nejdrive cviceni (diference pro vektory dat, srovnej vypocet vynosu z cen)

(********* differences ***********)
ListDifference[data_,n_Integer,s_?(#>5&)]:=Nest[diff[#,s]&,data,n]
diff[data_,s_]:=Drop[data,s]-Drop[data,-s]

Clear[Δ,f,j];
Δ[f_,1,n_]:=f[n]-f[n-1]
Δ[f_,j_,n_]:=Δ[f,j,n]=Δ[f,j-1,n]-Δ[f,j-1,n-1]
Map[Δ[f,#,n]&,Range[4]]//TableForm

Zabudovana funkce

Differences

Srovnej: Derivace vyssich radu

D[f[x],{x,4}]//FullForm

Derivative[-1][f][x]

Clear[x];
g[x_]:=(x^5)
(g^((4)))[x]
g''[x]
(∂_({x,4})) g[x]
Derivative[-1][g][z]

Dalsi priklad (diferencni rovnice, pokracovani)

Ryzi funkce je nepojmenovana funkce urcena k okamzitemu pouziti

(#1+Sin[#2]&)[x,y]

Function[{x,y},x+Sin[y]][x,y]
Function[{x,y},x+Sin[y]][z,v]

Zde dostaneme reseni pomoci ryzi funkce (Function).

solpure=RSolve[{y[x+2]==ay[x+1]+y[x],y[0]==0,y[1]==1},y,x]

y/.solpure[[1]]

(y/.solpure[[1]])[z]

Chceme-li reseni ve forme funkce x, pouzijeme kostrukci (predposledni argument zmenime na y[x])

solr=RSolve[{y[x+2]==ay[x+1]+y[x],y[0]==0,y[1]==1},y[x],x]

Pripomenme, jak vypreparovat reseni:

y[x]/.solr[[1]]

a definovat tuto funkci jako dale pouzitelnou (zkuste, co takova definice udela bez pouziti Evaluate).

yy[x_]:=Evaluate[(y[x]/.solr[[1]])]
yy[z]

(y[x]/.solr[[1]])/.x->z

MaxMemoryUsed[]

Priklady naivniho sifrovani

Julius Caesar

abeceda=CharacterRange["a","z"]

a===First[abeceda]

caesarcoderules=Thread[abeceda->RotateLeft[abeceda]]

caesardecoderules=Thread[abeceda->RotateRight[abeceda]]

Poznamka (uziti pri formovani posloupnosti vynosu z posloupnosti cen

ceny=Array[p,{5}]
rotceny=RotateLeft[ceny]
((rotceny-ceny)/ceny)//Most

p[1]-p[5]//FullForm
p[5]-p[1]//FullForm
y x//FullForm

text="dominika"
roztext=Characters[text]
codetext=roztext/.caesarcoderules
decodetext=codetext/.caesardecoderules
StringJoin[decodetext]

Lepsi zpusob sifrovani

Needs["Combinatorica`"]

SeedRandom[13];
permcoderules=Thread[abeceda->(perm=RandomPermutation[abeceda])]

permc=roztext/.permcoderules

permc/.Thread[perm->abeceda]
%//StringJoin

Data Report

Clear[datareport];
datareport[data_,noclass___]:=Module[{mi=Min[data],ma=Max[data],bin,h,k,le,mea,med,se},
k=If[noclass===Null,10,noclass];
bin=BinCounts[data,{mi,ma,h=((ma-mi)/k)}];Print["No. = ",le=Length[data],"\nMin = ",mi," Max = ",ma,"\nMean= ",mea=Mean[data]," Median = ",med=Median[data],"\nSE = ",se=StandardDeviation[data]," VaR 5% = ",Quantile[data,0.05`]," VaR 1% = ",Quantile[data,0.01`]];Histogram[data,k]]

data=RandomReal[{0,1},10000];
datareport[data]
datareport[data,30]

Kernel estimate

Jadrovy odhad - neparametricky odhad hustoty
x ... hodnota, ve ktere odhad pocitame
k ... jadro
h ... sirka okenka - nutno vyzkouset experimentalne; cim vetsi, tim hladsi krivku dostaneme

Clear[kernel];
kernel[x_,k_,h_,observ_]:=Module[{n=Length[observ],au},
		au=Map[k[(x-#)/h]&,observ];
		Total[au]/(n*h)
		             ]

k[x_]:=Exp[-x^2/2]/Sqrt[2*Pi] (* jadro = hustota N(0,1) *)

observ=Table[Random[]+2*Random[],{50}];

Plot[kernel[x,k,0.2,observ],{x,-1,4}]

Za cviceni sestrojte proceduru pomoci funkce Manipulate, kde budete manipulovat sirkou okenka.

Options and Binomial trees

In this Section we closely follow Chapter 18, Shaw, W.: Modelling Financial Derivatives with Mathematica.

Introduction

Let us consider a time interval (t,t+Δt). Under the Black-Scholes model the conditional distribution of (Subscript[S, t+Δt]) given (Subscript[S, t]) is log-normal LN(ln (Subscript[S, t])+(μ-(1/(σ^2)))Δt,(σ^2)Δt). In a risk neutral world μ=r, where r is the riskless rate and we have

E((Subscript[S, t+Δt])|(Subscript[S, t]))=(Subscript[S, t])(e^(r Δt))=:(Subscript[E, c])     and     var((Subscript[S, t+Δt])|(Subscript[S, t]))=((Subscript[S, t])^2)(e^(2r Δt))((e^((σ^2)Δt))-1)=:(Subscript[V, c]).

Under the discrete binomial random walk we suppose that the stock price (Subscript[S, t]) moves to one of two new values, (Subscript[S, t+Δt])=(Subscript[S, t])u and (Subscript[S, t+Δt])=(Subscript[S, t])d where d<1+r<u with probabilities p and q:=1-p, respectively. In this case we have

E((Subscript[S, t+Δt])|(Subscript[S, t]))=(Subscript[S, t])(p u + q d)=:(Subscript[E, b])    and     var((Subscript[S, t+Δt])|(Subscript[S, t]))=((Subscript[S, t])^2)p ((q(u-d))^2)=(Subscript[V, b]).

While in the continuous model we have two parameters r and (σ^2), ceteris paribus, in the discrete model we have three parameters, i.e., u, d, p. In practice we choose the parameters of the binomial model to match the parameters of the continuous one in terms of the moments:

(Subscript[E, b])=(Subscript[E, c])     and     (Subscript[V, b])=(Subscript[V, c]).

The reason is simple; usually we can estimate r and (σ^2) directly which is hardly possible for the binomial model. Since there is one more degree of freedom, we have to postulate a reasonable relationship between the parameters involved. The corresponding system of equations now reads:

Recall the functions call and put for the European options in a previous lecture:

Clear[Φ]
Φ[x_]:=(1/2) (1+Erf[(x/(Sqrt[2]))])
ϕ[x_]:=Exp[-x^2/2]/Sqrt[2*π]

d1=(Log[s/k]+(r+σ^2/2)*(T-t))/(σ*Sqrt[T-t]);
d2=(Log[s/k]+(r-σ^2/2)*(T-t))/(σ*Sqrt[T-t]);
call[s_,k_,σ_,T_,t_,r_]:=Module[{d1,d2},d1=(Log[s/k]+(r+σ^2/2)*(T-t))/(σ*Sqrt[T-t]);
d2=(Log[s/k]+(r-σ^2/2)*(T-t))/(σ*Sqrt[T-t]);
s*Φ[d1]-k*Exp[-r*(T-t)]*Φ[d2]]
put[s_,k_,σ_,T_,t_,r_]:=Module[{d1,d2},d1=(Log[s/k]+(r+σ^2/2)*(T-t))/(σ*Sqrt[T-t]);
d2=(Log[s/k]+(r-σ^2/2)*(T-t))/(σ*Sqrt[T-t]);
k*Exp[-r*(T-t)]*Φ[-d2]-s*Φ[-d1]]

Clear[p,u,d,r,σ]
eqs={p*u+(1-p)*d==Exp[r*Δt],p*(1-p)*(((u-d))^2)==Exp[2*r*Δt]*(Exp[(σ^2)Δt]-1)};

Jarrow-Rudd scheme

Under this scheme it is supposed that both up and down movements have the same probability.

Clear[p,u,d,r,σ]
p=(1/2);solJR=Solve[eqs,{u,d}]

Only the largest solution makes sense:

{uJR,dJR}={u,d}/.solJR[[2]]

This is an approximation:

approxJR={Series[uJR,{Δt,0,1}]//Normal//PowerExpand,
Series[dJR,{Δt,0,1}]//Normal//PowerExpand}

which is very close to frequently used

Exp[approxJR]/E

Cox-Ross-Rubinstein scheme

Under this scheme it is supposed that u=1/d:

Clear[p,u,d,r,σ]
u=(1/d);solCRR=Solve[eqs,{p,d}]//PowerExpand//FullSimplify

After a hand-made simplification (not too difficult) we get

Clear[CRRupBin];
CRRupBin[n_, r_, q_, σ_, T_] := Module[{f = Exp[-(r-q)*T/n], 
g = Exp[(r-q+σ^2)*T/n]},
1/2*(f+g + Sqrt[(f+g)^2-4])];
Clear[CRRdownBin];
CRRdownBin[n_, r_, q_, σ_, T_] := Module[{f = Exp[-(r-q)*T/n], 
g = Exp[(r-q+σ^2)*T/n]},
1/2*(f+g - Sqrt[(f+g)^2-4])];
Discount[n_, r_, T_] := N[Exp[-r T/n]];

Binomial tree for American options

Let us consider an n steps' evolution, time to expiry T and the time step Δt=T/n. Denote (i,j) the jth node at time iΔt, i=1,...,n, j=0,...,i. Let (Subscript[V, i,j]) be the value of the option (either CALL or PUT) at node (i,j). The price of the stock is (Subscript[S, 0])(u^j)(d^(i-j)) at that node. Since the value for an American PUT is max(K-(Subscript[S, T]),0) at expiry, the value is

(Subscript[V, n,j])=max(K-(Subscript[S, 0])(u^j)(d^(n-j)),0),     j=0,1,...,n.

With probability p the system moves from node (i,j) at time iΔt to (i+1,j+1) at time (i+1)Δt and with probability 1-p to (i+1,j). Without possibility of an early exercise we get

(Subscript[V, i,j])=(e^(-r Δt))(p (Subscript[V, i+1,j+1])+(1-p)(Subscript[V, i+1,j])),   i=0,...,n-1,   j=0,...,i.

If an early exercise is taken into account, the value must be compared with its intrinsic value. Thus we finally get

(Subscript[V, i,j])=(e^(-r Δt))max{K-(Subscript[S, 0])(u^j)(d^(i-j)),(e^(-r Δt))(p (Subscript[V, i+1,j+1])+(1-p)(Subscript[V, i+1,j]))},   i=0,...,n-1,   j=0,...,i.

Starting from level n we go backward through the tree to the desired value (Subscript[V, 0,0]).

Implementation

Here also a dividend yield q is considered.

$RecursionLimit

$RecursionLimit=5000;

afunc[n_, r_, q_, T_] := N[Exp[(r-q) T/n]];
Pfunc[up_,down_,a_,disc_]:=N[disc*(a-down)/(up-down)];
Qfunc[up_,down_,a_,disc_]:=N[disc-Pfunc[up,down,a,disc]];

Clear[CRRAmericanOption]
CRRAmericanOption[S_, n_, σ_, T_, r_, q_, exercise_Function] :=
Module[{u = CRRupBin[n,r,q,σ,T], d = CRRdownBin[n,r,q,σ,T],
a = afunc[n,r,q,T],
disc = Discount[n,r,T],
P,Q,TreeOp,value,delta,gamma,theta},
P = Pfunc[u,d,a,disc];
Q = Qfunc[u,d,a,disc];
TreeOp[node_,level_] := TreeOp[node,level] = 
If[level==n,exercise[S d^node u^(level-node)],
Max[{P,Q}.{TreeOp[node,level+1],TreeOp[node+1,level+1]},
exercise[S d^node u^(level-node)]]];
value=TreeOp[0,0];
delta = (TreeOp[0,1] - TreeOp[1,1])/((u-d)*S);
gamma = ((TreeOp[0,2] - TreeOp[1,2])/(u^2-1) - 
(TreeOp[1,2]-TreeOp[2,2])/(1-d^2))*2/(S^2 (u^2-d^2));
theta = (TreeOp[1,2] - TreeOp[0,0])/(2*T/n);
Clear[TreeOp];
{value, delta, gamma, theta}];

CRRAmericanPut[S_,K_,n_,r_,q_,σ_,T_]:=CRRAmericanOption[S,n,σ,T,r,q,Max[K-#,0]&];
CRRAmericanCall[S_,K_,n_,r_,q_,σ_,T_]:=CRRAmericanOption[S,n,σ,T,r,q,Max[#-K,0]&];

Recall the meaning of the parameters:

CRRAmericanCall[S_,K_,n_,r_,q_,σ_,T_]
call[s_,k_,σ_,T_,t_,r_]

CRRAmericanCall[9, 10,100, 0.1, 0, 0.2, 1]
call[9,10,0.2,1,0,0.1]

Look at the convergence:

Table[{m,CRRAmericanCall[9, 10,m, 0.1, 0, 0.2, 1][[1]]//Evaluate},{m,3,8}]

The following plot (number of steps equal to 80) shows that the convergence is not too convincing and that the recommended number of steps 30 may not be sufficient.

ListPlot[Table[{m,CRRAmericanCall[9, 10,m, 0.1, 0, 0.2, 1][[1]]//Evaluate},{m,3,80}],Joined->True]

The following Example comes from HULL, p.391, where the approximate values u=(e^(σ(Sqrt[Δt]))), d=1/u are used. His result (after 5 steps) is 4.49. A comparison with usual BS formula for an European CALL is also given:

CRRAmericanCall[S_,K_,n_,r_,q_,σ_,T_]
call[s_,k_,σ_,T_,t_,r_]

CRRAmericanPut[50,50,50,0.1,0,0.4,0.4167]
put[50,50,0.4,0.4167,0,0.1]

Demonstrations

http://demonstrations.wolfram.com/search.html?query=Black-Scholes

Sharpe ratio and Markowitz portfolio

Data

(* Needs["Statistics`"] *)
(* Needs["Statistics`MultiDescriptiveStatistics`"] *)
SetDirectory[NotebookDirectory[]]
d1=Import["Returns.dat","Table"]
d1[[{1,2,3,4}]]
d2=Rest[d1];
d3=Rest[d2];
datareturns=Map[Rest,d3];
expected=Mean[datareturns]
(cov=Covariance[datareturns])//MatrixForm
(Correlation[datareturns])//MatrixForm

{{5},{"Month","AMR","BS","UTX","WX","Z"},{1,11.594`,1.515`,-1.579`,-9.36`,-0.408`},{2,22.078`,5.299`,9.947`,7.578`,0.`},{3,-8.511`,-8.`,9.446`,-16.327`,0.694`},{4,-3.488`,0.395`,5.164`,-8.571`,-3.862`},{5,-9.639`,-9.764`,-3.571`,-18.243`,-10.949`},{6,-12.`,6.195`,-0.943`,-11.57`,-8.197`},{7,-10.606`,-0.417`,-7.619`,-8.441`,-6.786`},{8,-1.695`,-6.276`,11.856`,-8.22`,-13.725`},{9,-12.069`,-9.091`,-3.286`,-10.227`,-2.273`},{10,27.451`,15.`,11.186`,-3.797`,2.791`},{11,-21.538`,-4.465`,6.55`,-4.021`,-8.14`},{12,-19.608`,-6.132`,8.75`,12.676`,-5.063`},{13,48.78`,27.156`,4.215`,12.5`,35.2`},{14,19.672`,-2.372`,7.355`,15.271`,12.121`},{15,1.37`,15.19`,13.889`,18.`,-0.901`},{16,4.054`,15.751`,17.073`,-0.847`,18.548`},{17,-15.584`,-9.81`,1.042`,23.884`,-10.958`},{18,3.077`,4.27`,20.833`,4.196`,14.912`},{19,-4.478`,-6.485`,-7.543`,-4.027`,-5.058`},{20,4.688`,11.679`,3.263`,-6.333`,4.918`},{21,-11.94`,-7.285`,3.189`,-18.939`,-3.125`},{22,-1.695`,-3.214`,-7.506`,-2.856`,27.742`},{23,15.517`,0.369`,-3.58`,5.882`,16.026`},{24,5.97`,-1.866`,-7.25`,-0.926`,-2.762`},{25,11.268`,26.616`,15.903`,16.822`,13.295`}}

Cviceni: Prostudujte nasledujici alternativni pristup k extrakci potrebnych dat.

ret1=Import["E:\Dokumenty\VPFPM\VPFPM2007 2008\Returns.dat","Table"];
Short[ret1,5]
ret2=Rest[Rest[ret1]];
Short[ret2,5]
returns=Map[Rest,ret2];
Short[returns,5]

Markowitz Mean-Variance Portfolio

We assume that N assets have expected returns given by the N×1vector r and the covariance matrix of returns V. We will seek the optimal portfolio in the sense that the Sharpe's ratio (also called Sharpe's measure of portfolio) is maximum:

max(((x^T)r)/(Sqrt[(x^T) V x]))subject to conditionsx≥0and(∑_(n=1)^N)(Subscript[x, n])=1,

so that short sales are not allowed.

Example 1 (short sales allowed)

Clear[xx];
xx=Table[Unique[],{Length[expected]}]

crit1=((xx.expected)/(Sqrt[xx.cov.xx]))//Simplify

?Maxim*

?NMax*

sol11=NMaximize[crit1,xx]

xa=xx/.sol11[[2]]
xnorm=xa/Total[xa]
Total[xnorm]

val1=crit1/.{xx->xnorm//Thread}

First[val1]
val1[[1]]

Example 2 (short sales not allowed)

ineq=xx≥0//Thread
ineq=Apply[And,ineq]&&Total[xx]==1

sol2=NMaximize[{crit1,ineq},xx]

xx/.Last[sol2]

Example 2 (Revisited)

ineq2=xx≥0//Thread
ineq2=Apply[And,ineq2]&&Total[xx]≤1
sol2r=NMaximize[{crit1,ineq2},xx]
Total[xx/.sol2r[[2]]]

crit1/.{$16->1.9682812537604517`*^-9,$17->3.7308353824431636`*^-10,$18->(0.4244081277770084`+(1-0.49334391693596863`)),$19->7.572991168693882`*^-11,$20->0.06893578674186551`}

Example 3a (Portfolio with prescribed expected return and minimum risk, short sales not allowed)

crit3=xx.cov.xx//Simplify
cons3=Apply[And,Thread[xx≥0]]&&Total[xx]==1&&xx.expected==4

sol3=NMinimize[{crit3,cons3},xx]

Total[xx3=xx/.Last[sol3]]

xx3
expected

xx3.expected

Example 3b (inequality instead of equality)

crit3b=xx.cov.xx//Simplify;
cons3b=Apply[And,Thread[xx≥0]]&&Total[xx]≤ 1&&xx.expected==4;
sol3b=NMinimize[{crit3b,cons3b},xx]
Total[xx3b=xx/.Last[sol3b]]
xx3b.expected

Example 3c (inequalities instead of equalities)

crit3c=xx.cov.xx//Simplify;
cons3c=Apply[And,Thread[xx≥0]]&&Total[xx]≤ 1&&xx.expected≥4;
sol3c=NMinimize[{crit3c,cons3c},xx]
Total[xx3c=xx/.Last[sol3c]]
xx3c
xx3c.expected

Total[xx3d={1.289839352143948`*^-32,-4.880587180699747`*^-33,0.8613009345035485574`,3.132574895506314`*^-34,0.13860717880556636`}]
xx3d.expected
xx3d.cov.xx3d

Unfeasible requirement
(Neither the asset with maximum expected return reaches 4.3.)

crit34=xx.cov.xx//Simplify;
cons34=Apply[And,Thread[xx≥0]]&&Total[xx]==1&&xx.expected≥ 4.3;
sol3=NMinimize[{crit34,cons34},xx]

Example 4

We have 8 assets traded on Prague Stock Exchange with the expected returns

returns={-9.92488533927622*10^-04,5.72524775783378*10^-03,9.04818959590086*10^-03,-1.44968966764537*10^-03,-1.74045537023702*10^-02,5.70188394543127*10^-03,3.33476815091662*10^-03,2.15075479228363*10^-03};

and the covariance matrix of returns

covariancematrix={{1.35640293495777*10^-03,2.56330097904098*10^-04,3.33037161321956*10^-04,1.80223328810003*10^-04,2.81325992644196*10^-04,5.93263545509886*10^-04,7.13610711759175*10^-04,3.20932535820455*10^-04},{2.56330097904098*10^-04,2.35439634554186*10^-03,4.42407197078591*10^-04,6.57567906193998*10^-06,4.01068448708833*10^-04,5.93785356343409*10^-04,1.30961261626908*10^-04,2.93338399396577*10^-04},{3.33037161321956*10^-04,4.42407197078591*10^-04,1.57895651834543*10^-03,-1.16804549527263*10^-04,1.85336104268171*10^-04,.000176059157865,3.41400743910472*10^-04,1.86624552570071*10^-04},{1.80223328810003*10^-04,6.57567906193998*10^-06,-1.16804549527263*10^-04,6.16786555694974*10^-03,4.84828341822248*10^-04,8.50333422445824*10^-04,7.74803938060802*10^-05,4.39412705238888*10^-04},{2.81325992644196*10^-04,4.01068448708833*10^-04,1.85336104268171*10^-04,4.84828341822248*10^-04,9.45147795376653*10^-03,1.78331918975426*10^-03,5.0987274991634*10^-05,1.49823427374484*10^-03},{5.93263545509886*10^-04,5.93785356343409*10^-04,.000176059157865,8.50333422445824*10^-04,1.78331918975426*10^-03,7.38754498647722*10^-03,1.03084013803777*10^-03,9.34434830622224*10^-04},{7.13610711759175*10^-04,1.30961261626908*10^-04,3.41400743910472*10^-04,7.74803938060802*10^-05,5.0987274991634*10^-05,1.03084013803777*10^-03,2.1036055256453*10^-03,2.40818723671803*10^-04},{3.20932535820455*10^-04,2.93338399396577*10^-04,1.86624552570071*10^-04,4.39412705238888*10^-04,1.49823427374484*10^-03,9.34434830622224*10^-04,2.40818723671803*10^-04,2.77726103166242*10^-03}};

We see that the above matrix is badly conditioned:

Det[covariancematrix]

Next we make use of the procedure NMaximize implemented in the Standard Numerical Package Mathematica, version 4.2. First we create the vector of unknown proportions in the portfolio:

xxPS=Table[(Subscript[x, i]),{i,1,Length[returns]}]

The criterion is given by

criterion=(((xxPS.returns))/(Sqrt[xxPS.covariancematrix.xxPS]))//Simplify;

This loads the package: (Not necessary in version 5.2)

The constraints are

constraints1=Apply[And,Table[(Subscript[x, i])>=0,{i,1,Length[returns]}]];
constraints2=Sum[(Subscript[x, i]),{i,1,Length[returns]}]==1;
constraints=constraints1&&(constraints2)

Next we solve the problem using the option for the method:

soldif=NMaximize[{criterion,constraints},xxPS](*,Method->"DifferentialEvolution"]*)

soldif=NMaximize[{criterion,constraints},xxPS,Method->"DifferentialEvolution"]

Apply[Plus,xxPS/.soldif[[2]]]
Print["Expected return of the optimal portfolio = ",(xxPS/.soldif[[2]]).returns]

Looking on the criterion, we see that the solution is independent of a nonzero multiple. So we can suspend the second constraint and make the

constraints3=Apply[And,Table[(Subscript[x, i])>=0,{i,1,Length[returns]}]]

soldif=NMaximize[{criterion,constraints3},xxPS ]

We see that the optimal value is slightly higher than in the above case. The optimal portfolio together with its expected return (substantially higher than in the previous case) is simply obtained:

sol=xxPS/.soldif[[2]];
sol=sol/Apply[Plus,sol]
sol.returns

0.007859525763768638`

Lagrange multipliers

Example 3 (Markowitz, short sales allowed, μ prescribed expected return)

Clear[n,ones,xx,λ1,λ2,μ];
n=Length[expected];
ones=Table[1,{n}];
xx=Table[Unique[],{n}]
lagrange=(1/2)xx.cov.xx+λ1*(1-ones.xx)+λ2*(μ-expected.xx)//Simplify

Map[D[lagrange,#]&,xx]

expected
μ=3;
(* !!!! *)
eq1=Map[D[lagrange,#]&,xxx={xx,λ1,λ2}//Flatten]==0
Take[Flatten[xxx/.Solve[eq1,xxx]],Length[xx]]

Optimalizace

V optimalizacnich ulohach je mozne volit normy (Subscript[L, 2]), (Subscript[L, 1]), (Subscript[L, ∞]).

tutorial/ConstrainedOptimizationLinearProgramming

LinearProgramming

FindMinimum, FindMaximum, NMinimize, NMaximize, Minimize, and Maximize. LinearProgramming

Metody Monte Carlo

Pokus: dokaze clovek generovat nahodne cislice?

vera=;

(* Inovace *)
vera5=IntegerDigits[vera]

Clear[datareport];
datareport[data_,noclass___]:=Module[{mi=Min[data],ma=Max[data],bin,h,k,le,mea,med,se},
k=If[noclass===Null,10,noclass];
bin=BinCounts[data,{mi,ma,h=((ma-mi)/k)}];Print["No. = ",le=Length[data],"\nMin = ",mi," Max = ",ma,"\nMean= ",mea=Mean[data]," Median = ",med=Median[data],"\nSE = ",se=StandardDeviation[data]," VaR 5% = ",Quantile[data,0.05`]," VaR 1% = ",Quantile[data,0.01`]];Histogram[data,k]]

datareport[vera5//N]
(* Vsimnete si rozdilu mezi //N a bez numericke hodnoty! *)

Puvodne zvoleny postup. Ten sice funguje, ale je zbytecne slozity. Ten vyse uvedeny potrebuje jedinou instrukci.

vera1=ToString[vera];
vera2=Characters[vera1];
vera3=ToExpression[vera2];
(* Pozor, nutne iniciovat proceduru datareport! *)
datareport[vera3]

Test hypotezy, zda data maji diskretni rovnomerne rozdeleni 
DiscreteUniformDistribution[{0, 9}]. Pouziji vlastni balik obsahujici test dobre shody. Za cviceni si vytvorte analogickou proceduru.

AppendTo[$Path,"E:\Dokumenty\VPFPM\VPFPM2007 2008\Hurt"];Needs["StatJHpriv`"]

?chi*

chi2test[vera5,DiscreteUniformDistribution[{0,9}],10,0]

To, ze test neprojde, je patrne jiz z obrazku, pro statistika je tudiz zbytecne provadet test dobre shody. Nicmene zakaznici vyzaduji P-hodnotu i v techto pripadech a nelze jim to zazlivat. Dalsi testy nahodnosti budou nasledovat.

Nahodna cisla

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. Pro ilustraci uzijeme zabudovanou funkci Mod:

Mod[5,3]

Priklad 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]

Priklad 2    (Vetsi problem, rozdeleni dvojic. Stejny generator (Subscript[X, n+1])=5(Subscript[X, n])+1 mod 16.)   Hypoteticky kazda dvojice ma pst  (1/(16^2))=(1/256). Vzhledem k rekurenci vsak neni mozne generovat vice nez 16-1 dvojic; gnerovana cisla lezi na dvou systemech primek:

gen16shift=RotateLeft[gen16]

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

Priklad 3 (Generator kapesniho kalkulatoru TI-59). Nasledujici funkce generuje nahodna cisla smisenou kongruencni metodou delky n+1.

tutorial/RandomNumberGeneration

Clear[gen59];
gen59[n_]:=NestList[Mod[24298 # +9991,199017]&,0,n]

gen59[10]

Po sobe jdouci dvojice

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

gen592[199016]//ListPlot

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

Testovani nahodnosti, specielne generatoru nahodnych cisel

Integrate[Boole[x1<x2<x3],{x1,0,1},{x2,0,1},{x3,0,1}]

Generatory nahodnych cisel se testuji baterii statistickych testu, z nichz nektere uvadime. Podrobnosti HURT Simulacni metody, jeste podrobneji KNUTH The Art of Computer Programming.

1 Test dobre shody viz vyse

2 Test podle znamenek diferenci
Vhodny proti alternative trendu. Pro spojite iid X,Y je P(X≥Y)=(1/2). Pro diskretni priblizne take. Vezmeme-li  X,Y z rozdeleni DiscreteUniformDistribution[{0, m-1}], je

P(X≥Y)=(∑_(i=0)^(m-1))(∑_(j=0, i≥j)^(m-1))(1/(m^2))=(1/(m^2))(∑_(i=0)^(m-1))(∑_(j=0)^i)1=

spocitano jiz v systemu Mathematica

pxy=(1/(m^2))(∑_(i=0)^(m-1))(∑_(j=0)^i)1//Simplify
eC=(n-1)pxy

Necht testovana posloupnost je (Subscript[X, 1]),...,(Subscript[X, n]). Polozme (Subscript[Y, i])=I{(Subscript[X, i])≤(Subscript[X, i+1])}, i=1,...,n-1,

C=(∑_(i=1)^(n-1))(Subscript[Y, i]). Spocteme var C.

varY1=pxy-(pxy^2)//Simplify

Plati cov((Subscript[Y, 1]),(Subscript[Y, 2]))=E (Subscript[Y, 1])(Subscript[Y, 2]) - (((E (Subscript[Y, 1])))^2)  a E (Subscript[Y, 1])(Subscript[Y, 2])= P((Subscript[X, 1])≤(Subscript[X, 2])≤(Subscript[X, 3]))

pxyz=(1/(m^3))(∑_(i=0)^(m-1))(∑_(j=i)^(m-1))(∑_(l=j)^(m-1))1
covY1Y2=pxyz-(pxy^2)//Simplify

varC=(n-1)varY1+2(n-2)covY1Y2//Simplify

Nyni staci aplikovat CLV pro konecne zavisle nahodne veliciny.

rozdily=Rest[vera5]-Most[vera5];

Count[{1,2,3,4,0,0,-6,-7,I},_?NonNegative]

cc=Count[rozdily,_?NonNegative]

kriter=((cc-eC)/(Sqrt[varC]))/.{n->Length[vera5],m->9.}

P(|N(0,1)|>k)=1-Φ(k)+Φ(-k)=2(1-Φ(k)). Ale P-hodnotu muzeme jednodusse spocitat jako

1- Integrate[PDF[NormalDistribution[0,1],x],{x,-k,k}]/.k->kriter

Select[Partition[{1,2,3,4,0,0,-6},3,1],#[[1]]≤#[[2]]≥#[[3]]&]

3 Test extremalnich hodnot (Lokalni extremy)
Pro spojite iid X,Y, Z je P(X≤ Y≥Z)=(1/3). (Staci vzit vsechny permutace trojic (tech je 6) a zjistit, kolik z nich ma lokalni maximum.) Podobne pst lokalniho minima je (1/3). Pst lokalniho extremu v trojici je tedy (2/3). Vezmeme-li  X,Y, Z z rozdeleni DiscreteUniformDistribution[{0, m-1}] a dostaneme

P(X≤Y≥Z)=(1/(m^3))(∑_(i=0)^(m-1))(∑_(j=0,j≥i)^(m-1))(∑_(l=0,l≤j)^(m-1))1=

pstmax=(1/(m^3))(∑_(i=0)^(m-1))(∑_(j=i)^(m-1))(∑_(l=0)^j)1

pstextr=2pstmax (* maxima + minima *)

V nahodne posloupnosti delky n je tedy prumerne

eP=(n-2)(2/(m^3))(∑_(i=0)^(m-1))(∑_(j=i)^(m-1))(∑_(l=0)^j)1

lokalnich extremu. Rozptyl P je priblizne

varP=((16n-29)/90);

Partition[{x1,x2,x3,x4,x5,x6,x7},3,1]

Select[{{1,2,3},{1,2,1},{5,6,7},{3,4,1}},(#[[1]]≤#[[2]]≥#[[3]]&)]
Length[%]

Count[{{1,2,3},{1,2,1},{5,6,7},{3,4,1},{7,3,8}},_?(#[[1]]≤#[[2]]≥#[[3]]&)]

Count[{{1,2,3},{1,2,1},{5,6,7},{3,4,1},{7,3,8}},_?(Or[#[[1]]≤#[[2]]≥#[[3]],#[[1]]≥ #[[2]]≤ #[[3]]]&)]

Clear[extremalnibody];extremalnibody[data_,m_:10]:=Module[{par,pp,n,eP,varP},
	n=Length[data];
	par=Partition[data,3,1];			pp=Count[par,_?(Or[#[[1]]≤#[[2]]≥#[[3]],#[[1]]≥ #[[2]]≤ #[[3]]]&)];
	eP=(((1+m) (1+2 m) (-2+n))/(3 (m^2)));
varP=((16n-29)/90);
Print["n = ",n,"  P = ",pp,"  P-hodnota = ",
2 CDF[NormalDistribution[0.,1],-Abs[((pp-eP)/(Sqrt[varP]))]]]
]

extremalnibody[{1,2,3,4,3,5,7,6,8,9}]

extremalnibody[vera5]

Correlation[vera5,RotateLeft[vera5]]//N

Korelogram

ListPlot[Table[{k,Correlation[vera5,RotateLeft[vera5,k]]},{k,0,15}]//N,Joined->True,PlotRange->All]

Tally[{{4,3},{3,2},{4,3}}]

Metody redukce rozptylu I

f1[u_]:=u
f2[u_]:=(E^u)-1
du=UniformDistribution[{0,1}];

{intf1=(∫_0^1)f1[u]\[DifferentialD]u,intf2=(∫_0^1)f2[u]\[DifferentialD]u}//N

Plot[{f1[u],f2[u]},{u,0,1}]

Obycejna metoda Monte Carlo

Needs["HypothesisTesting`"]

SeedRandom[13]
Timing[randomuniform=RandomReal[{0,1},100000];]
(f1observ=f1[randomuniform])//Length
Mean[f1observ]
Variance[f1observ]
MeanCI[f1observ]

SeedRandom[13]
Timing[randomuniform=RandomReal[{0,1},10000000];]
(f1observ=f1[randomuniform])//Length
Mean[f1observ]
Variance[f1observ]
MeanCI[f1observ]

SeedRandom[13]
Timing[randomuniform=RandomReal[{0,1},10000000];]
(f2observ=f2[randomuniform])//Length
Mean[f2observ]
Variance[f2observ]
MeanCI[f2observ]

Korelacni koeficient mezi X a (X^2)

Clear[correlationcoefficient];
correlationcoefficient[f_,g_,distribution_]:=
	Module[{e1,e2,v1,v2,cov12},
		e1=ExpectedValue[f[u],distribution,u]//N;
		e2=ExpectedValue[g[u],distribution,u];
		v1=ExpectedValue[((f[u])^2),distribution,u]-(e1^2);
		v2=ExpectedValue[((g[u])^2),distribution,u]-(e2^2);
		cov12=ExpectedValue[f[u]g[u],distribution,u]-e1 e2;
		(cov12/(Sqrt[v1 v2]))//Simplify
		]

correlationcoefficient[Function[u,u],Function[u,(u^2)],UniformDistribution[{0,1}]]

e1=ExpectedValue[u,du,u]
e2=ExpectedValue[(u^2),du,u]
cov12=ExpectedValue[(u^3),du,u]-e1 e2
v1=Variance[du]
v2=ExpectedValue[(u^4),du,u]-(e2^2)
corr=(cov12/(Sqrt[v1 v2]))//Simplify
corr//N

e1=ExpectedValue[f1[u],du,u]
e2=ExpectedValue[f2[u],du,u]
cov12=ExpectedValue[f1[u] f2[u],du,u]-e1 e2
v1=Variance[du]
v2=ExpectedValue[((f2[u])^2),du,u]-(e2^2)
corr=(cov12/(Sqrt[v1 v2]))//Simplify
corr//N

Korelovany vyber f(u)=(e^u)-1, g(u)=u

Mean[fg=f2[randomuniform]-f1[randomuniform]]
Length[fg]
Variance[fg]
MeanCI[fg]+intf1

Konec 9. 4. 2008

Historicke lekce (doporuceno k prohlednuti)

Lekce 5 (2008)

Zabudovane algoritmy II

Miliardte prvocislo

Prime[(10^9)]

{PrimeQ[22801763489],PrimeQ[22801763489-1]}

Mersennova prvocisla

jsou prvocisla ve tvaru (Subscript[M, n])=(2^n)-1, n musi byt prvocislo

http ://mathworld.wolfram.com/news/2005-02-26/mersenne/

prvocisla=Prime[Range[30]]

{2,3,5,7,11,13,17,19,23,29,31,37,41,43,47,53,59,61,67,71,73,79,83,89,97,101,103,107,109,113}

(2^prvocisla)-1//Last

pole=PrimeQ[(2^prvocisla)-1]

Count[pole,True]

Select[(2^prvocisla)-1,PrimeQ]

Numerika

Solve[(x^6)+(x^5)+1==0,x]

NSolve[(x^6)+(x^5)+1==0,x]

Plot[(E^x)-Sin[x],{x,-15,1}]

NSolve[eq1= (E^x)-Sin[x]==0,x]

(tabsol=Table[{j,FindRoot[eq1,{x,j}]},{j,-15,1}])//TableForm

Z eq1 vypreparuji funkci

FullForm[eq1]
lhsminusrhs=Apply[List,eq1]

Transpose[tabsol][[2]]

x/.{x->5}

Map[(x^2)/.{x->#}&,Range[5]]
x

Map[lhsminusrhs[[1]]/.#&,Transpose[tabsol][[2]]]//Chop

Integrate[f2=Log[x+Sin[x]],{x,0,2}]

i3=NIntegrate[Log[x+Sin[x]],{x,0,2}]

∫((Sign[x])/(x^(1/3)))\[DifferentialD]x

Assuming[x∈Reals,∫((Sign[x])/(x^(1/3)))\[DifferentialD]x]

i33=(∫_0^∞)Exp[a x]\[DifferentialD]x

Rule-based programming

i33/.If[a_,b_,c___]->b

Integrate[Exp[a x],x]

Integrate[Exp[a x],{x,0,∞},Assumptions->a<0]

(∫_0^∞)Exp[(-1+I) x]\[DifferentialD]x

s1=Series[f2=Log[x+Sin[x]],{x,0,10}]

Series[Gamma[n],{n,∞,6}]

s1*s1

n1=Normal[s1]

i4=Integrate[n1,{x,0,2}]
i3-i4

i5=NIntegrate[n1,{x,0,2}]
i4-i5

i6=Integrate[n1,x]

i6/.x->0

i60=Limit[i6,x->0]
i62=i6/.x->2
i62-i60

Multiple integrals, beware of the integration order!!

i3=Module[{x,y},Integrate[x*y,{x,a,b},{y,0,x}]]

Clear[x,y];
(∫_a^b)(∫_0^x)x y\[DifferentialD]y\[DifferentialD]x

Booleovske funkce

Plot3D[Boole[x^2+ y^2≤1],{x,-10,10},{y,-10,10}]

Integrate[Boole[x^2+ y^2≤1],{x,-10,10},{y,-10,10}]

ineq=y^2-4x^2+4x^4≤0;
Integrate[Boole[ineq],{x,-Infinity,Infinity},{y,-Infinity,Infinity}]

Integration of an ugly function

Plot[Sin[200x],{x,0,10}]

i0 the exact value, i1 its numerical value, i2 the result of a numerical integration

i0=Integrate[Sin[200x],{x,0,10}]
i1=N[i0]
i2=NIntegrate[Sin[200x],{x,0,10}]
i1-i2

Options[NIntegrate]

Reseni soustav nelinearnich rovnic

Solve[eq3={Sin[x]==x-y,Cos[y]==x+y},{x,y}]

NSolve[eq3,{x,y}]

FindRoot[eq3,{x,10},{y,10}]

Za odmenu pekne obrazky

fun1[i_,j_]:=If[GCD[i,j]==1,1,0]
fun2[i_,j_]:=Which[
GCD[i,j]==1,1,
True,0]

array1[n_]:=Table[fun1[i,j],{i,n},{j,n}]
array2[n_]:=Table[fun2[i,j],{i,n},{j,n}]

array1[4]==array2[4]

ListDensityPlot[array1[5],InterpolationOrder->0]

ListDensityPlot[array1[50],InterpolationOrder->0]

ListDensityPlot[Table[Random[],{30},{30}]]

John McLoon principles

jonm@wolfram.co.uk

1 Function names should label tasks not algorithms.
2 Use a minimum of data structures: the underlying principle: everything is an expression
3 Automate as much as possible

A warying example on precision
Varující příklad týkající se přesnosti v systému Mathematica a v systému MS Excel.

4 Try in MS Excel: 0.2->C2. Then copy formula 11*C2-2 15 times. This results in:

0.20.20.20.20.20.20.20.20.20.20.20.2000050.2000510.2005570.2061320.2674570.9420288.36231289.98543987.839710864.24

5 The same in Mathematica is solved as follows:

Clear[f,x];
NestList[f,x,10]

fjo[x_]:=11*x-2
NestList[fjo,0.2,15]

0.2
BaseForm[0.2,2]
0.2`20
BaseForm[0.2`20,2]
BaseForm[0.2`30,2]

But (e.g.):

fjo30=NestList[fjo,0.2`30,30]

Max[Precision,fjo30]

6 Needs["Optimization`UnconstrainedProblems`"]

Needs["Optimization`UnconstrainedProblems`"]

?*FindMinimumPlot

7 Another calculation of a factorial

facn[n_]:=Apply[Times,Range[n]]
facn[5]

8 An open problem:  How is it possible to get the original matrix, i. e. Log, from MatrixExp?

Financni systemy pracujici pod systemem Mathematica

Time Series

Financial Essentials

Velmi jednoduchy, ale poucny.

UnRisk

V soucasne dobe asi nejlepsi. MFF by k nemu mohla mit pristup zdarma, ale na fakulte neni nikdo, kdo by mel kapacitu tento system provozovat.

http://www.wolfram.com/products/applications/unrisk/features.html

Derivatives Expert

Ukazka:

PresentValue[cashflow1,termstructure1,Compounding->{Exponential,2}] // ToN

Zaver

Vetsinou jedna se o otevrene systemy, coz umoznuje uzivateli adaptovat procedury dle sve potreby. Nejsou soucasti dodavky systemu Mathematica a cenove jsou bohuzel srovnatelne s cenou vlastniho systemu.

Find Instance

Je velmi uzitecna funkce. Neumime-li najit vsechna reseni, tato funkce nam poskytne aspon nektere. A nahodnym, nam neznamym mechanismem jich najde i vice, viz:

Clear[x,y];
FindInstance[x^2+y^2≤1,{x,y}]
fi=FindInstance[x^2+y^2≤1,{x,y},1000];
fi//Short

{x,y}/.First[fi]

ca=Map[{x,y}/.#&,fi];
Short[ca]
xx=Map[First,ca]//N;yy=Map[Last,ca]//N;
Length[xx]
{Mean[xx],Mean[yy]}
ϕϕ=Map[If[#[[1]]==0,(π/2)Sign[#[[2]]],ArcTan[((#[[2]])/(#[[1]]))]]&,{xx,yy}//Transpose];
rr=(Sqrt[(xx^2)+(yy^2)]);

ListPlot[ca,AspectRatio->Automatic]

Select[ca,(((#[[1]])^2)+((#[[2]])^2)==1)&]

crit=(((#[[1]])^2)+((#[[2]])^2)==1)&;
Count[ca,_?crit]

{Mean[rr],Mean[ϕϕ]}

Distribution of (r,ϕ)

Clear[x,y,xy,r,ϕ];

(X,Y) possesses the uniform distribution on the unit circle with the PDF f(x,y)=(π^(-1)) if (x^2)+(y^2)≤1 an 0 otherwise. After the transform

xy={r Cos[ϕ],r Sin[ϕ]}

jac=D[xy,{{r,ϕ},1}]
det=Det[jac]//Simplify

the density of (r,ϕ) becomes f(x,y)=r/π, 0≤r≤1, 0≤ϕ≤2π.

Cviceni

Analyzujte a popiste vsechny pouzite kroky v prikladu Find Instance.

Lekce 6 (2008)

FAQ

Poznamky k predchozimu vykladu

1. poznamka

V teto lekci nejsou nazvy In a Out (tj. ocislovani, viz Edit -> Preferences -> Evaluation -> Show In/Out Names)

Poznamka k praci s jednotkami, import atypickych formatu a jejich zpracovani

Aktualni kursy

http ://www.cnb.cz/cs/financni_trhy/devizovy_trh/kurzy_devizoveho _trhu/denni_kurz.jsp

-> Stazeni v textovem formatu (nazev souboru denni_kurz.txt)

Rozsireni $Path:

(* AppendTo[$Path,"E:\Dokumenty\VPFPM\VPFPM2007 2008"] *)
SetDirectory[NotebookDirectory[]]

Cteni souboru

kursytxt=Import["denni_kurz20090429.txt"];
kursytxt//FullForm

Nejprve zamenit desetinnou carku (vymozenost z roku 1941 za Protektoratu Boehmen und Maehren, kdy CSN prevzala normu DIN, kterou se dosud v ceske matematicke literature psala desetinna tecka)

kursydesetinnatecka=StringReplace[kursytxt,","->"."]

Dale oddelovac "|" nahradit carkami:

(kursyoddelovace=StringReplace[kursydesetinnatecka,"|"->","])//FullForm

Jeednotlive zaznamy (radky) jsou zde oddeleny (jako v UNIX) "\n". Vyse uvedeny jeden kompaktni retezec "kursyoddelovace" tedy rozdelime na retezce podle radku:

(kursyporadcich=StringSplit[kursyoddelovace,"\n"])//FullForm

Prvni dva radky jsou pro dalsi praci nezadouci: odstranime je:

(kursyporadcichbezhlavicky=Drop[kursyporadcich,2])//FullForm

Vytvorime matici pocet men × 5:

(kursymatice5=Map[StringSplit[#,","]&,kursyporadcichbezhlavicky])//MatrixForm

Kazdy zaznam jej jeden zaznam tvoreny 5 polozkami. Vytvorime tedy novy soubor ve tvaru matice pocet men × 2:

Potrebujeme jen posledni dva sloupce: (Pozor, i cislo ma tvar String!!!)

(kursymatice2string=kursymatice5[[All,{4,5}]])//FullForm

(kursymatice2=kursymatice2string//ToExpression)//FullForm

kursymatice2

((kursymatice2[[1]])/(kursymatice2[[2]]))

(kursyAUD=Map[((#)/(kursymatice2[[1]]))&,kursymatice2])//TableForm

HoldForm[(AUD/AUD)]

Select[kursyAUD,#[[2]]>.15&]

Select[kursyAUD,0.63092≤ #[[2]]≤ 0.63093&]

Select[kursyAUD,#[[1]]==RUB/AUD&]

CZK=1;

EUR=25.30 CZK
15 EUR

retez={"245,345","Nevim, jestli to dokazu.", "5,34"}

r1=StringReplace[retez,", "->"XXX"]

(r2=StringReplace[r1,","->"."])

r3=StringReplace[r2,"XXX"->", "]//FullForm

Funkce SortBy

SortBy[{{1,2,3},{2,3,1},{3,1,2}, {2,2}}, Last]

SortBy[{{1,2,3},{2,3,1},{3,1,2}, {2,2}}, Total]

SortBy[Range[10], Mod[#,3]&]

SortBy[{{1,2,3},{2,3,1},{3,1,2}, {2,2}}, Total]//Reverse

Diferencni rovnice

Pro reseni diferencnich rovnic (nebo jejich soustav) je zabudovana mohutna funkce RSolve. Syntaxe zadani je prakticky stejna jako u obycejnych rovnic, diferencialnich rovnic, &c. Jak jest obecne z kalkulu znamo, spojite pripady jsou vetsinou daleko snadneji zvladnutelne nez pripady diskretni. Mathematica poskytuje neuveritelne rozsahly aparat pro symbolicke i numericke reseni diferencnich rovnic. Ilustraci vsak zacneme na banalnim prikladu vypoctu Fibonacciovych cisel. Pomineme na chvili skutecnost, ze funkce pro vypocet Fibonacciovych cisel jest funkce zabudovana, a z metodickych duvodu uvedeme nekolik algoritmu jejich vypoctu. Take se pritom poucime, jak efektivne vyuzit prikazy pro rekursivni procedury.

Fibonacci A

Spatna definice

Remove[fibS];
fibS[1]=1;
fibS[2]=1;
fibS[n_]:=fibS[n-1]+fibS[n-2]

Table[fibS[j],{j,6}]
Map[fibS,Range[6]]

fibS[3.5]

fibS[270]

$RecursionLimit=5000;

fibS[270]

fibS[-3]

Trivialni definice:

Remove[fibA];
fibA[1]=1;
fibA[2]=1;
fibA[n_?IntegerQ]:=fibA[n-1]+fibA[n-2]/;Positive[n]

Tady uzivame test, zda n je cele cislo a navic, jestli je kladne.

{fibA[3],fibA[4],fibA[-6],fibA[3.4],fibA[10]}

Casova narocnost je velka (za cviceni zjistete řád a proč)

fibA[29]//Timing

(fibA[24]//Timing)[[1]]

ListPlot[Table[(fibA[j]//Timing)[[1]],{j,24}],Joined->True,PlotRange->All]

Konec 19. 3. 2008

Fibonacci B

Vyhodou je, kdyz si jiz vypocitana cisla muzeme zapamatovat. To Mathematica umoznuje nasledujici konstrukci:

$RecursionLimit=256;

Clear[fibB];
fibB[1]=fibB[2]=1;
fibB[n_?IntegerQ]:=fibB[n]=fibB[n-1]+fibB[n-2]

?fibB

Vsimneme si, jak se po vyvolani definice funkci A a B lisi:

{fibA[6],fibB[6]}

?fibA

?fibB

Vsimneme si tez casu potrebneho pro vypocet:

$RecursionLimit=500;

fibB[270]//Timing

?fibB

Ale i tato konstrukce selze pri vyvolani

fibB[300]

To lze odstranit zmenou hodnoty systemove promenne, uvidime, ze to bude nezbytne pri hodnoceni americkych opci:

$RecursionLimit

$RecursionLimit=5000;
fibB[500]//Timing

Na tomto miste je take mozne demonstrovat "listovatelnost" funkci, ktera jiz byla demonstrovana pri vyvolavani BuiltIn funkci. Pro uzivatelem definovane funkce pouzijeme nasledujici postup:

Nejprve zkusme aplikovat funkci fibA na list:

fibA[Range[6]]
Attributes[fibA]

Map[fibA,Range[6]]

Po zmene vlastnosti vsak jiz dosahneme kyzeneho vysledku:

SetAttributes[fibA,Listable];
fibA[Range[6]]
Attributes[fibA]

Fibonacci M

fibM[n_]:=First[MatrixPower[(({{1, 1}, {1, 0}})),n-2].{1,1}]/;n>2

Timing[fibM[1600];]

Fibonacci C (Reseni diferencnich rovnic)

Vcetne pocatecnich podminek se jedna o soustavu 3 rovnic.

reseniC=RSolve[{fibC[n]==fibC[n-1]+fibC[n-2],fibC[1]==fibC[2]==1},fibC[n],n]

To je obecne reseni. Neni vsak patrne, jak vypada v konkretnim pripade.

Pokusy o zjednoduseni

reseniC//Simplify

reseniC//FullSimplify

Lepe uz to ocividne nejde. Nicmene konkretni hodnoty dostaneme takto:

reseniC/.n->200

Je to sice spravne, ale neni to to, co bychom si prali videt. Pouzijeme tedy zjednoduseni:

reseniC/.n->200//FullSimplify

-(((-5+(Sqrt[5])) (-((((1/2)-((Sqrt[5])/2)))^200)+((((1/2)+((Sqrt[5])/2)))^200)))/(5 (-1+(Sqrt[5]))))//PowerExpand

&c.

Fibonacci D (zabudovana funkce)

Fibonacci[200]//Timing

FAQ

Diference vyssich radu

Nejdrive cviceni (diference pro vektory dat, srovnej vypocet vynosu z cen)

(********* differences ***********)
ListDifference[data_,n_Integer,s_?(#>5&)]:=Nest[diff[#,s]&,data,n]
diff[data_,s_]:=Drop[data,s]-Drop[data,-s]

Clear[Δ,f,j];
Δ[f_,1,n_]:=f[n]-f[n-1]
Δ[f_,j_,n_]:=Δ[f,j,n]=Δ[f,j-1,n]-Δ[f,j-1,n-1]
Map[Δ[f,#,n]&,Range[4]]//TableForm

Derivace vyssich radu

D[f[x],{x,4}]//FullForm

Derivative[-1][f][x]

Clear[x];
g[x_]:=(x^5)
(g^((4)))[x]
g''[x]
(∂_({x,4})) g[x]
Derivative[-1][g][z]

Import grafu ve formatu rtf aj.

Plot[Sin[(x^2)],{x,0,π}];

$ImportFormats

$Path//FullForm

Import["C:\\Documents and Settings\\kpms\\My Documents\\VPFPM Hurt\\doc1.xml"]

Import["C:\\Dokumenty\\VPFPM\\test.xml","XML"]

i1=Import["C:\Documents and Settings\kpms\My Documents\VPFPM Hurt\\GRAPH4.BMP","BMP"]

Show[i1]

Import["C:\\STAT\\EXAMPLES\\GRAPH1.WMF"]

Spolehliva reseni:

(i) Je mozne otevrit notebook: File ⟶ Open ⟶ File of type ⟶ Metafile (*.?MF). Pak je mozne graf oznacit, Edit ⟶ Save Selection As ⟶ napr. RTF, pak otevrit treba ve Wordu.

(ii) Obrazek je mozne primo prepastovat.

Grafy funkci definovanych po castech.

dominika[x_]:=Which[
x≤-1,(x^2),
-1<x≤1,1,
x>1,(x^2)]

Plot[dominika[x],{x,-2,2}]

p1=Plot[(x^2),{x,-3,-1},PlotRange->All];
p2=Plot[1,{x,-1,1}];
p3=Plot[(x^2),{x,1,3}];
Show[p1,p2,p3]

Dalsi priklad (diferencni rovnice, pokracovani)

Ryzi funkce je nepojmenovana funkce urcena k okamzitemu pouziti

(#1+Sin[#2]&)[x,y]

Function[{x,y},x+Sin[y]][x,y]
Function[{x,y},x+Sin[y]][z,v]

Zde dostaneme reseni pomoci ryzi funkce (Function).

solpure=RSolve[{y[x+2]==ay[x+1]+y[x],y[0]==0,y[1]==1},y,x]

y/.solpure[[1]]

(y/.solpure[[1]])[z]

Chceme-li reseni ve forme funkce x, pouzijeme kostrukci (predposledni argument zmenime na y[x])

solr=RSolve[{y[x+2]==ay[x+1]+y[x],y[0]==0,y[1]==1},y[x],x]

Pripomenme, jak vypreparovat reseni:

y[x]/.solr[[1]]

a definovat tuto funkci jako dale pouzitelnou (zkuste, co takova definice udela bez pouziti Evaluate).

yy[x_]:=Evaluate[(y[x]/.solr[[1]])]
yy[z]

(y[x]/.solr[[1]])/.x->z

MaxMemoryUsed[]

Priklady naivniho sifrovani

Julius Caesar

abeceda=CharacterRange["a","z"]

a===First[abeceda]

caesarcoderules=Thread[abeceda->RotateLeft[abeceda]]

caesardecoderules=Thread[abeceda->RotateRight[abeceda]]

Poznamka (uziti pri formovani posloupnosti vynosu z posloupnosti cen

ceny=Array[p,{5}]
rotceny=RotateLeft[ceny]
((rotceny-ceny)/ceny)//Most

p[1]-p[5]//FullForm
p[5]-p[1]//FullForm
y x//FullForm

text="dominika"
roztext=Characters[text]
codetext=roztext/.caesarcoderules
decodetext=codetext/.caesardecoderules
StringJoin[decodetext]

Lepsi zpusob sifrovani

Needs["Combinatorica`"]

SeedRandom[13];
permcoderules=Thread[abeceda->(perm=RandomPermutation[abeceda])]

permc=roztext/.permcoderules

permc/.Thread[perm->abeceda]
%//StringJoin

Standardni funkce

fr=FactorInteger[77542455900]

trfr=Transpose[fr]

mocniny=((First[trfr])^(Last[trfr]))
Apply[Times,mocniny]

Apply[Times,((Map[First,fr])^(Map[Last,fr]))]

Elementarni symbolicke vypocty

poly=Factor[(y^3)+3(y^2)+3y+1]

Expand[poly]

rac=(1/((x^2)-4))-(1/(x-2))
torac=Together[rac]
Numerator[torac]
Denominator[torac]
Denominator[torac]//Expand

Poznamka: Pri zjednodussovani je dobre pouzivat funkce Simplify, PowerExpand, FullSimplify.

poly3=(∏_(j=1)^5)(1+j (x^j))

rac2=(rac/poly3)
Expand[rac2]
ExpandAll[rac2]

poly3=(∏_(j=1)^15)(1+j (x^j));
vytvorfce=(1/poly3);

D[vytvorfce,{x,120}];

apavytvorfce=Apart[vytvorfce];

D[apavytvorfce,{x,200}];

{((x-2)/((x^2)-4))//Simplify,((x-2)/((x^2)-4))//Cancel}

{((x-2)/((x^2)-4))+((x-2)/(((x-2))^2))//Simplify,((x-2)/((x^2)-4))+((x-2)/(((x-2))^2))//Cancel}

Limit[(((((x+h))^3)-(x^3))/h),h->0]

Clear[f];
Limit[((f[x+h]-f[x])/h),h->0]//Simplify

Numericke reseni soustavy diferencialnich rovnic

Clear[x,y,t];
ndsol=NDSolve[{x'[t]==y[t],y'[t]==-4x[t]-0.1y[t],x[0]==1,y[0]==1},{x,y},{t,1,10}]

x/.ndsol

Plot[(x/.ndsol)[[1]][t],{t,1,10}]

ParametricPlot[{(x/.ndsol)[[1]][t],(y/.ndsol)[[1]][t]},{t,1,10}]

Uzivatelem definovane funkce

SeedRandom[13];
addrandomone[x_]=x+Random[]-0.5
addrandomtwo[x_]:=x+Random[]-0.5

Nest[Sin,x,3]
NestList[Sin,x,3]

NestList[addrandomtwo,0,10]

NestList[addrandomone,0,10]

walktwo=NestList[addrandomtwo,0,1000];
walkone=NestList[addrandomone,0,1000];

ListPlot[walktwo,Joined->True]
ListPlot[walkone,Joined->True]

Ryzi funkce (Pure functions)

Clear[g];
g[x_]:=(x^2)
g[2]

Function[w,(w^2)]

%[2]

(((#1)^2)&)[2]

(#1+#2^2&)

%[2,3]

%%[a,b]

SeedRandom[13];
NestList[addrandomtwo,0,3]
SeedRandom[13];
NestList[addrandomtwo[#]&,0,3]
SeedRandom[13];
NestList[#+Random[]-(1/2)&,0,3]

FoldList[Plus,0,{a,b,c}]

FoldList[f,x,{a,b,c}]

RandomReal[{-.5,.5},5]

walk3=FoldList[Plus,0,nahcisla=RandomReal[{-.5,.5},5000]];
ListPlot[walk3,Joined->True]

{Mean[nahcisla],Variance[nahcisla]}

RandomReal[{-.5,.5},{5,2}]

walkplane=FoldList[Plus,{0,0},RandomReal[{-.5,.5},{1000,2}]];
ListPlot[walkplane,Joined->True]
walkplane//Last

ListPlot3D[walkplane]

Show[{ListPlot[walkplane,Joined->True],ListPlot3D[walkplane]}]

Uzivateleme definovane funkce - pokracovani

Je dany objekt permutace?

{4,,3,3,2,1}//Sort

Clear[permQ];
permQ[x_?VectorQ]:=(Sort[x]==Range[Length[x]])
permQ[2]
permQ[{4,1,2,3}]
permQ[{4,,3,3,2,1}]
permQ[{4,3,3,2,1}]
permQ[{1,2,{3},4,5}]

Definice logaritmu

Clear[log];
log[x_*y_]:=log[x]+log[y]
log[((x_)^(n_))]:=n log[x]
log[a (b^3)/(c^5)]

Konec 26. 3. 2008

Lekce 7  (2008)

Komentare k predchozi lekci

Nahrazeni desetinnych carek teckami a jejich ponechani v textu

retez={"245,345","Nevim, jestli to dokazu.", "5,34"}
(* Vychazime z pravidel ceskeho pravopisu; v textu je za carkou mezera, pri psani cisel se mezera nepise *)
r1=StringReplace[retez,", "->"XXX"]
(r2=StringReplace[r1,","->"."])
(r3=StringReplace[r2,"XXX"->", "])//FullForm

Zatim jsou vsechny polozky listu retezce. Z numericky interpretovatelnych retezcu se numericke hodnoty dosahne aplikaci funkce ToExpression. Avsak

((First[r3])/(Last[r3]))

r3//ToExpression

If[ToExpression[r3[[2]]]===$Failed,r3[[2]],ToExpression[r3[[2]]]]//FullForm

(r4=Map[If[ToExpression[#]===$Failed,#,ToExpression[#]]&,r3])//FullForm
((First[r4])/(Last[r4]))//FullForm

%/2

%%/2

%%/2

%%%$/2

Cviceni. Prozkoumejte, jak je mozne vyuzit funkci 
?StringQ
"\!\(\*RowBox[{\"StringQ\", \"[\", StyleBox[\"expr\", \"TI\"], \"]\"}]\) gives True if \!\(\*StyleBox[\"expr\", \"TI\"]\) is a string, and False otherwise." »

PRIBOR

http ://www.cnb.cz/cs/financni_trhy/penezni_trh/pribor/denni.jsp

Chybove hlaseni pro definici funkci

Clear[prumer1];
prumer1[data_]:=Module[{delka},
Which[(delka=Length[data])==0,Print["Pocet pozorovani je nula."],
!VectorQ[data],Print["Data nejsou vektor."],
True,((Total[data])/delka)]
		]

prumer1[Range[5]]
prumer1[{1,2,3,{4},5}]
prumer1[{}]
prumer1[π]

tutorial/Messages

f::overflow="Factorial argument `1` too large."

Here is the function f.

f[x_]:=If[x>10,(Message[f::overflow,x];Infinity),x!]

When the argument of f is greater than 10, the message is generated.

{f[20],f[10],f[-3.5],f[-3]}

This switches off the message.

Off[f::overflow]

Now the message is no longer generated.

f[20]

FAQ

Poznamka k funkci Which

Clear[jh];
jh[x_]:=Which[x== 2,7777,
			  x=="graf",Plot[Sin[(y^2)],{y,0,π}],
			  True,"nebavim se s Tebou"]
jh[2]
jh["graf"]
jh["gras"]

Priklad nekonvergujici rady

s1=(∑_(n=1)^∞)(1/n)

castecnesoucty=(∑_(n=1)^k)(1/n)
Limit[castecnesoucty,k->∞]
Limit[HarmonicNumber[n],n->∞]

Stromy zivota viz Wolfram Demonstrations -> Trees

Stromova struktura viz Cominatorica; take

{a,b,{c,d,{e,f,{g},{{s,v},{r,t,g,x},c}}}}//TreeForm

{ja,{prvnitatinek,{druhamaminka,{baba,dedek}},{druhytatinek,prvnimaminka}}}//TreeForm

Vypreparovani polozek se sudymi indexy

Clear[vyber];
vyber[data_]:=Module[{delka,delka2,indexy},
	delka=Length[data];
	delka2=Floor[(delka/2)];
	indexy=2 Range[delka2];
	data[[indexy]]
]
vyber[((Range[11])^2)]

Lepsi reseni

Data Report

<<Calendar`(* ano *)
<<Histograms`

Date[]
{Take[Date[],3],Drop[Date[],-3],Date[][[Range[3]]]}
startday=Date[];
difdat=DaysBetween[{2019,2,7},startday]

{FullForm[_],FullForm[__],FullForm[___]}

hh[xx_,yy___]:=If[yy=== Null,xx,yy]
hh[2]
hh[3,4]

a==b

a===b

Clear[datareport];
datareport[data_,noclass___]:=Module[{mi=Min[data],ma=Max[data],bin,h,k,le,mea,med,se},
k=If[noclass===Null,10,noclass];
bin=BinCounts[data,{mi,ma,h=((ma-mi)/k)}];Print["No. = ",le=Length[data],"\nMin = ",mi," Max = ",ma,"\nMean= ",mea=Mean[data]," Median = ",med=Median[data],"\nSE = ",se=StandardDeviation[data]," VaR 5% = ",Quantile[data,0.05`]," VaR 1% = ",Quantile[data,0.01`]];Histogram[data,k]]

data=RandomReal[{0,1},10000];
datareport[data]
datareport[data,30]

Mnozinove funkce

?f

Join[{a,b,c,d},{f},{d,e,b,a}]

Union[{a,b,c,d,f,d,e,b,a}]

Pozor, Union setridi!!!

s1={3,2,4,5,2,1,3,3};s2={6,3,4,5,1,2,3,2,1};
s1⋃ s2
s1⋂s2
Complement[Range[10],s1]

Map a Apply - opakovani

Clear[f,g,a,b,c];
li={a,b,c};
soucin=a b c;
{Head[li],Head[soucin]}
{FullForm[li],FullForm[soucin]}

Map preskoci hlavicku druheho argumentu a applikuje prvni argument na to, co je za hlavickou druheho argumentu

Head[g[a,b,c]]
Map[f,g[a,b,c]]
Map[Plus,g[-a,b,c]]

Attributes[Plus]

Plus[Plus[a,b],c]

Plus[{a,b,c}]

Total[{a,b,c}]

Apply je funkce, ktera vymeni hlavicku druheho argumentu za prvni argument

Apply[f,g[a,b,c]]

Priklad vypoctu aritmetickeho prumeru s pojistkou na to, aby argument funkce byl vektor delky alespon jedna.

prumer[x_]:=Apply[Plus,x]/Length[x]/;VectorQ[x]&&Length[x]>0
{prumer[{a,b,c}],prumer[{{a},b,c}],prumer[a,b,c],prumer[2]}

li

soucin

Head[3 I/4]

Apply[Times,li]
Apply[List,soucin]
Apply[Plus,soucin]

Lekce 8  (2008)

FAQ

Podtrzitko zastupuje cokoli

VectorQ[{1,2}]

FullForm[_]

Vypreparovani polozek se sudymi indexy

Clear[vybersudych];
vybersudych[data_]:=Module[{delka,delka2,indexy},
	delka=Length[data];
	delka2=Floor[(delka/2)];
	indexy=2 Range[delka2];
	data[[indexy]]
]
vybersudych[vs=((Range[11])^2)]

Lepsi reseni (Analyzujte za cviceni!)

Map[Last,Partition[vs,2]]

Last[Transpose[Partition[vs,2]]]

Cviceni

Select[Range[20],Mod[#,3]==1&]

Select[Range[11],EvenQ]

With[{data=((Range[11])^2)},data[[Select[Range[Length[data]],EvenQ]]]]

(*******************************************)

BarChart

data=RandomInteger[{0,9},{10000,2}];
Take[data,10]

The following bar chart shows the marginal distribution of the first dimension.

Needs["BarCharts`"]

frqs=Sort[Tally[data[[All,1]]]];
BarChart[frqs[[All,2]]/10^4,BarLabels->frqs[[All,1]]]

Metody Monte Carlo

Pokus: dokaze clovek generovat nahodne cislice?

vera=;

(* Inovace *)
vera5=IntegerDigits[vera]

Clear[datareport];
datareport[data_,noclass___]:=Module[{mi=Min[data],ma=Max[data],bin,h,k,le,mea,med,se},
k=If[noclass===Null,10,noclass];
bin=BinCounts[data,{mi,ma,h=((ma-mi)/k)}];Print["No. = ",le=Length[data],"\nMin = ",mi," Max = ",ma,"\nMean= ",mea=Mean[data]," Median = ",med=Median[data],"\nSE = ",se=StandardDeviation[data]," VaR 5% = ",Quantile[data,0.05`]," VaR 1% = ",Quantile[data,0.01`]];Histogram[data,k]]

datareport[vera5//N]
(* Vsimnete si rozdilu mezi //N a bez numericke hodnoty! *)

Puvodne zvoleny postup. Ten sice funguje, ale je zbytecne slozity. Ten vyse uvedeny potrebuje jedinou instrukci.

vera1=ToString[vera];
vera2=Characters[vera1];
vera3=ToExpression[vera2];
(* Pozor, nutne iniciovat proceduru datareport! *)
datareport[vera3]

Test hypotezy, zda data maji diskretni rovnomerne rozdeleni 
DiscreteUniformDistribution[{0, 9}]. Pouziji vlastni balik obsahujici test dobre shody. Za cviceni si vytvorte analogickou proceduru.

AppendTo[$Path,"E:\Dokumenty\VPFPM\VPFPM2007 2008\Hurt"];Needs["StatJHpriv`"]

?chi*

chi2test[vera5,DiscreteUniformDistribution[{0,9}],10,0]

To, ze test neprojde, je patrne jiz z obrazku, pro statistika je tudiz zbytecne provadet test dobre shody. Nicmene zakaznici vyzaduji P-hodnotu i v techto pripadech a nelze jim to zazlivat. Dalsi testy nahodnosti budou nasledovat.

Nahodna cisla

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. Pro ilustraci uzijeme zabudovanou funkci Mod:

Mod[5,3]

Priklad 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]

Priklad 2    (Vetsi problem, rozdeleni dvojic. Stejny generator (Subscript[X, n+1])=5(Subscript[X, n])+1 mod 16.)   Hypoteticky kazda dvojice ma pst  (1/(16^2))=(1/256). Vzhledem k rekurenci vsak neni mozne generovat vice nez 16-1 dvojic; gnerovana cisla lezi na dvou systemech primek:

gen16shift=RotateLeft[gen16]

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

Priklad 3 (Generator kapesniho kalkulatoru TI-59). Nasledujici funkce generuje nahodna cisla smisenou kongruencni metodou delky n+1.

tutorial/RandomNumberGeneration

Clear[gen59];
gen59[n_]:=NestList[Mod[24298 # +9991,199017]&,0,n]

gen59[10]

Po sobe jdouci dvojice

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

gen592[199016]//ListPlot

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

Testovani nahodnosti, specielne generatoru nahodnych cisel

Integrate[Boole[x1<x2<x3],{x1,0,1},{x2,0,1},{x3,0,1}]

Generatory nahodnych cisel se testuji baterii statistickych testu, z nichz nektere uvadime. Podrobnosti HURT Simulacni metody, jeste podrobneji KNUTH The Art of Computer Programming.

1 Test dobre shody viz vyse

2 Test podle znamenek diferenci
Vhodny proti alternative trendu. Pro spojite iid X,Y je P(X≥Y)=(1/2). Pro diskretni priblizne take. Vezmeme-li  X,Y z rozdeleni DiscreteUniformDistribution[{0, m-1}], je

P(X≥Y)=(∑_(i=0)^(m-1))(∑_(j=0, i≥j)^(m-1))(1/(m^2))=(1/(m^2))(∑_(i=0)^(m-1))(∑_(j=0)^i)1=

spocitano jiz v systemu Mathematica

pxy=(1/(m^2))(∑_(i=0)^(m-1))(∑_(j=0)^i)1//Simplify
eC=(n-1)pxy

Necht testovana posloupnost je (Subscript[X, 1]),...,(Subscript[X, n]). Polozme (Subscript[Y, i])=I{(Subscript[X, i])≤(Subscript[X, i+1])}, i=1,...,n-1,

C=(∑_(i=1)^(n-1))(Subscript[Y, i]). Spocteme var C.

varY1=pxy-(pxy^2)//Simplify

Plati cov((Subscript[Y, 1]),(Subscript[Y, 2]))=E (Subscript[Y, 1])(Subscript[Y, 2]) - (((E (Subscript[Y, 1])))^2)  a E (Subscript[Y, 1])(Subscript[Y, 2])= P((Subscript[X, 1])≤(Subscript[X, 2])≤(Subscript[X, 3]))

pxyz=(1/(m^3))(∑_(i=0)^(m-1))(∑_(j=i)^(m-1))(∑_(l=j)^(m-1))1
covY1Y2=pxyz-(pxy^2)//Simplify

varC=(n-1)varY1+2(n-2)covY1Y2//Simplify

Nyni staci aplikovat CLV pro konecne zavisle nahodne veliciny.

rozdily=Rest[vera5]-Most[vera5];

Count[{1,2,3,4,0,0,-6,-7,I},_?NonNegative]

cc=Count[rozdily,_?NonNegative]

kriter=((cc-eC)/(Sqrt[varC]))/.{n->Length[vera5],m->9.}

P(|N(0,1)|>k)=1-Φ(k)+Φ(-k)=2(1-Φ(k)). Ale P-hodnotu muzeme jednodusse spocitat jako

1- Integrate[PDF[NormalDistribution[0,1],x],{x,-k,k}]/.k->kriter

Select[Partition[{1,2,3,4,0,0,-6},3,1],#[[1]]≤#[[2]]≥#[[3]]&]

3 Test extremalnich hodnot (Lokalni extremy)
Pro spojite iid X,Y, Z je P(X≤ Y≥Z)=(1/3). (Staci vzit vsechny permutace trojic (tech je 6) a zjistit, kolik z nich ma lokalni maximum.) Podobne pst lokalniho minima je (1/3). Pst lokalniho extremu v trojici je tedy (2/3). Vezmeme-li  X,Y, Z z rozdeleni DiscreteUniformDistribution[{0, m-1}] a dostaneme

P(X≤Y≥Z)=(1/(m^3))(∑_(i=0)^(m-1))(∑_(j=0,j≥i)^(m-1))(∑_(l=0,l≤j)^(m-1))1=

pstmax=(1/(m^3))(∑_(i=0)^(m-1))(∑_(j=i)^(m-1))(∑_(l=0)^j)1

pstextr=2pstmax (* maxima + minima *)

V nahodne posloupnosti delky n je tedy prumerne

eP=(n-2)(2/(m^3))(∑_(i=0)^(m-1))(∑_(j=i)^(m-1))(∑_(l=0)^j)1

lokalnich extremu. Rozptyl P je priblizne

varP=((16n-29)/90);

Partition[{x1,x2,x3,x4,x5,x6,x7},3,1]

Select[{{1,2,3},{1,2,1},{5,6,7},{3,4,1}},(#[[1]]≤#[[2]]≥#[[3]]&)]
Length[%]

Count[{{1,2,3},{1,2,1},{5,6,7},{3,4,1},{7,3,8}},_?(#[[1]]≤#[[2]]≥#[[3]]&)]

Count[{{1,2,3},{1,2,1},{5,6,7},{3,4,1},{7,3,8}},_?(Or[#[[1]]≤#[[2]]≥#[[3]],#[[1]]≥ #[[2]]≤ #[[3]]]&)]

Clear[extremalnibody];extremalnibody[data_,m_:10]:=Module[{par,pp,n,eP,varP},
	n=Length[data];
	par=Partition[data,3,1];			pp=Count[par,_?(Or[#[[1]]≤#[[2]]≥#[[3]],#[[1]]≥ #[[2]]≤ #[[3]]]&)];
	eP=(((1+m) (1+2 m) (-2+n))/(3 (m^2)));
varP=((16n-29)/90);
Print["n = ",n,"  P = ",pp,"  P-hodnota = ",
2 CDF[NormalDistribution[0.,1],-Abs[((pp-eP)/(Sqrt[varP]))]]]
]

extremalnibody[{1,2,3,4,3,5,7,6,8,9}]

extremalnibody[vera5]

Correlation[vera5,RotateLeft[vera5]]//N

Korelogram

ListPlot[Table[{k,Correlation[vera5,RotateLeft[vera5,k]]},{k,0,15}]//N,Joined->True,PlotRange->All]

Tally[{{4,3},{3,2},{4,3}}]

Metody redukce rozptylu I

f1[u_]:=u
f2[u_]:=(E^u)-1
du=UniformDistribution[{0,1}];

{intf1=(∫_0^1)f1[u]\[DifferentialD]u,intf2=(∫_0^1)f2[u]\[DifferentialD]u}//N

Plot[{f1[u],f2[u]},{u,0,1}]

Obycejna metoda Monte Carlo

Needs["HypothesisTesting`"]

SeedRandom[13]
Timing[randomuniform=RandomReal[{0,1},100000];]
(f1observ=f1[randomuniform])//Length
Mean[f1observ]
Variance[f1observ]
MeanCI[f1observ]

SeedRandom[13]
Timing[randomuniform=RandomReal[{0,1},10000000];]
(f1observ=f1[randomuniform])//Length
Mean[f1observ]
Variance[f1observ]
MeanCI[f1observ]

SeedRandom[13]
Timing[randomuniform=RandomReal[{0,1},10000000];]
(f2observ=f2[randomuniform])//Length
Mean[f2observ]
Variance[f2observ]
MeanCI[f2observ]

Korelacni koeficient mezi X a (X^2)

Clear[correlationcoefficient];
correlationcoefficient[f_,g_,distribution_]:=
	Module[{e1,e2,v1,v2,cov12},
		e1=ExpectedValue[f[u],distribution,u]//N;
		e2=ExpectedValue[g[u],distribution,u];
		v1=ExpectedValue[((f[u])^2),distribution,u]-(e1^2);
		v2=ExpectedValue[((g[u])^2),distribution,u]-(e2^2);
		cov12=ExpectedValue[f[u]g[u],distribution,u]-e1 e2;
		(cov12/(Sqrt[v1 v2]))//Simplify
		]

correlationcoefficient[Function[u,u],Function[u,(u^2)],UniformDistribution[{0,1}]]

e1=ExpectedValue[u,du,u]
e2=ExpectedValue[(u^2),du,u]
cov12=ExpectedValue[(u^3),du,u]-e1 e2
v1=Variance[du]
v2=ExpectedValue[(u^4),du,u]-(e2^2)
corr=(cov12/(Sqrt[v1 v2]))//Simplify
corr//N

e1=ExpectedValue[f1[u],du,u]
e2=ExpectedValue[f2[u],du,u]
cov12=ExpectedValue[f1[u] f2[u],du,u]-e1 e2
v1=Variance[du]
v2=ExpectedValue[((f2[u])^2),du,u]-(e2^2)
corr=(cov12/(Sqrt[v1 v2]))//Simplify
corr//N

Korelovany vyber f(u)=(e^u)-1, g(u)=u

Mean[fg=f2[randomuniform]-f1[randomuniform]]
Length[fg]
Variance[fg]
MeanCI[fg]+intf1

Konec 9. 4. 2008

Lekce 9 (2008)

FAQ

Maximize: Null: Notebook nebyl zaslan.

Video

Uzitecny odkaz - demonstrace uziti systemu Mathematica v ruznych oborech:  http://www.wolfram.com/products/mathematica/portraits/

Cyklus

Ukol: Realizovat radu simulacnich uloh pro ruzne pocatecni hodnoty generatoru nahodnych cisel. Ilustrace na obycejne metode MC. Prvni reseni neni sice uplne v duchu sytemu Mathematica, je spise v proceduralnim stylu (pouziva iterator), ale je prehlednejsi. Druhe reseni pouziva zabudovanou funkci Do. Prirozene jsou tez zabudovane funkce For, While. Vetsinou je zbytecne je pouzivat: vetsinu uloh vyresime vyuzitim daleko efektivnejsich funkci typu NestWhile, TakeWhile, LengthWhile FixedPoint, FixedPointList, FoldList apod. Viz guide/FunctionalIteration.

Needs["HypothesisTesting`"]
f1[u_]:=u
f2[u_]:=(E^u)-1
du=UniformDistribution[{0,1}];
Table[SeedRandom[13+i];
{Timing[randomuniform=RandomReal[{0,1},100000];]//First,
(f2observ=f2[randomuniform])//Length,
Mean[f2observ],
Variance[f2observ],
MeanCI[f2observ]},{i,5}]//TableForm

Nerad, ale prozrazuji druhe, proceduralni reseni pomoci skutecneho cyklu. Nejprve priklady jednoduche:

Do[Print[n^2],{n,4}]

Do[Print[n],{n,-3,5,2}]

Do[Print[x],{2}]

Do[Print[i],{i,10,8,-1}]

Do[Print[n],{n,x,x+3y,y}]

Do[Print[{i,j}],{i,4},{j,i-1}]

t=67;Do[Print[t];t=Floor[t/2],{3}]

Nyni nase reseni

Do[SeedRandom[13+i];
{Timing[randomuniform=RandomReal[{0,1},100000];]//First,
(f2observ=f2[randomuniform])//Length,
Mean[f2observ],
Variance[f2observ],
MeanCI[f2observ]}//Print,{i,5}]

Procedura pro korelacni koeficient mezi dvema funkcemi

Clear[correlationcoefficient];
correlationcoefficient[f_,g_,distribution_]:=
	Module[{e1,e2,v1,v2,cov12},
		e1=ExpectedValue[f[u],distribution,u]//N;
		e2=ExpectedValue[g[u],distribution,u];
		v1=ExpectedValue[((f[u])^2),distribution,u]-(e1^2);
		v2=ExpectedValue[((g[u])^2),distribution,u]-(e2^2);
		cov12=ExpectedValue[f[u]g[u],distribution,u]-e1 e2;
		(cov12/(Sqrt[v1 v2]))//Simplify
		]

correlationcoefficient[Function[u,u],Function[u,(u^2)],UniformDistribution[{0,1}]]

Generovani nahodnych cisel z jinych rozdeleni

Univerzalni je metoda inverzni transformace. Pro nektera rozdeleni je casto efektivnejsi pouzit specielni metody.

N(0,1) 
a) Box-Mullerova transformace (metoda presna, nekdy pomala, pro kongruencni smisenou metodou dava kuriosni vysledky)

bm1[x1_,x2_]:=(Sqrt[-2Log[x1]])Sin[2π x2]
bm2[x1_,x2_]:=(Sqrt[-2Log[x1]])Cos[2π x2]

Clear[bm];
bm[{x1_,x2_}]:={(Sqrt[-2Log[x1]])Sin[2π x2],(Sqrt[-2Log[x1]])Cos[2π x2]}

SeedRandom[13];
nv={RandomReal[{0,1},500000],RandomReal[{0,1},500000]}//Transpose;
nv[[{1,2}]]
norvar=Map[bm,nv]//Flatten;
datareport[norvar,100]

chi2test[norvar,CDF[NormalDistribution[0,1],#]&,10,0]

b) CLV: (Subscript[Y, 1]),...(Subscript[Y, 12]) iid R(0,1)  ⟹	(∑_(i=1)^12)(Subscript[Y, i]) -6 ∼ N(0,1) 
Presna hustota se ziska jako konvoluce viz Convolution of Uniforms III

Convolution of uniforms III

See also explicit formula in Renyi, A.: Theory of Probability

Clear[funi,con,x,z];
funi[x_]:=Which[0≤x≤1,1,
	True,0]
con[x_,1]=(Piecewise[{{0, x<0}, {1, 0≤x≤1}, {0, x>1}}])
con[x_,n_]:=con[x,n]=(∫_0^n)(con[z,n-1]/.z->x-t)funi[t]\[DifferentialD]t//Simplify

con[x,1]

c3=con[x,3]

c6=(∫_0^6)(c3/.x->x-t)(c3/.x->t)\[DifferentialD]t//Simplify

c12=(∫_0^12)(c6/.x->x-t)(c6/.x->t)\[DifferentialD]t//Simplify

c12input=(Piecewise[{{-((((-12+x))^11)/39916800), 11<x<12}, {((x^11)/39916800), 0<x≤1}, {(1/39916800)(-10+x) (74284944372-60673478460 x+22310165230 (x^2)-4863566720 (x^3)+696087844 (x^4)-68343584 (x^5)+4661680 (x^6)-218120 (x^7)+6700 (x^8)-122 (x^9)+(x^10)), x==11}, {-(1/19958400)(-9+x) (107974382202-102689692542 x+43546218192 (x^2)-10860918147 (x^3)+1766716557 (x^4)-196061157 (x^5)+15045759 (x^6)-788949 (x^7)+27069 (x^8)-549 (x^9)+5 (x^10)), x==10}, {(1/6652800)(6105755284-13251797328 x+13029593500 (x^2)-7654933440 (x^3)+2983069320 (x^4)-808780896 (x^5)+155499960 (x^6)-21178080 (x^7)+2000460 (x^8)-124740 (x^9)+4620 (x^10)-77 (x^11)), 5<x<6}, {(1/13305600)(12945728-47385096 x+78768800 (x^2)-78459480 (x^3)+51997440 (x^4)-24049872 (x^5)+7909440 (x^6)-1845360 (x^7)+298320 (x^8)-31680 (x^9)+1980 (x^10)-55 (x^11)), 3<x<4}, {(1/6652800)(211242138736-320976156732 x+220638695200 (x^2)-90533256660 (x^3)+24630254880 (x^4)-4664006424 (x^5)+627211200 (x^6)-59902920 (x^7)+3982440 (x^8)-175560 (x^9)+4620 (x^10)-55 (x^11)), 7<x<8}, {(1/39916800)(3919268323356-4517350959132 x+2357542443300 (x^2)-735516395460 (x^3)+152455299480 (x^4)-22049990424 (x^5)+2271293640 (x^6)-166664520 (x^7)+8539740 (x^8)-291060 (x^9)+5940 (x^10)-55 (x^11)), 9<x<10}, {(1/3326400)(3104094230-6387049746 x+6039682880 (x^2)-3434650065 (x^3)+1295420280 (x^4)-338227582 (x^5)+62109432 (x^6)-7998540 (x^7)+706860 (x^8)-40810 (x^9)+1386 (x^10)-21 (x^11)), x==6}, {(1/39916800)(12-132 x+660 (x^2)-1980 (x^3)+3960 (x^4)-5544 (x^5)+5544 (x^6)-3960 (x^7)+1980 (x^8)-660 (x^9)+132 (x^10)-11 (x^11)), 1<x<2}, {-(((((-1+x))^6) (-258-96 x-116 (x^2)+39 (x^3)-36 (x^4)+5 (x^5)))/19958400), x==2}, {(1/1663200)(1175112-5262477 x+9463905 (x^2)-9655470 (x^3)+6396720 (x^4)-2928464 (x^5)+948486 (x^6)-216480 (x^7)+33825 (x^8)-3410 (x^9)+198 (x^10)-5 (x^11)), x==4}, {(1/1663200)(41310443132-60669031731 x+39954760615 (x^2)-15573918360 (x^3)+3993010560 (x^4)-707401772 (x^5)+88420794 (x^6)-7803840 (x^7)+477015 (x^8)-19250 (x^9)+462 (x^10)-5 (x^11)), x==8}, {(1/39916800)(-2680731676644+2742649040868 x-1272457556700 (x^2)+353483604540 (x^3)-65344700520 (x^4)+8442009576 (x^5)-777906360 (x^6)+51135480 (x^7)-2350260 (x^8)+71940 (x^9)-1320 (x^10)+11 (x^11)), 10<x<11}, {(1/13305600)(-656721822976+816825452256 x-455962034440 (x^2)+150963005160 (x^3)-32976989760 (x^4)+4995810512 (x^5)-536115888 (x^6)+40790640 (x^7)-2158200 (x^8)+75680 (x^9)-1584 (x^10)+15 (x^11)), x==9}, {(1/13305600)(214184-474012 x+194480 (x^2)+462330 (x^3)-790680 (x^4)+630476 (x^5)-321552 (x^6)+113520 (x^7)-27720 (x^8)+4400 (x^9)-396 (x^10)+15 (x^11)), x==3}, {(1/6652800)(-44214794856+76868873004 x-60370376220 (x^2)+28168558110 (x^3)-8653010520 (x^4)+1834245028 (x^5)-273528024 (x^6)+28687560 (x^7)-2074380 (x^8)+98560 (x^9)-2772 (x^10)+35 (x^11)), x==7}, {(1/6652800)(-391056256+986665680 x-1157665740 (x^2)+827160180 (x^3)-396285120 (x^4)+132480040 (x^5)-31281096 (x^6)+5181000 (x^7)-586740 (x^8)+43120 (x^9)-1848 (x^10)+35 (x^11)), x==5}, {(1/13305600)(-994854930208+1306889097096 x-776748491200 (x^2)+275693192280 (x^3)-64929416640 (x^4)+10655224272 (x^5)-1243482240 (x^6)+103221360 (x^7)-5974320 (x^8)+229680 (x^9)-5280 (x^10)+55 (x^11)), 8<x<9}, {(1/6652800)(-339557216+927890172 x-1150094000 (x^2)+852879060 (x^3)-420055680 (x^4)+144094104 (x^5)-35075040 (x^6)+6046920 (x^7)-722040 (x^8)+56760 (x^9)-2640 (x^10)+55 (x^11)), 4<x<5}, {(1/39916800)(-135156+743292 x-1857900 (x^2)+2785860 (x^3)-2783880 (x^4)+1945944 (x^5)-970200 (x^6)+344520 (x^7)-85140 (x^8)+13860 (x^9)-1320 (x^10)+55 (x^11)), 2<x<3}, {(1/6652800)(-49764991340+89177904816 x-72328491620 (x^2)+35024109120 (x^3)-11243278200 (x^4)+2510700192 (x^5)-397746888 (x^6)+44684640 (x^7)-3488100 (x^8)+180180 (x^9)-5544 (x^10)+77 (x^11)), 6<x<7}}])

Plot[c12,{x,0,12},PlotPoints->1000]

Plot[{c12input,PDF[NormalDistribution[6,1],x]}//Evaluate,{x,0,12},PlotPoints->1000]

Plot[{c12input,PDF[NormalDistribution[6,1],x]}//Evaluate,{x,6-.1,6+.1},PlotPoints->1000]

AppendTo[$Path,"E:\Dokumenty\VPFPM\VPFPM2007 2008\Hurt"];Needs["StatJHpriv`"]

SeedRandom[13];
nv1=RandomReal[{0,1},12*1000000];
nv2=Partition[nv1,12];
normalvariates=Map[Total,nv2]-6;

datareport[normalvariates,100]

chi2test[normalvariates,CDF[NormalDistribution[0,1],#]&,10,0]

chi2test[normalvariates,CDF[NormalDistribution[0,1],#]&,100,0]

Metody redukce rozptylu II

?f2

Korelovany vyber f(u)=(e^u)-1, g(u)=u

Mean[fg=f2[randomuniform]-f1[randomuniform]]
Length[fg]
Variance[fg]
MeanCI[fg]+intf1

Variance[fg]
MeanCI[fg]+intf1

2Sqrt[0.04318586885151399`/10000]

0.004156241997358382`+ 0.7159502441735873`

ikor=Mean[f2[randomuniform]-f1[randomuniform]]+(1/2)

Stratifikovany (oblastni) vyber

Viz skripta Hurt, J.: Simulacni metody

Stejne velke oblasti

Nestejne velke oblasti

Madras: 3.1 The Basic Problem (p. 29)

Analyticky nespocita, Mathematica to, co neumi, vrati v puvodnim tvaru:

Clear[fmadras];
fmadras[x_]:=((E^(Sqrt[x]))/(3+Cos[x]))
(∫_0^1)fmadras[x]\[DifferentialD]x

Priblizne pomoci NIntegrate (N vzdy znamena pocitej numericky):

intmadras=NIntegrate[fmadras[x],{x,0,1}]

Najdeme ridici funkci pomoci asymptotickeho rozvoje fce fmadras, zde do pateho radu. Z takove funkce je pak snadne spocitat integral presne.

pfmadras=Plot[fmadras[x],{x,0,1},AxesOrigin->{0,0}]
sermadras=Series[fmadras[x],{x,0,5}]
psermadras=Plot[(nsermadras=Normal[sermadras])//Evaluate,{x,0,1},AxesOrigin->{0,0}];
Show[pfmadras,psermadras]

fmadras-nsermadras

Plot[fmadras[x]-nsermadras,{x,0,1},PlotRange->{Automatic,{-.0001,.0001}}]

intnsermadrasnum=NIntegrate[nsermadras,{x,0,1}]
intnsermadrasexact=Integrate[nsermadras,{x,0,1}]
intnsermadrasnum-intnsermadrasexact

intmadras-intnsermadras

Je jasne, ze pri volbe teto funkce bude mit metoda korelovaneho vyberu bude mit uspech

Twodimensional case

Clear[fmadrasxy]; (* Integral ani po 284 sec nespocita *)
fmadrasxy[x_,y_]:=((E^(Sqrt[x]))/(3+Cos[x]))((E^(Sqrt[y]))/(3+Cos[y]))((E^(Sqrt[x+y]))/(3+Cos[x+y]))
(∫_0^1)(∫_0^1)fmadrasxy[x,y]\[DifferentialD]x\[DifferentialD]y

NIntegrate[fmadrasxy[x,y],{x,0,1},{y,0,1}]

nsermadrasxy=Series[fmadrasxy[x,y],{x,0,3},{y,0,3}]//Simplify//Normal

Plot3D[fmadrasxy[x,y]-nsermadrasxy,{x,0.01,1},{y,0.01,1}]

h1=(∫_0^1)(∫_0^1)nsermadrasxy\[DifferentialD]y\[DifferentialD]x

h1//N

(∫_0^1)(∫_0^1)nsermadrasxy\[DifferentialD]x\[DifferentialD]y

Cviceni. Pouzijte rozvoj jako ridici funkci.

Antiteticke promenne

Zakladni myslenka: Pouzit dva nebo vice zaporne korelovanych odhadu. Nejjednodussi pripad:

Priklad 1

correlationcoefficient[Function[u,u],Function[u,1-u],UniformDistribution[{0,1}]]

correlationcoefficient[Function[u,u],Function[u,1-u],ExponentialDistribution[1]]

Ukazka, ze neplati pro "rychle se menici" funkce:

correlationcoefficient[Function[u,Sin[10u]],Function[u,Sin[10(1-u)]],UniformDistribution[{0,1}]]

correlationcoefficient[Function[u,Sin[2u]],Function[u,Sin[2-2u]],UniformDistribution[{0,1}]]

Needs["HypothesisTesting`"]
f1[u_]:=u
f2[u_]:=(E^u)-1
du=UniformDistribution[{0,1}];
SeedRandom[13]
Timing[randomuniform=RandomReal[{0,1},10000000];]
(f2observAnt1=f2[randomuniform])//Length
(f2observAnt2=f2[1-randomuniform])//Length
f2observAnt=(1/2)(f2observAnt1+f2observAnt2);
Mean[f2observAnt]
Variance[f2observAnt]
MeanCI[f2observAnt]

Regrese

Idea: Optimalne zkombinovat dva nebo vice odhadu pomoci v Gauss-Markovove modelu

Clear[yy];
yy[f_,u_]:={0.5(f[u]+f[1-u]),
0.25(f[0.5u]+f[0.5(1-u)]+f[0.5(1+u)]+f[1-0.5u])}

simulations[f_,n_]:=Map[yy[f,#]&,RandomReal[{0,1},n]]

(∫_0^1)Sin[5(x^2)]\[DifferentialD]x//N
s3=simulations[Sin[5((#)^2)]&,30000];
me=Mean[s3]
cov=Covariance[s3]
Det[cov]

βest=(1/({1,1}.(inv=Inverse[cov]).{1,1})){1,1}.inv.me

Total[inv]

βest=(1/(Total[Flatten[inv=Inverse[cov]]]))Total[inv].me

{1,1}.(({{a, b}, {c, d}}))

(({{a, b}, {c, d}}))//Total

regr[f_,n_,seed___]:=Module[
{sim,me,cov,inv,den1,β},
If[seed>0,SeedRandom[seed]];
sim=Map[{0.5(f[#]+f[1-#]),
0.25(f[0.5#]+f[0.5(1-#)]+f[0.5(1+#)]+f[1-0.5#])}&,
RandomReal[{0,1},n]];
me=Mean[sim];
cov=Covariance[sim];
inv=Inverse[cov];
den1=(1/(Total[Flatten[inv]]));
β=den1 Total[inv].me;
Print["β est = ",β,"  with  variance = ",(den1/n),"  and standard error = ",(Sqrt[(den1/n)])]
]

regr[Sin[5((#)^2)]&,100000]

regr[(E^(#))-1&,10000,13]

(∫_0^1)((E^x)-1)\[DifferentialD]x

%//N

Packages

Using packages

Context

User's package

Otevrit File -> New -> Package (.m)

BeginPackage["Hurt`Stat2`","Histograms`"]
(* usage messages for the exported functions and the context itself *)
datareport::usage = "datareport[data,noclass] returns the basic statistical characteristics and the histogram of the data. Data should be in the form of vector. Variable noclass should be a positive integer and is the number of the categories for the histogram. If missing, the default value is set to 10."
datareport::baddata = "Argument `1` should be a vector!"
datareport::badclasses = "Argument `2` should be a positive integer!"
Begin["`Private`"]    (* begin the private context (implementation part) *)
(* Clear[datareport]; *)
datareport[data_,noclass_:10]:=Module[{mi=Min[data],ma=Max[data],bin,h,k,le,mea,med,se},
k=noclass;
If[Not[VectorQ[data]],Message[datareport::baddata]];
If[Not[IntegerQ[k]&&k>0],Message[datareport::badclasses,k]];
bin=BinCounts[data,{mi,ma,h=((ma-mi)/k)}];Print["No. = ",le=Length[data],"\nMin = ",mi," Max = ",ma,"\nMean= ",mea=Mean[data]," Median = ",med=Median[data],"\nSE = ",se=StandardDeviation[data]," VaR 5% = ",Quantile[data,0.05`]," VaR 1% = ",Quantile[data,0.01`]];Histogram[data,HistogramCategories->k,(* BarStyle->{GrayLevel[0.95`]}, *)ApproximateIntervals->True]]
End[ ]         (* end the private context *)
EndPackage[ ]  (* end the package context *)

Ulozit  do stejne pojmenovaneho adresare s priponou *.m (package)

$Path

Quit[]

AppendTo[$Path,"E:\Dokumenty\VPFPM\VPFPM2007 2008\Hurt"];

Directory[]

Needs["Stat2`"]

?datareport

datareport[((Range[500])^2),15]

Vybrane financni aplikace I

Kalendar Actual/Actual

Needs["Calendar`"]

Datum = {YYYY,MM,DD}

04==4

dat1=Date[][[{1,2,3}]]
dattab=Map[DaysPlus[dat1,30*#]&,Range[12]]

DaysBetween[dat1,dattab//Last]
DaysBetween[dattab//Last,dat1]

Soucasna hodnota I

Diskontni faktory

(((1+r))^(-Range[5]+1))

Clear[pv1];
pv1[cf_?VectorQ,r_]:=(((1+r))^(-Range[Length[cf]]+1)).cf
cf1=.;
cf11={cf0,cf1,cf2,cf3}
pv1[cf11,r]

Promenne urokove miry (okamzite, spotove, z vynosove krivky)

rr={r1,r2,r3}
(((1+rr))^(-Range[Length[rr]]))

pv2[cf_,r_]:=Which[
Length[cf]=!=Length[r]+1,"Nestejne delky",
True,Print["Soucasna hodnota = "];First[cf]+Rest[cf].(((1+r))^(-Range[Length[r]]))]
pv2[cf11,rr]
pv2[cf11,{r1,r2}]

Konec 16. 4. 2008

Lekce 10 (2008)

Sharpen

FAQ

Mathematica pro zacatecniky (Antonin Slavik KDM)

verze 6.0.2.1

Rekordy

SeedRandom[13];
dataprorekordy=RandomReal[{0,1},100];
p1=ListPlot[dataprorekordy,Joined->True]
rekordyaux=FoldList[Max,First[dataprorekordy],dataprorekordy];
ListPlot[rekordyaux,Joined->True]
rekordy=Union[rekordyaux]

Position[{a,b,c,d,e},d]

posice=Map[Position[dataprorekordy,#]&,rekordy]//Flatten

p2=ListPlot[{posice,rekordy}//Transpose,PlotStyle->PointSize[0.03]];

Show[p1,p2]

dataprorekordy[[40]]
dataprorekordy[[posice]]

ListPlot[dataprorekordy,Ticks->{Automatic,Table[z,{z,0,1,0.1}]},Joined->True]

Demonstrations

http : // demonstrations.wolfram.com/search.html?query = Black + Scholes & limit = 20

Optimalizace

V optimalizacnich ulohach je mozne volit normy (Subscript[L, 2]), (Subscript[L, 1]), (Subscript[L, ∞]).

tutorial/ConstrainedOptimizationLinearProgramming #293106634

Vybrane financni aplikace II

Soucasna hodnota II

Needs["Calendar`"]

Platby v neekvidistatnich casech

Date[]
dat1=Date[][[{1,2,3}]]
dattab=Map[DaysPlus[dat1,30*#]&,Range[12]]

platby=Map[DaysBetween[dat1,#]&,dattab]-190

pv3[cf_,r_,startday_,datumy_]:=Module[{da=Map[DaysBetween[startday,#]&,datumy],diskontnifaktory},
diskontnifaktory=(1+r)^(-da/365);
cf.diskontnifaktory
	     ]

pv3[platby,r,dat1,dattab]

pv3[platby,0.03,dat1,dattab]

Plot[pv3[platby,r,dat1,dattab],{r,0,0.5}]

IRR

Pro vyse uvedeny priklad

NSolve[pv3[platby,r,dat1,dattab]==0,r]

FindRoot[pv3[platby,r,dat1,dattab]==0,{r,0.1},EvaluationMonitor->True]

Pro ekvidistantni casove intervaly

Clear[pv1];
pv1[cf_?VectorQ,r_]:=(((1+r))^(-Range[Length[cf]]+1)).cf
cf1=.;
cf11={cf0,cf1,cf2,cf3}
pv1[cf11,r]

Priklad

x=.;
NestList[Function[x,(x/(1+r))],c,4]

FoldList[Function[{x,r},(x/(1+r))],c,{(Subscript[r, 1]),(Subscript[r, 2]),(Subscript[r, 3])}]

(* Kuponova obligace, F=1000, MV=1000, c=0.1 *)
cf33={-1000, Table[100,{9}],1100}//Flatten
pv1[cf33,i]

pvbond=pv1[cf33//Flatten,i]//Simplify
Plot[pvbond,{i,0,0.2}]

IRR

sol=NSolve[pvbond,i]

vsechnykoreny=i/.sol

Clear[jh];
jh[x_]:=x/;x≥0&&x≤1.5;
jh[x_]:=-x/;x<0;
jh[x_]=Sin[x];
Plot[jh[x],{x,-3,3}]
jh[0.5]

{Cases[vsechnykoreny,z_Real/;z≥0]//First,
Cases[vsechnykoreny,z_Real/;z≥0][[1]]}

(* npv, duration, convexity,  IRR *)
Clear[npv,dur,con];
npv[cf_,i_]:=cf.(1/(1+i))^(Range[Length[cf]]-1)
dur[cf_,i_]:=(cf*(Range[Length[cf]]-1)).(1/(1+i))^(Range[Length[cf]]-1)/npv[cf,i]
con[cf_,i_]:=Module[{times=Range[Length[cf]]-1},(cf*times*(times+1)).(1/(1+i)^times)]/npv[cf,i]
altIRR[cf_]:=Module[{a},a=N[Solve[npv[cf,x]==0,x]];Cases[x/.a,z_Real/;z≥0]]

Vicenasobna simulace nahodne prochazky

ListPlot[NestList[#+Random[]-.5&,10,250],Joined->True]

SeedRandom[13];
nahodneprochazkygrafy=Table[ListPlot[NestList[#+Random[]-.5&,10,250],Joined->True,PlotRange->All],{100}];

Show[nahodneprochazkygrafy]

Table[ListPlot[NestList[#+Random[]-.5&,10,250],Joined->True,PlotRange->All],{100}]//Show

Nalezeni obalky

SeedRandom[13];
nahodneprochazky=Table[NestList[#+Random[]-.5&,10,250],{100}];

nahodneprochazky//First

Dimensions[nahodneprochazky]

Attributes[Max]

Max[{Range[5],Range[5]+10}]
Mean[{Range[5],Range[5]+10}]

obalkamax=Map[Max,nahodneprochazky//Transpose];
obalkamin=Map[Min,nahodneprochazky//Transpose];

pmax=ListPlot[obalkamax,PlotStyle->RGBColor[1,0,0],Joined->True,DisplayFunction->Identity];
pmin=ListPlot[obalkamin,PlotStyle->RGBColor[0,1,0],Joined->True,DisplayFunction->Identity];
Show[nahodneprochazkygrafy,pmax,pmin,DisplayFunction->$DisplayFunction]

MaxMemoryUsed[]

Statistika dosazenych cen (v case 251)

Map[Last,nahodneprochazky]//Length

konecneceny=Map[Last,nahodneprochazky];

Poznamka k lokalizaci promennych

$RecursionLimit
With[{$RecursionLimit=1000,nahodneprochazky=10},Print[$RecursionLimit,"   ",nahodneprochazky]]
$RecursionLimit

Clear[datareport];
Needs["Histograms`"];
(* Needs["Graphics`"]; *)
datareport[data_,noclass___]:=Module[{mi=Min[data],ma=Max[data],bin,h,k,le,mea,med,se},
k=If[noclass===Null,10,noclass];
bin=BinCounts[data,{mi,ma,h=((ma-mi)/k)}];Print["No. = ",le=Length[data],"\nMin = ",mi," Max = ",ma,"\nMean= ",mea=Mean[data]," Median = ",med=Median[data],"\nSE = ",se=StandardDeviation[data]," VaR 5% = ",Quantile[data,0.05`]," VaR 1% = ",Quantile[data,0.01`]];Histogram[data(*,HistogramCategories->k,(* BarStyle->{GrayLevel[0.95`]}, *)ApproximateIntervals->True*)]]

datareport[konecneceny]

(* Needs["Statistics`"] *)
Needs["Graphics`"]
LocationReport[konecneceny]
DispersionReport[konecneceny]
BarChart[konecneceny]
Histogram[konecneceny]
{Skewness[konecneceny],Kurtosis[konecneceny]}

(251./12)

Rozdeleni konecnych cen: kazda konecna cena je (∑_(n=1)^251)R(10-(1/2),10+(1/2)) nezavislych a z CLV ma rozdeleni N(10,(251/12)) .

BinCounts[konecneceny,{-10,30,5}]

(χ^2)-test dobre shody

Clear[chi2test];
chi2test[data_,k_,hypoF_]:=Module[{pom,cetnosti,n=Length[data],tot},
pom=Map[hypoF,data];
cetnosti=BinCounts[pom,{0,1,(1/k)}];
tot=Total[((((cetnosti-n/k))^2)/(n/k))];
Print["P-value = ",1.- CDF[ChiSquareDistribution[k-1],tot]]
]

chi2test[konecneceny,10,CDF[NormalDistribution[10,(Sqrt[(251/12)])],#]&]

Poznamka k mapovani

mat=RandomReal[{0,1},{10,10}];
(* radkova maxima *)
Map[Max,mat]
(* Sloupcova maxima *)
Map[Max,mat//Transpose]

Map[Sin,(({{a, b}, {c, d}}))]

Map[Max,(({{3, 5}, {4, 3}}))]

Max[(({{3, 5}, {4, 3}}))]

Sharpe ratio and Markowitz portfolio

Data

(* Needs["Statistics`"] *)
(* Needs["Statistics`MultiDescriptiveStatistics`"] *)
d1=Import["E:\Dokumenty\VPFPM\VPFPM2007 2008\Returns.dat","Table"];
d1[[{1,2,3,4}]]
d2=Rest[d1];
d3=Rest[d2];
datareturns=Map[Rest,d3];
expected=Mean[datareturns]
(cov=Covariance[datareturns])//MatrixForm
(Correlation[datareturns])//MatrixForm

Cviceni: Prostudujte nasledujici alternativni pristup k extrakci potrebnych dat.

ret1=Import["E:\Dokumenty\VPFPM\VPFPM2007 2008\Returns.dat","Table"];
Short[ret1,5]
ret2=Rest[Rest[ret1]];
Short[ret2,5]
returns=Map[Rest,ret2];
Short[returns,5]

Markowitz Mean-Variance Portfolio

We assume that N assets have expected returns given by the N×1vector r and the covariance matrix of returns V. We will seek the optimal portfolio in the sense that the Sharpe's ratio (also called Sharpe's measure of portfolio) is maximum:

max(((x^T)r)/(Sqrt[(x^T) V x]))subject to conditionsx≥0and(∑_(n=1)^N)(Subscript[x, n])=1,

so that short sales are not allowed.

Example 1 (short sales allowed)

Clear[xx];
xx=Table[Unique[],{Length[expected]}]

crit1=((xx.expected)/(Sqrt[xx.cov.xx]))//Simplify

?Maxim*

?NMax*

sol11=NMaximize[crit1,xx]

xa=xx/.sol11[[2]]
xnorm=xa/Total[xa]
Total[xnorm]

val1=crit1/.{xx->xnorm//Thread}

First[val1]
val1[[1]]

Example 2 (short sales not allowed)

ineq=xx≥0//Thread
ineq=Apply[And,ineq]&&Total[xx]==1

sol2=NMaximize[{crit1,ineq},xx]

xx/.Last[sol2]

Example 2 (Revisited)

ineq2=xx≥0//Thread
ineq2=Apply[And,ineq2]&&Total[xx]≤1
sol2r=NMaximize[{crit1,ineq2},xx]
Total[xx/.sol2r[[2]]]

crit1/.{$16->1.9682812537604517`*^-9,$17->3.7308353824431636`*^-10,$18->(0.4244081277770084`+(1-0.49334391693596863`)),$19->7.572991168693882`*^-11,$20->0.06893578674186551`}

Example 3a (Portfolio with prescribed expected return and minimum risk, short sales not allowed)

crit3=xx.cov.xx//Simplify
cons3=Apply[And,Thread[xx≥0]]&&Total[xx]==1&&xx.expected==4

sol3=NMinimize[{crit3,cons3},xx]

Total[xx3=xx/.Last[sol3]]

xx3
expected

xx3.expected

Example 3b (inequality instead of equality)

crit3b=xx.cov.xx//Simplify;
cons3b=Apply[And,Thread[xx≥0]]&&Total[xx]≤ 1&&xx.expected==4;
sol3b=NMinimize[{crit3b,cons3b},xx]
Total[xx3b=xx/.Last[sol3b]]
xx3b.expected

Example 3c (inequalities instead of equalities)

crit3c=xx.cov.xx//Simplify;
cons3c=Apply[And,Thread[xx≥0]]&&Total[xx]≤ 1&&xx.expected≥4;
sol3c=NMinimize[{crit3c,cons3c},xx]
Total[xx3c=xx/.Last[sol3c]]
xx3c
xx3c.expected

Total[xx3d={1.289839352143948`*^-32,-4.880587180699747`*^-33,0.8613009345035485574`,3.132574895506314`*^-34,0.13860717880556636`}]
xx3d.expected
xx3d.cov.xx3d

Unfeasible requirement
(Neither the asset with maximum expected return reaches 4.3.)

crit34=xx.cov.xx//Simplify;
cons34=Apply[And,Thread[xx≥0]]&&Total[xx]==1&&xx.expected≥ 4.3;
sol3=NMinimize[{crit34,cons34},xx]

Example 4

We have 8 assets traded on Prague Stock Exchange with the expected returns

returns={-9.92488533927622*10^-04,5.72524775783378*10^-03,9.04818959590086*10^-03,-1.44968966764537*10^-03,-1.74045537023702*10^-02,5.70188394543127*10^-03,3.33476815091662*10^-03,2.15075479228363*10^-03};

and the covariance matrix of returns

covariancematrix={{1.35640293495777*10^-03,2.56330097904098*10^-04,3.33037161321956*10^-04,1.80223328810003*10^-04,2.81325992644196*10^-04,5.93263545509886*10^-04,7.13610711759175*10^-04,3.20932535820455*10^-04},{2.56330097904098*10^-04,2.35439634554186*10^-03,4.42407197078591*10^-04,6.57567906193998*10^-06,4.01068448708833*10^-04,5.93785356343409*10^-04,1.30961261626908*10^-04,2.93338399396577*10^-04},{3.33037161321956*10^-04,4.42407197078591*10^-04,1.57895651834543*10^-03,-1.16804549527263*10^-04,1.85336104268171*10^-04,.000176059157865,3.41400743910472*10^-04,1.86624552570071*10^-04},{1.80223328810003*10^-04,6.57567906193998*10^-06,-1.16804549527263*10^-04,6.16786555694974*10^-03,4.84828341822248*10^-04,8.50333422445824*10^-04,7.74803938060802*10^-05,4.39412705238888*10^-04},{2.81325992644196*10^-04,4.01068448708833*10^-04,1.85336104268171*10^-04,4.84828341822248*10^-04,9.45147795376653*10^-03,1.78331918975426*10^-03,5.0987274991634*10^-05,1.49823427374484*10^-03},{5.93263545509886*10^-04,5.93785356343409*10^-04,.000176059157865,8.50333422445824*10^-04,1.78331918975426*10^-03,7.38754498647722*10^-03,1.03084013803777*10^-03,9.34434830622224*10^-04},{7.13610711759175*10^-04,1.30961261626908*10^-04,3.41400743910472*10^-04,7.74803938060802*10^-05,5.0987274991634*10^-05,1.03084013803777*10^-03,2.1036055256453*10^-03,2.40818723671803*10^-04},{3.20932535820455*10^-04,2.93338399396577*10^-04,1.86624552570071*10^-04,4.39412705238888*10^-04,1.49823427374484*10^-03,9.34434830622224*10^-04,2.40818723671803*10^-04,2.77726103166242*10^-03}};

We see that the above matrix is badly conditioned:

Det[covariancematrix]

Next we make use of the procedure NMaximize implemented in the Standard Numerical Package Mathematica, version 4.2. First we create the vector of unknown proportions in the portfolio:

xxPS=Table[(Subscript[x, i]),{i,1,Length[returns]}]

The criterion is given by

criterion=(((xxPS.returns))/(Sqrt[xxPS.covariancematrix.xxPS]))//Simplify;

This loads the package: (Not necessary in version 5.2)

The constraints are

constraints1=Apply[And,Table[(Subscript[x, i])>=0,{i,1,Length[returns]}]];
constraints2=Sum[(Subscript[x, i]),{i,1,Length[returns]}]==1;
constraints=constraints1&&(constraints2)

Next we solve the problem using the option for the method:

soldif=NMaximize[{criterion,constraints},xxPS](*,Method->"DifferentialEvolution"]*)

soldif=NMaximize[{criterion,constraints},xxPS,Method->"DifferentialEvolution"]

Apply[Plus,xxPS/.soldif[[2]]]
Print["Expected return of the optimal portfolio = ",(xxPS/.soldif[[2]]).returns]

Looking on the criterion, we see that the solution is independent of a nonzero multiple. So we can suspend the second constraint and make the

constraints3=Apply[And,Table[(Subscript[x, i])>=0,{i,1,Length[returns]}]]

soldif=NMaximize[{criterion,constraints3},xxPS ]

We see that the optimal value is slightly higher than in the above case. The optimal portfolio together with its expected return (substantially higher than in the previous case) is simply obtained:

sol=xxPS/.soldif[[2]];
sol=sol/Apply[Plus,sol]
sol.returns

0.007859525763768638`

Lagr. mult.

Example 3 (Markowitz, short sales allowed, μ prescribed expected return)

Clear[n,ones,xx,λ1,λ2,μ];
n=Length[expected];
ones=Table[1,{n}];
xx=Table[Unique[],{n}]
lagrange=(1/2)xx.cov.xx+λ1*(1-ones.xx)+λ2*(μ-expected.xx)//Simplify

Map[D[lagrange,#]&,xx]

expected
μ=3;
(* !!!! *)
eq1=Map[D[lagrange,#]&,xxx={xx,λ1,λ2}//Flatten]==0
Take[Flatten[xxx/.Solve[eq1,xxx]],Length[xx]]

Ruzne programovaci styly

Nasledujici funkce zjistuji cetnosti.

Proceduralni styl

SeedRandom[130];
data=Table[Random[Integer,{1,10}],{100}];

proc[x_]:=Module[{bins=Table[0,{10}]},
Do[bins[[x[[i]]]]++,{i,1,Length[x]}];bins]
pr=proc[data]
Total[pr]

Funkcionalni styl

Clear[fun];
fun[x_?VectorQ]:=fun[Map[{#,1}&,x]]
fun[{1,2,3,2}]

Poznamka: _ prave jeden argument, __ jeden a vice argumentu, ___ nula a vice argumentu, x___

fun[{a___,{b_,n_},c___,{b_,m_},d___}]:=fun[{a,{b,n+m},c,d}]

fun[{1,2,3,2}]

fun[a_List]:=Map[Last,a];
data
fun[data//Sort]

Mathematica styl

Table[Count[data,i],{i,10}]

Map[Count[data,#]&,Range[10]]

Tally[data]

Nahrazovaci pravidla (replacement rules)

{1,2,3,2,3,1,2,4}/. 3->"a"/.x_Integer->0

v1={1,2,3,2,3,1,2,4}/. 3->"a"/.x_Integer->0/."a"->1

Total[v1]

pravidla[data_List,n_]:=
Total[data/.n->"a"/.x_Integer->0/."a"->1]
pravidla[{1,2,3,2,3,1,2,4},2]

Map[pravidla[data,#]&,Range[10]]

Databaze

Filtrace

list2=Range[20];
Select[list2,#>17&]
Select[list2,9≤#≤12&]
Select[list2,#≤3&&#==2&]
Select[list2,And[#≤3,#==2]&]
Select[list2,#≥18||#==2&]

Vyber uzivajici kriterii

jmena1={"Petr","Monika","Tomas","Zuzana"};
prijmeni={"Svitekova","Sotona","Stastkova","Buncak","Valaskova"};
jmena=Prepend[jmena1,"Zuzana"]

(skupina=Transpose[{jmena,prijmeni}])//TableForm

zuzanatest[name_String]:=If[name=="Zuzana",True,False]
zuzanatest["Petr"]

{{"Zuzana","Svitekova"},{"Petr","Sotona"},{"Monika","Stastkova"},{"Tomas","Buncak"},{"Zuzana","Valaskova"}}

Map[StringJoin[#[[1]]," ",#[[2]]]&,skupina]//FullForm

Select[skupina,zuzanatest[#[[1]]]&]

Select[skupina,#[[2]]=="Buncak"&]//First

(* Priradime nahodne veky *)

veky=Table[Random[Integer,{18,26}],{Length[skupina]}]

Append[{"Zuzana","Svitekova"},12]

RowJoin[skupina,{veky}]

sk2=MapThread[Append,{skupina,veky}]

Select[sk2,#[[1]]=="Zuzana"&&#[[3]]>24&]//First

Associative Arrays JH

mat=(({{"names", "BC", "EFG", "HIJ", "KLM"}, {"a", 1, 2, 3, 4}, {"b", 5, 6, 7, 8}, {"c", 9, 10, 11, 12}, {"d", 13, 14, 15, 16}}));

Clear[data,firstrow];
firstrow=mat[[1]]
(* POZOR, z vystupu neni patrne, ze se jedna o promenne typu String ! Viz:*)

mat[[1]]//FullForm

fields=Map[data[#]&,firstrow]

matdrop=Drop[mat,1]//Transpose

MapThread[Set,{fields,matdrop}]

fields

data["KLM"]

Select[mat,#[[2]]>7&]
Select[fields,#[[4]]==16&]

??fields

fields//Dimensions

fields

{fields[[1,All]],fields[[All,1]]}

firstrow

Toto je obtiznejsi:

(symbols=ToExpression[firstrow])//FullForm
nnn=Length[symbols]
MapThread[Set,{symbols,Range[nnn]}]
?symbols
symbols

BC je nyni promenna:

BC

Financial Derivatives (Options)

Black-Scholes Formula

Notation: s ... current price of underlying asset at time t, k ... strike price, r ..  riskless interest rate, T ... expiry date, σ ... volatility of the underlying asset
call (put) ... the value of a CALL (PUT) option

Clear[Φ]
Φ[x_]:=(1/2) (1+Erf[(x/(Sqrt[2]))])
ϕ[x_]:=Exp[-x^2/2]/Sqrt[2*π]

d1=(Log[s/k]+(r+σ^2/2)*(T-t))/(σ*Sqrt[T-t]);
d2=(Log[s/k]+(r-σ^2/2)*(T-t))/(σ*Sqrt[T-t]);
call[s_,k_,σ_,T_,t_,r_]:=Module[{d1,d2},d1=(Log[s/k]+(r+σ^2/2)*(T-t))/(σ*Sqrt[T-t]);
d2=(Log[s/k]+(r-σ^2/2)*(T-t))/(σ*Sqrt[T-t]);
s*Φ[d1]-k*Exp[-r*(T-t)]*Φ[d2]]
put[s_,k_,σ_,T_,t_,r_]:=Module[{d1,d2},d1=(Log[s/k]+(r+σ^2/2)*(T-t))/(σ*Sqrt[T-t]);
d2=(Log[s/k]+(r-σ^2/2)*(T-t))/(σ*Sqrt[T-t]);
k*Exp[-r*(T-t)]*Φ[-d2]-s*Φ[-d1]]

?Put

Implied volatility: price ... the market price of an option at time t

impliedcall[s_,k_,T_,t_,r_,price_]:=sd/.FindRoot[call[s,k,sd,T,t,r]==price,{sd,0.2}]
impliedput[s_,k_,T_,t_,r_,price_]:=sd/.FindRoot[put[s,k,sd,T,t,r]==price,{sd,0.2}]

call[52.44,50,0.2,20/52,0,0.05]
put[49,50,0.2,20/52,0,0.05]

plcall=Plot[call[s,50,0.2,20/52,0,0.05],{s,40,100},PlotRange->All]
plput=Plot[put[s,50,0.2,20/52,0,0.05],{s,40,100},PlotRange->All]
plintr=Plot[s-50,{s,40,100},PlotStyle->{Hue[0.3]}]
Show[plcall,plput,plintr]

Plot[call[s,46,0.3,28/360,0,0.05],{s,40,53},PlotRange->All];
Plot[call[s,50,0.3,28/360,0,0.05],{s,40,53},PlotRange->All];

Plot[{call[(* s,k,σ,T,t,r *)s,50,0.2,22/52,0,0.05],put[s,50,0.2,20/52,0,0.05]},{s,0,70},PlotRange->All];

Actual data from the German stock exchange

impliedcall[70.72,k1=67.5,T=58/360,t=0,r=0.05,c1=6.31]
impliedcall[70.72,k2=70,58/360,0,0.05,c2=4.92]
impliedcall[70.72,k3=72.5,58/360,0,0.05,c3=3.37]
impliedput[70.72,67.5,58/360,0,0.05,p1=2.92]
impliedput[70.72,70,58/360,0,0.05,p2=4.08]
impliedput[70.72,72.5,58/360,0,0.05,p3=5.48]

{k1,T,r}

Definition of the object OPTION:
Real data: Underlying Deutsche Bank April 23, 1999 - Price EUR 52.44 
pos=+1 long, pos=-1 short

optionDB46c[pos_,S_]={pos,1,5.26,S,46};
optionDB48c[pos_,S_]={pos,1,3.75,S,48};
optionDB50c[pos_,S_]={pos,1,2.51,S,50};
optionDB46p[pos_,S_]={pos,-1,0.74,S,46};
optionDB48p[pos_,S_]={pos,-1,1.39,S,48};
optionDB50p[pos_,S_]={pos,-1,2.33,S,50};

optionDB46c[-1,s]

Portfolio of options

Butterfly spread (long CALL k1, long PUT k3, short CALL k2, short PUT k2)

valuepofportfolio[s_]:=-c1-p3+c2+p2+Max[s-k1,0]+Max[k3-s,0]-Max[s-k2,0]-Max[k2-s,0]

plbut1=Plot[valuepofportfolio[s],{s,55,85}];

valueBS[s_,t_]:=-c1-p3+c2+p2 +call[s,k1=67.5,0.2,T=58/360,t,r=0.05]+put[s,72.5,0.2,58/360,t,0.05]-call[s,k2=70,0.2,58/360,t,0.05]-put[s,70,0.2,58/360,t,0.05]

valueBS[s,t]//Simplify

plbut2=Plot[valueBS[s,30/360],{s,55,85}];

Show[plbut1,plbut2,PlotRange->All]

Plot[valueBS[s,(58-0.01)/360],{s,55,86},PlotRange->All];

tabBS=Table[Plot[valueBS[s,(t)/360],{s,55,86},PlotRange->All,DisplayFunction->Identity],{t,0,55,5}];

Show[tabBS,DisplayFunction->$DisplayFunction]

Sensitivities (The greeks)

Clear[s,k,σ,T,t,r];
call[s,k,σ,T,t,r]//Simplify

Delta

(Subscript[Δ, c])=D[call[s,k,σ,T,t,r],s]//Simplify

Plot[(Subscript[Δ, c])/.{k->50,σ->0.2,T->20/52,t->0,r->0.05},{s,20,70}]

rho=D[call[s,k,σ,T,t,r],r]//FullSimplify

Plot3D[rho/.{k->50,σ->0.2,T->20/52,t->0},{s,20,70},{r,0.02,0.07}]

Plot[Table[rho/.{k->50,σ->0.2,T->20/52,t->0,s->j},{j,20,70,10}]//Evaluate,{r,0.02,0.07}];

Vyzkousejte bez Evaluate a podivejte se do Helpu.

Binomial trees

In this Section we closely follow Chapter 18, Shaw, W.: Modelling Financial Derivatives with Mathematica.

Introduction

Let us consider a time interval (t,t+Δt). Under the Black-Scholes model the conditional distribution of (Subscript[S, t+Δt]) given (Subscript[S, t]) is log-normal LN(ln (Subscript[S, t])+(μ-(1/(σ^2)))Δt,(σ^2)Δt). In a risk neutral world μ=r, where r is the riskless rate and we have

E((Subscript[S, t+Δt])|(Subscript[S, t]))=(Subscript[S, t])(e^(r Δt))=:(Subscript[E, c])     and     var((Subscript[S, t+Δt])|(Subscript[S, t]))=((Subscript[S, t])^2)(e^(2r Δt))((e^((σ^2)Δt))-1)=:(Subscript[V, c]).

Under the discrete binomial random walk we suppose that the stock price (Subscript[S, t]) moves to one of two new values, (Subscript[S, t+Δt])=(Subscript[S, t])u and (Subscript[S, t+Δt])=(Subscript[S, t])d where d<1+r<u with probabilities p and q:=1-p, respectively. In this case we have

E((Subscript[S, t+Δt])|(Subscript[S, t]))=(Subscript[S, t])(p u + q d)=:(Subscript[E, b])    and     var((Subscript[S, t+Δt])|(Subscript[S, t]))=((Subscript[S, t])^2)p ((q(u-d))^2)=(Subscript[V, b]).

While in the continuous model we have two parameters r and (σ^2), ceteris paribus, in the discrete model we have three parameters, i.e., u, d, p. In practice we choose the parameters of the binomial model to match the parameters of the continuous one in terms of the moments:

(Subscript[E, b])=(Subscript[E, c])     and     (Subscript[V, b])=(Subscript[V, c]).

The reason is simple; usually we can estimate r and (σ^2) directly which is hardly possible for the binomial model. Since there is one more degree of freedom, we have to postulate a reasonable relationship between the parameters involved. The corresponding system of equations now reads:

Recall the functions call and put for the European options in a previous lecture:

Clear[Φ]
Φ[x_]:=(1/2) (1+Erf[(x/(Sqrt[2]))])
ϕ[x_]:=Exp[-x^2/2]/Sqrt[2*π]

d1=(Log[s/k]+(r+σ^2/2)*(T-t))/(σ*Sqrt[T-t]);
d2=(Log[s/k]+(r-σ^2/2)*(T-t))/(σ*Sqrt[T-t]);
call[s_,k_,σ_,T_,t_,r_]:=Module[{d1,d2},d1=(Log[s/k]+(r+σ^2/2)*(T-t))/(σ*Sqrt[T-t]);
d2=(Log[s/k]+(r-σ^2/2)*(T-t))/(σ*Sqrt[T-t]);
s*Φ[d1]-k*Exp[-r*(T-t)]*Φ[d2]]
put[s_,k_,σ_,T_,t_,r_]:=Module[{d1,d2},d1=(Log[s/k]+(r+σ^2/2)*(T-t))/(σ*Sqrt[T-t]);
d2=(Log[s/k]+(r-σ^2/2)*(T-t))/(σ*Sqrt[T-t]);
k*Exp[-r*(T-t)]*Φ[-d2]-s*Φ[-d1]]

Clear[p,u,d,r,σ]
eqs={p*u+(1-p)*d==Exp[r*Δt],p*(1-p)*(((u-d))^2)==Exp[2*r*Δt]*(Exp[(σ^2)Δt]-1)};

Jarrow-Rudd scheme

Under this scheme it is supposed that both up and down movements have the same probability.

Clear[p,u,d,r,σ]
p=(1/2);solJR=Solve[eqs,{u,d}]

Only the largest solution makes sense:

{uJR,dJR}={u,d}/.solJR[[2]]

This is an approximation:

approxJR={Series[uJR,{Δt,0,1}]//Normal//PowerExpand,
Series[dJR,{Δt,0,1}]//Normal//PowerExpand}

which is very close to frequently used

Exp[approxJR]/E

Cox-Ross-Rubinstein scheme

Under this scheme it is supposed that u=1/d:

Clear[p,u,d,r,σ]
u=(1/d);solCRR=Solve[eqs,{p,d}]//PowerExpand//FullSimplify

After a hand-made simplification (not too difficult) we get

Clear[CRRupBin];
CRRupBin[n_, r_, q_, σ_, T_] := Module[{f = Exp[-(r-q)*T/n], 
g = Exp[(r-q+σ^2)*T/n]},
1/2*(f+g + Sqrt[(f+g)^2-4])];
Clear[CRRdownBin];
CRRdownBin[n_, r_, q_, σ_, T_] := Module[{f = Exp[-(r-q)*T/n], 
g = Exp[(r-q+σ^2)*T/n]},
1/2*(f+g - Sqrt[(f+g)^2-4])];
Discount[n_, r_, T_] := N[Exp[-r T/n]];

Binomial tree for American options

Let us consider an n steps' evolution, time to expiry T and the time step Δt=T/n. Denote (i,j) the jth node at time iΔt, i=1,...,n, j=0,...,i. Let (Subscript[V, i,j]) be the value of the option (either CALL or PUT) at node (i,j). The price of the stock is (Subscript[S, 0])(u^j)(d^(i-j)) at that node. Since the value for an American PUT is max(K-(Subscript[S, T]),0) at expiry, the value is

(Subscript[V, n,j])=max(K-(Subscript[S, 0])(u^j)(d^(n-j)),0),     j=0,1,...,n.

With probability p the system moves from node (i,j) at time iΔt to (i+1,j+1) at time (i+1)Δt and with probability 1-p to (i+1,j). Without possibility of an early exercise we get

(Subscript[V, i,j])=(e^(-r Δt))(p (Subscript[V, i+1,j+1])+(1-p)(Subscript[V, i+1,j])),   i=0,...,n-1,   j=0,...,i.

If an early exercise is taken into account, the value must be compared with its intrinsic value. Thus we finally get

(Subscript[V, i,j])=(e^(-r Δt))max{K-(Subscript[S, 0])(u^j)(d^(i-j)),(e^(-r Δt))(p (Subscript[V, i+1,j+1])+(1-p)(Subscript[V, i+1,j]))},   i=0,...,n-1,   j=0,...,i.

Starting from level n we go backward through the tree to the desired value (Subscript[V, 0,0]).

Implementation

Here also a dividend yield q is considered.

$RecursionLimit

$RecursionLimit=5000;

afunc[n_, r_, q_, T_] := N[Exp[(r-q) T/n]];
Pfunc[up_,down_,a_,disc_]:=N[disc*(a-down)/(up-down)];
Qfunc[up_,down_,a_,disc_]:=N[disc-Pfunc[up,down,a,disc]];

Clear[CRRAmericanOption]
CRRAmericanOption[S_, n_, σ_, T_, r_, q_, exercise_Function] :=
Module[{u = CRRupBin[n,r,q,σ,T], d = CRRdownBin[n,r,q,σ,T],
a = afunc[n,r,q,T],
disc = Discount[n,r,T],
P,Q,TreeOp,value,delta,gamma,theta},
P = Pfunc[u,d,a,disc];
Q = Qfunc[u,d,a,disc];
TreeOp[node_,level_] := TreeOp[node,level] = 
If[level==n,exercise[S d^node u^(level-node)],
Max[{P,Q}.{TreeOp[node,level+1],TreeOp[node+1,level+1]},
exercise[S d^node u^(level-node)]]];
value=TreeOp[0,0];
delta = (TreeOp[0,1] - TreeOp[1,1])/((u-d)*S);
gamma = ((TreeOp[0,2] - TreeOp[1,2])/(u^2-1) - 
(TreeOp[1,2]-TreeOp[2,2])/(1-d^2))*2/(S^2 (u^2-d^2));
theta = (TreeOp[1,2] - TreeOp[0,0])/(2*T/n);
Clear[TreeOp];
{value, delta, gamma, theta}];

CRRAmericanPut[S_,K_,n_,r_,q_,σ_,T_]:=CRRAmericanOption[S,n,σ,T,r,q,Max[K-#,0]&];
CRRAmericanCall[S_,K_,n_,r_,q_,σ_,T_]:=CRRAmericanOption[S,n,σ,T,r,q,Max[#-K,0]&];

Recall the meaning of the parameters:

CRRAmericanCall[S_,K_,n_,r_,q_,σ_,T_]
call[s_,k_,σ_,T_,t_,r_]

CRRAmericanCall[9, 10,100, 0.1, 0, 0.2, 1]
call[9,10,0.2,1,0,0.1]

Look at the convergence:

Table[{m,CRRAmericanCall[9, 10,m, 0.1, 0, 0.2, 1][[1]]//Evaluate},{m,3,8}]

The following plot (number of steps equal to 80) shows that the convergence is not too convincing and that the recommended number of steps 30 may not be sufficient.

ListPlot[Table[{m,CRRAmericanCall[9, 10,m, 0.1, 0, 0.2, 1][[1]]//Evaluate},{m,3,80}],Joined->True]

The following Example comes from HULL, p.391, where the approximate values u=(e^(σ(Sqrt[Δt]))), d=1/u are used. His result (after 5 steps) is 4.49. A comparison with usual BS formula for an European CALL is also given:

CRRAmericanCall[S_,K_,n_,r_,q_,σ_,T_]
call[s_,k_,σ_,T_,t_,r_]

CRRAmericanPut[50,50,50,0.1,0,0.4,0.4167]
put[50,50,0.4,0.4167,0,0.1]

Recurssion (Fibonacci)

?Fibonacci

{Fibonacci[1],Fibonacci[2]}

fibA[1]=fibA[2]=1;
fibA[n_]:=fibA[n-1]+fibA[n-2]

(tabA=Table[{n,Timing[fibA[n]][[1]]/Second},{n,20,28}])//TableForm

ListPlot[tabA,PlotJoined->True]

?fibA

fibB[1]=fibB[2]=1;
fibB[n_]:=fibB[n]=fibB[n-1]+fibB[n-2]

(tabB=Table[{n,Timing[fibB[n]][[1]]/Second},{n,100,105}])//TableForm
ListPlot[tabB]

?fibB

fibC[n_]:=First[MatrixPower[(({{1, 1}, {1, 0}})),n-2].{{1},{1}}]

fibC[100]

Lekce 11 (2008)

FAQ

Kresleni dvou grafu

Needs["PlotLegends`"]

p1=Plot[Sin[(x^2)],{x,0,π}]
p2=Plot[Cos[x],{x,0,Pi} ,PlotLegend-> {"cos x"},PlotRange->  All ]

Show[p2,p1]

Show[Plot[Cos[x],{x,0,π} ,PlotLegend-> {"cos x"},PlotRange->  All ],Plot[Sin[(x^2)],{x,0,π}]]

Plot[{Cos[x],Sin[(x^2)]},{x,0,π},PlotRange->  All,PlotLegend-> {"cos x","haf"} ]

Plot[x^2,{x,0,3.5},PlotLegend->{"haf"}]

tab=Table[{x,Sin[x]},{x,0.,π,(π/6)}]

inttab=Interpolation[tab, InterpolationOrder->1]

inttab[2]

Plot[inttab[x],{x,0,π}]

Plot[{inttab[x],Sin[x]},{x,0,π},PlotLegend->{"HAf","Hafhaf"}]

Zobrazeni cisel s mantisou

FortranForm[π//N]

ScientificForm[{123450000.0,0.000012345,123.45}]

StringReplace["237E-10","E"->" 10^"]//ToExpression

AppendTo[$Path,"E:\Dokumenty\VPFPM\VPFPM2007 2008"]

Import["Haf.xls"]

Import["Haf.xls"]//Flatten//First

Solve a FindRoot

s1=Solve[(x^3)+(x^2)+x==1,x]
x/.s1
(x/.s1)//First

s2=Solve[Sum[(x^i)+i,{i,5}]==0,x]

s2//N

s3=NSolve[Sum[(x^i)+i,{i,5}]==0,x]

s4=NSolve[(E^x)==x,x]

(E^x)==x//FullForm

FullForm[a-b]

Clear[f,g];
Map[f,g[x]]
Apply[f,g[x]]

(s41=Apply[Plus,(E^x)==x])

s42=Apply[Plus[#1,-#2]&,(E^x)==x]

s42/.s4

s42/.s4//Chop//First

s5=Solve[eq1=((Sin[x+(E^x)])^(Sin[x+1]))==0,x]

s51=NSolve[eq1=((Sin[x+(E^x)])^(Sin[x+1]))==0,x]

s52=FindRoot[eq1=((Sin[x+(E^x)])^(Sin[x+1]))==0,{x,0.5}]

Viz tez Reduce

Plot[((Sin[x+(E^x)])^(Sin[x+1])),{x,0,π},PlotRange->All]

Table[{y,FindRoot[eq1=((Sin[x+(E^x)])^(Sin[x+1]))==0,{x,y}]},{y,0,π,(π/12)}]

Options[FindRoot]

Specifikace oboru pripustnych hodnot parametru

(∫_0^∞)(E^(-a (x^2)))\[DifferentialD]x

Assuming[a>0,(∫_0^∞)(E^(-a (x^2)))\[DifferentialD]x]

Integrate[(E^(-a (x^2))),{x,0,Infinity},Assumptions->a>0]

Assuming[Re[a]>0,(∫_0^∞)(E^(-a (x^2)))\[DifferentialD]x]

(∫_0^∞)(E^(-a (x^2)))\[DifferentialD]x/.a-> 1+I//FullSimplify

Graf s popisy os

Plot[x,{x,-2,2},AxesLabel->{"a","bb aa"}]

Options[Plot]

Obrazky v TeXu:

\vglue 4 mm
\epsfysize = 4 cm
\centerline {\epsfbox {cf1.eps}}
\vskip 2 mm
\centerline {{\sl Figure 1 :} A cash flow}
\vskip 2 mm

Doporucene ve Wolfram Demonstration Project

Monte Carlo
finance
Black Scholes

Lekce 12 (2008)

FAQ

Cestina v grafech

Plot[Sin[(x^2)],{x,0,π},AxesLabel->{"č ě ž" , "ýžřčč šč"}]

list=Range[5]
Do[list[[i]]=list[[i]]+1,{i,1,5}]
list

list+1