Metoda Monte Carlo ⬩ Adamátorovy zápisky

Pojištění

Pojďme si spočítat nějaký příklad, který se vyskytuje ve finanční matematice:

S_{α} ≔ \sum_{n = 1}^{\infty} \frac{1}{n^{1 + α}} .

Vyjádříme si:

S_{α} = \sum_{n = 1}^{\infty} \frac{1}{n \cdot (n + 1)} \cdot \frac{n + 1}{n^{α}} .

Z analýzy víme, že

\sum_{n = 1}^{\infty} \frac{1}{n \cdot (n + 1)} = 1

, takže si můžeme zavést náhodnou proměnnou

ξ_{α}

s diskrétním rozdělením

P [ξ_{α} = \frac{n + 1}{n^{α}}] ≔ \frac{1}{n \cdot (n + 1)} .

Potom bude platit

S_{α} = E ξ_{α}

. Můžeme tedy použít aproximaci

{\bar{ξ}}_{α} ≔ \frac{1}{N} \cdot \sum_{i = 1}^{N} ξ_{i},

σ^{2} = \frac{1}{N \cdot (N - 1)} \sum_{i = 1}^{N} {| {\bar{ξ}}_{α} - ξ_{i} |}^{2} .

Obecně

Mějme náhodnou veličinu

ξ

se střední hodnotou

E ξ

a rozptylem

D ξ = σ^{2} < \infty

. Tu aproximujeme jako

\bar{ξ} ≔ \frac{1}{N} \sum_{i = 1}^{N} ξ_{i},

E \bar{ξ} = μ,

D \bar{ξ} = \frac{σ^{2}}{N^{2}} .

Platí odhad:

P (| \bar{ξ} - μ | \geq ε) \leq \frac{σ^{2}}{N^{2} \cdot ε^{2}},

který si můžeme upravit jako

P (| \bar{ξ} - μ | < \frac{σ}{N \cdot \sqrt{η}}) \geq 1 - η .

Převod rovnoměrného rozdělení na diskrétní

Mějme diskrétní rozdělení, kde

P [x = x_{i}] = p_{i}

. Máme-li k dispozici názhodnou proměnnou

r \sim u (0, 1)

a chceme pomocí ní náhodně vybrat

x

, stačí najít

i

takové, že

\sum_{j = 1}^{i - 1} p_{j} \leq r < \sum_{j = 1}^{i} p_{j} .

Speciálně pro Poissonovo rozdělení

Po (λ)

existuje lepší algoritmus. Budeme postupně generovat náhodné proměnné

γ_{i} \sim R (0, 1)

a skončíme, jakmile

\sum_{i = 0}^{n} γ_{i} \geq exp (- λ)

. Potom

n \sim Po (λ)

Více exponenciálně rozdělených veličin

Chceme vygenerovat exponenciálně rozdělenou veličinu:

F_{X} (x) = 1 - exp (λ x), x \geq 0

Můžeme ji vygenerovat obecným způsobem, tedy vezmeme

γ \sim R (0, 1)

a spočteme si

x = - \frac{ln (1 - γ)}{λ}

.K tomu potřebujeme jeden výpočet logaritmu.Ovšem existuje způsob, jak vygenerovat několik nezávislých exponenciálně rozdělených veličin s jen jedním logaritmem. Mějme

γ_{1}, \dots, γ_{2 n - 1} \sim R (0, 1)

. Vezměme

γ_{n + 1}, \dots, γ_{2 n - 1}

a seřaďme je:

α_{0} = 0 \leq α_{1} = γ_{něco} \leq \dots \leq α_{n - 1} = γ_{něco} \leq α_{n} = 1

. Potom pro každé

i \in \hat{n}

platí (bez důkazu)

(α_{i} - α_{i - 1}) \cdot ln (γ_{1} \dots γ_{n}) \sim exp (- x) .

Gamma rozdělení

Chceme mít veličinu s rozdělením

