Винеровский хаос или Еще один способ подбросить монетку08.12.2017 09:18

Теория вероятности никогда не переставала меня удивлять, начиная ещё с того момента, как я впервые с ней столкнулся, и до сих пор. В разное время в разной степени меня настигали, назовём их «вау-эффекты», шоковые удары в мозжечок, от которых меня накрывало эффектом третьего ока, и мир навсегда переставал быть прежним.

Первый «вау-эффект» я испытал от Центральной предельной теоремы. Берем кучу случайных величин, устремляем их количество в бесконечность и получаем нормальное распределение. И совсем неважно как распределены эти величины, неважно, будь это подбрасывания монетки или капли дождя на стекле, вспышки на Солнце или остатки кофейной гущи, результат будет всегда один — их сумма всегда стремится к нормальности. Разве что, нужно потребовать их независимость и существование дисперсии (позднее я узнал, что существует теорема и для экстремальных тяжелохвостых распределений с бесконечной дисперсией). Тогда этот парадокс долго не давал мне заснуть.
В какой-то момент учебы в университете такие предметы как дискретная математика и функциональный анализ слились вместе и всплыли в теорвере под видом выражения «почти наверное». Стандартный пример: вы случайно выбираете число от 0 до 1. С какой вероятностью вы ткнёте в рациональное число (привет, функция Дирихле)? Спойлер: 0. Ноль, Карл! Бесконечное множество не имеет никакой силы, если оно счетно. У вас бесконечное число вариантов, но вы не выберете ни один из них. Вы не выберете 0, или 1, или ½, или ¼. Вы и не выберете 3/2.
Да-да, что выбрать ½, что выбрать 3/2, вероятность нулевая. Вот только в 3/2 вы не ткнёте точно, таковы условия, а в ½ вы не попадёте ну… «почти наверное». Концепция «почти всюду»/«почти наверное» забавляет математика, а обывателя заставляет крутить пальцем у виска. Многие ломают себе мозг в попытке классифицировать нули, но результат того стоит.
Третий по счёту, но не по силе, «вау-эффект» настиг уже на переходе в advanced level
— при чтении книг по стохастическим исчислениям. Причиной тому стала лемма Ито. Со времён школьной скамьи, когда нашим девственным глазам впервые показали производную, мы нисколько не сомневались в правильности вот такой вот формулы:
$dX^2=2X\cdot dX.$
И она верна. Вот только, если — это не случайный процесс. Адовая смесь из свойств нормального распределения и «почти наверное» доказывает, что в обратном случае эта формула в общем случае неверна. Томик мат.анализа с решениями обыкновенных дифференциальных уравнений теперь можно выкинуть в топку. Люди в теме тихо хихикают, остальные нетерпеливо листают статьи в Вики с исчислениями Ито.

Но совсем недавно я испытал и четвёртый, так называемый, «вау-эффект». Это не какой-то отдельно взятый факт, а целая теория, которую я собираюсь поведать в серии из нескольких статей. И если предыдущие финты теории вероятности вас уже не удивляют, то прошу милости под кат (знаю, вы и так уже здесь).

Полиномы Эрмита

Начнем с обыкновенной алгебры — определим «вероятностные» (они немного отличаются от «физических») полиномы Эрмита:

$H_n(x) = (-1)^ne^{x^2/2} \frac{d^n}{dx^n}(e^{-x^2/2}), \quad n \in \mathbb{N}_0.$

Значения первых полиномов: $H_0(x)=1, H_1(x) = x, H_2(x) = x^2-1, H_3(x)=x^3-3x, \dots$

Полиномы Эрмита обладают следующими свойствами:

$H_n'(x) = nH_{n-1}(x),$
$H_{n}(x)=xH_{n-1}(x)-(n-1)H_{n-2}(x).$

Последнее соотношение поможет нам в вычислении $n$ -ых полиномов Эрмита для заданного $x$ . Программироваться мы будем на Haskell, ибо он позволяет математикам выражаться на привычном им языке — Haskell чист, строг и прекрасен как сама математика.

hermite :: (Enum a, Num a) => a -> [a]
hermite x = s
    where s@(_:ts) = 1 : x : zipWith3 (\hn2 hn1 n1 -> x * hn1 - n1 * hn2) s ts [1..]

