VPFPM 2008

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

Úplný začátek

Podívejte se na Help ⟶ Documentation Center -> First Five Minutes

Pokročilejší přístup

Help ⟶ Záložka "Getting Started"´. Projít všechny varianty.
Užitečné a z matematického hlediska zajímavé informace jsou  v záložce "Demos".

Důkladný průvodce

Help ⟶  Záložka "Tour". Také zde je prospěšné projít všechny varianty.

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), případně Tour.nb. 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

Motivace, odkazy

www.wolfram.com

"C:\Program Files\Wolfram Research\Mathematica\5.2\Documentation\English\Tour"

4 pilire

Numericke vypocty

Vypocet pomoci Shift + Enter

1+1

(2^1000)

N[π,20]

Short[npi,10]

Sdruzovani vice vypoctu do jedne bunky

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

Symbolicke vypocty

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

e2=Factor[e1]

TeXForm[e2/e1]

Grafika

Plot[Sin[x+x^2],{x,0,π}]
Plot3D[Sin[x+y^2],{x,0,π},{y,0,π},ViewPoint->{-0.007, -1.850, 2.833}]

Options[Plot3D]

Programovani

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.]

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

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

hu1[x_]:=(x/2)

hu1[y]

?hu1

Remove[hu1]
hu1[y]
?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}]

hu4[x]=2x

{hu4[5],hu4[y],hu4[x]}

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]]

sdfghjklertyuikjjfdjfkhh=3;

sdfghjklertyuikjjfdjfkhh 3

MatrixForm[h3c=hil1[3]]

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

(hil15=hil1[5])//MatrixForm

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

?Global`*

Konec 20. 2. 2008

Lekce 2

Mathematica book v papirove podobe v knihovne

FAQ

Technicke zalezitosti

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

N znamena numerickou hodnotu:

{π,π//N,Log[2.],Log[2],N[Log[2]],N[Log[2],10]}

0.6931471805599453`
0.6931471805599453093744803656`10.

{Det[hil15],Det[hil15]//N,Det[hil15//N]}

{Det[h40=hil1[40]],N[Det[hil1[40]]],Det[hil1[40]//N]}

Det[hil1[40]] Det[Inverse[hil1[40]]]

Cas dole v rohu a beep: Edit⟶Preferences⟶Notebook Options⟶Evaluation Options⟶EvaluationCompletionAct⟶{"Beep", "ShowTiming"}

Inverse[h40][[30,30]]

{x*y,x y,xy}

aa=(({{a, b}, {c, d}}))

Esc[[Esc dava [[

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

Poznamka k exportu a importu:

$ImportFormats

$ImportFormats//Length

$ExportFormats//Length

Import["F:\\Dokumenty\\VPFPM\\Text.RTF"]

Technicke zalezitosti, vstup a vystup

To, ze probiha vypocet se pozna podle zlute bunky.
Zruseni vypoctu: Kernel⟶Abort Evaluation nebo Kernel⟶Quit Kernel

inv=Inverse[hil1[50]];

inv[[1,3]]

FullForm[inv]

?hil1

(* Strednik znamena potlaceny vystup, ale vypocet se provede. *)
inv3=Inverse[hil1[3]];

inv3

inv30=Inverse[hil1[30]];

inv30

Doplneni vyrazu: Ctrl+K

hil

IdentityMatrix[3]

Typy cisel: Integer, Rational (= (Integer/Integer)), Real, Complex

{2,2/3,2./3,π,1+2 I,1+2 I,((1+I)/(((2-I))^3))}

Rationalize[0.6666666666666666`]

Rationalize[2007.,0.6666666666666666]

Listovatelne funkce

Log[{2,3.,4}]

Map[Rationalize,{2007.,0.6666666666666666}]

(E^(2π I))

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

i2=Integrate[(E^(-c x^2))/2,{x,0,∞},Assumptions->c∈Reals]

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

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

Zavorky: [ ] argumenty funkci, [[ ]] elementy listu, { } listy, ( ) priorita algebraickych operaci

aa[1,1]

aa(1,1)

Map[MatrixForm,{aa {1,2},aa*{1,2}}]
aa .{1,2}

*  v maticovem poctu oznacuje Hadamarduv soucin (nasobeni po slozkach)