ξ_{α} \sim f_{α (x)} = \frac{x^{α - 1} e^{- x}}{Γ (α)} .

Toto rozdělení má zajímavou vlastnost, že

ξ_{α} + ξ_{β} \sim ξ_{α + β}

\begin{align} (f_{α} * f_{β}) (x) & = \frac{1}{Γ (α) \cdot Γ (β)} \int_{0}^{x} y^{α - 1} e^{- y} {(x - y)}^{β - 1} e^{- (x - y)} d y \\ = \frac{e^{- x}}{Γ (α) \cdot Γ (β)} \int_{0}^{x} y^{α - 1} {(x - y)}^{β - 1} d y \\ = | \begin{matrix} y = t x \\ d y = x d t \end{matrix} | = \frac{e^{- x} \cdot x^{α + β - 2}}{Γ (α) \cdot Γ (β)} \int_{0}^{1} t^{α - 1} {(1 - t)}^{β - 1} . \end{align}

Normální rozdělení

Z centrální limitní věty víme, že pokud budeme pro veličiny

X_{i} \sim R (0, 1)

sčítat

\frac{X - E X}{\sqrt{D X}}

, dostaneme normální rozdělení

N (0, 1)

. Stačí si tedy spočíst

E X_{i} = \frac{1}{2}

D X_{i} = \frac{1}{12}

. Ale tím samozřejmě dostaneme normální rozdělení jenom přibližně. K získání přesného normálního rozdělení můžeme použít legrační tvrzení, které jsme si odvodili na cvičení. Mějme náhodné veličiny

γ_{1}, γ_{2} \sim R (0, 1)

, potom z nich můžeme získat dvě nezávislé normálně rozdělené veličiny

η_{1} ≔ \sqrt{- 2 \cdot ln γ_{1}} \cdot cos 2 π γ_{2},

η_{2} ≔ \sqrt{- 2 \cdot ln γ_{1}} \cdot sin 2 π γ_{2} .

$χ$ -rozdělení

f_{χ^{2}} (x) = \frac{2^{- \frac{n}{2}}}{Γ (\frac{n}{2})} x^{\frac{n}{2} - 1} exp (- \frac{x}{2})

f_{χ} (x) = \frac{2^{- \frac{n}{2} + 1}}{Γ (\frac{n}{2})} x^{n - 1} exp (- \frac{x^{2}}{2})

Vícerozměrné náhodné veličiny

Jak generovat náhodnou veličinu s vícerozměrným rozdělením

F (x_{1}, \dots, x_{n})

? Pokud dokážeme spočítat všechny potřebné integrály, můžeme prostě generovat složky postupně s tím, že je vždy podmíníme již vygenerovanými složkami. Ale to většinou moc nefunguje. Další možnost je zamítací metoda, ale ta bude většinou strašně pomalá, protože je malá šance, že se trefíme. Například u dvacetirozměrné koule je to

2.5 \times 10^{- 8}

. Ale v určitých případech to jde lépe.

Kovarianční matice

Máme danou kovarianční matici

𝐂

a chceme vygenerovat vektor s rozdělením

η \sim N (μ, 𝐂)

, to znamená

𝐂 = E ((η - E η) {(η - E η)}^{𝖳}) .

Pro jednoduchost uvažejme

μ = 0

.Vynásobíme-li to nějakou maticí

A

, máme

𝐂 (A η) = E (A η {(A η)}^{𝖳}) = E (A η η^{𝖳} A^{𝖳}) = A E (η η^{𝖳}) A^{𝖳} = A A^{𝖳} = 𝐂 .

Výpočet integrálu metodou Monte Carlo

Mějme nezápornou omezenou integrabilní funkci

f : ⟨ a, b ⟩ \to ⟨ 0, c ⟩

. Chceme najít

I ≔ \int_{a}^{b} f

. První způsob by byl prostě střílet do obdélníku

⟨ a, b ⟩ \times ⟨ 0, c ⟩