Функция hermite принимает на вход параметр $x$ , а на выходе даёт бесконечный лист из $n$ -ых полиномов для $n=0,1,...$ Кто не знаком с концепцией ленивых вычислений, очень советую ознакомиться. Для тех же, кто эту концепцию знает, но ещё не до конца может в функциональное программирование: что здесь происходит? Представьте, что у нас уже есть бесконечный лист со всеми значениями Эрмитовых полиномов:

s = [1, x, x^2-1, x^3-3x, x^4-6x^2+3, ... ]

Хвост этого листа (без первого элемента):

ts = [x, x^2-1, x^3-3x, x^4-6x^2+3, ... ]

Вдогонку мы возьмем ещё лист с натуральными числами:

[1, 2, 3, ... ]

Функция zipWith3 комбинирует последние три листа, используя данный ему оператор:

x * [    x,  x^2-1,     x^3-3x, ... ]
  - [  1*1,    2*x,  3*(x^2-1), ... ]
  = [x^2-1, x^3-3x, x^4-6x^2+3, ... ]

Добавляем впереди 1 и x, и получаем полный набор Эрмитовых полиномов. Иными словами, мы достали лист со значениями полиномов, используя лист с этими значениями, то есть, лист, который мы и пытаемся достать. Поговаривают, что полное осознание красоты и мощи ФП сродни умению заглянуть себе в ухо.
Проверим: первые 6 значений для $x=1$ :

Prelude> take 6 (hermite 1)
[1,1,0,-2,-2,6]

Что мы и ожидали увидеть.

Гильбертово пространство

Двинемся немного в другую степь — вспомним определение пространства Гильберта. Говоря научным языком, это полное метрическое линейное пространство с заданным на нём скалярным произведением $\langle X, Y \rangle$ На этом пространстве каждому элементу соответствует вещественное число, именуемое нормой и равное

$\| X\|=\sqrt{\langle X,X\rangle}.$

Ничего сверхъестественного. Когда я пытаюсь представить Гильбертово пространство, я начинаю от простого и постепенно прихожу к сложному.

Самый простой пример — это пространство вещественных чисел: $H=\mathbb{R}$ . В таком случае скалярным произведением двух чисел и у нас будет
$\langle X,Y \rangle=XY.$
Затем я перехожу в Эвклидово пространство $H=\mathbb{R}^n$ . Теперь
$\langle X,Y\rangle=\sum_{i=0}^nX_iY_i.$
Это пространство можно расширить до пространства комплексных векторов: $H=\mathbb{C}^n$ , для которого скалярное произведение будет
$\langle X, Y\rangle = \sum_{i=0}^nX_i\overline{Y_i}$
(верхняя черта обозначает комплексное сопряжение).
Ну и наконец прихожу в пространство для взрослых, пространство с бесконечной размерностью. В нашем случае это будет пространство квадратично интегрируемых функций, заданных на некотором множестве $\Omega$ с заданной мерой $\mu$ . Мы будем его обозначать в виде $H=L^2(\Omega,\mu)$ . Скалярное произведение на нем задается следующим образом:
$\langle X, Y \rangle =\int_\Omega (X \cdot Y) d\mu.$
Обычно под множеством $\Omega$ подразумевается интервал , а под мерой $\mu$ — равномерная мера (мера Лебега), т.е. $d\mu=\mu(d\omega)=d\omega$ . И тогда скалярное произведение записывается в виде обыкновенного интеграла Лебега
$\int_a^bX(\omega)Y(\omega) d\omega.$
Если же мы думаем в терминах теории вероятностей, то $\Omega$ — это пространство элементарных событий, $X=X(\omega)$ и $Y=Y(\omega)$ — случайные величины, а $\mu$ — вероятностная мера. У каждой такой меры есть своя функция плотности распределения $\rho$ , которая может быть отличной от константы, тогда $d\mu=\rho(\omega)d\omega$ и скалярное произведение совпадает с математическим ожиданием:
$\langle X, Y \rangle=\int_\Omega X(\omega) Y(\omega) \rho(\omega)d\omega = \mathbb{E}[XY].$

Гауссовский процесс

Настало время внести в наши размышления элемент случайности. Пусть у нас имеется гильбертово пространство $H$ . Тогда мы назовем $\{W(h)\}_{h\in H}$ (изонормальным) Гауссовским процессом, если