aa*aa

aa.aa

Psani textu, sdruzovani bunek apod.

2D zapis z palety nebo Ctrl + ...

(a/b)

(c/d)

(∑_(i=1)^n)(I^6)

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

Mathematica vhodne interpretuje vektory a matice

(aa=(({{a, b}, {c, d}}))).(({{e}, {f}}))

{MatrixForm[aa],(({{a, b}, {c, d}}))}

Duvod

(({{e}, {f}}))//FullForm

(({{a, b}, {c, d}}))(({{e, f}, {g, h}}))//MatrixForm

CDF je zkratka pro Cumulative Distribution Function, tj. ditribucni funkci: F(x)=P{X<x}.

Psani vzorcu: (∫_a^b)(x^2)dx je co? Scitali jsme (∑_(i=1)^n)(i^8) a dostali jsme prekvapivy vysledek.

(∫_a^b)(x^2)\[DifferentialD]x

Hold[(∫_a^b)(x^2)\[DifferentialD]x]//TeXForm

\int_a^b x^2 \,
   dx

TeXForm[Psani vzorcu: (∫_a^b)(x^2)\[DifferentialD]x je co? Scitali jsme (∑_(i=1)^n)(i^8) a dostali jsme prekvapivy vysledek.]

(∫_a^b)(x^2)\[DifferentialD]x=sin x

(∫_a^b)(x^2)\[DifferentialD]x

Popisky k obrazkum

Plot[Sin[(x^2)],{x,0,2 π},AxesLabel->{"x","Hodnoty funkce"},PlotLabel->"Sin[(x^2)]",Ticks->{{0,(π/2),π,((3 π)/2),2 π},{0,0.5`,0.75`,1.`}},PlotRange->All]

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

?FresnelS

Clear[fF];
fF[y_]:=Integrate[Sin[(x^2)],{x,0,y}]

fF[2.π]

?Global`*

maxfF=FindMaximum[fF[x],{x,π/2}]

{maxfF[[2]],Last[maxfF]}

x/.Last[maxfF]

?x

y=x/.Last[maxfF]

?y

maxfFexact=Maximize[{fF[x],2π>x>0},x]

Needs["Graphics`"]
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`]},HeadCenter->0,HeadLength->0.03`],Text["Sin(x)",{3.5`,0.8`},{0,-1}],Arrow[{1,-0.75`},{(π/2),Cos[(π/2)]}],Text["Cos(x)",{1,-0.9`},{0,-1}]}]

Operace s maticemi

Maticovy soucin: . (tecka, dot)

(({{a, b}, {c, d}})).(({{e}, {f}}))//MatrixForm

(({{a, b}, {c, d}})).(({{e, f}, {g, h}}))//MatrixForm

{e,f}.(({{a, b}, {c, d}}))

(({{a, b}, {c, d}})).{e,f}

Transponovat lze pouze matice a ne vektory!!

Transpose[{e,f}]

Transpose[{{e,f}}(* matice 1×2 *)] (* matice 2×1 *)

{{e,f}}//MatrixForm

x=.
xx=Table[(Subscript[x, i]),{i,3}]

xx.xx

Transpose[xx]

Transpose[{xx}]//MatrixForm

Transpose[{xx}].{xx}//MatrixForm

Zopakujme

xx*xx
xx.xx
y=.
Outer[Times,xx,yy={(Subscript[y, 1]),(Subscript[y, 2]),(Subscript[y, 3])}]//MatrixForm
yyxx=Transpose[{yy}].{xx}

GridBox[yyxx,RowLines->True,ColumnLines->True]//FrameBox//DisplayForm

Outer[Plus,xx,yy]

Konec 27. 2. 2008

Lekce 3 (5. 3. 2008)

FAQ

Prehlednejsi vystupy maticovych operaci