a počítat, kolikrát se trefíme. Formálně bychom měli náhodnou proměnnou

χ ≔ [γ_{2} \leq γ_{1}], (γ_{1}, γ_{2}) \in R (⟨ a, b ⟩ \times ⟨ 0, c ⟩) .

Porom

E χ = \frac{I}{c \cdot (b - a)}

. Náš odhad bude vypadat jako

θ_{1} ≔ \frac{c \cdot (b - a)}{N} \sum_{i = 1}^{N} χ_{i} .

D θ_{1} = \frac{{(c \cdot (b - a))}^{2}}{N} \cdot D χ = \frac{{(c \cdot (b - a))}^{2}}{N} \cdot (\frac{I}{c \cdot (b - a)} - \frac{I^{2}}{c^{2} \cdot {(b - a)}^{2}}) = \frac{I}{N} (c \cdot (b - a) - I) .

Můžeme na to ovšem jít míň hloupě. Stačí si říct, že kdybychom měli náhodnou proměnnou

γ \sim R (a, b)

, potom

E f (γ) = I

. Náš odhad bude

θ_{2} ≔ \frac{b - a}{N} \sum_{i = 1}^{N} f (γ_{i}), γ_{i} \sim R (a, b) .

Potom

E θ_{2} = (b - a) E f (γ_{i}) = I,

D θ_{2} = \dots = \frac{b - a}{N} \cdot \int_{a}^{b} f^{2},

\frac{N}{b - a} \cdot (D θ_{1} - D θ_{2}) = c \cdot I - \int_{a}^{b} f^{2} = \int_{a}^{b} (c - f (x)) \cdot f (x) d x \geq 0 .

Takže za našich předpokladů je lepší druhá metoda. Navíc potřebujeme dvakrát méně náhodných veličin. Mohli bychom to ještě nějak vylepšit? Máme-li odhad

g (x) \leq f (x) \leq g (x) + d

, kde funkci

g

umíme přesně zintegrovat, můžeme použít metodu vyloučení hlavní části: metodou Monte Carlo si spočteme

\int_{a}^{b} (f - g)

a přičteme

\int_{a}^{b} g

. Také pokud víme, že nějaká část funkce bude výrazně chaotičtější než jiná, můžeme si vzít nějakou distribuční funkci

p

, náhodnou proměnnou

ξ \sim p

a počítat

E f (ξ_{i}) = \int_{a}^{b} p (x) f (x) d x,

θ_{3} = \frac{1}{N} \sum_{i = 1}^{n} \frac{f (ξ_{i})}{p (ξ_{i})} \cdot ξ_{i},

E θ_{3} = \dots = I,

D θ_{3} = \frac{1}{N} D \frac{f (ξ)}{p (ξ)} = \frac{1}{N} (E {(\frac{f (ξ)}{p (ξ)})}^{2} - I^{2}) = \frac{1}{N} (\int_{a}^{b} \frac{f^{2} (x)}{p (x)} d x - I^{2})

Virius zkoušel různými metodami spočítat

\int_{0}^{1} exp = e - 1

. Nejlepší možnost byla použít odhad

1 + x \leq exp x \leq e - 1 + x

a hustotu

p (x) = 1 + x

Vícerozměrné integrály

Chceme najít integrál

I ≔ \int_{M} f

, kde

M \subset ℝ^{d}

. Je-li

M

d

-rozměrný kvádr, můžeme to dělat stejně jako u jednorozměrného, ale pokud to bude nějaká hnusná množina, máme problém. Pro malé

d

můžeme použít zamítací metodu k vytvoření rovnoměrného rozdělení na

M

, ale jak už víme, pro velké

d

to funguje blbě, protože může být malá šance, že se trefíme. Často se dá využít nějaká symetrie.

Neomezená funkce