вектор из случайных величин $(W(h_1), \dots, W(h_n))$ распределен нормально с нулевым мат.ожиданием для любых $h_1, \dots h_n \in H$ , и
для $h,g \in H$
$\mathbb{E}[W(h)\cdot W(g)] = \langle h, g \rangle.$

По своей математической сути $W(h)$ — это отображение из одного гильбертова пространства в другое, из некоторого $H$ в $L^2(\Omega, \mathcal{F}, \mathbb{P})$ — вероятностное пространство из случайных величин с конечной дисперсией, заданного триплетом $\Omega$ (множество элементарных событий), $\mathcal{F}$ (сигма-алгебра) и $\mathbb{P}$ (вероятностная мера). Несложно показать, что это отображение линейно:

$W(ah+bg)=aW(h)+bW(g) \quad \forall a,b \in \mathbb{R}, \forall h,g \in H.$

(в смысле равенства «почти наверное», привет «вау-эффект» #2)

Пример. Пусть $H=L^2((0,\infty), \lambda)$ , где $\lambda$ — равномерная мера (Лебега). Скалярное произведение на нём

$\langle f, g\rangle=\int (f \cdot g) d\lambda.$

Пусть $h(t) = 1_{[0,t]}$ — единичная функция на интервале $[0,t]$ . Тогда $\| h\|^2 = \int 1_{[0,t]}ds = t$ и

$B(t)=W(1_{[0,t]}) \sim \mathcal{N}(0,t)$

не что иное, как Броуновское движение (или Винеровский процесс). Более того,

$\int_0^t f(s)dB(s) = W(1_{[0,t]} f)$

называется интегралом Ито от функции $f$ относительно $B$ .

Для того, чтобы реализовать Гауссовский процесс я воспользуюсь пакетами, которые благородные люди уже написали за нас.

import Data.Random.Distribution.Normal
import Numeric.LinearAlgebra.HMatrix as H

gaussianProcess :: Seed -> Int -> Int -> ((Int, Int) -> Double) -> Matrix Double
gaussianProcess seed nSamples dim dotProducts = gaussianSample seed nSamples mean cov_matrix
    where mean = vector (replicate dim 0)
          cov_matrix = H.sym $ matrix dim (map (\i -> dotProducts (i `quot` dim, i `rem` dim)) $ take (dim * dim) [0, 1..])

Функция gaussianProcess принимает параметр seed (стандартная штука для генераторов), nSamples — размер выборки, dim — размерность вектора $(h_1,…,h_n)^T$ , dotProducts — функцию, принимающую на вход $(i,j)$ , индекс матрицы ковариации и возвращающую соответствующее этому индексу скалярное произведение $\langle h_i,h_j\rangle$ . На выход gaussianProcess выдает вектор $(W(h_1), \dots, W(h_n))$ .

Уже подходит время объединить все полученные нами знания вместе. Но прежде, стоит упомянуть об одном полезном свойстве эрмитовых полиномов и нормального распределения в совокупности. Пусть $F(t,x)=\exp(tx-t^2/2).$ Тогда, используя разложение Тейлора,

$\begin{aligned} F(t,x)& = \exp(x^2/2-(x-t)^2/2) \\ &=\exp(x^2/2) \sum_{n=0}^\infty \frac{t^n}{n!t}\frac{d^n}{dt^n}\exp(-(x-t)^2/2)\bigg |_{t=0}\\ &= \sum_{n=0}^\infty \frac{(-1)^n}{n!} \exp(x^2/2) \frac{d^n}{dz^n} \exp(-z^2/2)\bigg |_{z=x}\\ &=\sum_{n=0}^\infty \frac{t^n}{n!} H_n(x). \end{aligned}$

Возьмем $X, Y \sim \mathcal{N}(0,1)$ — две стандартные нормально распределенные случайные величины. Через производящую функцию нормального распределения мы можем вытащить следующее соотношение:

$\mathbb{E}[F(s,X)\cdot F(t,Y)]=\exp(st\mathbb{E}[XY]).$

Берем $(n+m)$ -ую частную производную $\frac{\partial^{n+m}}{\partial s^n \partial t^m}$ , приравниваем $s=t=0$ по обе стороны уравнения сверху и получаем

$\mathbb{E}[H_n(X) \cdot H_m(Y)] = \begin{cases} n!(\mathbb{E}[XY])^n, & n=m, \\ 0, & n \neq m. \end{cases}$

О чем это нам говорит? Во-первых, мы получили норму $\| H_n(X) \|^2=n!$ для $X \sim \mathcal{N}(0,1)$ , а, во-вторых, мы теперь знаем, что разные эрмитовы полиномы от нормальных случайных величин ортогональны друг другу. Вот сейчас мы готовы к осознанию нечто большего.

Разложим пространство в хаос

Пусть $\mathcal{H}_n=\overline{\operatorname{span}}\Big \{ H_n(W(h)) \Big| \| h\| = 1 \Big\}$ — n-й Винеровский хаос. Тогда

$L^2(\Omega, \mathcal{F}, \mathbb{P}) = \bigoplus_{n=0}^\infty \mathcal{H}_n.$

Воу-воу, палехче! Давайте разложим эту теорему о разложении по кусочкам и переведём с математического на человеческий. Мы не будем сильно вдаваться в детали, а лишь интуитивно поясним о чем тут речь. Значок $\operatorname{span}(X)$ обозначает линейную оболочку подмножества $X$ гильбертова пространства $H$ — пересечение всех подпространств $H$ , содержащих $X$ . Говоря проще, это множество всех линейных комбинаций элементов из $X$ . Черта сверху над $\operatorname{span}$ обозначает замыкание множества. Если $\overline{\operatorname{span}}(X) = H$ , то $X$ называется полным множеством (грубо говоря,» $X$ плотно в $H$ »). Следовательно, $\overline{\operatorname{span}} \{ H_n(W(h)) | \| h\| = 1 \}$ — замыкание линейной оболочки полиномов Эрмита от Гауссовского процесса на единичной гиперсфере.

С нотацией вроде разобрались. Теперь о том, что такое Винеровский хаос. Идем от простого: $\mathcal{H}_0$ содержит все линейные комбинации Эрмитовых полиномов со степенью 0, то есть различные комбинации чисел $a\cdot 1$ , то есть всё пространство вещественных чисел. Следовательно, $\mathcal{H}_0=\mathbb{R}$ . Идем дальше. Несложно увидеть, что $\mathcal{H}_1= \{ W(h) | h \in H \}$ , то есть пространство, составленное из Гауссовских процессов. Получается, что все центрированные нормальные величины принадлежат $\mathcal{H}_1$ . Если мы добавим еще $\mathcal{H}_0$ , то к ним присоединятся и остальные нормальные случайные величины, чье математическое ожидание отлично от нуля. Дальнейшие множества $\mathcal{H}_n$ уже оперируют с n-ми степенями $W(h)$ .

Пример. Пусть $H=L^2((0,\infty), \lambda)$ и $X=B(t)^2$ — квадрат Броуновского движения. Тогда

$\begin{aligned} B(t)^2&=W(1_{[0,t]})^2 \\ & = \|1_{[0,t]} \|^2 \cdot W\bigg(\frac{1_{[0,t]}}{\|1_{[0,t]} \|}\bigg)^2 \\ & = t\cdot W\bigg(\frac{1_{[0,t]}}{\sqrt{t}}\bigg)^2\\ &=tH_2\bigg(W\bigg(\frac{1_{[0,t]}}{\sqrt{t}}\bigg)\bigg)+t. \end{aligned}$

Первое слагаемое принадлежит $\mathcal{H}_2$ , второе — $\mathcal{H}_0$ . Это и называется разложением в Винеровский хаос.

Мы показали ранее, что $\mathcal{H}_n \perp \mathcal{H}_m$ для $n \neq m$ . Теорема о разложении гласит о том, что эти множества не только ортогональны друг другу, но также формируют полную систему в $L^2(\Omega, \mathcal{F}, \mathbb{P})$ . Что это означает на практике? Это значит, что любая случайная величина $X$ с конечной дисперсией может быть аппроксимирована полиномиальной функцией от нормально распределенной случайной величины.

На самом деле

На самом деле, такое разложение полезно, если распределение $X$

в определенном смысле близко к распределению нормальному. Например, если мы имеем дело с Броуновским движением или с логнормальным распределением. В действительности, плотность нормального распределения

$\rho(x)=\frac{1}{\sqrt{2\pi}}e^{-x^2/2}$

очень схожа с определением полинома Эрмита

$H_n(x) = (-1)^ne^{x^2/2} \frac{d^n}{dx^n}(e^{-x^2/2}), \quad n \in \mathbb{N}_0.$

Если же распределение $X$

далеко от Гаусса, то можно попробовать и другие ортогональные полиномы. Например, плотность Гамма-распределения:

$\rho(x)=\frac{x^{n-1} e^{-x}}{\Gamma(n)}.$

Ничего не напоминает? Да это же полиномы Лагерра

$L_n(x) = \frac{e^x}{n!}\frac{d^n}{dx^n}(x^ne^{-x})$

Равномерному распределению соответствуют полиномы Лежандра, биномиальному распределению — полиномы Кравчука, и т.п. Теория, развивающая идею разложения вероятностного пространства на ортогональные полиномы именуется в англоязычной литературе как «Polynomial chaos expansion».

Пример. Давайте теперь возьмем $H=\mathbb{R}$ , функцию $f$ и зададим случайную величину $X$ , такую что

$X=f(\xi) \in L^2(\Omega, \mathcal{F}, \mathbb{P}),$

где $\xi = W(1) \sim \mathcal{N}(0,1)$ . По теореме о разложении мы можем представить её в виде взвешенной суммы из полиномов Эрмита

$f(\xi)=\sum_{n=0}^\infty f_nH_n(\xi),$

где коэффициенты задаются формулой

$f_n=\frac{1}{n!}\mathbb{E}[f(\xi) \cdot H_n(\xi)].$

Эти значения $f_n$ мы получили следующим образом:

$\begin{aligned} \mathbb{E}[f(\xi)\cdot H_n(\xi)] & = \langle f(\xi), H_n(\xi) \rangle \\ &= \langle \sum_{k=0}^\infty f_k H_k(\xi) , H_n(\xi) \rangle \\ &= \sum_{k=0}^\infty f_k \langle H_k(\xi), H_n(\xi) \rangle \\ & = f_n \| H_n(\xi) \|^2 = f_n n!. \end{aligned}$

Поздравляю! Теперь, если у вас есть функция от стандартной нормально распределенной случайной величины, вы сможете её разложить по базису из Эрмитовых полиномов. Например, подбрасывание честной монетки 0–1 мы можем представить в виде

$X=1_{(0,\infty)}(\xi).$

Немного поколдовав с математикой (подсчет несложных интегралов мы оставим читателю), мы получаем разложение:

$X= \frac{1}{2} + \frac{1}{\sqrt{2\pi}} \sum_{n=0}^\infty \frac{(-1)^n}{2^n(2n+1)n!} H_{2n+1}(\xi).$

Заметьте, что каждый второй элемент в разложении по базису равен нулю.

second (x:y:xs) = y : second xs
second _ = []

coinTossExpansion :: Int -> Double -> Double
coinTossExpansion n xi = sum (take n $ 0.5 : zipWith (*) fn (second $ hermite xi))
    where fn = 1.0 / (sqrt $ 2 * pi) : zipWith ( \fn1 k -> -fn1 * k / ((k + 1) * (k + 2)) ) fn [1, 3..]

Функция coinTossExpansion возвращает сумму, полученную разложением случайной монетки в винеровский хаос, для данного $\xi$ от $0$ до $n$ . На графике показана постепенная сходимость для выбранных случайным образом $\xi$ с возрастанием $n$ .

Судя по этому графику, где-то после $n \approx 100$ мы можем обрезать сумму, округлить и вернуть в качестве $X$ .

coinTossSequence :: Seed -> Int -> [Int]
coinTossSequence seed n = map (round.coinTossExpansion 100) (toList nvec)
    where nvec = (toColumns $ gaussianProcess seed n 1 (\(i,j)->1) ) !! 0

Проверим, как будет выглядеть последовательность из 20 подбрасываний.

Prelude> coinTossSequence 42 20
[0,0,1,0,0,0,1,1,0,1,0,0,0,1,0,1,1,1,0,1]

Теперь, когда вас попросят сгенерировать подбрасывания монетки, вы знаете, что им показать.

Ну, а без шуток, мы что-то посчитали и что-то разложили, а какой в этом всем толк, спросите вы. Не спешите чувствовать себя обманутыми. В последующих статьях мы покажем, как это разложение позволяет брать производную от случайной величины (в некотором смысле), расширим стохастическое интегрирование (и ваше сознание), и найдем всему этому практическое применение в машинном обучении.