Clear[x,y,xx,yy];
xx=Map[(Subscript[x, #])&,Range[3]]
Outer[Times,xx,yy={(Subscript[y, 1]),(Subscript[y, 2]),(Subscript[y, 3])}]//MatrixForm
yyxx=Transpose[{yy}].{xx}

MatrixForm[yyxx,TableSpacing->{6,4}]

Take

yyxx//MatrixForm[#,TableSpacing->{0,4}]&

yyxx//GridBox//DisplayForm

GridBox[yyxx,ColumnSpacings->5]//DisplayForm

GridBox[yyxx,RowLines->True,ColumnLines->True]//FrameBox//DisplayForm

Demonstrations

Integrated Data Source

FinancialData["GE"]

Financial Data (GE)

fGE=FinancialData["GE",{2007,1,31}]

Financial Data (cont)

DateListPlot[FinancialData["GE",{1999,1,31}]]

f1=FinancialData["EUR/USD",{2008,2,1}]

f2=Map[Last,f1]

fmore=Map[Last,FinancialData["EUR/USD",{2000,2,1}]];

ListPlot[fmore,Joined->True]

j=Joined->True;
ListPlot[fmore,j]

DateList[{2008,1,4}]

FinancialData["EUR/USD","1/1/08"]

Country Data

CountryData["SA","Flag"]

CountryData["CZ","Population"]

Jak zaridit, abychom nemuseli stale vypisovat Joined->True?

SetOptions[ListPlot,Joined->True]

Options[ListPlot]

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

j=Joined->True;
ListPlot[fmore,j]

listplot[data_]:=ListPlot[data,Joined->True]
listplot[fmore]

Prace s bunkami (cells, uzavirani, oznacovani, automatizace)

Ctrlu (underlined), Ctrli (italics), Ctrlb (bold)

Notace

Vyuziti F1

Inverse

Alternativni definice Hilbertovy matice (Cviceni) a ryzi funkce

Ryzi funkce (Pure function) je funkce bez nazvu:
Function[x, (x^2)] zkracene ((#)^2)&
Pro dva a vice argumentu
Function[{x,y},x-y]

Clear[hil6,hil5]; (* Nasledujici funkce jsou ekvivalentni *)
hil5[n_]:=(1/(Array[Function[{x,y},x+y-1],{n,n}]))
hil6[n_]:=Array[(1/(#1+#2-1))&,{n,n}]
(* # argument ryzi funkce (pure function) *)

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

((#)^2)&

Sin[#1] Cos[#2]&[2,3]
Sin[#1] Cos[#2]&[2.,3]

{hil5[3],hil6[3]}

Definice Hilbertovy matice v duchu systemu Mathematica

Outer[Plus,Range[3],Range[3]]-1//MatrixForm

Clear[hilbertouter];
hilbertouter[n_]:=(1/(Outer[Plus,Range[n],Range[n]]-1))
hilbertouter[5]//MatrixForm

Dalsi operace s maticemi

CharacterRange["A","Z"]//FullForm

aachar=CharacterRange["a","d"]
FullForm[aachar]
aa=aachar//ToExpression
FullForm[aa]
aa22=Partition[aa,2]
Partition[aa,3]
Flatten[aa22]

aachar.aachar
aachar.aachar//FullForm
aa.aa
aa.aa//FullForm

Cviceni

Clear[a,b,c,d];
aa=(({{a, b}, {c, d}}))

s1=aa//ToString
s2=Map[ToString,aa]
s3=Map[ToString,aa,{2}]
s4=Flatten[aa]//ToString
s5=Map[ToString,Flatten[aa]]
s6=s5//ToExpression

Map[FullForm,{s1,s2,s3,s4,s5,s6}]//TableForm

FullForm[s1]
FullForm[s2]
FullForm[s3]
FullForm[s4]
FullForm[s5]
FullForm[s6]

Konec cviceni

Vlastni cisla a vlastni hodnoty, pseudoinverse

aa=Partition[aa,2];

Eigenvalues[aa]

Eigenvectors[aa]//PowerExpand//Simplify

hil5[5]

Solve[Det[hil5[3]-λ IdentityMatrix[3]]==0,λ]//Simplify
λ/.%
%//N

{1.408318927123654`+4.934324553889585`*^-18 I,0.002687340355773592`+0.` I,0.12232706585390582`+7.894919286223336`*^-17 I}//Chop

Eigenvalues[hil5[3]]//N

Apply[Times,Eigenvalues[hil5[5]]]* Apply[Times,Eigenvalues[hil5[5]//Inverse]]//N

h7=Map[hil5[#]&,Range[7]];

eigen7=Map[Eigenvalues,h7//N]

Vypocet determinantu

detpresny=Map[Det,h7]

detnepresny=Map[Apply[Times,#]&,eigen7]

(* Kontrola presnosti, relativni presnost *)
((detpresny-detnepresny)/detpresny)

e11=Expand[Product[x^2-2 i,{i,7}]]
Solve[e11==0,x]

Inverse[aa]//MatrixForm

aa2=aa/.{b->a,d->c}

Inverse[aa2]

PseudoInverse[aa2]

mm={{1,2},{1,2}};
PseudoInverse[mm]//MatrixForm
mm.PseudoInverse[mm].mm
mm.PseudoInverse[mm].mm==mm

{2.==2,2.===2}
{a==b,a===b}

protoze neni mozne rozhodnout, jakych hodnot a a b nabyvaji:

?a

?b

Generovani pozitivne definitni matice, vlastni cisla a vektory

aa100=RandomReal[{0,1},{100,100}];
Eigenvalues[posdefmat=Transpose[aa100].aa100]

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

Head[g[a,b,c]]
Map[f,g[a,b,c]]
Map[Plus,g[a,b,c]]

Total[{a,b,c}]

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

Apply[f,g[a,b,c]]

Apply[Times,a+b+c]

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

Aplikace obou funkci je u vypoctu detnepresny.

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//Flatten//First

dt={{{a}},b,c}//Flatten

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

Konec 5. 3. 2008, chyby opraveny

Lekce 4 (12. 3. 2008)

Odkazy

Obecne:

http://www.wolfram.com

http://www.wolfram.cz

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

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

Velmi dukladna uvodni instruktaz, viz take dalsi zajimave soubory v adresari Demos, GetingStarted:

"C:\Program Files\Wolfram Research\Mathematica\5.2\Documentation\English\Tour\Tour.nb"

Dalsi zajimave zdroje pouceni

"C:\Program Files\Wolfram Research\Mathematica\5.2\AddOns\ExtraPackages"

Help⟶ Master Index ⟶ Built-in Functions ⟶ Advanced Documentation

Konecne je take vyhodne se podivat na nektere zdrojove programy ze Standard Add On Packages. Nektere jsou snadno srozumitelne, u nekterych je tomu naopak.

FAQ

Rozdil mezi = (prirazeni, set), := (odlozene prirazeni, set delayed).

r1=Random[]
r2:=Random[]
Table[r1,{5}]
Table[r2,{5}]

Table[Random[],{5}]

Table[RandomInteger[],{2000000}]//Mean//N

Cviceni

Clear[a,b,c,d];
aa=(({{a, b}, {c, d}}))

s1=aa//ToString
s2=Map[ToString,aa]
s3=Map[ToString,aa,{2}]
s4=Flatten[aa]//ToString
s5=Map[ToString,Flatten[aa]]
s6=s5//ToExpression

Map[FullForm,{s1,s2,s3,s4,s5,s6}]//TableForm

FullForm[s1]
FullForm[s2]
FullForm[s3]
FullForm[s4]
FullForm[s5]
FullForm[s6]

Konec cviceni

Options

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

GoldenRatio//N

Options[Plot]

Toto ocividne nedava kruznici:

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

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

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

Statistika

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[]

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]

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]

Pouziti funkce ListPlot

Reverse[{a,b}]

Map[f,{a,b,c}]
Map[f,{{a,b},{c,d}}]

fr5rev=Map[Reverse,fr5]

ListPlot[fr5rev,Joined->True]

Ovsem prehlednejsi je

ListPlot[fr5//Sort,Joined->True]

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,13},PlotRange->All]

Show[p1,p2]

FindMinimum[fit5,{x,9}]

Needs["LinearRegression`"]

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

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]

Nelinearni regrese I

Needs["NonlinearRegression`"]

This set of data come from real world. They are in the form {maturity (number of days), spot interest rate}.

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

int3=Interpolation[ir3];
int31=Interpolation[ir3,InterpolationOrder->1];
plotint3=Plot[int3[t],{t,7,10000}];
Show[plot3,plotint3]

Nelinearni regrese II

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.

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.

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

Packages in Mathematica 6.0.1

Since we were not able to find the complete list of packages, we present a (stupid) method for finding it.

$UserDocumentsDirectory

$Path

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

FileNames[]

SetDirectory[$UserDocumentsDirectory]

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]

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

Zabudovane algoritmy I

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

d1=D[i1,x]

d1==v1

d1===v1

d1//Simplify

x=5;
sol3=Solve[eq3={a (x^2)+y==x,b x+c (y^2)==y},{x,y}]

x=.;y=.;
sol3=Solve[eq3={a (x^2)+y==x,b x+c (y^2)==y},{x,y}];
sol3//Length

Plus[a,b,c]
Attributes[Plus]

Zkouska, zda druhy koren skutecne splnuje dany system rovnic

eq3/. sol3[[2]]//Simplify

Clear[x,y];
mat2=Table[(x^i)+(y^j),{i,3},{j,4}]
Dimensions[mat2]
MatrixRank[mat2]
NullSpace[mat2]

Diferencialni rovnice symbolicky

eqdif=y'[x]+x y[x]==x Exp[-(x^2)]
DSolve[eqdif,y[x],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]

Numericke reseni diferencialnich rovnic

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

Plot[ChebyshevT[5,x],{x,-1,1}]

Plot[Table[ChebyshevT[j,x],{j,7}],{x,-1,1}]

Plot[Evaluate[Table[ChebyshevT[j,x],{j,7}]],{x,-1,1}]

Soucty a souciny

?n

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

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

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

Konec 12. 3. 2008

Lekce 5 (19. 3. 2008)

Jiny stylesheet

Odkazy

Vytvoreni hyperlink na licencni cislo

Obecne :

http://www.wolfram.com

http://www.wolfram.cz

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

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

Konecne je take vyhodne se podivat na nektere zdrojove programy ze Standard Add On Packages. Nektere jsou snadno srozumitelne, u nekterych je tomu naopak.

FAQ

Ruzne

Slozeny graf

Clear[f1];
f1[x_]:=Which[
	x<0,2x,
	0≤x<2,x,
	True, x/2
			]
Plot[f1[x],{x,-1,5},Exclusions->{x==2}]

Clear[p1,p2,p3];
p1=Plot[2x,{x,-1,0},PlotRange->All];
p2=Plot[x,{x,0,2}];
p3=Plot[x/2,{x,2,5}];
Show[p1,p2,p3]

Options[Plot]

a) Units revisited

b) Jeste dalsi moznost vymazani obsahu promennych (symbolu)

xxx=5;x3=4;y=5;
Clear["Global`x*"]
{xxx,x3,y}
xxx=5;x3=4;y=5;
Clear["Global`xx*"]
{xxx,x3,y}

c) Lokalizace promennych

xxx=5;x3=100;y=5;
With[{xxx=x3},Print["(xxx^2) = ",xxx xxx]]
xxx

xxx=5;x3=100;y=5;
Module[{xxx=1,x3,y},
x3=xxx+5;
y=x3 x3;
{xxx,x3,y}
	]
{xxx,x3,y}

Module[{y},y]

Module[{y},y]

d) Tvorba vlastniho rozdeleni: Pomerne prozaicke, viz verze 5.2

E:\\Dokumenty\\VPFPM\\VPFPM2007 2008\\ContinuousDistributions.m

PDF[GammaDistribution[α,p],x]

ExpectedValue[Sin,GammaDistribution[α,p]]
ExpectedValue[Sin[x],GammaDistribution[α,p],x]

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 (26. 3. 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"]

Cteni souboru

kursytxt=Import["denni_kurz.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

Legenda grafu (revidovano pro versi 6.0.2)

<<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`]},HeadCenter->0,HeadLength->0.03`],Text["Sin(x)",{3.5`,0.8`},{0,-1}],Arrow[{1,-0.75`},{(π/2),Cos[(π/2)]}],Text["Cos(x)",{1,-0.9`},{0,-1}]}]

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

Demonstrations Insurance, Finance

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 (2. 4. 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},
Needs["Histograms`"];
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]]

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 (9. 4. 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]

(* Initialize Histograms` & Module datareport first!!! *)
Needs["Histograms`"]
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,HistogramCategories->k,(* BarStyle->{GrayLevel[0.95`]}, *)ApproximateIntervals->True]]

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

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 (16. 4. 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 po284 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 (23. 4. 2008)

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

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 (30. 4. 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