Co když máme funkci, která sice je integrabilní, ale má singularitu? Kdybychom to zkoušeli spočíst vzorkováním funkce, mohli bychom se náhodou trefit hodně blízko, což by nám naprosto rozhodilo výpočet. Místo toho můžeme zkusit na nějakém malém okolí singularity integrál odhadnout a počítat ho jenom na zbytku. Další možnost je vhodně zvolit

p

nebo provést nějakou substituci, aby singularita zmizela.

Integrály na neomezené množině

Co kdybychom chtěli spočíst například

I ≔ \int_{- \infty}^{\infty} e^{- x^{2}} d x

? Použijeme hustotu

ξ \sim p (x) ≔ \frac{1}{π \cdot (1 + x^{2})}

. Potom budeme mít

θ = \frac{π}{N} \cdot \sum_{i = 1}^{N} e^{- ξ_{i}^{2}} (1 + ξ_{i}^{2}), ξ \sim \frac{1}{π \cdot (1 + x^{2})},

D θ = \frac{π}{N} \int_{- \infty}^{\infty} e^{- 2 x^{2}} (1 + x^{2}) d x < \infty .

Takže tato metoda je použitelná, přestože integrujeme na nekonečné množině. Kdybychom ale tímhle způsobem zkusili spočíst

I ≔ \int_{- \infty}^{\infty} \frac{1}{1 + x^{2}}

s náhodnou proměnnou

ξ \sim N (0, 1)

, vyšlo by nám

D θ^{2} = \infty

Integrál přes kouli

Dá se ukázet, že pokud chceme mít rovnoměrné rozdělení v kouli, tak si stačí v nějakých upravených kulových souřadnicích vzít

φ \sim R (0, 2 π), μ \sim R (- 1, 1), r \sim ???

. Pokud bychom ale i

r

volili rovnoměrně, dostaneme integrál z

\frac{1}{r^{2}}

, což se může někdy hodit. Značí-li

U_{1}

jednotkovou kouli, zkusme například zintegrovat

\int_{U_{1}} \int_{U_{2}} f (| P - Q |) d P d Q .

Označíme-li

ρ ≔ | P - Q |

, pro

ρ \to 0

má

f (ρ) \cdot ρ^{2}

konečnou hustotu. Vzhledem k symetrii můžeme pro kulové souřadnice volit

φ_{P} = ϑ_{P} = 0, r_{P} \sim 3 r^{2}

. (Toto

r_{P}

můžeme vygenerovat tak, že zvolíme tři rovnoměrně náhodná čísla a vezmeme největší z nich.) Stejně tak můžeme volit

φ_{Q} = 0

. Zbylé dvě proměnné zvolíme jako

ϑ_{q} \sim R (- 1, 1), r_{Q} \sim R (0, l)

, kde

l ≔ μ_{Q} r_{P} + \sqrt{1 - r_{P} \cdot (1 - μ_{Q}^{2})}

. Potom

1 = \int_{U} \frac{1}{2 π} \frac{1}{2} \frac{1}{l} \frac{1}{ρ^{2}} ρ^{2} d φ d r d μ .

Zpět k pojištění

Na začátku jsme počítali sumáž

S_{α} ≔ \sum_{i = 1}^{\infty} \frac{1}{n^{1 + α}} .

To můžeme také počítat pomocí integrálu:

\int_{0}^{\infty} \frac{1}{{⌊ 1 + x ⌋}^{1 + α}} d x = 1 + \int_{1}^{\infty} \frac{1}{{⌊ 1 + x ⌋}^{1 + α}} d x .

To shora odhadneme

\frac{1}{x^{1 + α}}

a spočteme pomocí vzorkování plochy nebo vzorkování funkce.

Markovské procesy

Mějme proces s konečně mnoha stavy, kde stav, do kterého se dostaneme v dalším kroku, závisí jen na stavu, ve kterém jsem teď. Takový proces můžeme vyjádřit pomocí matice přechodu:

v_{k + 1} = (\begin{matrix} p_{1, 1} & \dots & p_{1, n} \\ ⋮ & ⋱ & ⋮ \\ p_{n, 1} & \dots & p_{n, n} \end{matrix}) \cdot v_{k} = P \cdot v_{k},

