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

в 4:10, , рубрики: haskell, вероятностное программирование, исчисления Маллявэна, математика, полиномы Эрмита, стохастический анализ
Винеровский хаос или Еще один способ подбросить монетку - 1

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

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

    $dX^2=2Xcdot dX.$

    И она верна. Вот только, если $X$ — это не случайный процесс. Адовая смесь из свойств нормального распределения и «почти наверное» доказывает, что в обратном случае эта формула в общем случае неверна. Томик мат.анализа с решениями обыкновенных дифференциальных уравнений теперь можно выкинуть в топку. Люди в теме тихо хихикают, остальные нетерпеливо листают статьи в Вики с исчислениями Ито.

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

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

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

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

Винеровский хаос или Еще один способ подбросить монетку - 6

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

  • $H_n'(x)=nH_{n-1}(x),$
  • $H_n(-x)=(-1)^nH_n(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,Xrangle}.$

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

  1. Самый простой пример — это пространство вещественных чисел: $H=mathbb{R}$. В таком случае скалярным произведением двух чисел $X$ и $Y$ у нас будет

    $langle X,Y rangle=XY.$

  2. Затем я перехожу в Эвклидово пространство $H=mathbb{R}^n$. Теперь

    $langle X,Yrangle=sum_{i=0}^nX_iY_i.$

    Это пространство можно расширить до пространства комплексных векторов: $H=mathbb{C}^n$, для которого скалярное произведение будет

    $langle X, Yrangle=sum_{i=0}^nX_ioverline{Y_i}$

    (верхняя черта обозначает комплексное сопряжение).

  3. Ну и наконец прихожу в пространство для взрослых, пространство с бесконечной размерностью. В нашем случае это будет пространство квадратично интегрируемых функций, заданных на некотором множестве $Omega$ с заданной мерой $mu$. Мы будем его обозначать в виде $H=L^2(Omega,mu)$. Скалярное произведение на нем задается следующим образом:

    $langle X, Y rangle=int_Omega (X cdot Y) dmu.$

    Обычно под множеством $Omega$ подразумевается интервал $[a,b]$, а под мерой $mu$ — равномерная мера (мера Лебега), т.е. $dmu=mu(domega)=domega$. И тогда скалярное произведение записывается в виде обыкновенного интеграла Лебега

    $int_a^bX(omega)Y(omega) domega.$

    Если же мы думаем в терминах теории вероятностей, то $Omega$ — это пространство элементарных событий, $X=X(omega)$ и $Y=Y(omega)$ — случайные величины, а $mu$ — вероятностная мера. У каждой такой меры есть своя функция плотности распределения $rho$, которая может быть отличной от константы, тогда $dmu=rho(omega)domega$ и скалярное произведение совпадает с математическим ожиданием:

    $langle X, Y rangle=int_Omega X(omega) Y(omega) rho(omega)domega=mathbb{E}[XY].$

Винеровский хаос или Еще один способ подбросить монетку - 42

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

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

  1. вектор из случайных величин $(W(h_1), dots, W(h_n))$ распределен нормально с нулевым мат.ожиданием для любых $h_1, dots h_n in H$, и
  2. для $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, grangle=int (f cdot g) dlambda.$

Пусть $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_jrangle$. На выход gaussianProcess выдает вектор $(W(h_1), dots, W(h_n))$.

Винеровский хаос или Еще один способ подбросить монетку - 70

Уже подходит время объединить все полученные нами знания вместе. Но прежде, стоит упомянуть об одном полезном свойстве эрмитовых полиномов и нормального распределения в совокупности. Пусть $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(stmathbb{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.$

Винеровский хаос или Еще один способ подбросить монетку - 83

Воу-воу, палехче! Давайте разложим эту теорему о разложении по кусочкам и переведём с математического на человеческий. Мы не будем сильно вдаваться в детали, а лишь интуитивно поясним о чем тут речь. Значок $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, то есть различные комбинации чисел $acdot 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 Wbigg(frac{1_{[0,t]}}{|1_{[0,t]} |}bigg)^2 \ &=tcdot Wbigg(frac{1_{[0,t]}}{sqrt{t}}bigg)^2\ &=tH_2bigg(Wbigg(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{2pi}}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{2pi}} 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$.

Винеровский хаос или Еще один способ подбросить монетку - 135

Судя по этому графику, где-то после $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]

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

Винеровский хаос или Еще один способ подбросить монетку - 138

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

Автор: The_Freeman

Источник

* - обязательные к заполнению поля

Главная   |  Архив новостей  |   Android  |   Google  |   Apple  |   Microsoft  |   Информационная безопасность  |   Веб – разработка
Публикации RSS  |  Комментарии RSS
https://ajax.googleapis.com/ajax/libs/jquery/3.4.1/jquery.min.js