kde

{(v_{k})}_{i}

značí pravděpodobnost, že v

k

-tém kroku jsme v

i

-tém stavu.Stav

i

je absorpční, pokud

p_{i, j} = δ_{i, j}

, tedy se z něj nejde dostat do jiného stavu. Stavy jsou sdružené, pokud existuje cesta z jednoho do druhého. Řetězec je zanikající, pokud každý stav je sdružený s nějakým absorpčním stavem. Řetězec je ergodický, pokud všechny stavy jsou sdružené.

Lineární rovnice

Mějme lineární rovnici

A x = b

. Tu můžeme přepsat do tvaru

x = B x + b

, kde

B ≔ I - A

. Je-li

ρ (B) < 1

, můžeme ji řešit metodou postupných aproximací:

x = \sum_{i = 0}^{\infty} B^{i} \cdot f,

x_{m} \approx b_{m} + \sum_{i = 1}^{n} B_{m, i} b_{i} + \sum_{i = 1}^{n} \sum_{j = 1}^{n} B_{m, i} B_{i, j} b_{j} + \dots

Je-li matice stochastická (všechny řádky mají součet

1

), můžeme použít metodu Monte Carlo. Chceme-li aproximovat například druhý člen, vygenerujeme náhodnou veličinu podle

m

-tého řádku, čímž dostaneme

i

. Poté vygenerujeme náhodnou veličinu podle

i

-tého řádku, čímž dostaneme

j

. Do výsledného součtu přihodíme

b_{j}

. Co když ale matice není stochastická? A jak se vypořádat s tím, že máme nekonečně mnoho sčítanců? Rozepíšeme si

B ≕ P \otimes F, b = p \otimes f

, kde

\forall i : \sum_{j = 1}^{n} P_{i, j} + p_{j} = 1

. Potom máme

x_{m} = p_{m} f_{m} + \sum_{i = 1}^{n} P_{m, i} p_{i} F_{m, i} f_{i} + \sum_{i = 1}^{n} \sum_{j = 1}^{n} P_{m, i} P_{i, j} p_{j} F_{m, i} F_{i, j} f_{j} .

Budeme to modelovat tak, že si přidáme jeden absorpční stav a vezmeme přechodovou matici

(\begin{matrix} P_{1, 1} & \dots & P_{1, n} & p_{1} \\ ⋮ & ⋱ & ⋮ & ⋮ \\ P_{n, 1} & \dots & P_{n, n} & p_{n} \\ 0 & \dots & 0 & 1 \end{matrix}) .

Budeme simulovat proces s touto přechodovou maticí. Začneme s jednotkovým vektorem

e_{m}

. Pokud se po

k + 1

krocích dostaneme do absorpčního stavu a trajektorie předtím je

i_{1}, \dots, i_{k}

, za hodnotu náhodné veličiny vezmeme

η ≔ (\prod_{j = 1}^{k} F_{i_{j - 1}, i_{j}}) \cdot f_{i_{k}}

. Střední hodnota potom bude hledané

x

. Ale jaká bude disperze?

D η = E η^{2} - {(E η)}^{2},

E η^{2} = p_{m} f_{m} + \sum_{i = 1}^{n} P_{m, i} p_{i} F_{m, i}^{2} f_{i}^{2} + \dots .

To znamená, že aby disperze byla konečná, musí konvergovat metoda postupných aproximací pro rovnici

x = C x + d

, kde

C ≔ F \otimes F \otimes P, d ≔ f \otimes f \otimes p

. Takže musíme

P, p

zvolit šikovně, na což neexistuje žádná obecná metoda.

Laplaceova rovnice

Mějme rovnici

\nabla^{2} u = 0

na oblasti

G

neboli

\partial_{x}^{2} u + \partial_{y}^{2} u = 0, u |_{\partial G} = f

.Při použití diferenční náhrady dostaneme

u_{i, j} = \frac{1}{4} (u_{i + 1, j} + u_{i - 1, j} + u_{i, j + 1} + u_{i, j - 1}) .

Když to interpretujeme jako střední hodnotu, tak to můžeme vzít jako náhodný proces, kde začneme v bodě

(i, j)

a budeme se náhodně pohybovat, až narazíme na hranici, a vezmeme hodnotu

f

v bodě, kam narazíme. Snadno vidíme, že k získání konečné disperze stačí, aby

f

byla omezená.Mimochodem, funkce

u

splňující vlastnost

\nabla^{2} u = 0

se nazývá harmonická. Pro ty existuje věta o průměrné hodnotě, která říká, že hodnotu v nějakém bodě můžeme spočíst jako průměr něčeho kolem něj. To je vlastně spojitá verze toho, co jsme si teď odvodili.

Co kdybychom měli rovnici s pravou stranou, $\nabla^{2} u = - g$ ? Její diskrétní verze je

u_{i, j} = \frac{1}{4} (u_{i + 1, j} + u_{i - 1, j} + u_{i, j + 1} + u_{i, j - 1}) + \frac{h^{2}}{4} \cdot g_{i, j} .

No to je jednoduché – stačí v každém kroku ještě připočíst

\frac{h^{2}}{4} \cdot g_{i, j}

. Jaká je střední hodnota počtu kroků, které provedeme? Dostáváme rekurentní rovnici

K_{i, j} = \frac{1}{4} (K_{i + 1, j} + K_{i - 1, j} + K_{i, j + 1} + K_{i, j - 1}) + 1 .

To jde vyřešit pomocí nějaké analýzy vlastních čísel, ale do toho nebudeme zabíhat.

Něčí rovnice

φ (x) = \int_{ℝ^{n}} k (y, x) \cdot φ (y) d y + f (x) .

Postupné aproximace:

φ (x) = \int_{ℝ^{n}} k (y, x) \cdot f (y) d y + f (x),

φ_{2} (x) = \int_{ℝ^{n}} \int_{ℝ^{n}} k (y, x_{1}) \cdot k (x_{1}, x) \cdot f (y) d y d x_{1} + φ_{1} (x) .

Nějaká magie s Lebesgueovými prostory, jádrem, sdruženými prostory a dalšími funkcionálními šílenostmi.

φ = 𝔎 φ + f

‖ 𝔎 ‖ = sup_{f \in ℒ_{+}} \frac{‖ 𝔎 f ‖}{‖ f ‖} \leq sup_{x \in ℝ^{n}} \int_{ℝ^{n}} | k (y, x) | d y

Testování prvočíselnosti

Je-li

n = d \cdot 2^{s} + 1

prvočíslo, potom

\forall a < n : a^{d} \equiv 1 (mod n)

. Pokud je složené, tak to pro nějaká

a

platí a pro nějaká ne. Těm, pro která to platí, říkáme lháři. Ukazuje se, že pravděpodobnost lháře je přibližně

\frac{1}{4}

. Takže pokud to náhodně zkusíme pro

k

různých hodnot

a

, máme šanci

{(\frac{1}{4})}^{k}

na falešné pozitivum. To je Millerův-Rabinův test, podobný Fermatovu testu prvočíselnosti, viz DIM1.

Metoda Monte Carlo

Pojištění

Obecně

Převod rovnoměrného rozdělení na diskrétní

Více exponenciálně rozdělených veličin

Gamma rozdělení

Normální rozdělení

χ-rozdělení

Vícerozměrné náhodné veličiny

Kovarianční matice

Výpočet integrálu metodou Monte Carlo

Vícerozměrné integrály

Neomezená funkce

Integrály na neomezené množině

Integrál přes kouli

Zpět k pojištění

Markovské procesy

Lineární rovnice

Laplaceova rovnice

Něčí rovnice

Testování prvočíselnosti

$χ$ -rozdělení